Â
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.
In this script, we finalise the integrated SDM database (iSDMdb) into an accessible format for use in management and decision-making. As mentioned before, the final data fields fall under seven categories (site identification, size, model uncertainty, restoration features, human impacts, additional suitability, locations of interest), based on the four types of assessment steps (SDM prediction, human impacts, novel preferences in novel spaces, and locations of inquiry). The raw data extractions were completed in Appendix S5.
Some main tasks for making the database more accessible include:
The final products from this script are as follows:
A polygon shapefile and CSV data table of sites with 34 data fields summarising the assessment outputs.
An interactive HTML map showing suitable site locations, areas of human impacts, and a menu of these 34 data field results for each suitable site.
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') # visualising maps
Change raster and tmap options.
# change raster memory options
rasterOptions(maxmemory = 2.5e10)
# change tmap plotting options
tmap_options(max.raster = c(plot = 879305280, view = 879305280))
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('#DDCC77','#CC3311')
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)
}
Suitable site polygons (minimum 35 females) from Appendix S5. This is the file named MMU_s123_poly_features_RAW.shp.
# load and ensure projections
sites <- get.shp('MMU_s123_poly_features_RAW', polydir = lay.dir)
## OGR data source with driver: ESRI Shapefile
## Source: "G:\R\NZSL_MSSDM_MCA_2019_NZ\layers", layer: "MMU_s123_poly_features_RAW"
## with 395 features
## It has 57 fields
# show column names
colnames(sites@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"
As many of the categorical features extracted were in numeric format, the categories have to be converted back to text. The original fields will not be included in the final outputs, as the values are not easily interpretable (see main text of current publication for more information). In most cases, numeric values are sorted by the alphabetical name or the numerical order of the ID column for the original dataset before it was rasterised. They also correspond with presence or absence of features. We recommend inspecting the data before back-converting to ensure proper matching, as we had done while making this script.
Numeric values are sorted by the FID of the original polygon:
1 - Northland Region; 2 - Auckland Region; 3 - Waikato Region; 4 - Bay of Plenty Region; 5 - Gisborne Region; 6 - Hawke’s Bay Region; 7 - Taranaki Region; 8 - Manawatu-Wanganui Region; 9 - Wellington Region; 10 - West Coast Region; 11 - Canterbury Region; 12 - Otago Region; 13 - Southland Region; 14 - Tasman Region; 15 - Nelson Region; 16 - Marlborough Region; 17 - Area Outside Region
# Change names
sites$region_nm <- mapvalues(sites$region,
from=c(1:17),
to=c('Northland Region',
'Auckland Region',
'Waikato Region',
'Bay of Plenty Region',
'Gisborne Region',
"Hawke's Bay Region",
'Taranaki Region',
'Manawatu-Wanganui Region',
'Wellington Region',
'West Coast Region',
'Canterbury Region',
'Otago Region',
'Southland Region',
'Tasman Region',
'Nelson Region',
'Marlborough Region',
'Area Outside Region'))
# get counts of sites per NZ region
count(sites$region_nm)
Numeric values are sorted by the FID of the original polygon:
1 - Eastern South Island; 2 - Central North Island; 3 - Northern South Island; 4 - Western South Island; 5 - Hauraki-Waikato-Taranaki; 6 - Lower North Island; 7 - Northern North Island; 8 - Southern South Island
# Change names
sites$DOC_reg_nm <- mapvalues(sites$DOC_region,
from=c(1:8),
to=c('Eastern South Island',
'Central North Island',
'Northern South Island',
'Western South Island',
'Hauraki-Waikato-Taranaki',
'Lower North Island',
'Northern North Island',
'Southern South Island'))
# get counts of sites per DOC region
count(sites$DOC_reg_nm)
Numeric values are sorted in alphabetical order of the variable names:
0 - cliff edges; 1 - coast distance; 2 - forest distance; 3 - grass distance; 4 - land cover; 5 - sand distance; 6 - slope; 7 - water distance
# Change names
sites$MOD_nm <- mapvalues(sites$MOD_md,
from=c(0:7),
to=c('cliff edges', #not present; not limiting per site
'coast distance', #not present; not limiting per site
'forest distance',
'grass distance',
'land cover',
'sand distance',
'slope',
'inland water distance'))
# get frequency of most limiting factors
count(sites$MOD_nm)
Repeat for S1, S2 and S3 limiting factors.
# make list of names to change
to_list <- c('cliff edges', #not present; not limiting per site
'coast distance', #not present; not limiting per site
'forest distance',
'grass distance',
'land cover',
'sand distance',
'slope',
'inland water distance')
# Change names
sites$S1_lim_nm <- mapvalues(sites$S1_lim_md, from=c(0:7), to=to_list)
sites$S2_lim_nm <- mapvalues(sites$S2_lim_md, from=c(0:7), to=to_list)
sites$S3_lim_nm <- mapvalues(sites$S3_lim_md, from=c(0:7), to=to_list)
# get frequency of most limiting factors
count(sites$S1_lim_nm); count(sites$S2_lim_nm); count(sites$S3_lim_nm)
1 - present; 0 - absent
# Change names
sites$fences <- mapvalues(sites$fence_pa,
from=c(0:1),
to=c('absent','present'))
# check
(sites$fence_pa)[200:215]; (sites$fences)[200:215]
## [1] 1 1 1 1 1 0 0 1 1 0 0 1 0 0 1 1
## [1] "present" "present" "present" "present" "present" "absent" "absent"
## [8] "present" "present" "absent" "absent" "present" "absent" "absent"
## [15] "present" "present"
# get summary
summary(as.factor(sites$fences))
## absent present
## 204 191
Overwriting and converting all area and distance features.
# area features
# divide by 1,000,000m
sites$area_km <- sites$area/1e6
sites$DOC_size <- sites$DOC_size/1e6
sites$grz_size <- sites$grz_size/1e6
# distance features
# divide by 1000m
sites$rd_sl_me <- sites$rd_sl_me/1e3
sites$rd_unsl_me <- sites$rd_unsl_me/1e3
sites$rd_sl_mi <- sites$rd_sl_mi/1e3
sites$rd_unsl_mi <- sites$rd_unsl_mi/1e3
sites$rd_sl_mx <- sites$rd_sl_mx/1e3
sites$rd_unsl_mx <- sites$rd_unsl_mx/1e3
sites$in_watr_me <- sites$in_watr_me/1e3
sites$in_watr_mi <- sites$in_watr_mi/1e3
sites$in_watr_mx <- sites$in_watr_mx/1e3
# get example summary
summary(sites$area_km)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3538 0.5678 1.0069 2.0208 2.1053 49.6362
# DOC percentages
sites$DOC_perc <- sites$DOC_perc*100
# suitability percentages
sites$S1_suit_pc <- sites$S1_suit_pc*100
sites$S2_suit_pc <- sites$S2_suit_pc*100
sites$S3_suit_pc <- sites$S3_suit_pc*100
# human impact percentages
sites$hum_im_pc <- sites$hum_im_pc*100
sites$grz_perc <- sites$grz_perc*100
# pine percentages
sites$pine_perc <- sites$pine_perc*100
Note that coefficient of variation was already calculated as a percentage in previous steps (Appendix S3).
# ignore any character strings to round all numbers in dataframe
sites@data <- sites@data %>%
mutate_if(is.numeric, round, digits=2)
We categorised MESS grid values into five extrapolation classes:
We then reclassified based on the mean MESS value.
# copy mean values to to new column
sites$MESS_class <- sites$MESS_me
# reclassify numerically
sites$MESS_class[sites$MESS_class >= 0] <- 0
sites$MESS_class[(sites$MESS_class >= -100) & (sites$MESS_class < 0)] <- 1
sites$MESS_class[(sites$MESS_class >= -500) & (sites$MESS_class < -100)] <- 2
sites$MESS_class[(sites$MESS_class >= -1000) & (sites$MESS_class < -500)] <- 3
sites$MESS_class[(sites$MESS_class >= -1600) & (sites$MESS_class < -1000)] <- 4
# change to words
sites$MESS_class <- mapvalues(sites$MESS_class,
from=c(0:4),
to=c('none','low','intermediate','strong','very strong'))
The class, very strong, was not present in this dataset.
In order to assist in differentiating or querying among the sites, we will assign unique identification codes to them. These will alphanumeric, using the first letter of each word in the region name (initials), followed by the sorted number of the site.
Get capital letters from each region name to use in the coding system.
# make a letter field for the region names
sites$DOC_reg_ID <- gsub("[^::A-Z::]","", sites$DOC_reg_nm)
# show names (there should be 8 total)
unique(sites$DOC_reg_ID)
## [1] "NNI" "HWT" "CNI" "LNI" "NSI" "WSI" "ESI" "SSI"
# check
(sites$DOC_region)[200:215];(sites$DOC_reg_nm)[200:215];(sites$DOC_reg_ID)[200:215]
## [1] 2 2 2 2 2 2 6 6 6 6 6 6 6 6 6 6
## [1] "Central North Island" "Central North Island" "Central North Island"
## [4] "Central North Island" "Central North Island" "Central North Island"
## [7] "Lower North Island" "Lower North Island" "Lower North Island"
## [10] "Lower North Island" "Lower North Island" "Lower North Island"
## [13] "Lower North Island" "Lower North Island" "Lower North Island"
## [16] "Lower North Island"
## [1] "CNI" "CNI" "CNI" "CNI" "CNI" "CNI" "LNI" "LNI" "LNI" "LNI" "LNI" "LNI"
## [13] "LNI" "LNI" "LNI" "LNI"
Also add a new field to sort by main island.
# rename by name pattern
sites$main_isld <- gsub(".*SI","South Island", sites$DOC_reg_ID)
sites$main_isld <- gsub("^HWT","North Island", sites$main_isld)
sites$main_isld <- gsub(".*NI","North Island", sites$main_isld)
Create unique ID numbers
# make a unique ID number field (site_ID); has to be done as dataframe
a <- ddply(sites@data, .(DOC_reg_nm), mutate,
site_ID = paste0(DOC_reg_ID,'-',seq_along(id)))
# make a copy of sites and then join dataframe by original ID
sites2 <- sites[1]
sites2@data <- join(sites2@data, a, by='id')
# show
(sites2$DOC_region)[200:215];(sites2$site_ID)[200:215]; (sites2$DOC_reg_nm)[200:215]
## [1] 2 2 2 2 2 2 6 6 6 6 6 6 6 6 6 6
## [1] "CNI-35" "CNI-36" "CNI-37" "CNI-38" "CNI-39" "CNI-40" "LNI-21" "LNI-22"
## [9] "LNI-23" "LNI-24" "LNI-25" "LNI-26" "LNI-27" "LNI-28" "LNI-29" "LNI-30"
## [1] "Central North Island" "Central North Island" "Central North Island"
## [4] "Central North Island" "Central North Island" "Central North Island"
## [7] "Lower North Island" "Lower North Island" "Lower North Island"
## [10] "Lower North Island" "Lower North Island" "Lower North Island"
## [13] "Lower North Island" "Lower North Island" "Lower North Island"
## [16] "Lower North Island"
We reorganise the sites and add another ID field. This time, the ID is only numeric.
# sort
sites2 <- sites2[order(sites2$main_isld,sites2$DOC_reg_ID),]
Overwrite ‘id’ field with newly-organised numbers (this will link to the FID field when opened in e.g. ArcGIS).
# change 'id' row name to 'RAW_ID' if needed
sites2$RAW_ID <- sites2$id
# make new 'id' field
sites2$id <- 1:nrow(sites2)
Add coordinates of the centrepoints of the site polygons.
# get XY coordinates of each
sites_pt <- coordinates(sites2)
# extract coordinates
sites_pt <- data.frame(id=c(1:395), X=sites_pt[,1], Y=sites_pt[,2])
# add to sites@data
sites2@data <- join(sites2@data, sites_pt, by='id')
# save as final polygon (for analysis)
shapefile(sites2, paste0(lay.dir,'\\NZSL_integrated_SDM_database_CLEAN_RAW.shp'),
overwrite=TRUE)
# save as CSV (without coords)
write.csv(sites2@data,
paste0(tab.dir,paste0('NZSL_integrated_SDM_database_CLEAN_RAW.csv')),
row.names = FALSE)
# add coordinates to data.frame
sitesxy <- as.data.frame(as(as(sites2,
"SpatialLinesDataFrame"),"SpatialPointsDataFrame"))
# save as CSV (with coords)
write.csv(sitesxy,
paste0(tab.dir,paste0('NZSL_integrated_SDM_database_CLEAN_RAW_XY.csv')),
row.names = FALSE)
# get column names
names(sites2)
## [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" "region_nm" "DOC_reg_nm" "MOD_nm"
## [61] "S1_lim_nm" "S2_lim_nm" "S3_lim_nm" "fences" "area_km"
## [66] "MESS_class" "DOC_reg_ID" "main_isld" "site_ID" "RAW_ID"
## [71] "X" "Y"
There still remain the raw data fields here. The following columns will be selected for the final iSDMdb database (see Table 1 for data field descriptions and interpretations):
Site identification features:
Size features:
Model uncertainty:
Restoration features:
Human impact features:
Additional suitability features:
Locations of inquiry features:
# copy polygon
sites_output <- sites2
# select columns
sites_output@data <- sites_output@data %>%
select(#site identification features
'site_ID','main_isld','region_nm','DOC_reg_nm',
#size features
'area_km','DOC_code','DOC_name','DOC_size','DOC_perc',
#habitat suitability model features
'S1_suit_pc','S2_suit_pc','S3_suit_pc',
#restoration features
'S1_lim_nm','S2_lim_nm','S3_lim_nm',
#model uncertainty
'S1_CV_me','S2_CV_me','S3_CV_me',
'MESS_class','MOD_nm',
#human impact features
'hum_im_pc', 'rd_sl_mi', 'rd_unsl_mi',
'fences', 'grz_perc',
#additional suitability and locations of interest
'in_watr_mi','in_watr_me','in_watr_mx',
'pine_perc','curr_name','hist_name',
#coordinate data
'X','Y',
'id' # this is only at the end for HTML purposes
)
# check
sites_output
## class : SpatialPolygonsDataFrame
## features : 395
## extent : 1167775, 2086175, 4758650, 6191975 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs
## variables : 34
## names : site_ID, main_isld, region_nm, DOC_reg_nm, area_km, DOC_code, DOC_name, DOC_size, DOC_perc, S1_suit_pc, S2_suit_pc, S3_suit_pc, S1_lim_nm, S2_lim_nm, S3_lim_nm, ...
## min values : CNI-1, North Island, Auckland Region, Central North Island, 0.35, 00001, Abel Tasman National Park; Abel Tasman Scenic Reserve; Local Purpose Reserve - Cemetery - Awaroa Inlet, 0, 0, 0.1, 0, 0.09, grass distance, grass distance, coast distance, ...
## max values : WSI-9, South Island, West Coast Region, Western South Island, 49.64, W15111, Whitiau Scenic Reserve; Whitiau Scientific Reserve; Kaitoke Marginal Strip, 8.77, 100, 100, 100, 100, sand distance, sand distance, slope, ...
Rename columns based on table above (we are doing all of them, to ensure the changes).
# get column names
#names(sites_output)
# rename columns
names(sites_output) <- c(#name features
'site_ID','main_isld','region','DOC_region',
#descriptive features
'area','DOC_code','DOC_name','DOC_size','DOC_pc',
#habitat suitability model features
'S1_area_pc','S2_area_pc','S3_area_pc',
#restoration features
'S1_limits','S2_limits','S3_limits',
#model uncertainty
'S1_uncrt','S2_uncrt','S3_uncrt',
'MESS_class','dissim_var',
#human impact features
'hum_im_pc','rd_sl_mi','rd_unsl_mi',
'fences', 'graze_pc',
#additional suitability and locations of interest
'in_watr_mi','in_watr_me','in_watr_mx',
'pine_pc','curr_NZSL','histr_NZSL',
#coordinate data
'X','Y',
'id' # this is only at the end for HTML purposes
)
# show new column names
names(sites_output)
## [1] "site_ID" "main_isld" "region" "DOC_region" "area"
## [6] "DOC_code" "DOC_name" "DOC_size" "DOC_pc" "S1_area_pc"
## [11] "S2_area_pc" "S3_area_pc" "S1_limits" "S2_limits" "S3_limits"
## [16] "S1_uncrt" "S2_uncrt" "S3_uncrt" "MESS_class" "dissim_var"
## [21] "hum_im_pc" "rd_sl_mi" "rd_unsl_mi" "fences" "graze_pc"
## [26] "in_watr_mi" "in_watr_me" "in_watr_mx" "pine_pc" "curr_NZSL"
## [31] "histr_NZSL" "X" "Y" "id"
# save as final polygon (for analysis)
shapefile(sites_output,
paste0(lay.dir,'\\NZSL_integrated_SDM_database.shp'),
overwrite=TRUE)
# save as CSV (without coords)
write.csv(sites_output@data,
paste0(tab.dir,
paste0('NZSL_integrated_SDM_database.csv')),
row.names = FALSE)
# add coordinates to data.frame
sitesxy <- as.data.frame(as(as(sites_output,
"SpatialLinesDataFrame"),"SpatialPointsDataFrame"))
# save as CSV (with coords)
write.csv(sitesxy, paste0(tab.dir,
paste0('NZSL_integrated_SDM_database_XY.csv')),
row.names = FALSE)
Load base map items.
# DOC operations regions polygon
DOC_region <- get.shp('DOC_op_regions_mainland')
# NZ polygons
NZ_polygon <- get.shp('NZ_polygon_dissolved')
Get NZSL location data (historic and current).
# get shapefiles
NZSL_current <- get.shp('current_NZSL_locations_9-20-2019',polydir = int.dir)
NZSL_historic <- get.shp('historic_NZSL_locations',polydir = int.dir)
# keep name columns only
NZSL_current <- (NZSL_current)[5]
NZSL_historic <- (NZSL_historic)[5]
# add column for location type
NZSL_current$loc_type <- '1994-2019 (Dec - Mar)'
NZSL_historic$loc_type <- 'Historic'
# combine
NZSL_locations <- rbind(NZSL_current,NZSL_historic)
# human impact raster
# prep - get raster and reclass 1s to NA's (1 = no impact)
hum_ras <- get.ras('MCA_thresh.tif')
hum_ras <- reclassify(hum_ras, rcl=cbind(1,NA))
writeRaster(hum_ras,paste0(int.dir,"\\MCDA_0s_only.tif"),options="COMPRESS=LZW",
overwrite=TRUE)
hum_ras <- get.ras("MCDA_0s_only.tif", ras_dir=int.dir)
# change from 25x25m to 100x100m resolution
hum_ras_100 <- aggregate(hum_ras, fact=4)
writeRaster(hum_ras_100,paste0(int.dir,"\\MCDA_0s_only_100m.tif"),options="COMPRESS=LZW",
overwrite=TRUE)
hum_ras_100 <- get.ras("MCDA_0s_only_100m.tif", ras_dir=int.dir)
# convert raster values to categorical (factors)
human_impact_areas <- ratify(hum_ras_100)
levels(human_impact_areas)[[1]]$Human_impact = c('areas of human impact')
human_impact_areas_ras <- human_impact_areas
# convert raster to polygons and save
hum_poly <- rasterToPolygons(human_impact_areas, dissolve=TRUE)
names(hum_poly) <- 'hum_imp'
hum_poly$hum_imp <- 'human impact areas'
shapefile(hum_poly,paste0(int.dir,'\\human_impact_areas.shp'),overwrite=TRUE)
# load human impact areas again
human_impact_areas <- get.shp('human_impact_areas',polydir = int.dir)
## OGR data source with driver: ESRI Shapefile
## Source: "G:\R\NZSL_MSSDM_MCA_2019_NZ\data\intermediate", layer: "human_impact_areas"
## with 1 features
## It has 1 fields
# change column name
human_impact_areas$hum_imp <- 'human impact areas'
Initiate tmap’s view mode, which is used to generate an interactive HTML map.
# change tmap to view mode
tmap_mode("view")
Prepare the base map, which will comprise of DOC operations regions, a map of New Zealand, and the human impact areas from the MCDA.
# DOC operations regions polygon
DOC_region$OBJECTID <- DOC_region$RegionName
DOC_region$Region <- DOC_region$OBJECTID
# make base layers
ops_map <- tm_shape(DOC_region) +
tm_polygons(col='Region', alpha = 0.2,
palette='white',
title='', legend.show=FALSE)
NZ_map <- tm_shape(NZ_polygon) +
tm_borders("grey27")
hum_map <- tm_shape(human_impact_areas) +
tm_fill(col='hum_imp',
palette='#EE6677',alpha = 0.7,
title = '') +
tm_legend(show=TRUE)
# make suitable sites polygon layers
suitable_sites <- sites_output
suit_map <- tm_shape(suitable_sites) +
tm_polygons(border.col='#000000', lwd=3, alpha = 0, border.alpha = 1,
title = 'suitable sites\nfor \u226535 females') +
tm_legend(show=TRUE)
# add points for current sighting locations
nzsl_pts <- tm_shape(NZSL_locations) +
tm_dots(col = 'loc_type', size = 0.1,
palette=c('#33BBEE','#0077BB'),
title = 'known female/pup locations')+ #legend titles at layer
tm_legend(show = TRUE)
Compile the layers to make the HTML map.
# legend for suitable sites
suit_leg <- tm_add_legend(type='fill',
col='black',
labels = 'suitable sites\nfor \u2265 35 females'
)
# view map
html_map <- ops_map +
NZ_map +
hum_map +
suit_map + suit_leg +
nzsl_pts +
tm_scale_bar(breaks = c(0,100,200), text.size = 1) +
tm_layout(frame = FALSE) +
tm_layout(title =
"Predicted suitable sites for NZ sea lion breeding colonies on mainland NZ",
title.position = c("left", "top"),
title.size = 1.5)
# view map in R-studio
#html_map
# save as stand-alone HTML file (in "view" mode)
tmap_save(html_map,
filename = paste0(dir,
"\\NZSL_integrated_SDM_database.html"),
selfcontained=FALSE)
As this is in HTML format, a screenshot of the map is shown here:
And a zoom-in of the map is shown here:
In this format, end-users can click on each suitable site and manually explore the 34 data field values/summaries for that site. Using the menu on the upper left, they can select open access world mapping layers (from ESRI and OpenStreetMap) for reference, and add or remove layers that we inputted here (corresponding to the right-hand side).
# save workspace to load later if needed
save.image("G:/R/final_database_NZ.RData")
This is the end of the script and the creation of the iSDMdb. See Appendix S7 for a CSV of the integrated SDM database and its data, and Appendix S8 for an example of how the predicted sites can be further evaluated.