inst/doc/custom_function_example.R

In steps: Spatially- and Temporally-Explicit Population Simulator

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  dpi = 150,
  fig.width = 7,
  out.width = "100%",
  cache = FALSE
)

## ---- message = FALSE---------------------------------------------------------
library(steps)
library(raster)
library(viridisLite)
library(future)

## ---- message = FALSE, eval = FALSE-------------------------------------------
#  
#  mortality <- function (mortality_layer, stages = NULL) {
#  
#    pop_dynamics <- function (landscape, timestep) {
#  
#        population_raster <- landscape$population
#        nstages <- raster::nlayers(population_raster)
#  
#        # get population as a matrix
#        idx <- which(!is.na(raster::getValues(population_raster[[1]])))
#        population_matrix <- raster::extract(population_raster, idx)
#  
#        mortality_prop <- raster::extract(landscape[[mortality_layer]][[timestep]], idx)
#  
#        if (is.null(stages)) stages <- seq_len(nstages)
#  
#        for (stage in stages) {
#  
#          population_matrix[ , stage] <- ceiling(population_matrix[ , stage] * mortality_prop)
#  
#        }
#  
#        # put back in the raster
#        population_raster[idx] <- population_matrix
#  
#        landscape$population <- population_raster
#  
#  
#      landscape
#  
#    }
#  
#    as.population_modification(pop_dynamics)
#  
#  }

## ---- message = FALSE, fig.align = "center"-----------------------------------
egk_pop

par(mar=c(0,0,0,0), oma=c(0,0,0,0))
spplot(egk_pop,
       col.regions = viridis(100),
       par.settings=list(strip.background = list(col = "white")))

## ---- message = FALSE---------------------------------------------------------
cellStats(egk_pop[[1]], sum)

## ---- message = FALSE, eval = FALSE-------------------------------------------
#  percent_mortality <- function (percentage, stages = NULL) {
#  
#    pop_dynamics <- function (landscape, timestep) { ...

## ---- message = FALSE, eval = FALSE-------------------------------------------
#  percent_mortality <- function (percentage, stages = NULL) {
#  
#    pop_dynamics <- function (landscape, timestep) {
#        population_raster <- landscape$population
#        nstages <- raster::nlayers(population_raster)
#  
#        # get population as a matrix
#        idx <- which(!is.na(raster::getValues(population_raster[[1]])))
#        population_matrix <- raster::extract(population_raster, idx) ...

## ---- message = FALSE, eval = FALSE-------------------------------------------
#  percent_mortality <- function (percentage, stages = NULL) {
#  
#    pop_dynamics <- function (landscape, timestep) {
#        population_raster <- landscape$population
#        nstages <- raster::nlayers(population_raster)
#  
#        # get population as a matrix
#        idx <- which(!is.na(raster::getValues(population_raster[[1]])))
#        population_matrix <- raster::extract(population_raster, idx)
#  
#        if (is.null(stages)) stages <- seq_len(nstages)
#        for (stage in stages) {
#          initial_pop <- sum(population_matrix[ , stage])
#          changing_pop <- sum(population_matrix[ , stage])
#          while (changing_pop > initial_pop * percentage) {
#            non_zero <- which(population_matrix[idx , stage] > 0)
#            i <- sample(non_zero, 1)
#            population_matrix[i , stage] <- population_matrix[i , stage] - 1
#            changing_pop <- sum(population_matrix[ , stage])
#          }
#        } ...

## ---- message = FALSE---------------------------------------------------------
percent_mortality <- function (percentage, stages = NULL) {
  
  pop_dynamics <- function (landscape, timestep) {
    population_raster <- landscape$population
    nstages <- raster::nlayers(population_raster)
    
    # get population as a matrix
    idx <- which(!is.na(raster::getValues(population_raster[[1]])))
    population_matrix <- raster::extract(population_raster, idx)
    
    if (is.null(stages)) stages <- seq_len(nstages)
    for (stage in stages) {
      initial_pop <- sum(population_matrix[ , stage])
      changing_pop <- sum(population_matrix[ , stage])
      while (changing_pop > initial_pop * (1 - percentage)) {
        non_zero <- which(population_matrix[ , stage] > 0)
        i <- sample(non_zero, 1)
        population_matrix[i , stage] <- population_matrix[i , stage] - 1
        changing_pop <- sum(population_matrix[ , stage])
      }
    }
    
    population_raster[idx] <- population_matrix
    landscape$population <- population_raster
    landscape
  }
} 

## ---- message = FALSE---------------------------------------------------------
ls <- landscape(population = egk_pop,
                suitability = NULL,
                carrying_capacity = NULL)

pd <- population_dynamics(change = NULL,
                          dispersal = NULL,
                          modification = percent_mortality(percentage = 0.1,
                                                           stages = 1),
                          density_dependence = NULL)

results <- simulation(landscape = ls,
                      population_dynamics = pd,
                      habitat_dynamics = NULL,
                      demo_stochasticity = "none",
                      timesteps = 3,
                      replicates = 1,
                      verbose = FALSE)

## ---- message = FALSE---------------------------------------------------------
cbind(c("t0", "t1", "t2", "t3"),
      c(cellStats(egk_pop[[1]], sum),
        cellStats(results[[1]][[1]][["population"]][[1]], sum),
        cellStats(results[[1]][[2]][["population"]][[1]], sum),
        cellStats(results[[1]][[3]][["population"]][[1]], sum)))


## ---- message = FALSE, fig.align = "center"-----------------------------------
pop_stack <- stack(egk_pop[[1]],
                   results[[1]][[1]][["population"]][[1]],
                   results[[1]][[2]][["population"]][[1]],
                   results[[1]][[3]][["population"]][[1]])

names(pop_stack) <- c("t0", "t1", "t2", "t3")

par(mar=c(0,0,0,0), oma=c(0,0,0,0))
spplot(pop_stack,
       col.regions = viridis(100),
       par.settings=list(strip.background = list(col = "white")))

## ---- message = FALSE---------------------------------------------------------
percent_mortality_hab <- function (percentage, threshold, stages = NULL) {
  
  pop_dynamics <- function (landscape, timestep) {
    population_raster <- landscape$population
    nstages <- raster::nlayers(population_raster)
    
    # get population as a matrix
    idx <- which(!is.na(raster::getValues(population_raster[[1]])))
    population_matrix <- raster::extract(population_raster, idx)
    
    # get suitability cell indexes
    suitable <- raster::extract(landscape$suitability, idx)
    
    if (is.null(stages)) stages <- seq_len(nstages)
    for (stage in stages) {
      initial_pop <- sum(population_matrix[ , stage])
      changing_pop <- sum(population_matrix[ , stage])
      
      highly_suitable <- which(suitable >= threshold) # check for suitable cells
      
      while (changing_pop > initial_pop * (1 - percentage)) {
        non_zero <- which(population_matrix[ , stage] > 0)
        i <- sample(intersect(highly_suitable, non_zero), 1) # change the sampling pool
        population_matrix[i , stage] <- population_matrix[i , stage] - 1
        changing_pop <- sum(population_matrix[ , stage])
      }
    }
    
    population_raster[idx] <- population_matrix
    landscape$population <- population_raster
    landscape
  }
} 

## ---- message = FALSE---------------------------------------------------------
ls <- landscape(population = egk_pop,
                suitability = egk_hab,
                carrying_capacity = NULL)

pd <- population_dynamics(change = NULL,
                          dispersal = NULL,
                          modification = percent_mortality_hab(percentage = .10,
                                                               threshold = 0.7,
                                                               stages = 1),
                          density_dependence = NULL)

results <- simulation(landscape = ls,
                      population_dynamics = pd,
                      habitat_dynamics = NULL,
                      demo_stochasticity = "none",
                      timesteps = 3,
                      replicates = 1,
                      verbose = FALSE)

## ---- message = FALSE---------------------------------------------------------
cbind(c("t0", "t1", "t2", "t3"),
      c(cellStats(egk_pop[[1]], sum),
        cellStats(results[[1]][[1]][["population"]][[1]], sum),
        cellStats(results[[1]][[2]][["population"]][[1]], sum),
        cellStats(results[[1]][[3]][["population"]][[1]], sum)))


## ---- message = FALSE, fig.align = "center"-----------------------------------
pop_stack <- stack(egk_pop[[1]],
                   results[[1]][[1]][["population"]][[1]],
                   results[[1]][[2]][["population"]][[1]],
                   results[[1]][[3]][["population"]][[1]])

names(pop_stack) <- c("t0", "t1", "t2", "t3")

par(mar=c(0,0,0,0), oma=c(0,0,0,0))
spplot(pop_stack,
       col.regions = viridis(100),
       par.settings=list(strip.background = list(col = "white")))

## ---- message = FALSE, fig.align = "center"-----------------------------------
sampling_mask <- egk_hab # copy the habitat layer to assume its properties
sampling_mask[!is.na(egk_hab)] <- 0 # set all non-NA values to zero

# get index values of the top-right corner (15 x 15 cells)
idx_corner <- sort(unlist(lapply(0:15, function (x) seq(21, 541, by = 36) + x)))

# set non-NA raster values in top-right corner to one
sampling_mask[intersect(idx_corner, which(!is.na(egk_hab[])))] <- 1

plot(sampling_mask, axes = FALSE, box = FALSE)

ls <- landscape(population = egk_pop,
                suitability = NULL,
                carrying_capacity = NULL,
                "sampling_mask" = sampling_mask) # name the object to reference in our function

## ---- message = FALSE---------------------------------------------------------
percent_mortality_ras <- function (percentage, masking_raster, threshold, stages = NULL) {
  
  pop_dynamics <- function (landscape, timestep) {
    population_raster <- landscape$population
    nstages <- raster::nlayers(population_raster)
    
    # get population as a matrix
    idx <- which(!is.na(raster::getValues(population_raster[[1]])))
    population_matrix <- raster::extract(population_raster, idx)
    
    # Note, here we now use the name of the raster ('masking_raster' parameter)
    # to extract the object from the landscape object
    suitable <- raster::extract(landscape[[masking_raster]], idx)
    
    if (is.null(stages)) stages <- seq_len(nstages)
    for (stage in stages) {
      initial_pop <- sum(population_matrix[ , stage])
      changing_pop <- sum(population_matrix[ , stage])
      
      highly_suitable <- which(suitable >= threshold) # check for suitable cells
      
      while (changing_pop > initial_pop * (1 - percentage)) {
        non_zero <- which(population_matrix[ , stage] > 0)
        i <- sample(intersect(highly_suitable, non_zero), 1) # change the sampling pool
        population_matrix[i , stage] <- population_matrix[i , stage] - 1
        changing_pop <- sum(population_matrix[ , stage])
      }
    }
    
    population_raster[idx] <- population_matrix
    landscape$population <- population_raster
    landscape
  }
} 

## ---- message = FALSE---------------------------------------------------------
pd <- population_dynamics(change = NULL,
                          dispersal = NULL,
                          modification = percent_mortality_ras(percentage = .10,
                                                               masking_raster = "sampling_mask",
                                                               threshold = 0.7,
                                                               stages = 1),
                          density_dependence = NULL)

results <- simulation(landscape = ls,
                      population_dynamics = pd,
                      habitat_dynamics = NULL,
                      demo_stochasticity = "none",
                      timesteps = 3,
                      replicates = 1,
                      verbose = FALSE)

## ---- message = FALSE---------------------------------------------------------
cbind(c("t0", "t1", "t2", "t3"),
      c(cellStats(egk_pop[[1]], sum),
        cellStats(results[[1]][[1]][["population"]][[1]], sum),
        cellStats(results[[1]][[2]][["population"]][[1]], sum),
        cellStats(results[[1]][[3]][["population"]][[1]], sum)))


## ---- message = FALSE, fig.align = "center"-----------------------------------
pop_stack <- stack(egk_pop[[1]],
                   results[[1]][[1]][["population"]][[1]],
                   results[[1]][[2]][["population"]][[1]],
                   results[[1]][[3]][["population"]][[1]])

names(pop_stack) <- c("t0", "t1", "t2", "t3")

par(mar=c(0,0,0,0), oma=c(0,0,0,0))
spplot(pop_stack,
       col.regions = viridis(100),
       par.settings=list(strip.background = list(col = "white")))

## ---- message = FALSE, eval = FALSE-------------------------------------------
#  modified_transition <- function(survival_layer = NULL,
#                                  fecundity_layer = NULL) {
#  
#    fun <- function (transition_array, landscape, timestep) {
#  
#      transition_matrix <- transition_array[, , 1]
#      idx <- which(transition_matrix != 0)
#      is_recruitment <- upper.tri(transition_matrix)[idx]
#  
#      array_length <- dim(transition_array)[3]
#  
#      cell_idx <- which(!is.na(raster::getValues(landscape$population[[1]])))
#  
#      if (is.null(survival_layer)) {
#        surv_mult <- rep(1, length(cell_idx))
#      } else {
#        if (raster::nlayers(landscape$suitability) > 1) {
#          surv_mult <- landscape[[survival_layer]][[timestep]][cell_idx]
#        } else {
#          surv_mult <- landscape[[survival_layer]][cell_idx]
#        }
#      }
#  
#      if (is.null(fecundity_layer)) {
#        fec_mult <- rep(1, length(cell_idx))
#      } else {
#        if (raster::nlayers(landscape$suitability) > 1) {
#          fec_mult <- landscape[[fecundity_layer]][[timestep]][cell_idx]
#        } else {
#          fec_mult <- landscape[[fecundity_layer]][cell_idx]
#        }
#      }
#  
#      for (i in seq_len(array_length)) {
#        transition_array[, , i][idx[!is_recruitment]] <- transition_array[, , i][idx[!is_recruitment]] * surv_mult[i]
#        transition_array[, , i][idx[is_recruitment]] <- transition_array[, , i][idx[is_recruitment]] * fec_mult[i]
#      }
#  
#      transition_array
#  
#    }
#  
#    as.transition_function(fun)
#  
#  }

## ---- message = FALSE---------------------------------------------------------
ls <- landscape(population = egk_pop,
                suitability = NULL,
                carrying_capacity = NULL)

pd <- population_dynamics(change = growth(transition_matrix = egk_mat),
                          dispersal = NULL,
                          modification = NULL,
                          density_dependence = NULL)

results <- simulation(landscape = ls,
                      population_dynamics = pd,
                      habitat_dynamics = NULL,
                      timesteps = 5,
                      replicates = 1,
                      verbose = FALSE)

## ---- message = FALSE, fig.align = "center"-----------------------------------
plot_pop_trend(results, stages = 1)

## ---- message = FALSE, fig.align = "center"-----------------------------------
plot_pop_spatial(results,
     stage = 1,
     timesteps = 1:5)

## ---- message = FALSE, echo = FALSE-------------------------------------------
print(egk_mat)

paste("Rmax =", abs(eigen(egk_mat)$values[1]))

## ---- message = FALSE, fig.align = "center", warning = FALSE------------------
adult_fec <- egk_hab # copy the habitat layer to assume its properties
adult_fec[] <- NA # set all values to NA

ncells <- ncell(adult_fec)

# assign random numbers between 0.1 and 0.5 to 80% of the cells
adult_fec[sample(seq_len(ncells), 0.8 * ncells, replace = FALSE)] <- runif(ncells, 0.1, 0.5)

plot(adult_fec, axes = FALSE, box = FALSE)

## ---- message = FALSE, eval = FALSE-------------------------------------------
#  define_fecundity <- function(fecundity_layer) {
#  
#    fun <- function (transition_array, landscape, timestep) { ...

## ---- message = FALSE, eval = FALSE-------------------------------------------
#  define_fecundity <- function(fecundity_layer) {
#  
#    fun <- function (transition_array, landscape, timestep) {
#      transition_matrix <- transition_array[, , 1]
#      idx <- which(transition_matrix != 0)
#      is_recruitment <- upper.tri(transition_matrix)[idx]
#      cell_idx <- which(!is.na(raster::getValues(landscape$population[[1]])
#  
#      # this length should be pre-defined by non-NA cells in the landscape
#      array_length <- dim(transition_array)[3]
#      ...

## ---- message = FALSE, eval = FALSE-------------------------------------------
#  define_fecundity <- function(fecundity_layer) {
#  
#    fun <- function (transition_array, landscape, timestep) {
#      transition_matrix <- transition_array[, , 1]
#      idx <- which(transition_matrix != 0)
#      is_recruitment <- upper.tri(transition_matrix)[idx]
#      cell_idx <- which(!is.na(raster::getValues(landscape$population[[1]])
#  
#      # this length should be pre-defined by non-NA cells in the landscape
#      array_length <- dim(transition_array)[3]
#  
#      fec_values <- landscape[[fecundity_layer]][cell_idx]
#  
#      for (i in seq_len(array_length)) {
#        new_fec <- fec_values[i]
#        if (!is.na(new_fec)) {
#          transition_array[, , i][tail(idx[is_recruitment], 1)] <- new_fec
#        }
#      } ...

## ---- message = FALSE---------------------------------------------------------
define_fecundity <- function(fecundity_layer) {
  
  fun <- function (transition_array, landscape, timestep) {
    transition_matrix <- transition_array[, , 1]
    idx <- which(transition_matrix != 0)
    is_recruitment <- upper.tri(transition_matrix)[idx]
    cell_idx <- which(!is.na(raster::getValues(landscape$population[[1]])))
    
    # this length should be pre-defined by non-NA cells in the landscape
    array_length <- dim(transition_array)[3]
          
    fec_values <- landscape[[fecundity_layer]][cell_idx]
    
    for (i in seq_len(array_length)) {
      new_fec <- fec_values[i]
      if (!is.na(new_fec)) {
        transition_array[, , i][tail(idx[is_recruitment], 1)] <- new_fec
      }
    }

    transition_array
  }
}

## ---- message = FALSE---------------------------------------------------------
ls <- landscape(population = egk_pop,
                suitability = NULL,
                carrying_capacity = NULL,
                "adult_fec" = adult_fec)

pd <- population_dynamics(
  change = growth(transition_matrix = egk_mat,
                  transition_function = define_fecundity(fecundity_layer = "adult_fec")),
  dispersal = NULL,
  modification = NULL,
  density_dependence = NULL)

results <- simulation(landscape = ls,
                      population_dynamics = pd,
                      habitat_dynamics = NULL,
                      timesteps = 5,
                      replicates = 1,
                      verbose = FALSE)

## ---- message = FALSE, fig.align = "center"-----------------------------------
plot_pop_trend(results, stages = 1)

## ---- message = FALSE, fig.align = "center"-----------------------------------
plot_pop_spatial(results,
     stage = 1,
     timesteps = 1:5)

## ---- message = FALSE, eval = FALSE-------------------------------------------
#  exponential_dispersal_kernel <- function (distance_decay = 0.5, normalize = FALSE) {
#    if (normalize) {
#      fun <- function (r) (1 / (2 * pi * distance_decay ^ 2)) * exp(-r / distance_decay)
#    } else {
#      fun <- function (r) exp(-r / distance_decay)
#    }
#    as.dispersal_function(fun)
#  }

## ---- message = FALSE---------------------------------------------------------
r <- seq(1, 1000, 10)
distance_decay = 100
plot(r, exp(-r / distance_decay), type = 'l')

## ---- message = FALSE---------------------------------------------------------
logistic_dispersal_kernel <- function (dist_shape, dist_slope) {
fun <- function (dist) 1 / (1 + exp((dist - dist_shape) / dist_slope))
fun
}

r <- seq(1, 1000, 10)
dist_shape = 500
dist_slope = 20
plot(r, 1 / (1 + exp((r - dist_shape) / dist_slope)), type = 'l')

## ---- message = FALSE---------------------------------------------------------
egk_pop_thin <- egk_pop
egk_pop_thin[sample(1:ncell(egk_pop), 0.99 * ncell(egk_pop))] <- 0

ls <- landscape(population = egk_pop_thin,
                suitability = NULL,
                carrying_capacity = NULL)

pd_exp_disp <- population_dynamics(
  change = NULL,
  dispersal = kernel_dispersal(dispersal_kernel = exponential_dispersal_kernel(distance_decay = 100),
                               max_distance = 1000,
                               arrival_probability = "none"),
  modification = NULL,
  density_dependence = NULL)

results_01 <- simulation(landscape = ls,
                      population_dynamics = pd_exp_disp,
                      habitat_dynamics = NULL,
                      timesteps = 10,
                      replicates = 1,
                      verbose = FALSE)

plot_pop_spatial(results_01, stage = 3, timesteps = c(1, 5, 10))

pd_log_disp <- population_dynamics(
  change = NULL,
  dispersal = kernel_dispersal(dispersal_kernel = logistic_dispersal_kernel(dist_shape = 500, dist_slope = 20),
                               max_distance = 1000,
                               arrival_probability = "none"),
  modification = NULL,
  density_dependence = NULL)

results_02 <- simulation(landscape = ls,
                      population_dynamics = pd_log_disp,
                      habitat_dynamics = NULL,
                      timesteps = 10,
                      replicates = 1,
                      verbose = FALSE)

plot_pop_spatial(results_02, stage = 3, timesteps = c(1, 5, 10))

## ---- message = FALSE, eval = FALSE-------------------------------------------
#  plan(multisession)
#  results_02 <- simulation(landscape = ls,
#                        population_dynamics = pd_log_disp,
#                        habitat_dynamics = NULL,
#                        timesteps = 10,
#                        replicates = 3,
#                        verbose = FALSE,
#                        future.globals = list("logistic_dispersal_kernel" = logistic_dispersal_kernel))