R/forest.transition.r

In nuaquilue/SBW: Spruce budworm outbreaks model

#' Forest transitions
#'
#' Determines the presence of species within a buffer around target cells.
#' 
#' @param land A \code{landscape} data frame with forest stand records in rows
#' @param target.cell.ids A vector of \code{cell.id} codes for those cells that may change stand composition
#' @param suitab A data frame with climatic suitability and soil suitability for potential colonizing species 
#' for each forest stand. It is the result of the function \code{suitability()}
#' @param params A list of default parameters generated by the function \code{default.params()} or 
#' a customized list of model parameters
#' @param tbls A list of default input tables as in \code{data(default.tables()} or 
#' a customized list of input tables
#' @param type.trans A character indicating the type of transition: \code{B} for burnt cells, 
#'  \code{O} for cells affected by an outbreak, \code{B} for clear-cut cells, and 
#'  \code{S} for cells undergoing natural succession.
#' 
#' @return A vector with the name of the new species / group of species
#' 
#' @export
#' 
#' @examples
#' data(landscape)
#' params = default.params()
#' suitab = suitability(landscape, params) 
#' forest.trans(landscape, landscape$cell.id[runif(10,1,nrow(landscape))], suitab, params, type.trans = "S")
#' 

forest.transition = function(land, target.cells, suitab, params, tbls, type.trans = "S"){   
  
  ## If target data.frame is empty
  if(length(target.cells)==0)
    return(numeric())
  
  ## Input tables
  spp.colonize.persist = tbls$spp.colonize.persist
  post.sbw.reg = tbls$post.sbw.reg
  forest.succ = tbls$forest.succ
    
  ## Tracking
  cat(ifelse(type.trans=="B", "Forest regeneration after fires", 
        ifelse(type.trans=="O", "Forest regeneration after outbreaks", 
          ifelse(type.trans=="C", "Forest regeneration after clear-cutting", 
            ifelse(type.trans=="S", "Seral succession", stop("Transition non-suported"))))), "\n")
  
  
  ## Probability of regeneration
  if(type.trans=="B")
    prob.reg = post.fire.reg
  if(type.trans=="O")
    prob.reg = post.sbw.reg
  if(type.trans=="C")
    prob.reg = post.harvest.reg
  if(type.trans=="S")
    prob.reg = forest.succ
  
  ## Compute buffer of neighours for target cells
  buffer = buffer.mig(land, target.cells, tbls)
  
  ## Keep the current species in case any potential species can colonize the site.
  ## In that case, the current species, persist.
  subland = filter(land, cell.id %in% target.cells)
  subland$spp = as.character(subland$spp)
  current.spp = subland$spp
  
  ## Join to the subland data frame the probability of transition to potential.spp (according to initial spp)
  ## Then join buffer results indicating whether the potential species is present in the surrounding neighborhood
  ## Finally join the climatic-soil suitability index (i.e. modifier) of the potential spp
  subland = left_join(subland, prob.reg, by="spp") %>%
             left_join(buffer, by=c("cell.id", "potential.spp")) %>%
             left_join(suitab, by=c("cell.id", "potential.spp")) 
  
  ## Reset suitability for 'other species' because the group includes many species and 
  ## it's assumed that the climate would be suitable for at least one of those species. 
  subland$suit.clim[subland$potential.spp=="OTH"] = 1
  subland$suit.clim[subland$potential.spp=="NonFor"] = 1
  subland$suit.soil[subland$potential.spp=="NonFor"] = 1
  
  ## Eufeuillement volontaire suite aux coupes, that is, after clear-cut for a % of EPN and SAB
  ## to transform to PET
  if(type.trans=="C" & params$enfeuil>0){  
    vec.enfeuil = filter(subland, spp %in% c("EPN", "SAB") & potential.spp=="PET") %>% select(ptrans)
    vec.enfeuil = vec.enfeuil + (runif(length(vec.enfeuil)) < params$enfeuil)*1000
    subland$ptrans[subland$spp %in% c("EPN", "SAB") & subland$potential.spp=="PET"] = unlist(vec.enfeuil)
  }
  
  ## Reburning case: If burnt stands are too young, probability of successful natural regeneration is lower
  if(type.trans=="B"){ 
    subland$ptrans[subland$spp=="EPN" & subland$potential.spp=="EPN" & subland$age<params$age.seed] = 
      subland$ptrans[subland$spp=="EPN" & subland$potential.spp=="EPN" & subland$age<params$age.seed] * (1-params$p.failure)
  }
  
  ## Stability criteria: if the species is present in the target location, then
  ## soil conditions are assumed to be optimal (not limiting)
      # levels(subland$spp) = levels(subland$potential.spp)
  subland$press.buffer = (subland$spp == subland$potential.spp) | (subland$press.buffer)
  subland$suit.soil[subland$spp == subland$potential.spp] = 1

  ## Species persistence when climatic conditions become unfavorable: when persistence is allowed (1), 
  ## there is a floor probability of self-replacement corresponding to sub-optimal conditions 
  ## (under the assumption that competition is more limiting than  physiological response to climate)
  ## First, find which spp are allowed to persist then, upgrade climatic suitability to suboptimal
  ## in case this is lower than suboptimal.
  spp.persist = spp.colonize.persist$spp[spp.colonize.persist$persist==1]
  subland$suit.clim[subland$spp %in% spp.persist & 
                     subland$spp == subland$potential.spp & subland$suit.clim<params$suboptimal] = params$suboptimal
  
  ## Determine the final succession / regeneration probability 
  subland$p = subland$ptrans * subland$press.buffer * pmin(subland$suit.clim, subland$suit.soil)
  
  ## Reshape the data frame, so we have a column for each potential species with 
  ## the corresponding transition probability (one row per target cell)
  ## Substitute dcast by "gather" or "spread" from tidyverse
  aux = reshape2::dcast(subland, formula = cell.id ~ potential.spp, value.var = "p")
  
  ## Function that returns spp id according to probability x
  select.spp = function(x){
    if(sum(x)==0)
      return(0)
    id.spp = sample(1:length(x), 1, replace=FALSE, prob=x)
    return(id.spp)
  }
  
  ## Now select a new spp according to these probabilities and assign the corresponing species name
  ## If after all filters, p for all potential.spp is 0, the current species remains
  spp.names = names(aux)[2:ncol(aux)]
  id.spp = apply(aux[,2:ncol(aux)], 1, select.spp)
  new.spp = numeric(length=length(id.spp))
  new.spp[id.spp!=0] = spp.names[id.spp[id.spp!=0]]
  new.spp[id.spp==0] = as.character(current.spp[id.spp==0])
  # pour les cellules deja dominées par other, revenir à la même chose
  new.spp[new.spp=="OTH" & current.spp %in% c("OTH.RES.N","OTH.RES.S","OTH.FEU.N","OTH.FEU.S")] = 
   current.spp[new.spp=="OTH" & current.spp %in% c("OTH.RES.N","OTH.RES.S","OTH.FEU.N","OTH.FEU.S")]   
  if(any(new.spp=="OTH")){
    new.spp[new.spp=="OTH"] = .select.others(land, unique(subland$cell.id)[new.spp=="OTH"])
  }
  
  ## Return the vector with the name of the new spp
  return(new.spp)
}


.select.others = function(land, target.cells.oth){
  
  target.cells.oth.xy = land[land$cell.id %in% target.cells.oth,c("x","y")]
  
  micro.OthCB = land[land$spp%in% c("OTH.RES.N"),]
  list.cell.buff = nn2(micro.OthCB[,c("x","y")], target.cells.oth.xy, k=40 , searchtype='priority')
  vec.OthCB = rowSums(list.cell.buff$nn.dists < 20000)  
  
  micro.OthCT = land[land$spp%in% c("OTH.RES.S"),]
  list.cell.buff = nn2(micro.OthCT[,c("x","y")], target.cells.oth.xy, k=40 , searchtype='priority')
  vec.OthCT = rowSums(list.cell.buff$nn.dists < 20000)
  
  micro.OthDB = land[land$spp%in% c("OTH.FEU.N"),]
  list.cell.buff = nn2(micro.OthDB[,c("x","y")], target.cells.oth.xy, k=40 , searchtype='priority')
  vec.OthDB = rowSums(list.cell.buff$nn.dists < 20000)
  
  micro.OthDT = land[land$spp%in% c("OTH.FEU.S","ERS","BOJ"),]
  list.cell.buff = nn2(micro.OthDT[,c("x","y")],target.cells.oth.xy, k=40 , searchtype='priority')
  vec.OthDT = rowSums(list.cell.buff$nn.dists < 20000)
  
  count.oth =   cbind(vec.OthCB,vec.OthCT,vec.OthDB,vec.OthDT)
  others = c("OTH.RES.N","OTH.RES.S","OTH.FEU.N","OTH.FEU.S")
  
  # Choose others according to the number of others in the neighborhood
  sample.oth = function(x){
    if(sum(x)==0) # if there aren't others in the neigh, assign either OthCB or OthDB,
      # that are everywhere and in equal proportion
      return(sample(c(1,3), 1, replace=F))
    else
      return(sample(1:4, 1, replace=F, prob=x))
  }
  spp = others[apply(count.oth, 1, sample.oth)]
  
  return(spp)
}