Spatially and temporally explicit PVA for Lowland Tapir in Magdalena
Diego J. Lizcano, German Forero, et al.
2021-03-02
Load packages
library(raster)
library(viridisLite)
library(future)
library(foreach)
library(xtable)
library(sf)
library(tmap)
library(tmaptools)
library(tidyverse)
library(popbio)
# library(velox)
library(steps)
library(gridExtra)
library(rasterVis)
###############################
# resolve namespace conflicts
###############################
select <- dplyr::select
map <- purrr::map
projection <- raster::projection
landscape <- steps::landscape
population_dynamics <- steps::population_dynamics
simulation <- steps::simulation
####################################
## Function to random harvest
####################################
ceiling_density_harvest_random <- function (harvest, stages = NULL) {
pop_dynamics <- function(landscape, timestep) {
population_raster <- landscape$population
idx <- which(!is.na(raster::getValues(population_raster[[1]])))
cc <- landscape$carrying_capacity
if (is.null(cc)) {
stop("carrying capacity must be specified in the landscape object to use ceiling_density",
call. = FALSE)
}
population_matrix <- raster::extract(population_raster, idx)
carrying_capacity <- raster::extract(cc, idx)
if (is.null(stages)) {
stages <- seq_len(ncol(population_matrix))
}
overpopulation <- as.vector(carrying_capacity)/rowSums(population_matrix[, stages, drop = FALSE])
overpopulation[is.nan(overpopulation)] <- 0
overpopulation <- pmin(overpopulation, 1)
population_matrix[, stages] <- sweep(population_matrix[, stages, drop = FALSE], 1, overpopulation, "*")
# population_matrix <- round_pop(population_matrix)
population_raster[idx] <- population_matrix
# landscape$population <- population_raster
# harvest - hunt
ncells <- length(population_raster[idx])
population_raster[sample(seq_len(ncells), harvest, replace =FALSE)] <- 0
landscape$population <- population_raster
landscape
}
#result <- as.population_density_dependence(pop_dynamics)
#attr(result, "density_dependence_stages") <- stages
#return(result)
}
round_pop <- function (population) {
population_min <- floor(population)
if (steps_stash$demo_stochasticity == "full") {
if (sum(population) == 0) return(population)
return(rmultinom_large_int(population)[, 1])
}
if (steps_stash$demo_stochasticity == "none") return(population)
}
steps_stash <- new.env()
flush_stash <- function() {
for (name in names(steps_stash)) {
steps_stash[[name]] <- NULL
}
}
# replace the values in the steps stash with those in this new one (used to pass the stash onto parallel workers)
replace_stash <- function(new_stash) {
# flush the old one
flush_stash()
for (name in names(new_stash)) {
steps_stash[[name]] <- new_stash[[name]]
}
}1st step: Lefkovitch matrix
First we setup a stage-based transition matrix, also known as a Lefkovitch matrix.
Our matrix is based in a simple model of three stages, were only adults reproduce.
* Where:
- S is survival (supervivencia).
- P is probability of remain on this stage. Probabilidad de permanecer en esta etapa en 1 año.
- F is contribution to next generation. Una hembra (0.5) produce 1 tapir cada dos años (0.5/2).
t_mat_ind <- matrix(c("P0","0" ,"F3", # juveniles y sub adultos sin cria. (1/2 h con cria)
"S0","P1" ,"0", # * 100 adultos hay 0.25 juveniles
"0","S1" ,"P2"), # 80 % adultos sobrevive cada año
nrow =3,
ncol =3,
byrow =TRUE)
xtable (t_mat_ind, caption = "Lefkovitch matrix")tapir_mat <- matrix(c(0.00,0.001,0.25, # juveniles y sub adultos sin cria. (1/2 h con cria)
0.50,0.2 ,0.00, # * 100 adultos hay 0.25 juveniles
0.00,0.75 ,0.90), # 90 % adultos sobrevive cada año
nrow =3,
ncol =3,
byrow =TRUE)
colnames(tapir_mat) <- rownames(tapir_mat) <- c('juvenile','subadult','adult')
tapir_mat_stoch <- matrix(c(0.00 ,0.05 ,0.05,
0.05, 0.05 ,0.00,
0.00 ,0.05 ,0.03),
nrow =3,
ncol =3,
byrow =TRUE)
colnames(tapir_mat_stoch) <- rownames(tapir_mat_stoch) <- c('juvenile','subadult','adult')
xtable (tapir_mat, caption = "transition matrix")The first matrix represents transition probabilities for each of the life stages in our tapir model. In this matrix, fecundities - or numbers of offspring per surviving and contributing stage - are in the first row. These are always represented as positive numbers. Survival probabilities are shown in the remaining rows and are expressed as positive numbers between zero and one.
Rationale:
Todos los juveniles pasan a sub-adulto ningun juvenil se queda alli (0), ningun juvenil tiene cria (0.0). la mitad de las hembras tienen cria (0.25)
La mortalidad de los juveniles es 25% (75% sobrevive (0.8)). 20% de los sub adultos se queda en la categoria un año mas (0.2),
cero, 75% de los sub adultos pasan a adultos, 90 % adultos sobrevive cada año.
The second matrix describes the uncertainty around these transition probabilities (as standard deviations) and is used to simulate environmental stochasticity - zeros indicating no stochasticity, i.e. there is no uncertainty about environmental fluctuations in our system. Note, there are three life-stages and the matrices are symmetrical in the number of columns and rows (three columns, three rows).
## [1] 1.1725
## F P G
## 0.7986654 0.7311043 0.1792457Fertility F in last column, stasis P on diagonal, and growth G in bottom-left triangle.
Read Spatial Data. WCS “Paisaje” and vegetation cover
Next, we read in spatial inputs to be used for the simulations - steps operates primarily on raster data. A habitat map based on MODIS vegetation cover with a pixel resolution of 1Km x 1Km is used as the landscape for the metapopulation of tapirs. This is considered a grid-based metapopulation structure because each cell represents a place where tapirs can move to and from - dependent upon its unique attributes. The habitat suitability layer describes the relative likelihood of the environment to support organisms (the tapir in this example) in each cell and contains values between 0 (inhabitable) and 1 (habitable).
We are using the limits of Sabanas Orientales de Arauca y Vichada provided by Carlos Rios (WCS) and a vegetation cover ma derived from MODIS. Specifically we are using the product MCD12Q1 form year 2017 (the same year the occupancy data was taken) with a pixel of 250 x 250 m. The vegetation cover was rescaled to 3000 x 3000 m each pixel counting the percent of forest and woody vegetation.
As reference we can see the potential sampling grid (1Km x 1Km) and Bita River un blue.
##### rad SHP ####
rio.magna <- st_read("D:/BoxFiles/Box Sync/CodigoR/WCS/shp/rio_epsg3116_magna_bog.shp")## Reading layer `rio_epsg3116_magna_bog' from data source `D:\BoxFiles\Box Sync\CodigoR\WCS\shp\rio_epsg3116_magna_bog.shp' using driver `ESRI Shapefile'
## Simple feature collection with 22 features and 17 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 1618438 ymin: 1127838 xmax: 1718072 ymax: 1182112
## projected CRS: MAGNA_SIRGAS_Colombia_Bogota_zone## Reading layer `CeldasMagdalena' from data source `D:\BoxFiles\Box Sync\CodigoR\WCS\shp\CeldasMagdalena.shp' using driver `ESRI Shapefile'
## Simple feature collection with 256 features and 4 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 963723.4 ymin: 1220595 xmax: 999723.4 ymax: 1247595
## projected CRS: MAGNA_SIRGAS_Colombia_Bogota_zone## Reading layer `CLC2017Llanos' from data source `D:\BoxFiles\Box Sync\CodigoR\WCS\shp\CLC2017Llanos.shp' using driver `ESRI Shapefile'
## Simple feature collection with 3092 features and 18 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 1615707 ymin: 1111484 xmax: 1737188 ymax: 1191307
## projected CRS: MAGNA-SIRGAS / Colombia Bogota zone## Reading layer `vPg_Paisajes_WCSCol_202009_MB' from data source `D:\BoxFiles\Box Sync\CodigoR\WCS\shp\vPg_Paisajes_WCSCol_202009_MB.shp' using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 6 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 613515.2 ymin: 464590.2 xmax: 1739856 ymax: 1470398
## projected CRS: MAGNA-SIRGAS / Colombia Bogota zonerio.utm <- st_transform(rio.magna, "+proj=utm +zone=18 +ellps=intl +towgs84=307,304,-318,0,0,0,0 +units=m +no_defs")
celdas.utm <- st_transform(celdas.magna, "+proj=utm +zone=18 +ellps=intl +towgs84=307,304,-318,0,0,0,0 +units=m +no_defs")
veg_cover.utm <- st_transform(cov.magna, "+proj=utm +zone=18 +ellps=intl +towgs84=307,304,-318,0,0,0,0 +units=m +no_defs")
paisajes.utm <- st_transform(paisajes.magna, "+proj=utm +zone=18 +ellps=intl +towgs84=307,304,-318,0,0,0,0 +units=m +no_defs")
magdalena <- paisajes.utm %>% filter (Landscape =="Middle Magdalena Valley")
hab <- raster("D:/BoxFiles/Box Sync/CodigoR/WCS/data/modis/modis_MODIS_MCD12Q1_2017.tif")
hab.utm <- projectRaster(hab,
crs = "+proj=utm +zone=18 +ellps=intl +towgs84=307,304,-318,0,0,0,0 +units=m +no_defs", method='ngb', res = 250)
# load burned areas 2017
burned_area <- list.files("D:/BoxFiles/Box Sync/CodigoR/WCS/data/modis/burned2017",
full.names = TRUE) %>% stack() %>% sum()
burned_area.utm <- projectRaster(burned_area,
crs = "+proj=utm +zone=18 +ellps=intl +towgs84=307,304,-318,0,0,0,0 +units=m +no_defs", method='ngb', res = 500)
lc_names <- tibble(landcover = 0:15,
lc_name = c("00_water", # 1
"01_evergreen_needleleaf", # 2
"02_evergreen_broadleaf", # 3
"03_deciduous_needleleaf", # 4
"04_deciduous_broadleaf", # 5
"05_mixed_forest", # 6
"06_closed_shrubland", # 7
"07_open_shrubland", # 8
"08_woody_savanna", # 9
"09_savanna", # 10
"10_grassland", # 11
"11_wetland", # 12
"12_cropland",
"13_urban",
"14_mosiac",
"15_barren"))
hab.mask <- mask(hab.utm, magdalena)
# vals<-unique(hab.mask)
# recl<-matrix(c(vals, c(0, rep(1, 6), 9, 1,1, 13)),ncol=2)
# recl
# hab.utm
# par(mar=c(0,0,0,0),oma=c(0,0,0,0))
# plot(hab.utm,box =FALSE,axes =FALSE,col = viridis(100),main ="Habitat Suitability")
##################################################
# function to count the % for each vegetation type
##################################################
# the factor count in 2, 4, 6, 8 celdas vecinas
# stored as a list
cov_pct <- lapply(unique(hab.mask), function(land_class) {
aggregate(hab.mask, fact=13, fun=function(vals, na.rm) {
sum(vals==land_class, na.rm=na.rm)/length(vals)
})
})
#################################################
#########################
# Suitability index
suitable <- cov_pct[[1]] + (cov_pct[[2]]*0.9) + cov_pct[[3]] + cov_pct[[4]] # + (0.5*cov_pct[[6]]) #+(0.1*cov_pct[[7]])
#########################
tapir_hab <- suitable
names(tapir_hab) <- "tapir_habitat_3k"
legend_title1 = expression("Modis Vegetation Cover (MCD12Q1_2017)")
legend_title2 = expression("Forest percent in 3x3 Km")
tm_shape(hab.mask) +
tm_raster(palette = "YlGn", title = legend_title1) +
tm_shape(celdas.utm) +
tm_polygons(alpha = 0.5) +
tm_shape(magdalena) +
tm_polygons(alpha = 0.25) +
tm_shape(rio.utm) +
tm_polygons (border.col = "blue", border.alpha = 0.5) 
tm_shape(tapir_hab) +
tm_raster(palette = "YlGn", title = legend_title2) + tm_layout(bg.color = "white") +
tm_shape(celdas.utm) +
tm_polygons(alpha = 0.5) +
tm_shape(magdalena) +
tm_polygons(alpha = 0.25) +
tm_shape(rio.utm) +
tm_polygons (border.col = "blue", border.alpha = 0.5) 
Carrying Capacity Map
The carrying capacity layer describes the total number of organisms (here, tapirs) that may occur in each cell (3000 x 3000 m) and contains either zeros or positive integer values (not a fractional number). The Tapir carrying capacity is directly based on the habitat suitability multiplied by 0.4, which is the density per 1 sqr Km according to:
Medici PhD thesis 2010.
That paper found a density for Tapirus terrestris of: 0.8±0.2 lowland tapirs/km2 and population size of 200±33 individuals.
tapir_k <- floor(suitable * (0.4 * 9) ) #3000 x 3000 = 9 pixeles ceiling
# par(mar=c(0,0,0,0),oma=c(0,0,0,0))
# plot(tapir_k,box =FALSE,axes =FALSE,col = viridis(100))
tm_shape(tapir_k) +
tm_raster(palette = "YlGnBu", title = "Maximum number of tapirs") +
tm_shape(celdas.utm) +
tm_polygons(alpha = 0.5) +
tm_shape(magdalena) +
tm_polygons(alpha = 0.25) 
Potential number of tapirs acording to carrying capacity map
maxnumber <- as.data.frame(tapir_k) %>%
group_by(layer) %>%
tally() %>%
mutate(area = n * res(tapir_k)[1] * res(tapir_k)[2])
maxnumber$tapirs <- maxnumber$layer * maxnumber$n
xtable (maxnumber, caption = "count of pixels values and tapirs")The maximum number of tapirs is: 4258.
Initial Population Size
Populations are represented as a stack of rasters that describes the total number of individuals that occur in each cell for each life-stage at the beginning of the simulations (initial populations). In the tapir example, there are three life-stages and thus three individual raster layers must be in the stack. The values are either zeros or positive integers.
For this example, we have not used equal numbers of tapirs in each cell for each life-stage. Instead we derived the initial population abundances by using stable state age distributions multiplied by the carrying capacity of the landscape, and then drawn from a multinomial distribution to return whole integers:
If we carried out field work we can use those numbers (much better). Usually the popular hunting models derived from Robinson and Redford assumed the estimated density is equivalent to carrying capacity.
In several Tapir Vortex studies the initial population size for all scenarios was half of their carrying capacity. https://www.researchgate.net/publication/254214001_How_many_lowland_tapirs_Tapirus_terrestris_are_needed_in_Atlantic_Forest_fragments_to_ensure_long-term_persistence Also: Five hundred iterations were run for each scenario and a time frame of 100 years was used (Medici et al. 2007).
A population is considered demographically viable when it had a ≤ 5% probability of extinction during a 100-year period (Soulé 1987).
#initial population abundances
stable_states <- abs( eigen(tapir_mat)$vectors[,1] / sum(eigen(tapir_mat)$vectors[,1]))
popN <- stack(replicate(ncol(tapir_mat), tapir_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(tapir_mat)
pop## class : RasterStack
## dimensions : 97, 56, 5432, 3 (nrow, ncol, ncell, nlayers)
## resolution : 3250, 3250 (x, y)
## extent : 502504.3, 684504.3, 665862, 981112 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +ellps=intl +units=m +no_defs
## names : juvenile, subadult, adult
## min values : 0, 0, 0
## max values : 1, 1, 3#####
tapir_pop <- pop
#names(tapir_pop) <- colnames(tapir_mat)
par(mar=c(0,0,0,0), oma=c(0,0,0,0))
spplot(tapir_pop, col.regions =c("gray", viridis(100)))
Initial stage abundances.
Building a model
Landscape object
All three of these spatial components (habitat suitability, carrying capacity, and initial stage abundances) are combined to define a “landscape” object, however, only the population abundance stack (tapir_pop) is required to complete a simulation. Note, all input rasters in the simulations need to have the same projection, extent, resolution, and positions of null (NA) values. The landscape object is modified at each timestep in a single simulation based on population and habitat dynamics.
Population Dynamics
Above are the only landscape data input requirements for a simple population simulation, however, we also need to specify population dynamics. Dynamic functions are used to modify populations or habitats in a landscape at each timestep in a simulation. They can be selected as ‘off-the-shelf’ functions - i.e. those included in the steps package or custom defined functions created by the user. Four types of population dynamics can be specified: change in population sizes resulting from transition matrices, population dispersal around the landscape, modifications to the population size and configuration, and density dependence. In its most basic form, only population change will be applied to the landscape. As shown in the following code, we use the growth() function, which uses the probabilities in the transition matrix to change the population, and leave all of the other parameters at their default values (NULL). This is analogous to specifying a simple life-stage population growth model.
Simulations
Now that we have constructed the landscape object and defined a dynamic object, we can run a single simulation (i.e replicates = 1 which is the default). We simulate changes to the tapir population over twenty timesteps. Runtime will depend on the complexity of the landscape object and the configuration of the dynamic object(s). The “verbose = FALSE” suppresses information printed to the console during a simulation; if omitted or set to TRUE (default) the simulation progress bars are shown.
tapir_results <- steps::simulation(landscape = tapir_landscape,
population_dynamics = tapir_pop_dynamics,
habitat_dynamics =NULL,
timesteps =100, # 100 años
replicates =1,
verbose =FALSE)
plot(tapir_results, stage =0) # view the total tapir population trajectory 
#population distribution over the landscape for a single life-stage (only timesteps one, ten, and twenty shown)
plot(tapir_results,
type ="raster",
stage =3,
timesteps = c(1,25,50,75,100),
panels = c(5,1)) 
### Animation: animate the the population rasters of subadults for timesteps one, ten and twenty (change these to show other timesteps of the simulation):
# plot(tapir_results,
# type ="raster",
# stage =2,
# timesteps = c(1,10,20),
# animate =TRUE)
# We can also generate multiple replicates of a simulation. For the tapir, we specify 50 replicates of a hundred timestep simulation. The replicates are run sequentially by default, however, to improve computation time, we have included the ability to run all three replicates in parallel - each on a different processor. Type ?plan for more information on how to operate in parallel. Note plan() is only required to be called once at the beginning of a script.
Many replicates in parallell
#########################
# define the replicates
replicas <- 25
#########################
# plan(multisession, workers =2) # This is how we specify to simulate
# replicates on separate processors in parallel
tapir_results1 <- steps::simulation(landscape = tapir_landscape,
population_dynamics = tapir_pop_dynamics,
habitat_dynamics =NULL,
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results1)
Stocasticity
Up to this point, the projected population changes in a multiple replicate simulation will all be similar to the population projection in the single replicate simulation. This is because we have not added any stochastic dynamics to the landscape. However, demographic stochasticity is on by default because it’s incorporated into the transition probabilities that the growth function uses - this behaviour can be disabled by setting ‘demo_stochasticity = “none”’ in the simulation parameters. Let’s try adding some globally-acting (applies to all cells in the landscape) environmental stochasticity to the growth function:
tapir_landscape1 <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity =NULL)
tapir_pop_dynamics1 <- population_dynamics(
change = growth(transition_matrix =tapir_mat,
global_stochasticity =tapir_mat_stoch),
dispersal =NULL,
modification =NULL,
density_dependence =NULL)
tapir_results2 <- steps::simulation(landscape =tapir_landscape1,
population_dynamics =tapir_pop_dynamics1,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results2)
# cellStats(tapir_results2[[1]][[1]][[1]], sum)
#extinction probability using popbio
logN <- log(grizzly$N[-1]/grizzly$N[-39])
countCDFxt(mu=mean(logN), sig2=var(logN), nt=38, tq=38, Nc=99, Ne=20)
We can also specify how survival and fecundity values are influenced by spatial layers in the landscape object (e.g. habitat suitability, carrying capacity, etc.). We modify the survival and fecundity values in the transition matrix at each time step according to a spatial layer by specifying a transition function within the growth function. This is done by multiplying values in the spatial layers by the survival and fecundity values in the transition matrix. This produces new survival and fecundity values for each grid cell in the landscape. Note, if either the survival or fecundity layer are set to NULL (default), the respective values in the transition matrix will remain unchanged. A user may add any spatial layer(s) to the landscape object to modify the transition matrices; here, we specify the habitat suitability layer in the landscape for the survival only:
# plan(sequential)
tapir_landscape2 <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity =NULL)
tapir_pop_dynamics2 <- population_dynamics (
change = growth(transition_matrix = tapir_mat,
global_stochasticity = tapir_mat_stoch,
transition_function = modified_transition (survival_layer ="suitability",
fecundity_layer =NULL)),
dispersal =NULL,
modification =NULL,
density_dependence =NULL)
tapir_results3 <- simulation(landscape =tapir_landscape2,
population_dynamics =tapir_pop_dynamics2,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results3)
Base line model
Using carrying capacity on stage 3 (Adults)
The population will grow according to the transition matrix until it reaches a ceiling carrying capacity, at which point individuals will die to keep the population at the carrying capacity of the landscape. Here we specify that it is based on the number of adults only (“stages = 3” in the ceiling_density function).
tapir_landscape <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity =tapir_k)
tapir_pop_dynamics <- population_dynamics(
change = growth(transition_matrix = tapir_mat,
global_stochasticity =tapir_mat_stoch),
dispersal =NULL,
modification =NULL,
density_dependence = ceiling_density(stages =3))
tapir_results4 <- simulation (landscape =tapir_landscape,
population_dynamics =tapir_pop_dynamics,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results4)
plot(tapir_results4[10],type ="raster",stage =3,timesteps = c(1,20,30,40,50,60,70,80,90,100),panels = c(5,2))
Harvesting (Hunting Tapirs)
We simulated the hunting of 10, 25, 50, 75, 100, 250, and 500 tapirs annually. To perform the simulation we applied a mortality layer in each timestep. The mortality layer consists of values from 0 to 1 and modifies the population by multiplying the population of a cell by the value of the corresponding cell in a mortality layer. For example, a cell with ten individuals before the mortality function is applied, and corresponding mortality layer cell with a value of 0.1, would have one individuals remaining after modification.
The hunting occurs randomly.
Hunting 10 tapirs per year
# hunting <- tapir_k
# hunting[hunting == 0] <- NA
# hunting[hunting >= 1] <- 0
#
# hunt10 <- stack()
# # hunt10[sampleRandom(hunting,10,cells=TRUE, na.rm=T)[,'cell']]=0.1
# for (i in 1:100) {
# newraster <- hunting
# newraster[sampleRandom(hunting,10,cells=TRUE, na.rm=T)[,'cell']]=1
# hunt10 <- stack(hunt10, newraster)
# }
#
# plot(hunt10[[50]],col=rainbow(3), main="hunting step 50")
##############
#hunting_mortal <- mortality(mortality_layer = "hunt", stages = NULL)
tapir_landscape <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity = tapir_k
#"hunt"= hunt10,
)
tapir_pop_dynamics <- population_dynamics(
change = growth(transition_matrix = tapir_mat,
global_stochasticity =tapir_mat_stoch),
dispersal = NULL,
modification = ceiling_density_harvest_random(harvest = 10, stages = c(1,2,3)),
density_dependence = NULL
# dynamics_order = c("change", "modification", "density_dependence","dispersal")
)
tapir_results5 <- simulation (landscape =tapir_landscape,
population_dynamics =tapir_pop_dynamics,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results5)
plot(tapir_results5[10],type ="raster",stage =3,timesteps = c(1,20,30,40,50,60,70,80,90,100),panels = c(5,2))
Hunting 25 tapirs per year
tapir_landscape <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity = tapir_k
#"hunt"= hunt10,
)
tapir_pop_dynamics <- population_dynamics(
change = growth(transition_matrix = tapir_mat,
global_stochasticity =tapir_mat_stoch),
dispersal = NULL,
modification = ceiling_density_harvest_random(harvest = 25, stages = c(1,2,3)),
density_dependence = NULL
# dynamics_order = c("change", "modification", "density_dependence","dispersal")
)
tapir_results6 <- simulation (landscape =tapir_landscape,
population_dynamics =tapir_pop_dynamics,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results6)
plot(tapir_results6[10],type ="raster",stage =3,timesteps = c(1,20,30,40,50,60,70,80,90,100),panels = c(5,2))
Hunting 50 tapirs per year
tapir_landscape <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity = tapir_k
#"hunt"= hunt10,
)
tapir_pop_dynamics <- population_dynamics(
change = growth(transition_matrix = tapir_mat,
global_stochasticity =tapir_mat_stoch),
dispersal = NULL,
modification = ceiling_density_harvest_random(harvest = 50, stages = c(1,2,3)),
density_dependence = NULL
# dynamics_order = c("change", "modification", "density_dependence","dispersal")
)
tapir_results7 <- simulation (landscape =tapir_landscape,
population_dynamics =tapir_pop_dynamics,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results7)
plot(tapir_results7[10],type ="raster",stage =3,timesteps = c(1,20,30,40,50,60,70,80,90,100),panels = c(5,2))
Hunting 75 tapirs per year
tapir_landscape <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity = tapir_k
#"hunt"= hunt10,
)
tapir_pop_dynamics <- population_dynamics(
change = growth(transition_matrix = tapir_mat,
global_stochasticity =tapir_mat_stoch),
dispersal = NULL,
modification = ceiling_density_harvest_random(harvest = 75, stages = c(1,2,3)),
density_dependence = NULL
# dynamics_order = c("change", "modification", "density_dependence","dispersal")
)
tapir_results8 <- simulation (landscape =tapir_landscape,
population_dynamics =tapir_pop_dynamics,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results8)
plot(tapir_results8[10],type ="raster",stage =3,timesteps = c(1,20,30,40,50,60,70,80,90,100),panels = c(5,2))
Hunting 100 tapirs per year
tapir_landscape <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity = tapir_k
#"hunt"= hunt10,
)
tapir_pop_dynamics <- population_dynamics(
change = growth(transition_matrix = tapir_mat,
global_stochasticity =tapir_mat_stoch),
dispersal = NULL,
modification = ceiling_density_harvest_random(harvest = 100, stages = c(1,2,3)),
density_dependence = NULL
# dynamics_order = c("change", "modification", "density_dependence","dispersal")
)
tapir_results9 <- simulation (landscape =tapir_landscape,
population_dynamics =tapir_pop_dynamics,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results9)
plot(tapir_results9[10],type ="raster",stage =3,timesteps = c(1,20,30,40,50,60,70,80,90,100),panels = c(5,2))
Hunting 150 tapirs per year
tapir_landscape <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity = tapir_k
#"hunt"= hunt10,
)
tapir_pop_dynamics <- population_dynamics(
change = growth(transition_matrix = tapir_mat,
global_stochasticity =tapir_mat_stoch),
dispersal = NULL,
modification = ceiling_density_harvest_random(harvest = 150, stages = c(1,2,3)),
density_dependence = NULL
# dynamics_order = c("change", "modification", "density_dependence","dispersal")
)
tapir_results10 <- simulation (landscape =tapir_landscape,
population_dynamics =tapir_pop_dynamics,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results10)
plot(tapir_results10[10],type ="raster",stage =3,timesteps = c(1,20,30,40,50,60,70,80,90,100),panels = c(5,2))
Hunting 250 tapirs per year
tapir_landscape <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity = tapir_k
#"hunt"= hunt10,
)
tapir_pop_dynamics <- population_dynamics(
change = growth(transition_matrix = tapir_mat,
global_stochasticity =tapir_mat_stoch),
dispersal = NULL,
modification = ceiling_density_harvest_random(harvest = 250, stages = c(1,2,3)),
density_dependence = NULL
# dynamics_order = c("change", "modification", "density_dependence","dispersal")
)
tapir_results11 <- simulation (landscape =tapir_landscape,
population_dynamics =tapir_pop_dynamics,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results11)
plot(tapir_results11[10],type ="raster",stage =3,timesteps = c(1,20,30,40,50,60,70,80,90,100),panels = c(5,2))
Hunting 500 tapirs per year
tapir_landscape <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity = tapir_k
#"hunt"= hunt10,
)
tapir_pop_dynamics <- population_dynamics(
change = growth(transition_matrix = tapir_mat,
global_stochasticity =tapir_mat_stoch),
dispersal = NULL,
modification = ceiling_density_harvest_random(harvest = 500, stages = c(1,2,3)),
density_dependence = NULL
# dynamics_order = c("change", "modification", "density_dependence","dispersal")
)
tapir_results12 <- simulation (landscape =tapir_landscape,
population_dynamics =tapir_pop_dynamics,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results12)
plot(tapir_results12[10],type ="raster",stage =3,timesteps = c(1,20,30,40,50,60,70,80,90,100),panels = c(5,2))
Dispersal dinamycs
We can either specify assisted species movements (i.e. translocations or reintroductions) or natural species movements (dispersal) in the model simulations.
Translocations or reintroductions
To use this part we need to define the source and destiny as a raster (few pixeles)
To characterise a translocation, we specify two raster layers that map the locations of where we will take individuals from and where we will place individuals in the landscape. Cell values (integers) specify how many individuals will be added/removed. Origin and destination spatial layers are stored in the landscape object and referenced by name. We use a built-in translocation function to modify the population(s) in a simulation at each specified timestep. The function requires the name of origin and destination layers stored in the landscape object (“origins” and “destinations”), the life-stages to modify, and the timesteps at which to modify the population. In this example, we provide the tapir origins (“tapir_origins”) and destinations (“tapir_destinations”) layers, a life-stage of “3” (only adults moved), and timesteps one, five, ten and fifteen - meaning each translocation will occur every five years in each simulation replicate. Note, if the number of individuals is not available in the landscape when the translocation is applied, a warning will be generated and only the maximum available individuals will be used.
tapir_landscape <- landscape(population =tapir_pop,
suitability =NULL,
carrying_capacity =NULL,
"origins"=tapir_origins,
"destinations"=tapir_destinations)
tapir_pop_dynamics <- population_dynamics(
change = growth(transition_matrix =tapir_mat,
global_stochasticity =tapir_mat_stoch),
dispersal =NULL,
modification = translocation(origins_layer ="origins",
destinations_layer ="destinations",
stages =3,
effect_timesteps = c(1,5,10,15)),
density_dependence =NULL)
tapir_results <- simulation(landscape =tapir_landscape,
population_dynamics =tapir_pop_dynamics,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =20,
replicates =replicas,
verbose =FALSE)
plot(tapir_results)Species movements (dispersal)
We can also add dispersal to allow tapirs to move throughout the landscape. There are several built-in dispersal functions and most of them utilize habitat suitability, carrying capacity, or both to determine arrival probabilities see package help for more information by typing ?population_dispersal_functions at the R prompt. Thus, we will need to provide the appropriate habitat layer in the initial landscape. Below, we use kernel-based dispersal (with a maximum dispersal distance of 10000 meters (10 Km)) and habitat suitability to determine the arrival probabilities:
For individual-based movements, we also provide a cellular automata-based dispersal function. Because this is an individual-based model, a dispersal kernel is not required. Specifying a single number for the “max_cells” parameter will indicate the maximum number of cells movements (traversed around the landscape), for each dispersing individual, for all life-stages. A vector of max_cells (equal in length to the number of life-stages) can also be specified. In the example below, no juveniles (stage 1) disperse, subadults (stage 2) only disperse up to 5 cells (15 km), and adults (stage 3) disperse up to 8 cells (24 km):
tapir_landscape <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity =tapir_k)
tapir_pop_dynamics <- population_dynamics(
change = growth(transition_matrix =tapir_mat,
global_stochasticity=tapir_mat_stoch),
dispersal = cellular_automata_dispersal(
max_cells = c(0, 5, 8)), # kernel_dispersal(arrival_probability ="none",
# max_distance =15000,
# dispersal_kernel = exponential_dispersal_kernel(distance_decay =9000)),
density_dependence = NULL, #ceiling_density(stages =3))
modification = ceiling_density_harvest_random(harvest = 100, stages = c(1,2,3)) )
tapir_results15 <- simulation(landscape =tapir_landscape,
population_dynamics =tapir_pop_dynamics,
habitat_dynamics =NULL,
demo_stochasticity ="none",
timesteps =100,
replicates =replicas,
verbose =FALSE)
plot(tapir_results15)
# Let’s have a look at how the adult population is changing in the first replicate of the simulation:
plot(tapir_results15[10],type ="raster",stage =3,timesteps = c(1,20,30,40,50,60,70,80,90,100),panels = c(5,2))
Habitat stocasticity (fire, deforestation)
To characterise habitat disturbance, we can use a series of raster layers that map the locations and severity of disturbances. Each disturbance raster will be multiplied with the original habitat suitability layer and thus contain values that represent appropriate ranges - zero representing the maximum disturbance and one representing no disturbance. Note, the number of disturbance layers must match the intended number of timesteps in a single simulation. There is an existing spatial dataset of fires in the steps package (“tapir_fire”) which we store in the landscape object. We use a pre-defined fire_effects function to modify the habitat suitability in the simulation at each timestep. The function requires the name of fire layers stored in the landscape object, and an effect time which specifies the number of timesteps that each disturbance layer acts on the habitat suitability (in diminishing intensity). The effect of fires can be cumulative, so if there are two fires in the same cell in succession (i.e. in a shorter space of time than the effect time) then the habitat suitability will equal suitability multiplied by both the first and second disturbance effects. In this example, we provide the tapir fire disturbance input name (“fires”), and an effect time of five - meaning each fire layer will affect the habitat suitability (in linear decreasing intensity based on the regeneration function) for five timesteps in each simulation replicate. All functions that act on the habitat must be passed in as a list in the simulation call in the order to be executed:
Use fires to model population
# Fire (stored in the landscape object and called "fires") acts on the landscape for
# five years with an exponentially decaying intensity.
## Not run:
regen <- function (time) {-exp(time)}
plot(1:5, regen(1:5), type = "l")
fire <- fire_effects(fire_layers = "fires",
effect_time = 5,
regeneration_function = regen)
ls <- landscape(population = tapir_pop,
suitability = tapir_hab,
"fires" = tapir_fire)
pd <- population_dynamics(change = growth(tapir_mat))
sim <- simulation (landscape = ls,
population_dynamics = pd,
habitat_dynamics = list(fire),
timesteps = 20)
plot(sim)
plot(sim, object = "suitability", type = "raster", timesteps = 1:9)
plot(sim, object = "population", type = "raster", timesteps = 1:9)
plot(sim, object = "carrying_capacity", type = "raster", timesteps = 1:9)fires in a more elaborated way
let’s have a look at how the adult population is changing with fires occurring in the landscape:
tapir_landscape <- landscape(population =tapir_pop,
suitability =tapir_hab,
carrying_capacity =tapir_k,
"fires"=tapir_fire)
tapir_pop_dynamics <- population_dynamics(change =
growth(transition_matrix =tapir_mat,
global_stochasticity =tapir_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))
tapir_results <- simulation(landscape =tapir_landscape,
population_dynamics =tapir_pop_dynamics,
habitat_dynamics = list(
fire_effects(fire_layers ="fires",
effect_time =5,
regeneration_function = function (time) {-time})),
timesteps =20,
replicates =replicas,
verbose =FALSE)
plot(tapir_results)
# it is possible to view the habitat suitability throughout a single simulation
plot(tapir_results[1],object ="suitability",timesteps = c(1,10,20),panels = c(3,1))
# And now let’s have a look at how the adult population is changing with fires occurring in the landscape:
plot(tapir_results[1],type ="raster",stage =3,timesteps = c(1,10,20),panels = c(3,1))Info R session
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 14393)
##
## Matrix products: default
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rasterVis_0.49 latticeExtra_0.6-29 lattice_0.20-41
## [4] gridExtra_2.3 steps_1.2.0 popbio_2.7
## [7] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.4
## [10] purrr_0.3.4 readr_1.4.0 tidyr_1.1.2
## [13] tibble_3.0.6 ggplot2_3.3.3 tidyverse_1.3.0
## [16] tmaptools_3.1-1 tmap_3.3 sf_0.9-7
## [19] xtable_1.8-4 foreach_1.5.1 future_1.21.0
## [22] viridisLite_0.3.0 raster_3.4-5 sp_1.4-5
##
## loaded via a namespace (and not attached):
## [1] fs_1.5.0 lubridate_1.7.9.2 RColorBrewer_1.1-2 httr_1.4.2
## [5] tools_4.0.3 backports_1.2.1 bslib_0.2.4 rgdal_1.5-23
## [9] R6_2.5.0 KernSmooth_2.23-18 DBI_1.1.1 colorspace_2.0-0
## [13] withr_2.4.1 tidyselect_1.1.0 leaflet_2.0.4.1 compiler_4.0.3
## [17] leafem_0.1.3 cli_2.3.0 rvest_0.3.6 xml2_1.3.2
## [21] sass_0.3.1 scales_1.1.1 hexbin_1.28.2 classInt_0.4-3
## [25] digest_0.6.27 rmarkdown_2.7 jpeg_0.1-8.1 base64enc_0.1-3
## [29] dichromat_2.0-0 pkgconfig_2.0.3 htmltools_0.5.1.1 parallelly_1.23.0
## [33] highr_0.8 dbplyr_2.1.0 htmlwidgets_1.5.3 rlang_0.4.10
## [37] readxl_1.3.1 rstudioapi_0.13 jquerylib_0.1.3 generics_0.1.0
## [41] zoo_1.8-8 jsonlite_1.7.2 crosstalk_1.1.1 magrittr_2.0.1
## [45] Rcpp_1.0.6 munsell_0.5.0 abind_1.4-5 lifecycle_1.0.0
## [49] stringi_1.5.3 leafsync_0.1.0 yaml_2.2.1 grid_4.0.3
## [53] parallel_4.0.3 listenv_0.8.0 crayon_1.4.1 stars_0.5-1
## [57] haven_2.3.1 hms_1.0.0 knitr_1.31 pillar_1.4.7
## [61] codetools_0.2-18 reprex_1.0.0 XML_3.99-0.5 glue_1.4.2
## [65] evaluate_0.14 memuse_4.1-0 modelr_0.1.8 png_0.1-7
## [69] vctrs_0.3.6 cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
## [73] xfun_0.21 lwgeom_0.2-5 broom_0.7.5 e1071_1.7-4
## [77] class_7.3-18 iterators_1.0.13 units_0.6-7 globals_0.14.0
## [81] ellipsis_0.3.1