Â
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 summarised the outputs from the SDM prediction, human impacts, novel preferences in novel spaces and locations of inquiry assessments into seven categories:
Here, we summarised the output from the assessments and additional data layers by using the suitable sites polygons from the multi-state SDM (see Appendix S2) for spatial overlays over the other outputs. Within each site polygon and for each output layer, we extracted the mean, minimum, and/or maximum values, the mode (if the output data layer is categorical), percent coverage, presence or absence of a feature, and/or lists of place names.
We created several custom functions to extract and summarise these data. These are shown and explained in section 2.4 below. These functions work for raster and polygon data, and can be adapted for use in other studies creating integrated SDM databases. Throughout this appendix are examples of using these functions, and how to troubleshoot when warnings arise (e.g. section 6.4.3).
The outputted data fields extracted in this script will remain in their raw format. See Appendix S6 for how to make the finalised version of the database; the simplified data fields are further simplified for use by end-users (managers, rangers, decision-makers, and so on).
The final products from this script are as follows:
Polygon and point shapefiles of sites with 57 data fields summarising the assessment outputs in raw format.
Data table (CSV) with 57 data fields summarising the assessment outputs in raw format.
The script presented here was done using R (version 4.0.5; R Core Team 2021) and its packages.
# libraries
library("raster") # raster data
library("rgdal") # input/output; projections; reading ASCII files
library("tidyverse") # plots and working with complex tables
library("RColorBrewer") # colours for graphics
library("plyr") # sort data
library("reshape") # sort data
library("reshape2") # sort data
library("rgeos") # polygon and point conversions
library('tmap') # visualizing maps
# modify max memory from 10gb to 25gb (this is a 32GB computer)
rasterOptions(maxmemory = 2.5e10)
rasterOptions()
## format : raster
## datatype : FLT4S
## overwrite : FALSE
## progress : none
## timer : FALSE
## chunksize : 1e+08
## maxmemory : 2.5e+10
## memfrac : 0.6
## tmpdir : G:/R/temp\RtmpuQdrJa/raster/
## tmptime : 168
## setfileext : TRUE
## tolerance : 0.1
## standardnames : TRUE
## warn depracat.: TRUE
## header : none
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')
# Intermediate data and layers
dir.create(paste0(dir,'data\\intermediate'),recursive=TRUE)
int.dir <- paste0(dir,'data\\intermediate') # no '\\' here for read/write OGR
# Intermediate layers that have been clipped
dir.create(paste0(int.dir,'\\clipped'),recursive=TRUE)
clp.dir <- paste0(int.dir,'\\clipped') # no '\\' here for read/write OGR
# Final layers
dir.create(paste0(dir,'layers'),recursive=TRUE)
lay.dir <- paste0(dir,'layers') # no '\\' here for read/write OGR
# 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\\')
Colourblind-friendly colours for plotting.
# colours for state values (0-111)
state_cols <-c('#BBBBBB','#DC050C','#882E72','#1965B0','#6195CF','#4EB265',
'#90C987','#F7CB45','#EE8026')
# colours for thresholds or binary plots
thresh_cols <- c('#88CCEE','#CC3311')
# colours for variables
var_cols <-c('#88CCEE','#44AA99','#117733','#999933','#DDCC77','#CC6677',
'#882255','#AA4499')
Calculate the mode of a categorical feature.
# calculate the mode of a categorical feature
calc.mode <- function(x) {
uniq_vals <- unique(x)
uniq_vals[which.max(tabulate(match(x, uniq_vals)))]
}
Extract the name of an inputted variable to convert into a string and use for naming files that are saved.
# turn variable name into a string (for use in file naming)
get.name <- function(x){
y <- deparse(substitute(x))
return(y)
}
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(polydir, paste0(p_name))
# project and return point/polygon
pol <- spTransform(pol, CRS(paste0(prj)))
return(pol)
}
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)
}
Extract text fields from polygons. This works if there are no intersection issues (i.e., no more than one data polygon is found within the site polygon). If there are multiple polygons within a site and you want to have a list of all possible text fields, use the function get.text.more()
below.
Also, see section 6.4.3 for an example of how to handle warnings for duplicates.
# extract text fields from polygons.
# save table/polygon as intermediate but return POLYGON only.
get.text <- function(r, column_names,
# defaults
poly=MMU, poly_keep=c('id','area'),
polydir=int.dir, tabledir=tab.dir){
# extract the attribute(s)/column(s) to keep
col_nm <- column_names
keep <- c(col_nm) #column names to keep
ras <- r[,(names(r) %in% keep)]
# keep only 'id' and 'area' fields for the polygon (default)
poly <- poly[,(names(poly) %in% poly_keep)]
# extract names
inter_poly <- raster::intersect(ras, poly)
# check for duplicates
if (length(inter_poly$id) == length(poly$id)){
print('No duplicates')
} else if (length(inter_poly$id) > length(poly$id)) {
print('WARNING:\nDuplicates! Be sure to summarize/concatenate by ID afterwards.')
} else
print('WARNING:\nSome polygons were missed! Double-check!')
# save as intermediate polygon
shapefile(inter_poly,
paste0(polydir,'\\',paste0(deparse(substitute(r)),'_names.shp')),
overwrite=TRUE)
# save as CSV
tab <- data.frame(inter_poly)
write.csv(tab, paste0(tabledir,
paste0(deparse(substitute(r)),'_text_fields.csv')),
row.names = FALSE)
# return new polygon
return(inter_poly)
}
Sometimes the function get.text()
gives warnings due to intersection issues. When there are multiple polygons of interest within a site and you want a data field to have a list of results (semi-colon separated), then use get.text.more()
. Unlike get.text()
, only the table is returned, not the polygon. This table then has to be concatenated (see section 6.4.4) and appended to the site polygons (see section 7.1).
Also, see section 6.4.3 for an example of how to handle warnings for duplicates.
# extract text fields from polygons.
# save table/polygon as intermediate but return TABLE only.
get.text.more <- function(r, column_names,
# defaults
poly=MMU, poly_keep=c('id','area'),
polydir=int.dir, tabledir=tab.dir){
# extract the attribute(s)/column(s) to keep
col_nm <- column_names
keep <- c(col_nm) #column names to keep
ras <- r[,(names(r) %in% keep)]
# keep only 'id' and 'area' fields for the polygon (default)
poly <- poly[,(names(poly) %in% poly_keep)]
# convert to sf objects
require('sf')
poly_sf <- st_as_sf(poly)
ras_sf <- st_as_sf(ras)
# intersect polygons with points, keeping the information from both
inter_poly <- st_intersection(poly_sf, ras_sf)
# check for duplicates
if (length(inter_poly$id) == length(poly$id)){
print('No duplicates')
} else if (length(inter_poly$id) > length(poly$id)) {
print('WARNING:\nDuplicates! Be sure to summarize/concatenate by ID afterwards.')
} else
print('WARNING:\nSome polygons were missed! Double-check!')
# transform into a 'data.frame' by removing the geometry
st_geometry(inter_poly) <- NULL
# save as CSV
write.csv(inter_poly, paste0(tabledir,
paste0(deparse(substitute(r)),'_text_fields.csv')),
row.names = FALSE)
# print not to append to MMU table
print('WARNING: no polygons were saved. Still need to append table by ID.')
# return new polygon
return(inter_poly)
}
Extract mean, minimum and maximum values from raster layer or raster stack.
# calculate mean/min/max, save table/polygon as intermediate but return TABLE only.
get.mmm <- function(r,
# defaults
poly=MMU, polydir=int.dir,tabledir=tab.dir){
# extract mean, min, max for each polygon
r_means <- raster::extract(r, poly, fun=mean, df=TRUE)
r_min <- raster::extract(r, poly, fun=min, df=TRUE)
r_max <- raster::extract(r, poly, fun=max, df=TRUE)
# rename column names to match variables
colnames(r_means) <- c('ID', paste0(names(r),'_me'))
colnames(r_min) <- c('ID', paste0(names(r),'_mi'))
colnames(r_max) <- c('ID', paste0(names(r),'_mx'))
# combine columns by ID
tab <- join_all(list(r_means,r_min,r_max), by='ID', type='left')
# change 'ID' to 'id' to match shapefile
colnames(tab)[1] <- 'id'
# save as CSV (raw values--not rounded) and reload
write.csv(tab, paste0(tabledir,
paste0(deparse(substitute(r)),'_mean_min_max.csv')),
row.names = FALSE)
tab <- read.csv(paste0(tabledir,
paste0(deparse(substitute(r)),'_mean_min_max.csv')),
header=TRUE)
# round values (4 decimal places)
tab <- tab %>% mutate_at(vars(-1), funs(round(., digits=4)))
# append to polygon
pol <- merge(poly, tab, by='id')
# save as intermediate polygon
shapefile(pol, paste0(polydir,
'\\',paste0(deparse(substitute(r)),'_mean_min_max.shp')),
overwrite=TRUE)
# return table
return(tab)
}
Extract the mode of raster layer or raster stack values (typically for categorical data).
# calculate mode, save table/polygon as intermediate but return TABLE only.
get.mode <- function(r,
# defaults
poly=MMU,polydir=int.dir,tabledir=tab.dir){
# extract all values within polygons
r_vals <- raster::extract(r,poly)
if (nlayers(r) == 1){
# calculate mode (most frequent factor) for each polygon
mode_vals <- unlist(lapply(r_vals,
function(x) if (!is.null(x)) calc.mode(x) else NA ))
mode_vals <- data.frame(id=1:length(MMU),mode_vals)
} else {
# calculate mode (most frequent factor) for each polygon
list <- list()
#for-loop throught all polygons and get mode of each raster layer column
for (i in seq(1, length(poly), by=1)) {
list[[i]] <- unlist(apply(
r_vals[[i]],
2, # by column
function(x) if (!is.null(x)) calc.mode(x) else NA ))
}
# gather and rename columns
mode_vals <- data.frame(id=1:length(poly),do.call(rbind,list))
} #end of if-else loop
# rename columns
colnames(mode_vals) <- c('id', paste0(names(r),'_md'))
# save as CSV and reload
write.csv(mode_vals, paste0(tabledir,
paste0(deparse(substitute(r)),'_mode.csv')),
row.names = FALSE)
mode_vals <- read.csv(paste0(tabledir,
paste0(deparse(substitute(r)),'_mode.csv')),
header=TRUE)
# append to polygon and save as shapefile
pol <- merge(poly, mode_vals, by='id')
shapefile(pol, paste0(polydir,'\\',paste0(deparse(substitute(r)),'_mode.shp')),
overwrite=TRUE)
# return table
return(mode_vals)
}
Extract the presence/absence (1/0) of a feature in the raster layer or raster stack (binary maps required).
# calculate max (1 if present), save table/polygon as intermediate but return TABLE only.
get.pa <- function(r,
# defaults
poly=MMU, polydir=int.dir,tabledir=tab.dir){
# extract max for each polygon (if present, it would be a 1)
r_max <- raster::extract(r, poly, fun=max, df=TRUE)
# rename column names to match variables
colnames(r_max) <- c('ID', paste0(names(r),'_pa'))
# combine columns by ID
tab <- join_all(list(r_max), by='ID', type='left')
# change 'ID' to 'id' to match shapefile
colnames(tab)[1] <- 'id'
# save as CSV (raw values--not rounded) and reload
write.csv(tab, paste0(tabledir,
paste0(deparse(substitute(r)),'_presence.csv')),
row.names = FALSE)
tab <- read.csv(paste0(tabledir,
paste0(deparse(substitute(r)),'_presence.csv')),
header=TRUE)
# round values (4 decimal places)
tab <- tab %>% mutate_at(vars(-1), funs(round(., digits=4)))
# append to polygon
pol <- merge(poly, tab, by='id')
# save as intermediate polygon
shapefile(pol, paste0(polydir,'\\',paste0(deparse(substitute(r)),'_presence.shp')),
overwrite=TRUE)
# return table
return(tab)
}
Calculate the percent of 0’s (unsuitabile pixels) in raster layer or raster stack.
# calculate % of 0's, save table/polygon as intermediate but return TABLE only.
get.perc0 <- function(r,
# defaults
poly=MMU,polydir=int.dir,tabledir=tab.dir){
# extract all values within polygons
r_vals <- raster::extract(r,poly)
if (nlayers(r) == 1){
# calculate mode (most frequent factor) for each polygon
perc_vals <- unlist(lapply(r_vals,
function(x) if (!is.null(x)) (1-(sum(x,
na.rm=TRUE)/length(x))) else NA ))
perc_vals <- data.frame(id=1:length(MMU),perc_vals)
} else {
# calculate mode (most frequent factor) for each polygon
list <- list()
#for-loop throught all polygons and get mode of each raster layer column
for (i in seq(1, length(poly), by=1)) {
list[[i]] <- unlist(apply(
r_vals[[i]],
2, # by column
function(x) if (!is.null(x)) (1-(sum(x,
na.rm=TRUE)/length(x))) else NA ))
}
# gather
perc_vals <- data.frame(id=1:length(poly),do.call(rbind,list))
} #end of if-else loop
# rename columns
colnames(perc_vals) <- c('id', paste0(names(r),'_pc'))
# save as CSV and reload
write.csv(perc_vals, paste0(tabledir,
paste0(deparse(substitute(r)),'_percent.csv')),
row.names = FALSE)
perc_vals <- read.csv(paste0(tabledir,
paste0(deparse(substitute(r)),'_percent.csv')),
header=TRUE)
# append to polygon and save as shapefile
pol <- merge(poly, perc_vals, by='id')
shapefile(pol, paste0(polydir,'\\',paste0(deparse(substitute(r)),'_percent.shp')),
overwrite=TRUE)
# return table
return(perc_vals)
}
Calculate the percent of 1’s (suitable pixels) in raster layer or raster stack.
# calculate % of 1's, save table/polygon as intermediate but return TABLE only.
get.perc1 <- function(r,
# defaults
poly=MMU,polydir=int.dir,tabledir=tab.dir){
# extract all values within polygons
r_vals <- raster::extract(r,poly)
if (nlayers(r) == 1){
# calculate mode (most frequent factor) for each polygon
perc_vals <- unlist(lapply(r_vals,
function(x) if (!is.null(x)) (sum(x,
na.rm=TRUE)/length(x)) else NA ))
perc_vals <- data.frame(id=1:length(MMU),perc_vals)
} else {
# calculate mode (most frequent factor) for each polygon
list <- list()
#for-loop throught all polygons and get mode of each raster layer column
for (i in seq(1, length(poly), by=1)) {
list[[i]] <- unlist(apply(
r_vals[[i]],
2, # by column
function(x) if (!is.null(x)) (sum(x,
na.rm=TRUE)/length(x)) else NA ))
}
# gather
perc_vals <- data.frame(id=1:length(poly),do.call(rbind,list))
} #end of if-else loop
# rename columns
colnames(perc_vals) <- c('id', paste0(names(r),'_pc'))
# save as CSV
write.csv(perc_vals, paste0(tabledir,
paste0(deparse(substitute(r)),'_percent.csv')),
row.names = FALSE)
perc_vals <- read.csv(paste0(tabledir,
paste0(deparse(substitute(r)),'_percent.csv')),
header=TRUE)
# append to polygon and save as shapefile
pol <- merge(poly, perc_vals, by='id')
shapefile(pol, paste0(polydir,'\\',paste0(deparse(substitute(r)),'_percent.shp')),
overwrite=TRUE)
# return table
return(perc_vals)
}
Calculate the percent of a certain value (categorical feature) in a raster or raster stack.
# calculate % of a value (feat_val), save table/polygon as intermediate but return TABLE only.
get.percF <- function(r, feat_val,
# defaults
poly=MMU,polydir=int.dir,tabledir=tab.dir){
# extract all values within polygons
r_vals <- raster::extract(r,poly)
if (nlayers(r) == 1){
# calculate mode (most frequent factor) for each polygon
perc_vals <- unlist(lapply(r_vals,
function(x) if (!is.null(x)) ((sum(x,
na.rm=TRUE)/feat_val[1])/length(x)) else NA ))
perc_vals <- data.frame(id=1:length(MMU),perc_vals)
} else {
# calculate mode (most frequent factor) for each polygon
list <- list()
#for-loop throught all polygons and get mode of each raster layer column
for (i in seq(1, length(poly), by=1)) {
list[[i]] <- unlist(apply(
r_vals[[i]],
2, # by column
function(x) if (!is.null(x)) ((sum(x,
na.rm=TRUE))/length(x)) else NA ))
}
# gather
perc_vals <- data.frame(id=1:length(poly),do.call(rbind,list))
# divide by feature number to convert into proportion
for (n in 1:nlayers(r)){
perc_vals[n+1] <- perc_vals[n+1]/feat_val[n]
}
} #end of if-else loop
# rename columns
colnames(perc_vals) <- c('id', paste0(names(r),'_pc'))
# save as CSV
write.csv(perc_vals, paste0(tabledir,
paste0(deparse(substitute(r)),'_percent.csv')),
row.names = FALSE)
perc_vals <- read.csv(paste0(tabledir,
paste0(deparse(substitute(r)),'_percent.csv')),
header=TRUE)
# append to polygon and save as shapefile
pol <- merge(poly, perc_vals, by='id')
shapefile(pol, paste0(polydir,'\\',paste0(deparse(substitute(r)),'_percent.shp')),
overwrite=TRUE)
# return table
return(perc_vals)
}
Load the MMU (minimum mapping unit) multi-state SDM polygons created in Appendix S2.
# load and ensure projections
MMU <- get.shp('MMU_S123_poly_NZ', polydir = lay.dir)
MMU_pt <- get.shp('MMU_S123_point_NZ', polydir = lay.dir)
Overwrite ID row to make sequential.
# get field names
names(MMU)
# change 'Id' row name to 'MMU_ID' if needed
names(MMU)[2] <- 'MMU_ID'
# make new 'id' field
names(MMU)[1] <- 'id'
MMU@data$id <- 1:nrow(MMU)
# subset data fields to what is needed
MMU <- MMU[c(1:2,5)]
# change name of 3rd field
names(MMU)[3] <- 'area'
Show new field names.
# show fields
MMU
## class : SpatialPolygonsDataFrame
## features : 425
## extent : 1158250, 2089475, 4758650, 6191975 (xmin, xmax, ymin, ymax)
## crs : +init=epsg:2193 +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 3
## names : id, MMU_ID, area
## min values : 1, 1, 353750
## max values : 425, 99, 49636250
# NZ polygon
NZ_poly_shp <- get.shp('NZ_polygon_dissolved')
# NZ region names
region <- get.shp('Region_redo') #this one has a 500-m buffer
# DOC operations regions
op_region <- get.shp('DOC_op_regions')
# DOC conservation areas
cons_area <- get.shp('DOC_cons_areas')
# Grasslands for grazing
grazing <- get.shp('Grazing_grasslands')
# planted pine forests
pine <- get.shp('planted_pine')
# historic NZSL pupping location buffers (10km)
hist_nzsl <- get.shp('historic_NZSL_locations_10km_buffer', polydir = int.dir)
# current female NZSL/pup sighting location buffers (10km)
curr_nzsl <- get.shp('current_NZSL_locations_10km_buffer', polydir = int.dir)
Data from Maxent SDM
# mean suitability scores from Maxent (S1-3)
S1_sdm <- get.ras('S1_projection.tif')
S2_sdm <- get.ras('S2_projection.tif')
S3_sdm <- get.ras('S3_projection.tif')
# raster stack
suitability <- stack(S1_sdm, S2_sdm, S3_sdm)
names(suitability) <- c('S1_sdm', 'S2_sdm', 'S3_sdm')
# thresholds of suitability (MaxSSS) from Maxent (S1-3)
S1_suit <- get.ras('MaxSSS_S1.tif', ras_dir = int.dir)
S2_suit <- get.ras('MaxSSS_S2.tif', ras_dir = int.dir)
S3_suit <- get.ras('MaxSSS_S3.tif', ras_dir = int.dir)
# raster stack
suitability_thresh <- stack(S1_suit, S2_suit, S3_suit)
names(suitability_thresh) <- c('S1_suit', 'S2_suit', 'S3_suit')
# mode of the limiting factors from Maxent runs (S1-3)
S1_lim <- get.ras('S1_limiting_mode.tif')
S2_lim <- get.ras('S2_limiting_mode.tif')
S3_lim <- get.ras('S3_limiting_mode.tif')
# raster stack
limits <- stack(S1_lim, S2_lim, S3_lim)
names(limits) <- c('S1_limits', 'S2_limits', 'S3_limits')
# MESS grid (indicator of extrapolation error)
MESS <- get.ras('MESS.tif')
# MOD grid (shows most dissimilar variable from training area)
MOD <- get.ras('MoD.tif')
# coefficient of variation from Maxent runs (S1-3)
S1_CV <- get.ras('S1_CV.tif')
S2_CV <- get.ras('S2_CV.tif')
S3_CV <- get.ras('S3_CV.tif')
# raster stack
coeff_var <- stack(S1_CV, S2_CV, S3_CV)
Human impact features
# Multi-criteria decision analysis (MCDA) scores and threshold
MCDA <- get.ras('MCA.tif')
MCDA_thresh <- get.ras('MCA_thresh.tif')
# 3D road distances
rds_seal <- get.ras('sealed_roads_3D.tif')
rds_unseal <- get.ras('unsealed_roads_3D.tif')
# raster stack
roads <- stack(rds_seal, rds_unseal)
names(roads) <- c('roads_sealed','roads_unsealed')
# Presence of fences (barriers that can prohibit access to whole area)
fences <- get.ras('fences_pres_abs.tif')
Additional habitat features of interest
# inland water body distance (Euclidean)
in_water <- get.ras('inland_water.tif')
Plot
# plot and save
nzshp <- list("sp.polygons", NZ_poly_shp, fill = "transparent",
col = '#555555', first = TRUE, lwd=0.5)
plt <- spplot(suitability,
main=c("SDM results for each state"),
layout=c(3,1), #cols, rows
maxpixels=500000,
col.regions=c(rev(brewer.pal(11,"Spectral"))),
cuts=10,
colorkey=list(labels=list(at=seq(0,1,by=0.1)),
space="right"),
par.settings = list(panel.background=list(col="black")),
sp.layout=nzshp
)
# save
png(width=2400,height=1500,
paste0(fig.dir,"example_SDM_suitability.png"),res=300)
plt
dev.off()
# plot here
plt
Plot
# plot
nzshp <- list("sp.polygons", NZ_poly_shp, fill = "transparent",
col = '#555555', first = TRUE, lwd=0.5)
plt <- spplot(suitability_thresh,
main=c("SDM suitability thresholds"),
maxpixels=5000000,
col.regions=c(thresh_cols),
cuts=1,
colorkey=list(space="right",
breaks=list(0.5),tick.number=1,
labels=list(labels=c("unsuitable","suitable"),cex = 1)),
par.settings = list(panel.background=list(col="black")),
sp.layout=nzshp
)
# save
png(width = 2400, height = 1200,
paste0(fig.dir,"example_suitability_thresholds.png"),res=300)
print(plt)
dev.off()
# plot here
plt
Plot
# plot
var_labs <- c('cliff','coast\ndist.','forest\ndist.',
'grass\ndist.','land\ncover','sand\ndist.',
'slope','water\ndist.',' ')
plt <- spplot(limits,
main=c('Most frequent limiting factors (mode)'),
maxpixels=500000,
cuts=7,
col.regions=var_cols,
colorkey=list(space='bottom',
breaks=list(0,1,2,3,4,5,6,7,8),
labels=list(labels=var_labs,
at=seq(0,8,by=1),
cex = .8)),
par.settings = list(panel.background=list(col='black'))
)
# save
png(width=2400,height=1500,
paste0(fig.dir,"example_limiting_factors.png"),res=300)
plt
dev.off()
# plot here
plt
Plot
# plot
plt <- spplot(MESS,
main=c('MESS grid'),
maxpixels=500000,
col.regions=c(rev(brewer.pal(10,'Spectral'))),
cuts=9,
par.settings = list(panel.background=list(col='black'))
)
# save
png(width=1800,height=1800,
paste0(fig.dir,"example_MESS_grid.png"),res=300)
plt
dev.off()
# plot here
plt
Plot
# plot
var_labs <- c('cliff','coast\ndist.','forest\ndist.',
'grass\ndist.','land\ncover','sand\ndist.',
'slope','water\ndist.',' ')
plt <- spplot(MOD,
main=c('MOD grid'),
maxpixels=500000,
cuts=7,
col.regions=var_cols,
colorkey=list(space='bottom',
breaks=list(0,1,2,3,4,5,6,7,8),
labels=list(at=c(0,1,1.8,2.8,3.6,4.4,5.3,6.2,7),
labels=var_labs,
cex=0.8)),
par.settings = list(panel.background=list(col='black'))
)
# save
png(width=1800,height=1800,
paste0(fig.dir,"example_MoD_grid.png"),res=300)
plt
dev.off()
# plot here
plt
Plot
# plot
plt <- spplot(coeff_var,
main=c('Coefficient of variation'),
maxpixels=500000,
col.regions=c(rev(brewer.pal(10,'Spectral'))),
cuts=9,
par.settings = list(panel.background=list(col='black'))
)
# save
png(width=2400,height=1500,
paste0(fig.dir,"example_CV.png"),res=300)
plt
dev.off()
# plot here
plt
Plot
# plot
nzshp <- list("sp.polygons", NZ_poly_shp, fill = "transparent",
col = '#555555', first = TRUE, lwd=0.5)
m1 <- spplot(MCDA,
main=c("MCDA of human impacts"),
maxpixels=500000,
col.regions=c(rev(brewer.pal(11,"Spectral"))),
cuts=10,colorkey=list(space="bottom"),
par.settings = list(panel.background=list(col="black")),
sp.layout=nzshp
)
m2 <- spplot(MCDA_thresh,
main=c("MCDA suitability threshold"),
maxpixels=500000,
col.regions=c(thresh_cols[2]),
cuts=0,
colorkey=list(space="bottom",
breaks=list(0.5),tick.number=1,
labels=list(labels=c("suitable"),cex = 1)),
par.settings = list(panel.background=list(col="black")),
sp.layout=nzshp
)
# save
png(width=2400,height=1500, paste0(fig.dir,"example_MCDA.png"),res=300)
print(m1, position = c(0,0,.5,1),more = T)
print(m2, position = c(.465,0,.965,1))
dev.off()
# plot here
print(m1, position = c(0,0,.5,1),more = T)
print(m2, position = c(.465,0,.965,1))
Plot
# plot
nzshp <- list("sp.polygons", NZ_poly_shp, fill = "transparent",
col = '#555555', first = TRUE, lwd=0.5)
plt <- spplot(roads,
main=c("3D road distances in (meters)"),
layout=c(2,1), #cols, rows
maxpixels=500000,
col.regions=c(rev(brewer.pal(8,"Spectral"))),
cuts=7,
par.settings = list(panel.background=list(col="black")),
sp.layout=nzshp
)
# save
png(width=3600,height=1800,
paste0(fig.dir,"example_road_distances.png"),res=300)
plt
dev.off()
# plot here
plt
Plot
# plot
nzshp <- list("sp.polygons", NZ_poly_shp, fill = "transparent",
col = '#555555', first = TRUE, lwd=0.5)
plt <- spplot(fences,
main=c("Presence/absence of fences"),
maxpixels=500000,
col.regions=c(thresh_cols),
cuts=1,
colorkey=list(space="bottom",
breaks=list(0,1),tick.number=1,
labels=list(labels=c("absent","present"),cex = 1)),
par.settings = list(panel.background=list(col="black")),
sp.layout=nzshp
)
# save
png(width=1800,height=1800,
paste0(fig.dir,"example_fences.png"),res=300)
plt
dev.off()
# plot here
plt
We will use the most recent data from 2016 (‘SUBID_2016’). Dairy and Non-Dairy are ID numbers 503 and 502, which were already extracted and prepared.
Plot
# plot and save
plt <- spplot(grazing, "SUBID_2016",
main = "Grasslands for grazing (dairy and non-dairy)",
col.regions = c("#EE7733", "#009988"),
col = "transparent",
par.settings = list(panel.background=list(col="black")),
sp.layout=nzshp)
png(width=1800,height=1800,
paste0(fig.dir,"example_NZ_grazing_grasslands.png"),res=300)
plt
dev.off()
plt
We will use the most recent data from 2016 (‘SUBID_2016’). Commercial planted pine forests polygons are from ID number 201, which were already extracted and prepared.
Plot
# plot and save
plt <- spplot(pine, "SUBID_2016",
main = "Planted pine forests",
col.regions = "#009988",
col = "transparent",
par.settings = list(panel.background=list(col="black")),
sp.layout=nzshp)
png(width=1800,height=1800,
paste0(fig.dir,"example_NZ_planted_pine_forests.png"),res=300)
plt
dev.off()
plt
Plot
# plot
plt <- spplot(in_water,
main=c("Distance from inland water bodies"),
maxpixels=500000,
col.regions=c(rev(brewer.pal(11,"Spectral"))),
cuts=7,colorkey=list(space="bottom"),
par.settings = list(panel.background=list(col="black")),
sp.layout=nzshp
)
# save
png(width=1800,height=1800,
paste0(fig.dir,"example_inland_water_dist.png"),res=300)
plt
dev.off()
# plot here
plt
Plot
# plot and save
plt <- spplot(region, "REGC2019_1",
par.settings = list(panel.background=list(col="#8DCBE4")))
png(width=1800,height=1800,
paste0(fig.dir,"example_NZ_regions.png"),res=300)
plt
dev.off()
plt
Plot
# plot and save
plt <- spplot(op_region,"RegionName", #fill='gray96',
par.settings = list(panel.background=list(col="#8DCBE4")))
png(width=1800,height=1800,
paste0(fig.dir,"example_DOC_op_regions.png"),res=300)
plt
dev.off()
plt
Plot
# plot and save
plt1 <- spplot(cons_area, "Name", main = "DOC conservation areas (names)",
col = "transparent",
par.settings = list(panel.background=list(col="#8DCBE4")))
plt2 <- spplot(cons_area, "Type", main = "DOC conservation areas (type)",
col = "transparent",
par.settings = list(panel.background=list(col="#8DCBE4")),
sp.layout = list(list(NZ_poly_shp, fill="#BBBBBB",lwd=0.5, first=TRUE)))
plt3 <- spplot(cons_area, "Conservati", main = "DOC conservation areas (coded)",
col = "transparent",
par.settings = list(panel.background=list(col="#8DCBE4")))
# save
png(width=2000,height=3000,
paste0(fig.dir,"example_cons_area_names.png"),res=300)
plt1
dev.off()
png(width=2000,height=3000,
paste0(fig.dir,"example_cons_area_type.png"),res=300)
plt2
dev.off()
png(width=2000,height=3000,
paste0(fig.dir,"example_cons_area_codes.png"),res=300)
plt3
dev.off()
# plot example here (others have too much categorical data to see)
plt2
Plot
# plot and save
plt <- spplot(hist_nzsl, "name",
main = "Historic NZSL breeding sites (10km buffers)",
par.settings = list(panel.background=list(col="#8DCBE4")),
sp.layout = list(list(NZ_poly_shp, fill="#BBBBBB", first=TRUE)))
png(width=1800,height=1800,
paste0(fig.dir,"example_NZ_historic_breeding_sites_buffers.png"),res=300)
plt
dev.off()
plt
Plot
# plot and save
plt <- spplot(curr_nzsl, "name",
main = "Current female/pup sightings (10km buffers)",
par.settings = list(panel.background=list(col="#8DCBE4")),
sp.layout = list(list(NZ_poly_shp, fill="#BBBBBB", first=TRUE)))
png(width=1800,height=1800,
paste0(fig.dir,"example_NZ_Current_female_pup_sightings_buffers.png"),res=300)
plt
dev.off()
plt
Rasterize polygons to make a mask (snapped to sdm_1
).
First, create the template from sdm_1
.
# rasterize by ID and make non-overlaying polygons 'NA'
MMU_ras <- rasterize(MMU, suitability[[1]], 'id', background=NA)
# save
writeRaster(MMU_ras,paste0(clp.dir,'\\MMU_ras.tif'),options="COMPRESS=LZW",overwrite=TRUE)
Mask over the raster layers and save them.
# overwrite and mask the rasters and stacks to suitable sites
# save to clp.dir (clipped data directory) and compress.
MESS <- mask(MESS,MMU_ras)
writeRaster(MESS,paste0(clp.dir,'\\MESS.tif'),options="COMPRESS=LZW",
overwrite=TRUE)
MOD <- mask(MOD,MMU_ras)
writeRaster(MOD,paste0(clp.dir,'\\MOD.tif'),options="COMPRESS=LZW",
overwrite=TRUE)
MCDA <- mask(MCDA,MMU_ras)
writeRaster(MCDA,paste0(clp.dir,'\\MCDA.tif'),options="COMPRESS=LZW",
overwrite=TRUE)
MCDA_thresh <- mask(MCDA_thresh,MMU_ras)
writeRaster(MCDA_thresh,paste0(clp.dir,'\\MCDA_thresh.tif'),options="COMPRESS=LZW",
overwrite=TRUE)
roads <- mask(roads,MMU_ras)
writeRaster(roads,paste0(clp.dir,'\\',names(roads),'.tif'),
bylayer=TRUE,options="COMPRESS=LZW",overwrite=TRUE)
fences <- mask(fences,MMU_ras)
writeRaster(fences,paste0(clp.dir,'\\fences.tif'),options="COMPRESS=LZW",
overwrite=TRUE)
in_water <- mask(in_water,MMU_ras)
writeRaster(in_water,paste0(clp.dir,'\\in_water.tif'),options="COMPRESS=LZW",
overwrite=TRUE)
suitability <- mask(suitability, MMU_ras)
writeRaster(suitability,paste0(clp.dir,'\\',names(suitability),'.tif'),
bylayer=TRUE,options="COMPRESS=LZW",overwrite=TRUE)
suitability_thresh <- mask(suitability_thresh, MMU_ras)
writeRaster(suitability_thresh,paste0(clp.dir,'\\',names(suitability_thresh),'.tif'),
bylayer=TRUE,options="COMPRESS=LZW",overwrite=TRUE)
limits <- mask(limits,MMU_ras)
writeRaster(limits,paste0(clp.dir,'\\',names(limits),'.tif'),
bylayer=TRUE,options="COMPRESS=LZW",overwrite=TRUE)
coeff_var <- mask(coeff_var,MMU_ras)
writeRaster(coeff_var,paste0(clp.dir,'\\',names(coeff_var),'.tif'),
bylayer=TRUE,options="COMPRESS=LZW",overwrite=TRUE)
Create a table with the mean, minimum, and maximum SDM suitability scores per state and per polygon.
# rename to abbreviate
names(suitability) <- c('S1_sdm', 'S2_sdm', 'S3_sdm')
# calculate mean, minimum and maximum for each polygon
sdm_mmm <- get.mmm(suitability)
# rename to abbreviate
names(suitability_thresh) <- c('S1_suit', 'S2_suit', 'S3_suit')
# calculate mean, minimum and maximum for each polygon
suit_area <- get.percF(suitability_thresh, feat_val = c(1,10,100))
Create a table with the mean, minimum and maximum CV values per polygon.
# calculate mean, minimum and maximum for each polygon
CV_mmm <- get.mmm(coeff_var)
Create a table with the mode (most frequent) limiting factors per state SDM and per polygon.
# rename to abbreviate
names(limits) <- c('S1_lim', 'S2_lim', 'S3_lim')
# calculate mode for each polygon
lim_mode <- get.mode(limits)
Create a table with the mean, minimum and maximum MESS scores per polygon.
# calculate mean, minimum and maximum for each polygon
mess_mmm <- get.mmm(MESS)
Create a table with the mode (most frequent) limiting factors per polygon, which corresponds to the MESS grid.
# calculate mode for each polygon
MoD_mode <- get.mode(MOD)
Create a table with the mean, minimum and maximum MCDA scores per polygon.
# rename to abbreviate
names(MCDA) <- 'MCDA'
# calculate mean, minimum and maximum for each polygon
MCDA_mmm <- get.mmm(MCDA)
Get the proportion of area that is affected by human impacts (unsuitable MCDA scores; below the threshold).
# rename to abbreviate
names(MCDA_thresh) <- 'hum_im'
hum_im <- MCDA_thresh
# calculate mode for each polygon
hum_area <- get.perc0(hum_im)
Create a table with the mean, minimum and maximum road distances (3D) per polygon.
# rename to abbreviate
names(roads) <- c('rd_sl','rd_unsl')
# calculate mean, minimum and maximum for each polygon
rds_mmm <- get.mmm(roads)
# rename to abbreviate
names(fences) <- 'fence'
fence <- fences
# calc presence absence
fnc_pa <- get.pa(fence)
Get names of polygons, and calculate the area and proportion of MMU that has grazing grasslands.
# quick buffer of 0 to avoid ID error
grazing <- gBuffer(grazing, byid=TRUE, width=0)
# rename to abbreviate fields
names(grazing@data)[names(grazing@data)=='SUBID_2016'] <- 'grz_name'
# get names
col_nm <- c('grz_name')
grz_names <- get.text(grazing,col_nm)
# calculate area and summarize per ID. Include percent cover of these fields.
grz_names$grz_size <- round(area(grz_names))
# concatenate multiple fields
grz_tab <- ddply(grz_names@data, .(id, area), summarize,
#grz_name=paste(grz_name,collapse="; "),
grz_size=round(sum(grz_size)),
grz_perc=round((sum(grz_size)/unique(area)),digits = 4))
# new merged table to save (no data files will be NA)
grz_names <- sp::merge(MMU,grz_tab, by=c('id','area'), all.x=TRUE)
# see first few rows
head(grz_names); tail(grz_names)
# show fields with NA
#grz_names@data[!complete.cases(grz_names@data),]
# save as CSV
write.csv(grz_names, paste0(tab.dir,paste0('grazing_names_sizes.csv')),
row.names = FALSE)
# save as shapefile
shapefile(grz_names, paste0(int.dir,'\\grazing_names_sizes.shp'),overwrite=TRUE)
Get names of polygons, and calculate the area and proportion of MMU that has planted pine forests.
# rename to abbreviate fields
names(pine@data)[names(pine@data)=='SUBID_2016'] <- 'pine_name'
# get names
col_nm <- c('pine_name')
pine_names <- get.text(pine,col_nm)
# calculate area and summarize per ID. Include percent cover of these fields.
pine_names$pine_size <- round(area(pine_names))
# concatenate multiple fields
pine_tab <- ddply(pine_names@data, .(id, area), summarize,
pine_size=round(sum(pine_size)),
pine_perc=round((sum(pine_size)/unique(area)),digits = 4))
# new merged table to save (no data files will be NA)
pine_names <- sp::merge(MMU,pine_tab, by=c('id','area'), all.x=TRUE)
# see first few rows
head(pine_names); tail(pine_names)
# show fields with NA
#pine_names@data[complete.cases(pine_names@data),]
# save as CSV
write.csv(pine_names, paste0(tab.dir,paste0('pine_names_sizes.csv')),
row.names = FALSE)
# save as shapefile
shapefile(pine_names, paste0(int.dir,'\\pine_names_sizes.shp'),overwrite=TRUE)
Create a table with the mean, minimum, and maximum inland water distance (Euclidean) per polygon.
# rename to abbreviate
names(in_water) <- c('in_watr')
# calculate mean, minimum and maximum for each polygon
water_mmm <- get.mmm(in_water)
Get name of country region and append to polygons.
Because some borders between regions will be shared at some sites, we prepared region names as a raster and then used get.mode
to extract the region name where the majority of a site is located.
# load float raster of the region
reg_ras <- get.ras('Region_FLOAT.tif')
# mask to the suitable sites (MMUs)
reg_ras <- mask(reg_ras,MMU_ras)
# save
writeRaster(reg_ras,paste0(clp.dir,'\\Region.tif'),
options="COMPRESS=LZW",overwrite=TRUE)
Create a table with the mode (most frequent) region name per polygon (majority will then be the region name; accounts for sites at the borderlines).
# rename to abbreviate field
names(reg_ras) <- 'region'
region <- reg_ras
# calculate mode for each polygon
reg_names <- get.mode(region)
# rename column
colnames(reg_names) <- c('id', 'region')
# Manually edit one of the entries
#reg_names[reg_names$region > 17,]
reg_names[!complete.cases(reg_names),][2] <- 4 #ID#129
Get names of DOC operation regions and append to polygon. Similar to above, we do this with a raster, since shared borders can cause some problems with get.text()
.
# load float raster
op_region_ras <- get.ras('DOC_op_regions.tif')
# mask to the suitable sites (MMUs)
op_region_ras <- mask(op_region_ras,MMU_ras)
# save
writeRaster(op_region_ras,paste0(clp.dir,'\\DOC_op_regions.tif'),
options="COMPRESS=LZW",overwrite=TRUE)
Create a table with the mode (most frequent) region name per polygon (majority will then be the region name; this accounts for sites at the borderlines).
# rename to abbreviate field
names(op_region_ras) <- 'DOC_region'
op_region <- op_region_ras
# calculate mode for each polygon
op_names <- get.mode(op_region)
# rename column
colnames(op_names) <- c('id', 'DOC_region')
Get names of conservation areas, and calculate the area and proportion of MMU that is a conservation area.
# rename to abbreviate fields
names(cons_area@data)[names(cons_area@data)=='Name'] <- 'DOC_name'
names(cons_area@data)[names(cons_area@data)=='Conservati'] <- 'DOC_code'
names(cons_area@data)[names(cons_area@data)=='Conservati'] <- 'DOC_type'
# get names
col_nm <- c('DOC_name','DOC_code','DOC_type')
cons_names <- get.text(cons_area,col_nm)
# calculate area and summarize per ID. Include percent cover of these fields.
cons_names$DOC_size <- round(area(cons_names))
# get total area by aggregating duplicates of ID
DOC_ttl <- aggregate(DOC_size~id, data=cons_names, FUN=sum)
colnames(DOC_ttl) <- c('id','DOC_ttl')
sp::merge(cons_names,DOC_ttl, by='id')
# calculate proportion of total polygon size that has that feature
cons_names$DOC_perc <- round(cons_names$DOC_size/cons_names$area, digits = 4)
# save as CSV
write.csv(cons_names, paste0(tab.dir,paste0('cons_area_names_sizes.csv')),
row.names = FALSE)
# save as shapefile
shapefile(cons_names, paste0(int.dir,'\\cons_area_names_sizes.shp'),overwrite=TRUE)
[1] "WARNING:\nDuplicates! Be sure to summarize/concatenate by ID afterwards."
There are duplicates, so modifications need to be made before finalising the shapefile. It also cuts out polygons where there are no data (no polygons). We need to concatenate the multiple conservation area names, and then merge it once again with the MMU polygon, to account for the sites that have NAs (i.e. are not within protected areas). We use ddply()
to make lists of names, and also make some calculations on coverage.
# retrieve files again from 'get.text' command
cons_names <- get.shp('cons_area_names', polydir = int.dir)
cons_names_txt <- read.csv(paste0(tab.dir,paste0('cons_area_text_fields.csv')),
header = TRUE)
# calculate area and summarize per ID. Include percent cover of these fields.
cons_names$DOC_size <- area(cons_names)
# concatenate multiple fields
cons_tab <- ddply(cons_names@data, .(id, area), summarize,
#list of conservation area codes
DOC_code=paste(DOC_code,collapse="; "),
# list of conservation area names
DOC_name=paste(DOC_name,collapse="; "),
# size of the conservation area within the site
DOC_size=round(sum(DOC_size)),
# percent coverage of the conservation area within the site
DOC_perc=round((sum(DOC_size)/unique(area)),digits = 4))
# new merged table to save (no data files will be NA)
cons_names <- sp::merge(MMU,cons_tab, by=c('id','area'), all.x=TRUE)
# see first few rows
head(cons_names); tail(cons_names)
# show fields with NA
#cons_names@data[!complete.cases(cons_names@data),]
Save and overwrite the previous files that were made from get.text.more()
.
# save as CSV
write.csv(cons_names, paste0(tab.dir,paste0('cons_area_names_sizes.csv')),
row.names = FALSE)
# save as shapefile
shapefile(cons_names, paste0(int.dir,'\\cons_area_names_sizes.shp'),overwrite=TRUE)
List names of known historic pupping locations within 10km of a site.
# rename to abbreviate fields
names(hist_nzsl@data)[names(hist_nzsl@data)=='name'] <- 'hist_name'
# get names
col_nm <- c('hist_name')
hist_names <- get.text.more(hist_nzsl,col_nm)
# concatenate multiple fields
hist_tab <- ddply(hist_names, .(id), summarize,
hist_name=paste(hist_name,collapse="; "))
# new merged table to save (no data files will be NA)
hist_names <- sp::merge(MMU,hist_tab, by=c('id'), all.x=TRUE)
# see first few rows
head(hist_names); tail(hist_names)
# show completed fields
hist_tab
Save.
# save as CSV
write.csv(hist_names, paste0(tab.dir,paste0('hist_nzsl_names_sizes.csv')),
row.names = FALSE)
# save as shapefile
shapefile(hist_names, paste0(int.dir,'\\hist_nzsl_names_sizes.shp'),overwrite=TRUE)
Get names of current breeding locations within 10km of a site.
# rename to abbreviate fields
names(curr_nzsl@data)[names(curr_nzsl@data)=='name'] <- 'curr_name'
# get names
col_nm <- c('curr_name')
curr_names <- get.text.more(curr_nzsl,col_nm)
# concatenate multiple fields
curr_tab <- ddply(curr_names, .(id), summarize,
curr_name=paste(curr_name,collapse="; "))
# new merged table to save (no data files will be NA)
curr_names <- sp::merge(MMU,curr_tab, by=c('id'), all.x=TRUE)
# see first few rows
#head(curr_names); tail(curr_names)
# show completed fields
#curr_tab
Save
# save as CSV
write.csv(curr_names, paste0(tab.dir,paste0('curr_nzsl_names_sizes.csv')),
row.names = FALSE)
# save as shapefile
shapefile(curr_names,
paste0(int.dir,'\\curr_nzsl_names_sizes.shp'),overwrite=TRUE)
# Polygon to dataframe
reg_names <- data.frame(reg_names[-3]) #remove 'area' field (repeated)
op_names <- data.frame(op_names[-3]) #remove 'area' field (repeated)
cons_names <- data.frame(cons_names[-c(2,3)]) #remove 'area' and 'MMU_ID' fields
grz_names <- data.frame(grz_names[-c(2,3)]) #remove 'area' and 'MMU_ID' fields
pine_names <- data.frame(pine_names[-c(2,3)]) #remove 'area' and 'MMU_ID' fields
curr_names <- data.frame(curr_names[-c(2,3)]) #remove 'area' and 'MMU_ID' fields
hist_names <- data.frame(hist_names[-c(2,3)]) #remove 'area' and 'MMU_ID' fields
# combine all tables by ID
features <- Reduce(function(x, y) merge(x, y, by='id', all=TRUE),
list(
#name items first
reg_names,op_names,cons_names,
#environmental suitability items
sdm_mmm,suit_area,lim_mode,mess_mmm,
MoD_mode, CV_mmm,
#human impact items
MCDA_mmm,hum_area,rds_mmm,fnc_pa,grz_names,
#other habitat selection features
water_mmm, pine_names, curr_names, hist_names))
# view here
summary(features)
## id region DOC_region DOC_code
## Min. : 1 Min. : 1.000 Min. :1.000 Length:425
## 1st Qu.:107 1st Qu.: 3.000 1st Qu.:3.000 Class :character
## Median :213 Median : 7.000 Median :5.000 Mode :character
## Mean :213 Mean : 7.162 Mean :4.871
## 3rd Qu.:319 3rd Qu.:11.000 3rd Qu.:7.000
## Max. :425 Max. :16.000 Max. :8.000
##
## DOC_name DOC_size DOC_perc S1_sdm_me
## Length:425 Min. : 0 Min. :0.0000 Min. :0.0018
## Class :character 1st Qu.: 32236 1st Qu.:0.0259 1st Qu.:0.1141
## Mode :character Median : 142288 Median :0.1216 Median :0.1717
## Mean : 498179 Mean :0.2412 Mean :0.1623
## 3rd Qu.: 496177 3rd Qu.:0.3861 3rd Qu.:0.2141
## Max. :8768803 Max. :1.0000 Max. :0.3410
## NA's :159 NA's :159
## S2_sdm_me S3_sdm_me S1_sdm_mi S2_sdm_mi
## Min. :0.0023 Min. :0.0013 Min. :0.00000 Min. :0.00180
## 1st Qu.:0.0230 1st Qu.:0.1512 1st Qu.:0.00010 1st Qu.:0.00280
## Median :0.0978 Median :0.2558 Median :0.00030 Median :0.00450
## Mean :0.1129 Mean :0.2485 Mean :0.03604 Mean :0.00668
## 3rd Qu.:0.1821 3rd Qu.:0.3447 3rd Qu.:0.09360 3rd Qu.:0.00750
## Max. :0.4356 Max. :0.5984 Max. :0.21920 Max. :0.16410
##
## S3_sdm_mi S1_sdm_mx S2_sdm_mx S3_sdm_mx
## Min. :0.00100 Min. :0.1291 Min. :0.0034 Min. :0.0038
## 1st Qu.:0.00400 1st Qu.:0.3003 1st Qu.:0.0535 1st Qu.:0.4572
## Median :0.01290 Median :0.4809 Median :0.4591 Median :0.7942
## Mean :0.02247 Mean :0.4631 Mean :0.3819 Mean :0.6817
## 3rd Qu.:0.02780 3rd Qu.:0.5945 3rd Qu.:0.6223 3rd Qu.:0.9174
## Max. :0.20170 Max. :0.6874 Max. :0.7384 Max. :0.9824
##
## S1_suit_pc S2_suit_pc S3_suit_pc S1_lim_md
## Min. :0.001046 Min. :0.0000 Min. :0.0000 Min. :3.000
## 1st Qu.:0.343053 1st Qu.:0.0000 1st Qu.:0.4840 1st Qu.:3.000
## Median :0.612348 Median :0.2172 Median :0.7834 Median :3.000
## Mean :0.619212 Mean :0.2596 Mean :0.6627 Mean :3.861
## 3rd Qu.:0.998075 3rd Qu.:0.4474 3rd Qu.:0.9121 3rd Qu.:5.000
## Max. :1.000000 Max. :1.0000 Max. :1.0000 Max. :5.000
##
## S2_lim_md S3_lim_md MESS_me MESS_mi
## Min. :3.000 Min. :1.000 Min. :-538.930 Min. :-596.8150
## 1st Qu.:3.000 1st Qu.:3.000 1st Qu.: -36.520 1st Qu.:-110.9590
## Median :3.000 Median :3.000 Median : -17.154 Median : -80.0000
## Mean :3.504 Mean :3.325 Mean : -35.106 Mean : -87.0913
## 3rd Qu.:3.000 3rd Qu.:3.000 3rd Qu.: -3.079 3rd Qu.: -80.0000
## Max. :7.000 Max. :7.000 Max. : 3.438 Max. : 0.2387
##
## MESS_mx MOD_md S1_CV_me S2_CV_me
## Min. :-475.1720 Min. :2.000 Min. :0.4927 Min. :0.1278
## 1st Qu.: 0.3406 1st Qu.:4.000 1st Qu.:0.7018 1st Qu.:0.3797
## Median : 3.9428 Median :5.000 Median :1.0825 Median :0.4730
## Mean : -9.0145 Mean :5.205 Mean :1.0836 Mean :0.4764
## 3rd Qu.: 8.7376 3rd Qu.:6.000 3rd Qu.:1.3484 3rd Qu.:0.5673
## Max. : 15.6397 Max. :7.000 Max. :2.1955 Max. :0.7937
##
## S3_CV_me S1_CV_mi S2_CV_mi S3_CV_mi
## Min. :0.1988 Min. :0.0365 Min. :0.0286 Min. :0.0128
## 1st Qu.:0.3820 1st Qu.:0.1599 1st Qu.:0.0713 1st Qu.:0.0861
## Median :0.5439 Median :0.3097 Median :0.1363 Median :0.1356
## Mean :0.6550 Mean :0.3041 Mean :0.1893 Mean :0.2853
## 3rd Qu.:0.7968 3rd Qu.:0.4466 3rd Qu.:0.2866 3rd Qu.:0.4566
## Max. :1.9154 Max. :0.7110 Max. :0.6470 Max. :1.4889
##
## S1_CV_mx S2_CV_mx S3_CV_mx MCDA_me
## Min. :0.7457 Min. :0.4467 Min. :0.5461 Min. :0.0385
## 1st Qu.:0.9697 1st Qu.:0.8597 1st Qu.:0.9645 1st Qu.:0.5824
## Median :2.7121 Median :0.9002 Median :1.2343 Median :0.7123
## Mean :2.2749 Mean :0.8887 Mean :1.3172 Mean :0.7036
## 3rd Qu.:3.0084 3rd Qu.:0.9297 3rd Qu.:1.6795 3rd Qu.:0.9237
## Max. :3.2483 Max. :1.0752 Max. :2.1383 Max. :1.0000
##
## MCDA_mi MCDA_mx hum_im_pc rd_sl_me
## Min. :0.0000 Min. :0.1254 Min. :0.0000 Min. : 49.7
## 1st Qu.:0.2862 1st Qu.:0.6958 1st Qu.:0.0000 1st Qu.: 509.4
## Median :0.6550 Median :0.8326 Median :0.0000 Median : 1572.8
## Mean :0.5757 Mean :0.8157 Mean :0.1552 Mean : 16655.3
## 3rd Qu.:0.7790 3rd Qu.:1.0000 3rd Qu.:0.1673 3rd Qu.: 5037.4
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :337226.0
##
## rd_unsl_me rd_sl_mi rd_unsl_mi rd_sl_mx
## Min. : 98.55 Min. : 0.0 Min. : 0.00 Min. : 176.9
## 1st Qu.: 651.63 1st Qu.: 0.0 1st Qu.: 0.00 1st Qu.: 1423.9
## Median : 1097.68 Median : 135.4 Median : 75.44 Median : 2950.9
## Mean : 6415.80 Mean : 15346.9 Mean : 5419.48 Mean : 18300.8
## 3rd Qu.: 1954.90 3rd Qu.: 2908.8 3rd Qu.: 746.79 3rd Qu.: 7604.0
## Max. :149631.00 Max. :337226.0 Max. :149631.00 Max. :337226.0
##
## rd_unsl_mx fence_pa grz_size grz_perc
## Min. : 332.9 Min. :0.0000 Min. : 1 Min. :0.00000
## 1st Qu.: 1487.7 1st Qu.:0.0000 1st Qu.: 49891 1st Qu.:0.05618
## Median : 2355.9 Median :0.0000 Median : 151267 Median :0.15255
## Mean : 7749.4 Mean :0.4659 Mean : 342724 Mean :0.20356
## 3rd Qu.: 3861.9 3rd Qu.:1.0000 3rd Qu.: 413279 3rd Qu.:0.33345
## Max. :149631.0 Max. :1.0000 Max. :3750549 Max. :0.91720
## NA's :51 NA's :51
## in_watr_me in_watr_mi in_watr_mx pine_size
## Min. : 92.22 Min. : 0.0 Min. : 237.2 Min. : 3
## 1st Qu.: 505.42 1st Qu.: 0.0 1st Qu.: 1360.1 1st Qu.: 5640
## Median : 895.27 Median : 25.0 Median : 2269.9 Median : 12867
## Mean : 1677.31 Mean : 718.3 Mean : 3140.2 Mean : 42427
## 3rd Qu.: 1695.79 3rd Qu.: 246.2 3rd Qu.: 3893.9 3rd Qu.: 57694
## Max. :24100.36 Max. :23354.0 Max. :24554.3 Max. :157452
## NA's :418
## pine_perc curr_name hist_name
## Min. :0.0000 Length:425 Length:425
## 1st Qu.:0.0008 Class :character Class :character
## Median :0.0044 Mode :character Mode :character
## Mean :0.0366
## 3rd Qu.:0.0233
## Max. :0.2038
## NA's :418
# join table polygon layer
MMU2 <- MMU
MMU2@data <- join(MMU2@data,features,by='id')
# show column names
colnames(MMU2@data)
## [1] "id" "MMU_ID" "area" "region" "DOC_region"
## [6] "DOC_code" "DOC_name" "DOC_size" "DOC_perc" "S1_sdm_me"
## [11] "S2_sdm_me" "S3_sdm_me" "S1_sdm_mi" "S2_sdm_mi" "S3_sdm_mi"
## [16] "S1_sdm_mx" "S2_sdm_mx" "S3_sdm_mx" "S1_suit_pc" "S2_suit_pc"
## [21] "S3_suit_pc" "S1_lim_md" "S2_lim_md" "S3_lim_md" "MESS_me"
## [26] "MESS_mi" "MESS_mx" "MOD_md" "S1_CV_me" "S2_CV_me"
## [31] "S3_CV_me" "S1_CV_mi" "S2_CV_mi" "S3_CV_mi" "S1_CV_mx"
## [36] "S2_CV_mx" "S3_CV_mx" "MCDA_me" "MCDA_mi" "MCDA_mx"
## [41] "hum_im_pc" "rd_sl_me" "rd_unsl_me" "rd_sl_mi" "rd_unsl_mi"
## [46] "rd_sl_mx" "rd_unsl_mx" "fence_pa" "grz_size" "grz_perc"
## [51] "in_watr_me" "in_watr_mi" "in_watr_mx" "pine_size" "pine_perc"
## [56] "curr_name" "hist_name"
The multi-state SDM gave an output of 425 sites. These will be further reduced if there are sites that are too large and only suitable for the first behavioural state.
First, determine how many sites will be removed.
# make a subset of both S2 and S3 being unsuitable
MMU3 <- MMU2[((MMU2$S2_suit_pc==0)&(MMU2$S3_suit_pc==0)),]
# show how many will be removed
dim(MMU3)
## [1] 30 57
30 sites will be removed.
# Remove these polygons
MMU3 <- MMU2[!((MMU2$S2_suit_pc==0)&(MMU2$S3_suit_pc==0)),]
This is the final polygon.
# final polygon details
MMU3
## class : SpatialPolygonsDataFrame
## features : 395
## extent : 1167775, 2086175, 4758650, 6191975 (xmin, xmax, ymin, ymax)
## crs : +init=epsg:2193 +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 57
## names : id, MMU_ID, area, region, DOC_region, DOC_code, DOC_name, DOC_size, DOC_perc, S1_sdm_me, S2_sdm_me, S3_sdm_me, S1_sdm_mi, S2_sdm_mi, S3_sdm_mi, ...
## min values : 1, 1, 353750, 1, 1, 00001, Abel Tasman National Park; Abel Tasman Scenic Reserve; Local Purpose Reserve - Cemetery - Awaroa Inlet, 0, 0, 0.0018, 0.0041, 0.0194, 0, 0.0018, 0.001, ...
## max values : 425, 99, 49636250, 16, 8, W15111, Whitiau Scenic Reserve; Whitiau Scientific Reserve; Kaitoke Marginal Strip, 8768803, 1, 0.341, 0.4356, 0.5984, 0.2192, 0.1641, 0.2017, ...
There are 395 sites as the final output.
# save as final polygon (for analysis)
shapefile(MMU3, paste0(lay.dir,'\\MMU_s123_poly_features_RAW_shape.shp'),
overwrite=TRUE)
# Save
temp_dir <- file.path(paste0("G:/R/NZSL_MSSDM_MCA_2019_NZ/layers"))
writeOGR(obj=MMU3, dsn=temp_dir, layer="MMU_s123_poly_features_RAW",
driver="ESRI Shapefile", overwrite_layer = TRUE)
# save as CSV (without coords)
write.csv(MMU3@data, paste0(tab.dir,paste0('MMU_s123_poly_features_RAW.csv')),
row.names = FALSE)
# add coordinates to data.frame
MMU3xy <- as.data.frame(as(as(MMU3, "SpatialLinesDataFrame"),"SpatialPointsDataFrame"))
# save as CSV (with coords)
write.csv(MMU3xy, paste0(tab.dir,paste0('MMU_s123_poly_features_RAW_XY.csv')),
row.names = FALSE)
Save only the data (no spatial features).
# remove sites again
features2 <- features[!((features$S2_suit_pc==0)&(features$S3_suit_pc==0)),]
# save
write.csv(features2, paste0(tab.dir,'MMUs123_features_TABLE_RAW.csv'),
row.names = FALSE)
Plot the percent area in each site that is suitable for each of the three behavioural states.
# DOC operations regions polygon
DOC_region_shp <- get.shp('DOC_op_regions_mainland')
# NZ polygons
NZ_poly_shp <- get.shp('NZ_polygon_dissolved')
# quick plot of regions (test)
#qtm(DOC_region_shp, fill = "RegionName", fill.pallete = "RdYlGn")
# make base layers
ops_map <- tm_shape(DOC_region_shp) +
tm_polygons(col ='grey97', #alpha = 0.5,
title = 'DOC operation regions'
#palette = "YlOrRd")
#palette = "RdYlGn"
#fill='grey99'
) +
tm_borders(col ="grey27",alpha = 0.5) +
tm_legend(show=FALSE)
NZ_map <- tm_shape(NZ_poly_shp) +
#tm_polygons("OBJECTID", fill='white') +
tm_borders("grey27")
base <- ops_map + NZ_map
# state plots
S1_pts <- base +
tm_shape(MMU3) +
tm_dots(col = "S1_suit_pc", size = 0.075,
palette = 'RdYlGn') +
tm_layout(title="S1",
title.position = c('left','top'),
title.size = 1.5) +
tm_legend(show = FALSE)
S2_pts <- base +
tm_shape(MMU3) +
tm_dots(col = "S2_suit_pc", size = 0.075,
palette = 'RdYlGn') +
tm_layout(title="S2",
title.position = c('left','top'),
title.size = 1.5) +
tm_legend(show = FALSE)
S3_pts <- base +
tm_shape(MMU3) +
tm_dots(col = "S3_suit_pc",
title="% area of\nsuitability",
size = 0.075,
palette ="RdYlGn") +
tm_layout(title="S3",
title.position = c('left','top'),
title.size = 1.5) +
tm_scale_bar(breaks = c(0,100,200), size = .75) +
tm_legend(show = TRUE)
# multi-plot
s123_suit <- tmap_arrange(S1_pts,S2_pts,S3_pts)
# save an image ("plot" mode)
tmap_save(s123_suit, filename = paste0(fig.dir,"MMU_points_S123_suitability_pc.png"),
height = 6, width=10, units = 'in', dpi=600)
# view here
s123_suit
Plot all the locations, with the DOC operation regions in the background.
# DOC operations regions polygon
DOC_region_shp <- get.shp('DOC_op_regions_mainland')
# NZ polygons
NZ_poly_shp <- get.shp('NZ_polygon_dissolved')
# make base layers
ops_map <- tm_shape(DOC_region_shp) +
tm_polygons("RegionName", alpha = 0.5,
title = 'DOC operation regions'
)
NZ_map <- tm_shape(NZ_poly_shp) +
tm_borders("grey27")
base <- ops_map + NZ_map
# plot all points
#symbol<-"\u2265"
mmu_pts <- base +
tm_shape(MMU3) +
tm_dots(col = 'black', size = 0.1) +
tm_scale_bar(breaks = c(0,100,200), size = 1) +
tm_layout(frame = FALSE) +
tm_layout(title=bquote("suitable sites\nfor \u226535 females"),
title.position = c(0,.7),
title.size = 1.5) +
tm_legend(show = TRUE,
position = c(0.7,0.1),
legend.text.size = 1,
legend.title.size = 2)
## save an image ("plot" mode)
tmap_save(mmu_pts, filename = paste0(fig.dir,"MMU_points_FINAL.png"),
height = 10, width=8, units = 'in', dpi=600)
# show here
mmu_pts
Plot the percent area in each site that is suitable for each of the three behavioural states.
# DOC operations regions polygon
DOC_region_shp <- get.shp('DOC_op_regions_mainland')
# NZ polygons
NZ_poly_shp <- get.shp('NZ_polygon_dissolved')
# make base layers
ops_map <- tm_shape(DOC_region_shp) +
tm_polygons(col ='grey97',
title = 'DOC operation regions') +
tm_borders(col ="grey27",alpha = 0.5) +
tm_legend(show=FALSE)
NZ_map <- tm_shape(NZ_poly_shp) +
#tm_polygons("OBJECTID", fill='white') +
tm_borders("grey27")
base <- ops_map + NZ_map
# plot per state
S1_cvs <- base +
tm_shape(MMU3) +
tm_dots(col = "S1_CV_me", size = 0.075,
palette = 'RdYlGn') +
tm_layout(title="S1",
title.position = c('left','top'),
title.size = 1.5) +
tm_legend(show = FALSE)
S2_cvs <- base +
tm_shape(MMU3) +
tm_dots(col = "S2_CV_me", size = 0.075,
palette = 'RdYlGn') +
tm_layout(title="S2",
title.position = c('left','top'),
title.size = 1.5) +
tm_legend(show = FALSE)
S3_cvs <- base +
tm_shape(MMU3) +
tm_dots(col = "S3_CV_me",
title="Coefficient of\nvariation (%)",
size = 0.075,
palette ="RdYlGn") +
tm_layout(title="S3",
title.position = c('left','top'),
title.size = 1.5) +
tm_scale_bar(breaks = c(0,100,200), size = .75) +
tm_legend(show = TRUE)
# multi-plot
s123_cvplt <- tmap_arrange(S1_cvs,S2_cvs,S3_cvs)
# save an image ("plot" mode)
tmap_save(s123_cvplt, filename = paste0(fig.dir,"MMU_points_S123_CV_pc.png"),
height = 6, width=10, units = 'in', dpi=600)
# view here
s123_cvplt
# save workspace to load later if needed
save.image("G:/R/features_NZ_8-7-2021.RData")
This concludes the script.
Next, we simplify the extracted features to create an integrated SDM database. See Appendix S6.