---
title: "Appendix S5. Integrated SDM database (raw data extractions)"
author: 'Corresponding author: Veronica F. Frans (e-mail: verofrans[-at-]gmail.com)'
date: "September 8, 2021"
output:
  html_document:
    keep_tex: yes
    toc: yes
    toc_depth: 4
    toc_float:
      collapsed: false
    df_print: paged
  pdf_document:
    keep_tex: yes
    toc: yes
    toc_depth: 4
    fig_caption: no
header-includes:
- \usepackage{pdflscape}
- \newcommand{\blandscape}{\begin{landscape}}
- \newcommand{\elandscape}{\end{landscape}}
- \usepackage{fancyhdr}
- \pagestyle{fancy}
- \fancyhead[CO,CE]{Supporting Information}
- \fancyhead[LO,LE]{Frans et al. 2021 MEE}
- \fancyhead[RO,RE]{Integrated SDM Database}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = TRUE, cache.comments = TRUE,
                      warning = FALSE, message = FALSE)
```

&nbsp;

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

```{r,results='hide',warning=FALSE,message=FALSE}    
# 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

```{r}
# modify max memory from 10gb to 25gb (this is a 32GB computer)
rasterOptions(maxmemory = 2.5e10)
rasterOptions()
```

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

```{r, warning=FALSE,results='hide'}

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

```{r}
# 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. 

```{r}
# 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. 

```{r}
# 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).

```{r}
# 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).

```{r}
# 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.

```{r}
# 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.

```{r}
# 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.

```{r}
# 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).

```{r}
# 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).

```{r}
# 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.

```{r}
# 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.

```{r}

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

```{r}

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

```{r, results='hide'}
# 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.

```{r}
# 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.

```{r}
# show fields
  MMU
```

## 3.2 Shapefiles for feature extractions

```{r, message=FALSE, results='hide', warning=FALSE}
# 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

```{r, message=FALSE, results='hide', warning=FALSE}
# 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

```{r, message=FALSE, results='hide', warning=FALSE}
# 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

```{r, message=FALSE, results='hide', warning=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE, fig.width=9, fig.height=4}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r coeff_plot,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}

# 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

```{r,cache=TRUE,cache.comments=FALSE}

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

```{r,cache=TRUE,cache.comments=FALSE}
# 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.

```{r,cache=TRUE,cache.comments=FALSE}
# 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.

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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.

```{r cv_mmm,cache=TRUE,cache.comments=FALSE}
# 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.

```{r,cache=TRUE,cache.comments=FALSE}
# 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.

```{r,cache=TRUE,cache.comments=FALSE}
# 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. 

```{r,cache=TRUE,cache.comments=FALSE}
# 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.

```{r,cache=TRUE,cache.comments=FALSE}
# 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).

```{r,cache=TRUE,cache.comments=FALSE}
# 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.

```{r,cache=TRUE,cache.comments=FALSE}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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.

```{r,cache=TRUE,cache.comments=FALSE}
# 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)
```

```{r, cache=TRUE}
# see first few rows
  head(grz_names); tail(grz_names)
  
# show fields with NA
  #grz_names@data[!complete.cases(grz_names@data),]
  
```

```{r}
# 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.

```{r,cache=TRUE,cache.comments=FALSE}
# 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)
```

```{r, cache=TRUE}
# see first few rows
  head(pine_names); tail(pine_names)
  
# show fields with NA
  #pine_names@data[complete.cases(pine_names@data),]
```

```{r}
# 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.

```{r,cache=TRUE,cache.comments=FALSE}
# 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.

```{r}
# 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).

```{r,cache=TRUE,cache.comments=FALSE}
# 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()`.

```{r}
# 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).

```{r,cache=TRUE,cache.comments=FALSE}
# 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.

```{r,cache=TRUE,cache.comments=FALSE}
# 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.

```{r}
# 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)
```

```{r}
# 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()`.

```{r}
# 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.

```{r}

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

```{r}
# see first few rows
  head(hist_names); tail(hist_names)
```

```{r}
# show completed fields
  hist_tab
```

Save.

```{r}
# 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.

```{r}
# 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

```{r}
# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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))
```

```{r summ_feat,cache=TRUE,cache.comments=FALSE}
# view here
  summary(features)
```

## 7.2 Join table to the MMU polygon layer

```{r,cache=TRUE,cache.comments=FALSE}
# join table polygon layer
  MMU2 <- MMU
  MMU2@data <- join(MMU2@data,features,by='id')
```

```{r mmu_cols}
# show column names
  colnames(MMU2@data)
```

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

```{r}
# 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)
```

30 sites will be removed.

```{r}
# Remove these polygons
  MMU3 <- MMU2[!((MMU2$S2_suit_pc==0)&(MMU2$S3_suit_pc==0)),]
```

This is the final polygon.

```{r mmu3_details}
# final polygon details
  MMU3
```

There are 395 sites as the final output.

## 7.4 Save as shapefiles and CSVs.

```{r}
# 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).

```{r}
# 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.

```{r, fig.height=6, fig.width=8}
# 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.

```{r, fig.height=10, fig.width=8}
# 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.

```{r uncert_plt, fig.height=6, fig.width=8}

# 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

```{r,cache=TRUE,cache.comments=FALSE}
# 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**.
