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

We summarised the outputs from the SDM prediction, human impacts, novel preferences in novel spaces and locations of inquiry assessments into seven categories:

  • site identification
  • size
  • model uncertainty
  • restoration features
  • human impacts
  • additional suitability
  • locations of interest

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:

  1. Polygon and point shapefiles of sites with 57 data fields summarising the assessment outputs in raw format.

  2. Data table (CSV) with 57 data fields summarising the assessment outputs in raw format.

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')         # visualizing maps

2.1.1 Options

# 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

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('#88CCEE','#CC3311')
  
# colours for variables
  var_cols <-c('#88CCEE','#44AA99','#117733','#999933','#DDCC77','#CC6677',
               '#882255','#AA4499')

2.4 Functions

2.4.1 Calculate mode

Calc.mode()

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

2.4.2 Get variable names

get.name()

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

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

2.4.4 Extract polygon data from polygons

get.text()

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

get.text.more()

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

2.4.5 Extract polygon data from rasters

get.mmm()

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

get.mode()

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

get.pa()

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

get.perc0()

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

get.perc1()

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

get.percF()

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

3. Load data

3.1 Suitable site polygons (minimum 35 females)

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

3.2 Shapefiles for feature extractions

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

3.3 Rasters for feature extractions

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

4. Preview of feature datasets

4.1 SDM Prediction results

4.1.1 SDM suitability scores

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

4.1.2 Suitability thresholds

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

4.1.3 Mode of limiting factors across runs

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

4.1.4 MESS grid

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

4.1.5 MoD grid

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

4.1.6 Coefficient of variation

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

4.2 Human impact assessment results

4.2.1 MCDA of impacts from roads and residential areas

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

4.2.2 Road distances

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

4.2.3 Presence/absence of fences

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

4.2.4 Grazing grasslands

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

4.3 Novel preferences in novel spaces data

4.3.1 Commercial planted pine forests

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

4.3.2 Euclidean distance from inland water bodies

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

4.4. Locations of inquiry data

4.4.1 New Zealand mainland region names

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

4.4.2 DOC operation regions

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

4.4.3 DOC conservation areas

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

4.4.4 Known historic pupping sites

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

4.4.5 Known current NZSL female/pup locations in the breeding season

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

5. Raster data prep for feature extractions

5.1 Clip rasters to MMU polygon shapes

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)

6. Data feature extractions to suitable site polygons (MMU for >=35 females)

6.1 SDM prediction data

6.1.1 Habitat suitability scores

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)

6.1.2 Percent area of suitable sites by threshold per state

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

6.1.3 Coefficient of variation

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)

6.1.4 Mode of limiting factors across runs

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)

6.1.5 MESS grid

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)

6.1.6 MoD grid

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)

6.2 Human impact data

6.2.1 Proportion of areas impacted by roads and residential areas (MCDA)

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)

6.2.2 Road distances

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)

6.2.3 Presence/absence of fences

# rename to abbreviate
  names(fences) <- 'fence'
  fence <- fences

# calc presence absence
  fnc_pa <- get.pa(fence)

6.2.4 Area of grazing grasslands (dairy and non-dairy) within sites

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)

6.3 Novel preference in novel space features

6.3.1 Commercial planted pine forest coverage

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)

6.3.2 Inland water body distances

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)

6.4 Locations of interest features

6.4.1 New Zealand region names

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

6.4.2 DOC operation regions

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

6.4.3 DOC conservation areas

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

Handling duplicate warnings

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)

6.4.4 Sites near historic pupping locations

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)

6.4.5 Current (1990-2019) breeding locations

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)

7. Table cleanup and appending all features to shapefiles

7.1 Cleanup and append to make one whole table

# 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

7.2 Join table to the MMU polygon layer

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

7.3 Remove polygons of S1 suitability only (finalising suitable sites)

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.

7.4 Save as shapefiles and CSVs.

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

8. Example plots of the data summaries

8.1 Percent suitability for each state

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

8.2 Suitable Sites across DOC operation regions

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

8.3 Percent uncertainty for each state

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

9. Save workspace

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