## ----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--------------------------------------------------------- egk_mat <- matrix(c(0.00,0.55,0.85, 0.50,0.20,0.00, 0.00,0.55,0.85), nrow = 3, ncol = 3, byrow = TRUE) colnames(egk_mat) <- rownames(egk_mat) <- c('juvenile','subadult','adult') egk_mat_stoch <- matrix(c(0.00,0.05,0.05, 0.05,0.05,0.00, 0.00,0.01,0.01), nrow = 3, ncol = 3, byrow = TRUE) colnames(egk_mat_stoch) <- rownames(egk_mat_stoch) <- c('juvenile','subadult','adult') ## ---- message = FALSE, fig.align = "center"----------------------------------- egk_hab par(mar=c(0,0,0,0), oma=c(0,0,0,0)) plot(egk_hab, box = FALSE, axes = FALSE, col = viridis(100), main = "Habitat Suitability") ## ---- message = FALSE, fig.align = "center"----------------------------------- egk_k par(mar=c(0,0,0,0), oma=c(0,0,0,0)) plot(egk_k, box = FALSE, axes = FALSE, col = viridis(100)) ## ---- 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)) ## ---- message = FALSE--------------------------------------------------------- library(foreach) stable_states <- abs( eigen(egk_mat)$vectors[,1] / sum(eigen(egk_mat)$vectors[,1])) popN <- stack(replicate(ncol(egk_mat), egk_k)) * stable_states idx <- which(!is.na(getValues(popN[[1]]))) pop <- stack( foreach(i = 1:nlayers(popN)) %do% { max_pop <- ceiling(cellStats(popN[[i]], max, na.rm = T)) pop_values <- popN[[i]][idx] popN[[i]][idx] <- rbinom(prob = (pop_values/max_pop), size = max_pop, n = length(pop_values)) popN[[i]] }) names(pop) <- colnames(egk_mat) pop ## ---- message = FALSE--------------------------------------------------------- egk_landscape <- landscape(population = egk_pop, suitability = NULL, carrying_capacity = NULL) ## ---- message = FALSE--------------------------------------------------------- egk_pop_dynamics <- population_dynamics(change = growth(transition_matrix = egk_mat), dispersal = NULL, modification = NULL, density_dependence = NULL) ## ---- message = FALSE, results = 'hide', progress = FALSE-------------------- egk_results <- simulation(landscape = egk_landscape, population_dynamics = egk_pop_dynamics, habitat_dynamics = NULL, timesteps = 20, replicates = 1, verbose = FALSE) ## ---- message = FALSE, fig.align = "center"----------------------------------- plot_pop_trend(egk_results) ## ---- message = FALSE, fig.align = "center", fig.width = 3, fig.height = 4, out.width = "33%"---- plot_pop_trend(egk_results, stage = 0) ## ---- message = FALSE, fig.align = "center", fig.width = 3, fig.height = 4, out.width = "33%"---- plot_pop_trend(egk_results, stage = 2) ## ---- message = FALSE, fig.align = "center", fig.show='hold'------------------ plot_pop_spatial(egk_results, stage = 2, timesteps = c(1, 10, 20)) ## ---- message = FALSE, eval = FALSE------------------------------------------- # plot(egk_results) ## ---- message = FALSE, progress = FALSE, eval = FALSE------------------------- # plan(multisession, workers = 3) # This is how we specify to simulate # # replicates on separate processors in parallel # # egk_results <- simulation(landscape = egk_landscape, # population_dynamics = egk_pop_dynamics, # habitat_dynamics = NULL, # timesteps = 20, # replicates = 3, # verbose = FALSE) # ## ---- message = FALSE, fig.align = "center"----------------------------------- plot_pop_spatial(egk_results, replicate = 1, timesteps = c(1, 10, 20)) ## ---- message = FALSE, progress = FALSE, fig.align = "center"----------------- egk_pop_dynamics <- population_dynamics(change = growth(transition_matrix = egk_mat, global_stochasticity = egk_mat_stoch), dispersal = NULL, modification = NULL, density_dependence = NULL) egk_results <- simulation(landscape = egk_landscape, population_dynamics = egk_pop_dynamics, habitat_dynamics = NULL, demo_stochasticity = "none", timesteps = 20, replicates = 3, verbose = FALSE) plot_pop_trend(egk_results) ## ---- message = FALSE, progress = FALSE, fig.align = "center"----------------- plan(sequential) egk_landscape <- landscape(population = egk_pop, suitability = egk_hab, carrying_capacity = NULL) egk_pop_dynamics <- population_dynamics( change = growth(transition_matrix = egk_mat, global_stochasticity = egk_mat_stoch, transition_function = modified_transition(survival_layer = "suitability", fecundity_layer = NULL)), dispersal = NULL, modification = NULL, density_dependence = NULL) egk_results <- simulation(landscape = egk_landscape, population_dynamics = egk_pop_dynamics, habitat_dynamics = NULL, timesteps = 20, replicates = 3, verbose = FALSE) plot_pop_trend(egk_results) ## ---- message = FALSE, progress = FALSE, fig.align = "center"----------------- deterministic_transitions <- function(spatial_object) { fun <- function (transition_array, landscape, timestep) { #### This assumes that the spatial layers in the raster stack are named with #### a particular convention and will not work for all cases. # get metrics and constructor info current_timestep <- sprintf("%02i", timestep) cell_idx <- which(!is.na(raster::getValues(landscape$population[[1]]))) n_cells <- length(which(!is.na(raster::getValues(landscape$population[[1]])))) # get names of rasters and subset by timestep names_all <- names(landscape[[spatial_object]]) names_timestep <- grep(paste0("^.*", current_timestep, "\\_"), names_all, value = TRUE) # populate array: for (name in names_timestep) { r <- as.integer(substr(substr(name, nchar(name) - 3, nchar(name)), 1, 2)) c <- as.integer(substr(substr(name, nchar(name) - 1, nchar(name)), 1, 2)) transition_array[r, c, ] <- landscape[[spatial_object]][[name]][cell_idx] } #### Return array with required dimensions transition_array } } egk_landscape <- landscape(population = egk_pop, suitability = egk_hab, carrying_capacity = NULL, "sf_mods" = egk_sf) egk_pop_dynamics <- population_dynamics( change = growth(transition_matrix = egk_mat, transition_function = deterministic_transitions(spatial_object = "sf_mods")), dispersal = NULL, modification = NULL, density_dependence = NULL) egk_results <- simulation(landscape = egk_landscape, population_dynamics = egk_pop_dynamics, habitat_dynamics = NULL, timesteps = 20, verbose = FALSE) plot_pop_trend(egk_results) ## ---- message = FALSE, progress = FALSE, fig.align = "center", eval = FALSE---- # egk_landscape <- landscape(population = egk_pop, # suitability = egk_hab, # carrying_capacity = egk_k) # # egk_pop_dynamics <- population_dynamics(change = growth(transition_matrix = egk_mat, # global_stochasticity = egk_mat_stoch), # dispersal = NULL, # modification = NULL, # density_dependence = ceiling_density(stages = 3)) # # egk_results <- simulation(landscape = egk_landscape, # population_dynamics = egk_pop_dynamics, # habitat_dynamics = NULL, # timesteps = 20, # replicates = 3, # verbose = FALSE) ## ---- message = FALSE, progress = FALSE, fig.align = "center", eval = FALSE---- # egk_landscape <- landscape(population = egk_pop, # suitability = NULL, # carrying_capacity = NULL, # "origins" = egk_origins, # "destinations" = egk_destinations) # # egk_pop_dynamics <- population_dynamics( # change = growth(transition_matrix = egk_mat, # global_stochasticity = egk_mat_stoch), # dispersal = NULL, # modification = translocation(origins_layer = "origins", # destinations_layer = "destinations", # stages = 3, # effect_timesteps = c(1, 5, 10, 15)), # density_dependence = NULL) # # egk_results <- simulation(landscape = egk_landscape, # population_dynamics = egk_pop_dynamics, # habitat_dynamics = NULL, # timesteps = 20, # replicates = 3, # verbose = FALSE) ## ---- warning = FALSE, error = FALSE, message = FALSE, progress = FALSE, fig.align = "center"---- egk_landscape <- landscape(population = egk_pop, suitability = egk_hab, carrying_capacity = egk_k) egk_pop_dynamics <- population_dynamics( change = growth(transition_matrix = egk_mat, global_stochasticity = egk_mat_stoch), dispersal = kernel_dispersal( arrival_probability = "suitability", max_distance = 5000, dispersal_kernel = exponential_dispersal_kernel(distance_decay = 5000) ), density_dependence = ceiling_density(stages = 3)) egk_results <- simulation(landscape = egk_landscape, population_dynamics = egk_pop_dynamics, habitat_dynamics = NULL, timesteps = 20, replicates = 3, verbose = FALSE) ## ---- message = FALSE, fig.align = "center"----------------------------------- plot_pop_spatial(egk_results, timesteps = c(1, 10, 20)) ## ---- message = FALSE, progress = FALSE, fig.align = "center", eval = FALSE---- # egk_landscape <- landscape(population = egk_pop, # suitability = egk_hab, # carrying_capacity = egk_k) # # egk_pop_dynamics <- population_dynamics( # change = growth(transition_matrix = egk_mat, # global_stochasticity = egk_mat_stoch), # dispersal = kernel_dispersal( # arrival_probability = "suitability", # max_distance = 5000, # dispersal_kernel = exponential_dispersal_kernel(distance_decay = 5000), # dispersal_proportion = density_dependence_dispersing() # ), # modification = NULL, # density_dependence = ceiling_density(stages = 3)) # # egk_results <- simulation(landscape = egk_landscape, # population_dynamics = egk_pop_dynamics, # habitat_dynamics = NULL, # timesteps = 20, # replicates = 3, # verbose = FALSE) ## ---- message = FALSE, progress = FALSE, fig.align = "center", eval = FALSE---- # egk_landscape <- landscape(population = egk_pop, # suitability = egk_hab, # carrying_capacity = egk_k) # # egk_pop_dynamics <- population_dynamics( # change = growth(transition_matrix = egk_mat, # global_stochasticity = egk_mat_stoch), # dispersal = kernel_dispersal(arrival_probability = "suitability", # max_distance = 5000, # dispersal_kernel = function (r) exp(-r / 2000), # dispersal_proportion = density_dependence_dispersing()), # modification = NULL, # density_dependence = ceiling_density(stages = 3)) # # egk_results <- simulation(landscape = egk_landscape, # population_dynamics = egk_pop_dynamics, # habitat_dynamics = NULL, # timesteps = 20, # replicates = 3, # verbose = FALSE) ## ---- message = FALSE, progress = FALSE, fig.align = "center", eval = FALSE---- # egk_landscape <- landscape(population = egk_pop, # suitability = egk_hab, # carrying_capacity = egk_k) # # egk_pop_dynamics <- population_dynamics( # change = growth(transition_matrix = egk_mat, # global_stochasticity = egk_mat_stoch), # dispersal = cellular_automata_dispersal( # max_cells = c(0, 10, 20) # ), # modification = NULL, # density_dependence = ceiling_density(stages = 3)) # # egk_results <- simulation(landscape = egk_landscape, # population_dynamics = egk_pop_dynamics, # habitat_dynamics = NULL, # timesteps = 20, # replicates = 3, # verbose = FALSE) ## ---- message = FALSE, results='hide'----------------------------------------- egk_landscape <- landscape(population = egk_pop, suitability = egk_hab, carrying_capacity = egk_k, "fires" = egk_fire) egk_pop_dynamics <- population_dynamics( change = growth(transition_matrix = egk_mat, global_stochasticity = egk_mat_stoch), dispersal = kernel_dispersal( arrival_probability = "suitability", max_distance = 5000, dispersal_kernel = exponential_dispersal_kernel(distance_decay = 5000) ), modification = NULL, density_dependence = ceiling_density(stages = 3)) egk_results <- simulation(landscape = egk_landscape, population_dynamics = egk_pop_dynamics, habitat_dynamics = list( fire_effects(fire_layers = "fires", effect_time = 5, regeneration_function = function (time) {-time}) ), timesteps = 20, replicates = 3, verbose = FALSE) plot_pop_trend(egk_results) ## ---- message = FALSE, fig.align = "center"----------------------------------- plot_pop_spatial(egk_results, stage = 3, timesteps = c(1, 10, 20)) ## ---- message = FALSE, fig.align = "center"----------------------------------- plot_hab_spatial(egk_results, timesteps = c(1, 10, 20)) ## ---- message = FALSE, results='hide', eval = FALSE--------------------------- # carrying_cap_fun <- function (landscape, timestep) { # # fun <- function(suitability) { # 75 - round(75 * dlogis(suitability, scale = 0.25)) # } # # suit <- landscape$suitability # if (raster::nlayers(suit) > 1) { # suit <- suit[[timestep]] # } # # raster::calc(suit, fun) # # } # # suit_seq <- seq(0, 1, 0.1) # plot(suit_seq, 75 - round(75 * dlogis(suit_seq, scale = 0.25)), type = 'l') # # # egk_landscape <- landscape(population = egk_pop, # suitability = egk_hab, # carrying_capacity = carrying_cap_fun, # fires = egk_fire) # # egk_pop_dynamics <- population_dynamics( # change = growth(transition_matrix = egk_mat, # global_stochasticity = egk_mat_stoch), # dispersal = kernel_dispersal( # arrival_probability = "suitability", # max_distance = 5000, # dispersal_kernel = exponential_dispersal_kernel(distance_decay = 5000) # ), # modification = NULL, # density_dependence = ceiling_density(stages = 3)) # # egk_results <- simulation(landscape = egk_landscape, # population_dynamics = egk_pop_dynamics, # habitat_dynamics = list(fire_effects(fire_layers = "fires", # effect_time = 5, # regeneration_function = # function (time) {-time})), # timesteps = 20, # replicates = 3, # verbose = FALSE) ## ---- message = FALSE, results = 'hide', eval = FALSE------------------------- # egk_landscape <- landscape(population = egk_pop, # suitability = egk_hab, # carrying_capacity = carrying_cap_fun, # "fires" = egk_fire, # "roads" = egk_road) # # egk_pop_dynamics <- population_dynamics( # change = growth(transition_matrix = egk_mat, # global_stochasticity = egk_mat_stoch), # dispersal = kernel_dispersal( # arrival_probability = "suitability", # max_distance = 5000, # dispersal_kernel = exponential_dispersal_kernel(distance_decay = 5000) # ), # modification = NULL, # density_dependence = ceiling_density(stages = 3)) # # egk_results <- simulation(landscape = egk_landscape, # population_dynamics = egk_pop_dynamics, # habitat_dynamics = list( # fire_effects(fire_layers = "fires", # effect_time = 5, # regeneration_function = function (time) {-time}), # disturbance(disturbance_layers = "roads", # effect_time = 1) # ), # timesteps = 20, # replicates = 3, # verbose = FALSE) ## ---- message = FALSE, fig.align = "center"----------------------------------- egk_results[[2]][[5]][[1]][[1]] ## ---- message = FALSE, fig.align = "center"----------------------------------- egk_results[[2]][[3]][[2]] ## ---- message = FALSE, fig.align = "center"----------------------------------- egk_results[[2]][[3]][["suitability"]] ## ---- message = FALSE, fig.align = "center"----------------------------------- extract_spatial(egk_results, replicate = 2, timestep = 3, landscape_object = "suitability") ## ---- message = FALSE, fig.align = "center"----------------------------------- par(mar = c(0.1, 0.1, 0.1, 0.1)) plot(extract_spatial(egk_results), box = FALSE, axes = FALSE, col = viridis(100))