This is supporting information for the Methods in Ecology and Evolution article entitled, Integrated SDM database: enhancing the relevance and utility of species distribution models in conservation management, by Veronica F. Frans*, Amélie A. Augé, Jim Fyfe, Yuqian Zhang, Nathan McNally, Hendrik Edelhoff, Niko Balkenhol, & Jan O. Engler. Please contact the corresponding author (*) for any inquiries.
We used the New Zealand sea lion (Phocarctos hookeri; NZSL) Maxent model training background information data (.lambda files) from Frans et al. (2018) and predicted them on mainland New Zealand and Stewart Island (see Appendix S1 for demonstration of the model training procedure). In order to determine suitable sites for breeding colonies of a minimum of 35 females, we used the multi-state species distribution model (SDM) framework (Frans et al. 2018). Information on the original model training and the multi-state SDM framework are provided in Frans et al. (2018) and its supplementary materials. We also demonstrate it in Appendix S1 of this current publication (Frans et al. 2021).
The final products of the multi-state SDM are:
a raster map indicating the suitability (or preference) of each 25 X 25 m pixel for each of the species’ three behavioural states of inland movement during the breeding season (breeding–S1; transition–S2; dispersion–S3);
a polygon shapefile, highlighting suitable locations of a minimum area size for \(\geq\) 35 females.
The resulting polygon shapefile is also essential for creating the basic, simple, and spatial structure of the iSDMdb. Note that the iSDMdb is compatible with other SDM algorithms and/or other methods to create polygons of suitable sites.
The script presented here was done using R (version 3.6.0; R Core Team 2019) and its packages.
library("raster") # raster data
library("rgdal") # input/output; projections; reading ASCII files
library("maptools") # mapping
library("classInt") # changing colour intervals on a plot
library("gridExtra") # grid graphics
library("grid") # creating graphical objects
library("lattice") # required for grid graphics
library("RColorBrewer") # colours for graphics
library("reshape") # sort data
library("reshape2") # sort data
library("igraph") # clumping contiguous pixels
library("rgeos") # polygon and point conversions
Change raster options to maximise memory.
# view default raster options
rasterOptions()
## format : raster
## datatype : FLT4S
## overwrite : FALSE
## progress : none
## timer : FALSE
## chunksize : 1e+08
## maxmemory : 5e+09
## memfrac : 0.6
## tmpdir : G:/R/temp\RtmpEtZWtg/raster/
## tmptime : 168
## setfileext : TRUE
## tolerance : 0.1
## standardnames : TRUE
## warn depracat.: TRUE
## header : none
# modify max memory from 10gb to 30gb (this is a 32GB RAM computer)
rasterOptions(maxmemory = 3e10)
We create and use multiple folders to organise our inputs and outputs. Similar directories can be made anywhere by only changing the root directory object, dir
.
# Root directory
dir <- c('G:\\R\\NZSL_MSSDM_MCA_2019_NZ\\')
# Data directory
dat.dir <- paste0(dir,'data\\')
# Maxent.jar file directory
maxent.dir <- c(paste0(dir,'maxent'))
# Maxent cache files (training and projection environmental variables)
proj.cache <- c(paste0(dat.dir,'maxent.cache\\projection\\'))
train.cache <- c(paste0(dat.dir,'maxent.cache\\training\\'))
# Intermediate data and layers layers
dir.create(paste0(dir,'data\\intermediate'),recursive=TRUE)
int.dir <- paste0(dir,'data\\intermediate\\')
# Final layers
dir.create(paste0(dir,'layers'),recursive=TRUE)
lay.dir <- paste0(dir,'layers\\')
# Final tables
dir.create(paste0(dir,'tables'),recursive=TRUE)
tab.dir <- paste0(dir,'tables\\')
# Final figures
dir.create(paste0(dir,'figures'),recursive=TRUE)
fig.dir <- paste0(dir,'figures\\')
# colours for the multiple states
state_cols <-c('#BBBBBB','#DC050C','#882E72','#1965B0','#6195CF','#4EB265','#90C987',
'#F7CB45','#EE8026')
Load rasters and project. The default projection is set to NZDG2000/New Zealand Transverse Mercator 2000 (EPSG: 2193).
# load raster and project
get.ras <- function(r_name,
# defaults
ras_dir=dat.dir, prj=c("+init=epsg:2193")){
# read raster
r <- raster(paste0(ras_dir,"\\",r_name))
# project and return raster
crs(r) <- CRS(paste0(prj))
return(r)
}
Load points or polygons and project. The default projection is set to NZDG2000/New Zealand Transverse Mercator 2000 (EPSG: 2193).
# load point/polygon and project
get.shp <- function(p_name,
# defaults
polydir=dat.dir, prj=c("+init=epsg:2193")){
# read shapefile
pol <- readOGR(paste0(polydir,p_name))
# project and return point/polygon
pol <- spTransform(pol, CRS(paste0(prj)))
return(pol)
}
Defining a minimum area for n individuals (minimum mapping unit; MMU), where suitable habitats for each behavioural state are within x units of distance. This will be used to search for suitable areas of a minimum size, which will contain no unsuitable patches (contiguous suitable values within x). More information on MMU is provided in Frans et al. (2018).
We calculate the MMU for n individuals as:
\[MMU_S = ((a/n_1)*n_2)/r\]
where \(a\) and \(n_1\) represent the number of individuals (\(n_1\)) per area (\(a\)), multiplied by the minimum number of individuals defined by the user (\(n_2\)), and all divided by the resolution of the raster layer (\(r\)) in the same units as the initial area, \(a\) (e.g., meters, kilometers).
For the NZSL, the MMU is calculated for 35 breeding females. The resolution is 25 m, or, 625 m. These will be used in the default settings.
# calculate minimum area sizes
# input number (density) of individuals per given area (area_size)
mmu.calc <- function(density, area_size,
# defaults
min_indiv=35, res=625){
# calculate MMU
mmu <- ((area_size/density)*min_indiv)/res
if (mmu <= 1){
return(1)
} else {
# round to the nearest number of pixels
return(round(mmu, digits=0))
}}
We projected the training results to the entire New Zealand mainland by accessing the ‘density.project’ command in the Maxent.jar file for Maxent v.3.3.3k (Phillips & Dudik 2008). The model was trained for 100 runs for each state, so there was a total of 300 corresponding predictions. The training run data were originally stored in folders for each state and each run under the following pattern: ~\m_output\(state number)\(run number)\
(see Frans et al. 2018 and Appendix S1 of this current publication for more details). We then calculated the mean of these predictions for each state to get three raster maps.
First, we make the predictions on the mainland.
# predictions and mean
for(state in 1:3){
# output directory
maxent.out <- paste0(dir,'m_output\\',state)
# prediction
for(run in 1:100){
# Directory for each iteration's output
runs.out <- paste0(maxent.out,"\\",run,"\\")
# Maxent.jar command line
{
command <- paste0("java -mx1024m -cp ",maxent.dir,
"\\maxent.jar density.Project ",runs.out,
"*.lambdas ",proj.cache," ",runs.out,
"projection_S",state,"_",run,
".asc -t cliff_edges -t landcover pictures=true")
system(sprintf(command))
}
}
# mean prediction per state
# empty vector
maxent_rasters <- list()
# For-loop
for (run in 1:100){
# Directory
runs.out <- paste0(maxent.out,"\\",run,"\\")
# Search for ".asc" files and stack
filenames <- list.files(path=paste0(runs.out),
pattern=paste0("projection_S",state,"_",
run,".asc"),full.names=TRUE)
maxent_rasters[[run]] <- filenames
}
# stack
maxent_100_rasters <- do.call(stack, maxent_rasters)
# calculate mean
mean_proj <- calc(x=maxent_100_rasters, fun=mean)
# save
writeRaster(mean_proj,paste0(dat.dir,"\\S",state,"_projection.tif"),
options="COMPRESS=LZW",overwrite=TRUE)
}
We then plot the results. These three mean maps will be used in the multi-state SDM framework.
# Empty vector
s123 <- list()
# For-loop (S1-3):
for(state in 1:3){
raster_names <- list.files(path=paste0(dat.dir),
pattern=paste0("\\S",state,"_projection.tif"),
full.names=TRUE)
# Add to list
s123[[state]] <- raster_names
}
# stack
s123_sdm <- do.call(stack,s123)
# Set up plots (S1-3)
png(width = 2400, height = 1800,
paste0(fig.dir,"/SDMs_S1-3.png"),res=300)
sdmplot <- spplot(s123_sdm,
main=c("SDM results for each state"),
layout=c(3,1), #(ncols, nrows)
maxpixels=500000,
col.regions=c(rev(brewer.pal(11,"Spectral"))),
cuts=10,
colorkey=list(labels=list(at=seq(0,1,by=0.1)), #set up labels
space="right"),
par.settings = list(panel.background=list(col="black"))
)
sdmplot
dev.off()
Plot here.
# Display figure
sdmplot
Read table of threshold values from previous Maxent training results (also obtained from Frans et al. 2018 training step, and shown in Appendix S1).
# read CSV summary table
thresholds <- read.csv(paste0(dat.dir,"thresholds.csv"))
# Pull table to evaluate threshold values
thresholds_table <- acast(thresholds, thresholds[,3]~thresholds[,1])
# show table
knitr::kable(thresholds_table)
S1 | S2 | S3 | |
---|---|---|---|
Balance.training.omission..predicted.area.and.threshold.value | 0.0022 | 0.0130 | 0.0072 |
Equal.test.sensitivity.and.specificity | 0.5990 | 0.2320 | 0.1944 |
Equal.training.sensitivity.and.specificity | 0.1285 | 0.2146 | 0.1561 |
Equate.entropy.of.thresholded.and.original.distributions | 0.0514 | 0.1056 | 0.0984 |
Fixed.cumulative.value.1 | 0.0003 | 0.0216 | 0.0044 |
Fixed.cumulative.value.10 | 0.3395 | 0.4031 | 0.2303 |
Fixed.cumulative.value.5 | 0.0509 | 0.1908 | 0.0852 |
Maximum.test.sensitivity.plus.specificity | 0.5990 | 0.2085 | 0.1944 |
Maximum.training.sensitivity.plus.specificity | 0.1285 | 0.1641 | 0.1240 |
Minimum.training.presence | 0.1285 | 0.1641 | 0.1240 |
X10.percentile.training.presence | 0.4350 | 0.4380 | 0.2999 |
Use the Maximum Training Sensitivity Plus Specificity (MaxSSS) threshold to create a binary dataset of presence and absence, (or, suitability and unsuitability).
# Cast to make a new table, sorted by thresholds and states
cast_thresh <- cast(thresholds, state~threshold)
# Subset MaxSSS threshold
sub_thresh <- cast_thresh[,c(grep("Maximum.training",names(cast_thresh)))]
sub_thresh <- data.frame(sub_thresh)
sub_thresh
Reclassify the probability of presence for each state map to its respective threshold. 0 is for any values below the threshold (i.e. absence) and recl_val is for values greater than or equal to the threshold (i.e. presence). recl_val is a unique value for each state that increases by one digit per subsequent state in the for-loop: S1 = 1, S2 = 10, S3 = 100.
# For-loop to reclassify by threshold
for (state in 1:3){
# Extract threshold value
recl_thresh <- sub_thresh$sub_thresh[state]
# Pull up state predictions
s_proj <- raster(paste0(dat.dir,"S",state,"_projection.tif"))
# Set unique value for each state
recl_val <- ((10^state)/10)
# Reclassify to 0 and 1
s_proj_recl <- calc(s_proj,
fun=function(x){x[x>=recl_thresh]<-recl_val;
return(x)})
s_proj_recl <- calc(s_proj_recl,
fun=function(x){x[x<recl_thresh]<-0;
return(x)})
# Save
writeRaster(s_proj_recl,
paste0(int.dir,"MaxSSS_S",state,".tif"),
options="COMPRESS=LZW",overwrite=TRUE)
}
# Read new files and stack
# Empty vector
s123_rasters <- list()
# For-loop (S1-3):
for(state in 1:3){
raster_names <- list.files(path=paste0(int.dir),
pattern=paste0("\\MaxSSS_S",state,".tif"),
full.names=TRUE)
# Add to list
s123_rasters[[state]] <- raster_names
}
# Threshold raster stack
thresh_123_stack <- do.call(stack,s123_rasters)
# Ensure projection
crs(thresh_123_stack) <- CRS("+init=epsg:2193")
View intermediate results.
# Set up plots (S1-3)
thresh_S1 <- spplot(thresh_123_stack$MaxSSS_S1,
main=c("S1 sites"),
width = 1, maxpixels=50000000,
col.regions=colorRampPalette(c('#EEDD88','#CC3311')),
cuts=1,
colorkey=list(space="bottom",tick.number=1,
labels=list(labels=c("0","1"),cex=1)),
par.settings = list(panel.background=list(col="black"))
)
thresh_S2 <- spplot(thresh_123_stack$MaxSSS_S2,
main=c("S2 sites"),
width = 1, maxpixels=50000000,
col.regions=colorRampPalette(c('#DDCC77','#CC3311')),
cuts=1,
colorkey=list(space="bottom",tick.number=1,
labels=list(labels=c("0","10"),cex=1)),
par.settings = list(panel.background=list(col="black"))
)
thresh_S3 <- spplot(thresh_123_stack$MaxSSS_S3,
main=c("S3 sites"),
width = 1, maxpixels=50000000,
col.regions=colorRampPalette(c('#DDCC77','#CC3311')),
cuts=1,
colorkey=list(space="bottom",tick.number=4,
labels=list(labels=c("0","","100"),cex=1)),
par.settings = list(panel.background=list(col="black"))
)
# Save multiplot image to folder
png(width=2400,height=1500, paste0(fig.dir,"thresholds_MaxSSS_S1-3.png"),res=300)
print(thresh_S1, position = c(.012,0,.35,1),more = T)
print(thresh_S2, position = c(.312,0,.65,1),more = T)
print(thresh_S3, position = c(.612,0,.95,1))
dev.off()
Plot here.
# Display figure here
print(thresh_S1, position = c(.018,0,.35,1),more = T)
print(thresh_S2, position = c(.318,0,.65,1),more = T)
print(thresh_S3, position = c(.618,0,.95,1))
We then generate a map of overall suitability across states by adding the threshold layers together, indicating suitability for one, two, or all three behavioural states. Pixel values will be as follows (from Frans et al. 2018):
S1 | S2 | S3 | Sum | State Suitability |
---|---|---|---|---|
0 | 0 | 0 | 0 | None |
1 | 0 | 0 | 1 | S1 only |
0 | 10 | 0 | 10 | S2 only |
0 | 0 | 100 | 100 | S3 only |
1 | 10 | 0 | 11 | S1 and S2 |
1 | 0 | 100 | 101 | S1 and S3 |
0 | 10 | 100 | 110 | S2 and S3 |
1 | 10 | 100 | 111 | All states |
# Calculate sum
sum_s123 <- calc(x=thresh_123_stack, fun=sum)
# Save to file (this is a final layer)
writeRaster(sum_s123,paste0(lay.dir,"suitablility_all_states.tif"),
options="COMPRESS=LZW",overwrite=TRUE)
View image.
# read raster
sum_s123 <- get.ras("suitablility_all_states.tif", ras_dir=lay.dir)
# Use "cut()" to make 7 classifications (warning: this turns 0 into NA)
sum_s123_cut <- cut(sum_s123, c(0,1,10,11,100,101,110,111))
# Reclassify all values to 0 to include 0s in the calculation
sum_s123_0s <- calc(sum_s123, fun=function(x){x[x>=0] <-0;return(x)})
# Unify rasters to create classified layer that includes values of 0
sum_s123_stack <- stack(sum_s123_cut,sum_s123_0s)
sum_s123_cut <- max(sum_s123_stack,na.rm=TRUE)
# Set up spplot
sum_s123_plot <- spplot(sum_s123_cut,
main=c("Multi-state SDM: suitability for all states"),
width = 1, maxpixels=5000000,
cuts=7,
col.regions=c(state_cols),
colorkey=list(space="bottom",
breaks=list(0,1,2,3,4,5,6,7),
labels=list(labels=
c("0","1","10","11","100","101","110","111"),
cex = 1)),
par.settings = list(panel.background=list(col="black"))
)
# Save image
png(width=3600,height=3600, paste0(fig.dir,"suitablility_all_states.png"),res=600)
sum_s123_plot
dev.off()
Plot here.
# Display figure
sum_s123_plot
This is the final output of the multi-state SDM framework. The suitability across states shows locations suitable for one, two or all three behavioural states, according to the MaxSSS threshold calculated from the SDMs.
We next define a minimum area for 35 individuals, using the MMU function to calculate the minimum number of contiguous suitable pixels to use to extract these areas.
Note that for this species, using contiguous suitable pixels to define suitable sites was the best method. However, other methods for creating polygons can be used. For example, a moving window analysis could be used for species that can have non-contiguous suitable areas found within a reachable distance (see range module example in Frans et al. 2018).
For the NZSL, the MMU is calculated for 35 breeding females. At S1, female NZSLs are at a density of 85 females per 100 m. At S2, where densities are at 30 females per 100 m. S3 has densities of 1 female per 10,000 m (Augé et al. 2009).
Calculate MMUs for each state.
# calculate and append
S1 <- mmu.calc(85,100) #85 females per 100m2 in S1
S2 <- mmu.calc(30,100) #30 females per 100m2 in S2
S3 <- mmu.calc(1,10000) #1 female per 10000m2 in S3
mmu_table <- as.data.frame(rbind(S1,S2,S3))
colnames(mmu_table) <- "MMU"
# display table
knitr::kable(mmu_table)
MMU | |
---|---|
S1 | 1 |
S2 | 1 |
S3 | 560 |
Reload the MaxSSS layers and aggregate contiguous pixels for each state. Clumps greater than or equal to the size of the MMU are extracted and saved as new rasters.
# Create an empty vector
clmp_layers <- list()
# For-loop
for (state in 1:3)
{
# Extract threshold raster
recl_thresh <- subset(thresh_123_stack,paste0("MaxSSS_S",state))
# Reclassify 0 to NULL (NA), save as intermediate and reload
s_thresh_recl <- reclassify(recl_thresh, rcl=cbind(0,NA))
writeRaster(s_thresh_recl,paste0(int.dir,"s_thresh_recl_S",state,".tif"),
options="COMPRESS=LZW",overwrite=TRUE)
s_thresh_recl <- get.ras(paste0("s_thresh_recl_S",state,".tif"), ras_dir=int.dir)
# Extract unique class values in the raster
thr_val <- unique(s_thresh_recl)
# Pull up the MMU for each state
mmu_val <- mmu_table$MMU[state]
# skip if MMU is 1
if (mmu_val > 1){
# for-loop to detect contiguous groups of suitable pixels (clumps)
# and extract clumps of a certain cell count
for (i in thr_val)
{
# clump class raster
clmp_thresh <- clump(s_thresh_recl)
# save as intermediate raster, write over and reread
writeRaster(clmp_thresh,paste0(int.dir,"clump_S",state,".tif"),
options="COMPRESS=LZW",overwrite=TRUE)
clmp_thresh <- get.ras(paste0("clump_S",state,".tif"), ras_dir=int.dir)
# calculate frequency of each clump/patch
clmp_freq <- as.data.frame(freq(clmp_thresh))
# store unique clump IDs where cell count is < x (MMU)
remove_ID <- clmp_freq$value[which(clmp_freq$count < mmu_val)]
# assign NA to all clumps whose ID's have frequency x
# (this updates the raster layer; unwanted areas deleted)
s_thresh_recl[clmp_thresh %in% remove_ID] <- NA
}
}
# Save
writeRaster(s_thresh_recl,paste0(int.dir,"MMU_S",state,".tif"),
options="COMPRESS=LZW",overwrite=TRUE)
}
Take the sum of S1-3 to get values only for areas of a minimum size for all behavioural states combined.
# Read new files and stack
# Empty vector
s123_rasters <- list()
# For-loop (S1-3):
for(state in 1:3)
{
raster_names <- list.files(path=paste0(int.dir),
pattern=paste0("MMU_S",state,".tif"),
full.names=TRUE)
# Add to list
s123_rasters[[state]] <- raster_names
}
# Threshold raster stack of partitions (S1-3)
MMU_s123_stack <- do.call(stack,s123_rasters)
# Sum
MMU_sum <- calc(MMU_s123_stack,fun=sum,na.rm=TRUE)
# save intermediate raster
writeRaster(MMU_sum,paste0(int.dir,"MMU_123sum_w_NAs.tif"),
options="COMPRESS=LZW",overwrite=TRUE)
# Mask to study area by reclassifying all values to 1 and multiplying
# (use any raster as template. Here, we used 'recl_thres' from above)
# save for use later
rcl_df <- c(0, 100, 1)
rcl_mtrx <- matrix(rcl_df,ncol=3,byrow = TRUE)
all_1s <- reclassify(recl_thresh,rcl_mtrx)
all_1s <- reclassify(all_1s, rcl=cbind(0,1)) #don't forget the 0's!
writeRaster(all_1s,paste0(int.dir,"all_1s_mask.tif"),
options="COMPRESS=LZW",overwrite=TRUE)
MMU_sum2 <- prod(MMU_sum, all_1s, na.rm=TRUE)
# Save intermediate raster
writeRaster(MMU_sum2,paste0(int.dir,"S123_MMUs_sum_all.tif"),
options="COMPRESS=LZW",overwrite=TRUE)
Since suitable sites of all three behavioural states combined should be of a total MMU size of 562 cells, we reuse the clump() function to extract them. This is the final raster output for the multi-state SDM. It will only have suitable sites of this size, while all others are NULL (NA).
# Reclassify all values of 0 to NULL (NA) from MMU_sum
MMU_sum_recl <- reclassify(MMU_sum, rcl=cbind(0,NA),
filename = tempfile(fileext = ".tif"),
options = "COMPRESS=LZW")
# save intermediate raster
writeRaster(MMU_sum_recl,paste0(int.dir,"S123_MMU_sum_recl.tif"),
options="COMPRESS=LZW",overwrite=TRUE)
# Add total size of MMU across all 3 states
mmu_val <- sum(mmu_table$MMU[1:3])
# Detect contiguous groups of suitable pixels and extract by a certain cell count
# clump class raster
clmp_mmu <- clump(MMU_sum_recl,filename = tempfile(fileext=".tif"),
options="COMPRESS=LZW")
# calculate frequency of each clump/patch
clmp_freq <- as.data.frame(freq(clmp_mmu))
# store unique clump IDs where cell count is < x (MMU)
remove_ID <- clmp_freq$value[which(clmp_freq$count < mmu_val)]
# assign NA to all clumps whose ID's have frequency x (unwanted areas deleted)
MMU_sum_recl[clmp_mmu %in% remove_ID] <- NA
# Reclassify NA to 0
MMU_sum_recl <- reclassify(MMU_sum_recl, rcl=cbind(NA,0))
# Mask to study area by reclassifying all values to 1 and multiplying
MMU_sum_recl <- overlay(MMU_sum_recl, all_1s, fun=function(a,b){return(a*b)})
# Save and clear from workspace to save memory
writeRaster(MMU_sum_recl,paste0(lay.dir,"S123_MMUs_35_females.tif"),
options="COMPRESS=LZW",overwrite=TRUE)
Plot.
# reload
MMU_sum_recl <- get.ras(paste0("S123_MMUs_35_females.tif"), ras_dir=lay.dir)
# Use "cut()" to make 7 classifications (warning: this turns 0 into NA)
MMU_s123_cut <- cut(MMU_sum_recl, c(0,1,10,11,100,101,110,111))
# Reclassify all values to 0 to include 0s in the calculation
MMU_s123_0s <- calc(MMU_sum_recl, fun=function(x){x[x>=0] <-0;return(x)})
MMU_s123_NAs <- calc(MMU_sum_recl, fun=function(x){x[x==NA] <-0;return(x)})
# Unify rasters to create classified layer that includes values of 0
MMU_s123_stack <- stack(MMU_s123_cut,MMU_s123_0s)
MMU_s123_cut2 <- max(MMU_s123_stack,na.rm=TRUE)
# spplot
sum_MMU_plot <- spplot(MMU_s123_cut2,
main=c("Suitable sites for minimum of 35 females (S1-3)"),
width = 1,maxpixels=50000000,
cuts=7,
col.regions=c(state_cols),
colorkey=list(space="bottom",
tick.number=7,
breaks=list(0,1,2,3,4,5,6,7),
labels=list(labels=c("0","1","10","11","100",
"101","110","111",""),
cex = 1)),
par.settings = list(panel.background=list(col="black"))
)
# Save
png(width=1800,height=1800,
paste0(fig.dir,"Multi-state_SDM_MMU_35_females.png"),res=300)
print(sum_MMU_plot)
dev.off()
Plot here.
# View graphic
sum_MMU_plot
As shown here, the sites that were smaller than the MMUs were deleted.
Next, we read the overall suitability raster layer again (S123_MMUs_35_females.tif
) and reclassify it, where MMU clusters are valued to 1, and all other pixels are reclassified to NA (NULL).
# Read raster again
MMU_sum_recl <- raster(paste0(lay.dir,"S123_MMUs_35_females.tif"))
# Reclassify, where MMU clusters are valued at 1, and all other pixels are NA
# Reclassify all values of 0 to NULL (NA) from MMU_sum
MMU_s123_NA <- reclassify(MMU_sum_recl, rcl=cbind(0,NA),
filename = tempfile(fileext = ".tif"),
options = "COMPRESS=LZW")
# save intermediate raster
writeRaster(MMU_s123_NA,paste0(int.dir,"MMU_S123.tif"),
options="COMPRESS=LZW",overwrite=TRUE)
MMU_s123_NA <- calc(MMU_s123_NA, fun=function(x){x[x>0] <-1;return(x)})
# save intermediate raster
writeRaster(MMU_s123_NA,paste0(int.dir,"MMU_S123_all_1s.tif"),
options="COMPRESS=LZW",overwrite=TRUE)
We create a polygon that will be used to extract values from all clusters, making sure to first remove any polygons smaller than the minimum area size for all three states (562 pixels; 351,250 m).
# Convert raster to polygon and convert from multipart to singlepart (currently, features = 1)
MMU_poly <- rasterToPolygons(MMU_s123_NA, na.rm=TRUE, digits = 12, dissolve=TRUE)
MMU_poly <- disaggregate(MMU_poly)
# Project polygons to NZDG2000
proj4string(MMU_poly) <- CRS("+init=epsg:2193")
# Calculate the area of each individual polygon and round to nearest integer
MMU_poly$area <- round(area(MMU_poly))
# Extract polygons >= 351250 m^2 in size.
MMU_s123_poly_FINAL <- subset(MMU_poly, area >= 351250)
# Add a column assigning a unique ID number to each polygon
# (later used for polygon centroids)
MMU_s123_poly_FINAL$id <- 1:nrow(MMU_s123_poly_FINAL)
For the NZSL, suitable sites are areas of a minimum size that contain a location for at least S1, followed by S2 and/or S3. Hence, we prioritise values of 1, 11 and 111, and reclassify a new raster, where values 1, 11 and 111 are converted to 1, and 0 is NA.
# reclassify
MMU_s123_1s <- calc(MMU_sum_recl, fun=function(x){x[x<=0] <-NA;return(x)})
MMU_s123_1s <- calc(MMU_s123_1s, fun=function(x){x[x==11] <-1;return(x)})
MMU_s123_1s <- calc(MMU_s123_1s, fun=function(x){x[x==111] <-1;return(x)})
# save intermediate raster
writeRaster(MMU_s123_1s,paste0(int.dir,"MMU_S123_1_11_111.tif"),
options="COMPRESS=LZW",overwrite=TRUE)
Next, we add a column to the polygons, where the minimum value is extracted from the MMU rasters. We then extract all polygons with a minimum value of “1”. This is the final polygon used that highlights suitable sites for a minimum of 35 individuals across all three behavioural states.
# Add a column where minimum MMU raster value within each polygon is extracted
MMU_s123_poly_FINAL$min_val <- as.numeric(extract(MMU_s123_1s, MMU_s123_poly_FINAL,
fun=min, na.rm=TRUE))
# Make a subset
MMU_poly_1s_NZ <- subset(MMU_s123_poly_FINAL, min_val == 1)
# Save
temp_dir <- file.path(paste0("G:/R/NZSL_MSSDM_MCA_2019_NZ/layers"))
writeOGR(obj=MMU_poly_1s_NZ, dsn=temp_dir, layer="MMU_s123_poly_NZ",
driver="ESRI Shapefile", overwrite_layer = TRUE)
Also save the intermediate layer of all polygons, MMU_poly
(unfiltered by size), in case they’re needed.
# Save intermediate polygon as shapefile
temp_dir <- file.path(paste0("G:/R/NZSL_MSSDM_MCA_2019_NZ/data/intermediate"))
writeOGR(obj=MMU_poly, dsn=temp_dir, layer="MMU_S123_poly_UNFILTERED",
driver="ESRI Shapefile",overwrite_layer=TRUE)
Plot.
# Upload New Zealand coastline polygon
New_Zealand <-readShapePoly(paste0(dat.dir,"Region.shp"))
# Save image
png(width=1600,height=1800,
paste0(fig.dir,"Multi-state_SDM_MMU_35_females_polygons.png"),res=300)
{plot(New_Zealand,
col="gray96",bg="#8DCBE4",
axes=TRUE,
xlim=c(1000000,2005000))
plot(MMU_poly_1s_NZ,col="#E65518",cex=.5, border="#E65518", add=TRUE)
title("Suitable Sites for a minimum of 35 females")}
dev.off()
Plot here.
# plot here
{plot(New_Zealand,
col="gray96",bg="#8DCBE4",
axes=TRUE,
xlim=c(1000000,2005000))
plot(MMU_poly_1s_NZ,col="#E65518",cex=.5, border="#E65518", add=TRUE)
title("Suitable Sites for minimum of 35 females (multi-state SDM)")}
We convert the extracted polygons from above into polygon centroids.
# Polygon centroids by "id"
MMU_points_NZ = gCentroid(MMU_poly_1s, byid=TRUE)
# Save shapefile
shapefile(MMU_points_NZ, paste0(lay.dir,'MMU_S123_point_NZ.shp'),overwrite=TRUE)
Plot.
# Upload New Zealand coastline polygon
New_Zealand <-readShapePoly(paste0(dat.dir,"Region.shp"))
# Save image
png(width=1600,height=1800,
paste0(fig.dir,"Multi-state_SDM_MMU_35_females_points.png"),res=300)
{plot(New_Zealand,
col="gray96",bg="#8DCBE4",
axes=TRUE,
xlim=c(1000000,2005000))
points(MMU_points_NZ,col="#E65518",pch=20,cex=.4)
title("Suitable Sites for a minimum of 35 females")}
dev.off()
Plot here.
# plot here
{plot(New_Zealand,
col="gray96",bg="#8DCBE4",
axes=TRUE,
xlim=c(1000000,2005000))
points(MMU_points_NZ,col="#E65518",pch=20,cex=.4)
title("Suitable Sites for minimum of 35 females (multi-state SDM)")}
# save workspace to load later if needed
save.image("G:/R/mssdm_NZ.RData")
This concludes the multi-state SDM.
Next, we conducted CV, MESS, MoD and limiting factor evaluations. See Appendix S3.
Augé, A. A., Chilvers, B. L., Moore, A., Mathieu, R., & Robertson, B. C. (2009). Aggregation and dispersion of female New Zealand sea lions at the Sandy Bay breeding colony, Auckland Islands: how unusual is their spatial behaviour? Behaviour, 146(9), 1287–1311. https://doi.org/10.1163/15683909X427687
Frans, V.F., Augé, A.A., Edelhoff, H., Erasmi, S., Balkenhol, N. & Engler, J.O. (2018) Quantifying apart what belongs together: A multi-state species distribution modelling framework for species using distinct habitats. Methods in Ecology and Evolution, 9 (1), 98-108. https://doi.org/10.1111/2041-210X.12847
Phillips, S.J. & Dudik, M. (2008) Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography, 31 (2), 161-175. https://doi.org/10.1111/j.0906-7590.2008.5203.x