R/cohort.establish.r

In nuaquilue/medLDM: Mediterranean Landscape Dynamic Model

#' Cohort establishment 
#'
#' Determines the tree species after forest dieback
#' 
#' @param land A \code{landscape} data frame with forest stand records and land-cover types in rows
#' @param clim A data frame with minimum temperature (in ºC), maximum temperature (in ºC), 
#' and precipitation (in mm) per location
#' @param sdm A data frame with categorical species distribution model (0 or 1) for all tree species
#' (in columns) for all the locations in the study area. It is generated wiht the function \code{sdm.sqi}
#' @param params A list of model parameters, default ones are generated by the function \code{default.params()}  
#' 
#' @return A data frame with the new species and sqi for after drought-induced mortality
#' 
#' @export
#' 
#' @examples
#' data(landscape)
#' data(clim)
#' land = dplyr::left_join(landscape, sdm.sqi(landscape, clim), by="cell.id")
#' killed.cells = drought(land, 10, 1)
#' land$tsdist[land$cell.id %in% killed.cells] = 0
#' land$typdist[land$cell.id %in% killed.cells] = "drght"
#' params = default.params()
#' x = cohort.establish(land, clim, params)
#' 

cohort.establish = function(land, clim, sdm, params){
  
  ## Tracking
  cat("Cohort establishment", "\n") 

  ## Num of neighbours in a circular neighbourhood according to radius (radius is in pixels)
  ## Assume that the neighbourhood is a star, with the maximum number of pixels in the
  ## east-west or north-south direction is 2*radius + 1 (1 is the center cell).
  ## The num of pixels is sequentially: 3+1*2, 5+3*2+1*2, 7+5*2+3*2+1*2, ...
  nneigh = seq(3,201,2) + cumsum(seq(1,200,2)*2)

  ## Coordinates of killed cells, then find their closest neighbours 
  ## Important to sort killed.cells by cell.id, so the order of the list of neigh matches
  killed.cells = filter(land, tsdist==0, typdist=="drght") %>% left_join(coord, by = "cell.id") 
  killed.cells = killed.cells[order(killed.cells$cell.id),] 
  neighs = nn2(coord[,-1], select(killed.cells,x,y),  searchtype="priority", k=nneigh[params$colon.rad])
  
  ## Retrive the species of the nneigh neighbors 
  neigh.spp = matrix(land$spp[neighs$nn.idx[,-1]], nrow=nrow(killed.cells), ncol=nneigh[params$colon.rad]-1)*
              ifelse(matrix(neighs$nn.dists[,-1], nrow=nrow(killed.cells), ncol=nneigh[params$colon.rad]-1)<=params$colon.rad*100,1,0) 
  
  ## Check if a species is present (TRUE) or absent (FALSE) in the neigbor
  neigh.spp = t(apply(neigh.spp, 1, .count.spp))>=1 
  neigh.spp = as.data.frame(neigh.spp)
  names(neigh.spp) = paste0("X", 1:ncol(neigh.spp))
  
  ## Combine both data frames, and assume that shrubs are always presence in the neighbourhood
  killed.cells = cbind(killed.cells, neigh.spp)
  killed.cells$X14 = T
  
  ## Look up cells killed by drought, add sqi data, then add the sencondary species
  ## (according to dominant spp and sqi), then add sdm of all tree species and finally
  ## add the number of forest spp in the neighbourhood
  killed.cells = killed.cells %>% left_join(secondary.spp, by = c("spp", "sqi")) %>% 
    left_join(sdm, by = "cell.id") 
  
  ## Select spp among available
  new.cohort = data.frame(cell.id = killed.cells$cell.id,
                          spp = apply(select(killed.cells, phalepensis:shrub) * 
                                        select(killed.cells, spp1:spp14) * 
                                          select(killed.cells, X1:X14), 1, .select.cohort))
  
  ## Join climatic and orographic variables to compute sq and then sqi
  new.cohort = new.cohort %>% left_join(clim, by = "cell.id") %>% left_join(orography, by = "cell.id") %>% 
               mutate(spp=spp.x) %>% select(-spp.x, -spp.y) %>% 
               left_join(sq.tree, by = "spp") %>%
               mutate(aux = c0+c_mnan*tmin+c2_mnan*tmin*tmin+c_plan*precip+c2_plan*precip*precip+
                      c_aspect*ifelse(aspect!=1,0,1)+c_slope*slope.stand) %>%
               mutate(sq=1/(1+exp(-1*aux))) %>% mutate(sqi=ifelse(sq<=p50, 1, ifelse(sq<=p90, 2, 3))) 
  shrub.cells = new.cohort$cell.id[new.cohort$spp==14]
  if(length(shrub.cells)>0){
    sqi.shrub = filter(clim, cell.id %in% shrub.cells) %>% 
                mutate(aux.brolla=sq.shrub$c0_brolla+sq.shrub$c_temp_brolla*tmin+sq.shrub$c_temp2_brolla*tmin*tmin+sq.shrub$c_precip_brolla*precip+sq.shrub$c_precip2_brolla*precip*precip,
                       aux.maquia=sq.shrub$c0_maquia+sq.shrub$c_temp_maquia*tmin+sq.shrub$c_temp2_maquia*tmin*tmin+sq.shrub$c_precip_maquia*precip+sq.shrub$c_precip2_maquia*precip*precip,
                       aux.boix=sq.shrub$c0_boix+sq.shrub$c_temp_boix*tmin+sq.shrub$c_temp2_boix*tmin*tmin+sq.shrub$c_precip_boix*precip+sq.shrub$c_precip2_boix*precip*precip,
                       sq.brolla=1/(1+exp(-1*aux.brolla)), sq.maquia=1/(1+exp(-1*aux.maquia)), sq.boix=1/(1+exp(-1*aux.boix))) %>% 
                mutate(sqest.brolla=sq.brolla/max(sq.brolla), sqest.maquia=sq.maquia/max(sq.maquia), sqest.boix=sq.boix/max(sq.boix),
                       sqi=ifelse(sqest.brolla>=sqest.maquia & sqest.brolla>=sqest.boix, 1,
                            ifelse(sqest.maquia>=sqest.brolla & sqest.maquia>=sqest.boix, 2,
                               ifelse(sqest.boix>=sqest.brolla & sqest.boix>=sqest.maquia, 3, 0))))
    new.cohort$sqi[new.cohort$spp==14] = sqi.shrub$sqi  
  }
  
  return(select(new.cohort, cell.id, spp, sqi))
}