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.

1. Methods summary

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:

  • assigning names to numbered categorical data (from extracting from rasters).
  • converting proportions into percentages.
  • converting distance or area units from m or m to km or km, respectively.
  • classifying uncertainty in predictions (low or high).
  • classifying raw MESS values into extrapolation classes (from none to high).
  • creating an interactive HTML map.

The final products from this script are as follows:

  1. A polygon shapefile and CSV data table of sites with 34 data fields summarising the assessment outputs.

  2. 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.

2. R Setup

The script presented here was done using R (version 4.0.5; R Core Team 2021) and its packages.

2.1 Libraries

# 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

2.1.1 Options

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)) 

2.2 Directories

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\\')

2.3 Colours

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')

2.4 Functions

2.4.1 Load spatial data and ensure projection

get.shp()

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)
  }

get.ras()

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)
  }

3. Load data

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"

4. Reclassifying categorical features

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.

4.1 NZ region names

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)

4.2 DOC operation regions

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)

4.3 MOD and limiting factor grid values

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)

4.4 Presence of Fences

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

5. Unit conversions and rounding values

5.1 Convert from m to km

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

5.2 Convert proportions to percentages

# 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).

5.3 Round all numbers to two decimal places

# ignore any character strings to round all numbers in dataframe
  sites@data <- sites@data %>% 
                  mutate_if(is.numeric, round, digits=2)

6. Classifying impact of extrapolation (MESS grid values) at a site

We categorised MESS grid values into five extrapolation classes:

  • -1600 to -1000: very strong
  • -1000 to -500: strong
  • -500 to -100: intermediate
  • -100 to 0: low
  • 0 to 100: none

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.

7. Assigning unique ID’s to each site by DOC operations region

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"

8. Sorting sites by island, region and ID

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)

9. XY centrepoint coordinates

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')

10. Save raw dataset for use in the analysis and summaries

# 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)  

11. Final database for management/decision-making

# 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:

  • id
  • site_ID
  • DOC_reg_nm (will rename to DOC_region)
  • region_nm (will rename to region)
  • main_isld
  • X
  • Y

Size features:

  • area_km (will rename to area)
  • S1_suit_pc (will rename to S1_area_pc)
  • S2_suit_pc (will rename to S2_area_pc)
  • S3_suit_pc (will rename to S3_area_pc)

Model uncertainty:

  • S1_CV_me (will rename to S1_uncrt)
  • S2_CV_me (will rename to S2_uncrt)
  • S3_CV_me (will rename to S3_uncrt)
  • MESS_class
  • MOD_nm (will rename to dissim_var)

Restoration features:

  • S1_lim_nm (will rename to S1_limit)
  • S2_lim_nm (will rename to S2_limit)
  • S3_lim_nm (will rename to S3_limit)

Human impact features:

  • hum_im_pc
  • rd_sl_mi
  • rd_unsl_mi
  • fences
  • grz_perc (will rename to graze_pc)

Additional suitability features:

  • in_watr_mi
  • in_watr_me
  • in_watr_mx
  • pine_perc (will rename to pine_pc)

Locations of inquiry features:

  • curr_name (will rename to curr_NZSL)
  • hist_name (will rename to histr_NZSL)
  • DOC_code
  • DOC_name
  • DOC_size
  • DOC_perc (will rename to DOC_pc)

11.1 Select columns

# 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"

11.2 Save as final output layer

# 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)

12. Interactive database map (HTML) for managers/decision-makers

12.1 Load data to visualise

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)

12.2 Human impact raster

# 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'

12.3 HTML map

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).

13. Save workspace

# 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.