In Sz-Tim/gbPopMod: Mechanistic species distribution model for glossy buckthorn (*Frangula alnus*)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This is a wrapper for running the glossy buckthorn cellular automata model. It sets up the environment, allows the assignment of all necessary parameters, runs the simulation, and visualizes the results.

The land cover within each cell is compositional rather than the more common hard classification. As a result, all land cover-specific parameters are computed proportionally to the land cover in a cell (i.e., lc.prop %*% K, where lc.prop is a matrix with columns for land cover and a row for each cell, and K is a vector of carrying capacities with one per land cover type).

if(!require(pacman)) {install.packages("pacman"); require(pacman)}
p_load("tidyverse", "magrittr", "stringr", "parallel", "here", "devtools",
       "pbapply", "fastmatch", "gganimate"); theme_set(theme_bw())
library(gbPopMod)
data(lc.rct)
### lc.rct was produced from the raw GRANIT output as follows:
# grnt <- read_csv(here::here("data", "grnt_all.csv")) %>%
#   mutate(x=as.numeric(as.factor(.$left)),
#          y=as.numeric(factor(.$top, levels=rev(levels(factor(.$top))))),
#          x_y=paste(x, y, sep="_"))
# lc.rct <- as.tibble(expand.grid(x=1:max(grnt$x),
#                                 y=1:max(grnt$y))) %>%
#   mutate(x_y=paste(x, y, sep="_")) %>%
#   mutate(OpI=grnt$OpI[match(.$x_y, grnt$x_y)], 
#          Oth=grnt$Oth[match(.$x_y, grnt$x_y)], 
#          Dec=grnt$Dec[match(.$x_y, grnt$x_y)], 
#          Evg=grnt$Evg[match(.$x_y, grnt$x_y)], 
#          WP=grnt$WP[match(.$x_y, grnt$x_y)], 
#          Mxd=grnt$Mxd[match(.$x_y, grnt$x_y)],
#          inbd=!is.na(match(.$x_y, grnt$x_y)))
# lc.rct[is.na(lc.rct)] <- 0

##---
## set parameters
##---

  set.seed(225)
  n.sim <- 4; prl.sim <- FALSE
  g.p <- set_g_p(tmax=50, lc.r=50, lc.c=50, n.ldd=1, n.cores=4)
  control.p <- set_control_p()


##--
## initialize landscape
##--

  # land cover
  lc.df <- lc.rct %>% 
    filter(y >= (max(lc.rct$y) - g.p$lc.r) & x <= g.p$lc.c) %>%
    mutate(id=row_number(), 
           id.in=min_rank(na_if(inbd*id, 0)))
  ngrid <- nrow(lc.df)
  ncell <- sum(lc.df$inbd)

  # populations
  # indexed 1:ngrid including rows for out-of-bounds cells
  N.init <- pop_init(ngrid, g.p, lc.df)

  # dispersal probabilities
  # indexed 1:ncell with only inbound cells: sdd.pr[,,,i] = ith inbound cell
  # row i = dplyr::filter(lc.df, inbds)[i,] = lc.df[id.in==i,]
  # cell indexes in sdd.pr[,,2,] correspond with lc.df$id
  # this extra complication avoids calculating many, many unnecessary cells
  sdd.pr <- sdd_set_probs(ncell, lc.df, g.p)

if(prl.sim) {
  require(doParallel)
  p.c <- makeCluster(g.p$n.cores)
  registerDoParallel(p.c)
  out.p <- foreach(i=1:n.sim) %dopar% {
    gbPopMod::run_sim(ngrid, ncell, g.p, lc.df, sdd.pr, 
                      N.init, control.p, verbose=F)
  }
  stopCluster(p.c)
  out.ad <- map(out.p, ~.$N[,,max(g.p$age.f)]) %>% unlist %>% 
    array(., dim=c(ngrid, g.p$tmax+1, n.sim))
  out.sb <- map(out.p, ~.$N.sb) %>% unlist %>% 
    array(., dim=c(ngrid, g.p$tmax+1, n.sim))
} else {
  out.ad <- array(dim=c(ngrid, g.p$tmax+1, n.sim))
  out.sb <- array(dim=c(ngrid, g.p$tmax+1, n.sim))
  for(s in 1:n.sim) {
    out <- run_sim(ngrid, ncell, g.p, lc.df, sdd.pr, N.init, 
                   control.p, verbose=F)
    out.ad[,,s] <- out$N[,,max(g.p$age.f)]
    out.sb[,,s] <- out$N.sb
    rm(out)
    cat("Finished simulation", s, "of", n.sim, "\n")
  }
}

ad.mn <- apply(out.ad, 1:2, mean)
N.out <- cbind(lc.df, ad.mn) %>% as.tibble %>%
  gather(year, ad_Ab, (1:ncol(ad.mn)) + ncol(lc.df)) %>%
  mutate(ad_sd=c(apply(out.ad, 1:2, sd)),
         ad_pP=c(apply(out.ad > 0, 1:2, mean)),
         sb_Ab=c(apply(log(out.sb+1), 1:2, mean)),
         sb_sd=c(apply(log(out.sb+1), 1:2, sd)),
         sb_pP=c(apply(out.sb > 0, 1:2, mean)),
         ad_L5=c(ad.mn > 0))
N.out$ad_L5 <- N.out$ad_L5 + c(ad.mn > 5)
N.out$year <- str_pad(N.out$year, 3, "left", "0")

# prepare plot paths & filenames
p.wd <- paste0("out/", ncell, "_t", g.p$tmax, "_", g.p$edges, "/")
if(!dir.exists(here::here("out"))) dir.create(here::here("out"))
if(!dir.exists(here::here(p.wd))) dir.create(here::here(p.wd))

# save plots
save_pars(p.wd, g.p, control.p)
make_plots_final_t(p.wd, g.p, filter(N.out, year==max(N.out$year)))
make_plots_gifs(p.wd, g.p, N.out)
make_plots_lc(p.wd, lc.df)

Sz-Tim/gbPopMod documentation built on Dec. 7, 2020, 1:07 p.m.