---
title: 'Appendix S2. Multi-state species distribution model prediction and suitable sites'
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)
```

&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 used the New Zealand sea lion (*Phocarctos hookeri*; NZSL) Maxent model training background information data (*.lambda* files) from Frans et al. (2018) and predicted them on mainland New Zealand and Stewart Island (**see Appendix S1** for demonstration of the model training procedure). In order to determine suitable sites for breeding colonies of a minimum of 35 females, we used the multi-state species distribution model (SDM) framework (Frans et al. 2018). Information on the original model training and the multi-state SDM framework are provided in Frans et al. (2018) and its supplementary materials. We also demonstrate it in **Appendix S1** of this current publication (Frans et al. 2021).

The final products of the multi-state SDM are:

(1) a raster map indicating the suitability (or preference) of each 25 X 25 m pixel for each of the species' three behavioural states of inland movement during the breeding season (breeding--S1; transition--S2; dispersion--S3);

(2) a polygon shapefile, highlighting suitable locations of a minimum area size for $\geq$ 35 females.

The resulting polygon shapefile is also essential for creating the basic, simple, and spatial structure of the *iSDMdb*. Note that the *iSDMdb* is compatible with other SDM algorithms and/or other methods to create polygons of suitable sites.

# 2. R Setup

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

## 2.1 Libraries

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

  library("raster")       # raster data
  library("rgdal")        # input/output; projections; reading ASCII files
  library("maptools")     # mapping
  library("classInt")     # changing colour intervals on a plot
  library("gridExtra")    # grid graphics
  library("grid")         # creating graphical objects
  library("lattice")      # required for grid graphics
  library("RColorBrewer") # colours for graphics
  library("reshape")      # sort data
  library("reshape2")     # sort data
  library("igraph")       # clumping contiguous pixels
  library("rgeos")        # polygon and point conversions

``` 

### 2.1.1 Options

Change raster options to maximise memory.

```{r}
# view default raster options
  rasterOptions()
```

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

## 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', cache=TRUE}

# Root directory
  dir <- c('G:\\R\\NZSL_MSSDM_MCA_2019_NZ\\')

# Data directory
  dat.dir <- paste0(dir,'data\\')

# Maxent.jar file directory
  maxent.dir <- c(paste0(dir,'maxent'))
  
# Maxent cache files (training and projection environmental variables)
  proj.cache <- c(paste0(dat.dir,'maxent.cache\\projection\\'))
  train.cache <- c(paste0(dat.dir,'maxent.cache\\training\\'))
        
# Intermediate data and layers layers
  dir.create(paste0(dir,'data\\intermediate'),recursive=TRUE)
  int.dir <- paste0(dir,'data\\intermediate\\')
  
# Final layers
  dir.create(paste0(dir,'layers'),recursive=TRUE)
  lay.dir <- paste0(dir,'layers\\')
  
# 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 Custom colours

```{r}
# colours for the multiple states
  state_cols <-c('#BBBBBB','#DC050C','#882E72','#1965B0','#6195CF','#4EB265','#90C987',
                 '#F7CB45','#EE8026')
```

## 2.4 Custom functions

### 2.4.1 Load raster and ensure projection

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.2 Load points or polygons and project

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(paste0(polydir,p_name))

    # project and return point/polygon
      pol <- spTransform(pol, CRS(paste0(prj)))
      return(pol)
  }
```

### 2.4.3 Calculate minimum mapping units (sites for minimum of *n* individuals)

Defining a minimum area for *n* individuals (minimum mapping unit; MMU), where suitable habitats for each behavioural state are within *x* units of distance. This will be used to search for suitable areas of a minimum size, which will contain no unsuitable patches (contiguous suitable values within *x*). More information on MMU is provided in Frans et al. (2018).
 
We calculate the MMU for *n* individuals as:

$$MMU_S = ((a/n_1)*n_2)/r$$

where $a$ and $n_1$ represent the number of individuals ($n_1$) per area ($a$), multiplied by the minimum number of individuals defined by the user ($n_2$), and all divided by the resolution of the raster layer ($r$) in the same units as the initial area, $a$ (e.g., meters, kilometers).

For the NZSL, the MMU is calculated for 35 breeding females. The resolution is 25 m, or, 625 m\textsuperscript{2}. These will be used in the default settings.

```{r, mmu_function}
# calculate minimum area sizes
# input number (density) of individuals per given area (area_size)
  mmu.calc <- function(density, area_size,
                       # defaults
                       min_indiv=35, res=625){
    # calculate MMU
      mmu <- ((area_size/density)*min_indiv)/res
      if (mmu <= 1){
        return(1)
      } else {
        # round to the nearest number of pixels
        return(round(mmu, digits=0))
      }}
```

# 3. SDM prediction and results

We projected the training results to the entire New Zealand mainland by accessing the *'density.project'* command in the *Maxent.jar* file for Maxent v.3.3.3k (Phillips & Dudik 2008). The model was trained for 100 runs for each state, so there was a total of 300 corresponding predictions. The training run data were originally stored in folders for each state and each run under the following pattern: `~\m_output\(state number)\(run number)\` (see Frans et al. 2018 and Appendix S1 of this current publication for more details). We then calculated the mean of these predictions for each state to get three raster maps. 

First, we make the predictions on the mainland.

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

# predictions and mean
  for(state in 1:3){ 

      # output directory
        maxent.out <- paste0(dir,'m_output\\',state)
        
      # prediction
        for(run in 1:100){
            # Directory for each iteration's output
              runs.out <- paste0(maxent.out,"\\",run,"\\")
            
            # Maxent.jar command line
              {                   
               command <- paste0("java -mx1024m -cp ",maxent.dir,
                                 "\\maxent.jar density.Project ",runs.out,
                                 "*.lambdas ",proj.cache," ",runs.out,
                                 "projection_S",state,"_",run,
                                 ".asc -t cliff_edges -t landcover pictures=true")
               system(sprintf(command))
               }
        }
        
      # mean prediction per state
        
        # empty vector
          maxent_rasters <- list()
        
        # For-loop
          for (run in 1:100){
               # Directory
                 runs.out <- paste0(maxent.out,"\\",run,"\\")
                 
               # Search for ".asc" files and stack
                 filenames <- list.files(path=paste0(runs.out),
                                         pattern=paste0("projection_S",state,"_",
                                                        run,".asc"),full.names=TRUE)  
                 maxent_rasters[[run]] <- filenames
               }
          # stack
            maxent_100_rasters <- do.call(stack, maxent_rasters)
            
          # calculate mean
            mean_proj <- calc(x=maxent_100_rasters, fun=mean)
  
          # save
            writeRaster(mean_proj,paste0(dat.dir,"\\S",state,"_projection.tif"),
                        options="COMPRESS=LZW",overwrite=TRUE)  
  }

```

We then plot the results. These three mean maps will be used in the multi-state SDM framework.

```{r, cache=TRUE}

# Empty vector
  s123 <- list()

# For-loop (S1-3): 
  for(state in 1:3){
                    raster_names <- list.files(path=paste0(dat.dir),
                                       pattern=paste0("\\S",state,"_projection.tif"),
                                       full.names=TRUE)
                    # Add to list
                    s123[[state]] <- raster_names
  }

# stack
  s123_sdm <- do.call(stack,s123)

# Set up plots (S1-3)
  png(width = 2400, height = 1800,
      paste0(fig.dir,"/SDMs_S1-3.png"),res=300)
      sdmplot <- spplot(s123_sdm,
             main=c("SDM results for each state"),
             layout=c(3,1), #(ncols, nrows)
             maxpixels=500000,
             col.regions=c(rev(brewer.pal(11,"Spectral"))),
             cuts=10,
             colorkey=list(labels=list(at=seq(0,1,by=0.1)),   #set up labels
                           space="right"),
             par.settings = list(panel.background=list(col="black"))
             )
  sdmplot
  dev.off()
  
```

Plot here.

```{r}
# Display figure
  sdmplot
```

# 4. Multi-state SDM: Overall suitability across states

## 4.1 Statistical thresholds and reclassification

Read table of threshold values from previous Maxent training results (also obtained from Frans et al. 2018 training step, and shown in Appendix S1).

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

# read CSV summary table
  thresholds <- read.csv(paste0(dat.dir,"thresholds.csv"))

# Pull table to evaluate threshold values
  thresholds_table <- acast(thresholds, thresholds[,3]~thresholds[,1])

# show table
  knitr::kable(thresholds_table)

```

Use the Maximum Training Sensitivity Plus Specificity (MaxSSS) threshold to create a binary dataset of presence and absence, (or, suitability and unsuitability).

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

# Cast to make a new table, sorted by thresholds and states
  cast_thresh <- cast(thresholds, state~threshold)

# Subset MaxSSS threshold
  sub_thresh <- cast_thresh[,c(grep("Maximum.training",names(cast_thresh)))]
  sub_thresh <- data.frame(sub_thresh)
  sub_thresh

```

Reclassify the probability of presence for each state map to its respective threshold. 0 is for any values below the threshold (i.e. absence) and __*recl_val*__ is for values greater than or equal to the threshold (i.e. presence). __*recl_val*__ is a unique value for each state that increases by one digit per subsequent state in the for-loop: S1 = 1, S2 = 10, S3 = 100.

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

# For-loop to reclassify by threshold
  for (state in 1:3){         
                     # Extract threshold value
                       recl_thresh <- sub_thresh$sub_thresh[state]
                      
                     # Pull up state predictions
                       s_proj <- raster(paste0(dat.dir,"S",state,"_projection.tif"))

                     # Set unique value for each state
                       recl_val <- ((10^state)/10) 
                     
                     # Reclassify to 0 and 1
                       s_proj_recl <- calc(s_proj,
                                           fun=function(x){x[x>=recl_thresh]<-recl_val;
                                           return(x)})
                       s_proj_recl <- calc(s_proj_recl,
                                           fun=function(x){x[x<recl_thresh]<-0;
                                           return(x)})                         
                      
                     # Save 
                       writeRaster(s_proj_recl,
                                   paste0(int.dir,"MaxSSS_S",state,".tif"),
                                   options="COMPRESS=LZW",overwrite=TRUE)      
                    }

# Read new files and stack

  # Empty vector
    s123_rasters <- list()

  # For-loop (S1-3): 
    for(state in 1:3){
                      raster_names <- list.files(path=paste0(int.dir),
                                         pattern=paste0("\\MaxSSS_S",state,".tif"),
                                         full.names=TRUE)
                      
                      # Add to list
                      s123_rasters[[state]] <- raster_names
                     }
                  
  # Threshold raster stack
    thresh_123_stack <- do.call(stack,s123_rasters)
    
  # Ensure projection
    crs(thresh_123_stack) <- CRS("+init=epsg:2193")
```

View intermediate results.

```{r,warning=FALSE,cache=TRUE,cache.comments=FALSE,message=FALSE,fig.width=9, fig.height=4}

# Set up plots (S1-3)
  thresh_S1 <- spplot(thresh_123_stack$MaxSSS_S1,
                      main=c("S1 sites"),
                      width = 1, maxpixels=50000000,
                      col.regions=colorRampPalette(c('#EEDD88','#CC3311')),
                      cuts=1,
                      colorkey=list(space="bottom",tick.number=1,
                                    labels=list(labels=c("0","1"),cex=1)),
                      par.settings = list(panel.background=list(col="black"))
                     )
  
  thresh_S2 <- spplot(thresh_123_stack$MaxSSS_S2,
                      main=c("S2 sites"),
                      width = 1, maxpixels=50000000,
                      col.regions=colorRampPalette(c('#DDCC77','#CC3311')),
                      cuts=1,
                      colorkey=list(space="bottom",tick.number=1,
                                    labels=list(labels=c("0","10"),cex=1)),
                      par.settings = list(panel.background=list(col="black"))
                      )
    
  thresh_S3 <- spplot(thresh_123_stack$MaxSSS_S3,
                      main=c("S3 sites"),
                      width = 1, maxpixels=50000000,
                      col.regions=colorRampPalette(c('#DDCC77','#CC3311')),
                      cuts=1,
                      colorkey=list(space="bottom",tick.number=4,
                                    labels=list(labels=c("0","","100"),cex=1)),
                      par.settings = list(panel.background=list(col="black"))
                    )

# Save multiplot image to folder
  png(width=2400,height=1500, paste0(fig.dir,"thresholds_MaxSSS_S1-3.png"),res=300)
  print(thresh_S1, position = c(.012,0,.35,1),more = T)
  print(thresh_S2, position = c(.312,0,.65,1),more = T)
  print(thresh_S3, position = c(.612,0,.95,1))
  dev.off()
```

Plot here.

```{r, fig.width=9, fig.height=5}
# Display figure here
  print(thresh_S1, position = c(.018,0,.35,1),more = T)
  print(thresh_S2, position = c(.318,0,.65,1),more = T)
  print(thresh_S3, position = c(.618,0,.95,1))
```


## 4.2 Calculating overall suitability across states

We then generate a map of overall suitability across states by adding the threshold layers together, indicating suitability for one, two, or all three behavioural states. Pixel values will be as follows (from Frans et al. 2018):

&nbsp;

S1 | S2 | S3 | Sum | State Suitability
-- | -- | -- | --- | -----------------
0 | 0 | 0 | 0 | None
1 | 0 | 0 | 1 | S1 only
0 | 10 | 0 | 10 | S2 only
0 | 0 | 100 | 100 | S3 only
1 | 10 | 0 | 11 | S1 and S2
1 | 0 | 100 | 101 | S1 and S3
0 | 10 | 100 | 110 | S2 and S3
1 | 10 | 100 | 111 | All states

&nbsp;

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

# Calculate sum
  sum_s123 <- calc(x=thresh_123_stack, fun=sum)

# Save to file (this is a final layer)
  writeRaster(sum_s123,paste0(lay.dir,"suitablility_all_states.tif"),
              options="COMPRESS=LZW",overwrite=TRUE) 

```

View image.

```{r, warning=FALSE, cache=TRUE, cache.comments=FALSE, results='hide'}
# read raster
  sum_s123 <- get.ras("suitablility_all_states.tif", ras_dir=lay.dir)

# Use "cut()" to make 7 classifications (warning: this turns 0 into NA)
  sum_s123_cut <- cut(sum_s123, c(0,1,10,11,100,101,110,111))  

# Reclassify all values to 0 to include 0s in the calculation
  sum_s123_0s <- calc(sum_s123, fun=function(x){x[x>=0] <-0;return(x)})

# Unify rasters to create classified layer that includes values of 0
  sum_s123_stack <- stack(sum_s123_cut,sum_s123_0s)
  sum_s123_cut <- max(sum_s123_stack,na.rm=TRUE)

# Set up spplot
  sum_s123_plot <- spplot(sum_s123_cut,
                          main=c("Multi-state SDM: suitability for all states"),
                          width = 1, maxpixels=5000000,
                          cuts=7,
                          col.regions=c(state_cols),
                          colorkey=list(space="bottom",
                                        breaks=list(0,1,2,3,4,5,6,7),
                                        labels=list(labels=
                                          c("0","1","10","11","100","101","110","111"),
                                            cex = 1)),
                          par.settings = list(panel.background=list(col="black"))
                          )

# Save image
  png(width=3600,height=3600, paste0(fig.dir,"suitablility_all_states.png"),res=600)
  sum_s123_plot
  dev.off()
```

Plot here.

```{r, fig.width=9, fig.height=9}
# Display figure
  sum_s123_plot
```


This is the final output of the multi-state SDM framework. The suitability across states shows locations suitable for one, two or all three behavioural states, according to the MaxSSS threshold calculated from the SDMs.

# 5. Suitable sites for minimum of 35 females (minimum mapping units; MMU)

We next define a minimum area for 35 individuals, using the MMU function to calculate the minimum number of contiguous suitable pixels to use to extract these areas.

Note that for this species, using contiguous suitable pixels to define suitable sites was the best method. However, other methods for creating polygons can be used. For example, a moving window analysis could be used for species that can have non-contiguous suitable areas found within a reachable distance (see *range module* example in Frans et al. 2018).
 
## 5.1 Behavioural state minimum area requirements       

For the NZSL, the MMU is calculated for 35 breeding females. At S1, female NZSLs are at a density of 85 females per 100 m\textsuperscript{2}. At S2, where densities are at 30 females per 100 m\textsuperscript{2}. S3 has densities of 1 female per 10,000 m\textsuperscript{2} (Augé et al. 2009).

Calculate MMUs for each state.

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

# calculate and append
  S1 <- mmu.calc(85,100)  #85 females per 100m2 in S1
  S2 <- mmu.calc(30,100)  #30 females per 100m2 in S2
  S3 <- mmu.calc(1,10000) #1 female per 10000m2 in S3
  mmu_table <- as.data.frame(rbind(S1,S2,S3))
  colnames(mmu_table) <- "MMU"
  
# display table
  knitr::kable(mmu_table)

```

## 5.2 Aggregation of adjacent cells

Reload the MaxSSS layers and aggregate contiguous pixels for each state. Clumps greater than or equal to the size of the MMU are extracted and saved as new rasters.

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

# Create an empty vector
  clmp_layers <- list()

# For-loop
  for (state in 1:3)
         {         
          # Extract threshold raster
            recl_thresh <- subset(thresh_123_stack,paste0("MaxSSS_S",state))
            
          # Reclassify 0 to NULL (NA), save as intermediate and reload
            s_thresh_recl <- reclassify(recl_thresh, rcl=cbind(0,NA))
            writeRaster(s_thresh_recl,paste0(int.dir,"s_thresh_recl_S",state,".tif"),
                        options="COMPRESS=LZW",overwrite=TRUE) 
            s_thresh_recl <- get.ras(paste0("s_thresh_recl_S",state,".tif"), ras_dir=int.dir)
          
          # Extract unique class values in the raster
            thr_val <- unique(s_thresh_recl)
          
          # Pull up the MMU for each state
            mmu_val <- mmu_table$MMU[state]
            
          # skip if MMU is 1
            if (mmu_val > 1){
                              
          # for-loop to detect contiguous groups of suitable pixels (clumps)
          # and extract clumps of a certain cell count
            for (i in thr_val) 
                   {
                    # clump class raster
                      clmp_thresh <- clump(s_thresh_recl)
                      
                    # save as intermediate raster, write over and reread
                      writeRaster(clmp_thresh,paste0(int.dir,"clump_S",state,".tif"),
                                  options="COMPRESS=LZW",overwrite=TRUE) 
                      clmp_thresh <- get.ras(paste0("clump_S",state,".tif"), ras_dir=int.dir)
                      
                    # calculate frequency of each clump/patch
                      clmp_freq <- as.data.frame(freq(clmp_thresh))
                    
                    # store unique clump IDs where cell count is < x (MMU)
                      remove_ID <- clmp_freq$value[which(clmp_freq$count < mmu_val)]
                    
                    # assign NA to all clumps whose ID's have frequency x 
                    # (this updates the raster layer; unwanted areas deleted)
                      s_thresh_recl[clmp_thresh %in% remove_ID] <- NA
            }
            }                                
          # Save   
            writeRaster(s_thresh_recl,paste0(int.dir,"MMU_S",state,".tif"),
                        options="COMPRESS=LZW",overwrite=TRUE)                                
  }
```

Take the sum of S1-3 to get values only for areas of a minimum size for all behavioural states combined.

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

# Read new files and stack 
        
  # Empty vector
    s123_rasters <- list()
  
  # For-loop (S1-3): 
    for(state in 1:3)
         {
          raster_names <- list.files(path=paste0(int.dir),
                                     pattern=paste0("MMU_S",state,".tif"),
                                     full.names=TRUE)
          
          # Add to list
          s123_rasters[[state]] <- raster_names
         }
  
  # Threshold raster stack of partitions (S1-3)
    MMU_s123_stack <- do.call(stack,s123_rasters)

# Sum
  MMU_sum <- calc(MMU_s123_stack,fun=sum,na.rm=TRUE)
  
# save intermediate raster
  writeRaster(MMU_sum,paste0(int.dir,"MMU_123sum_w_NAs.tif"),
              options="COMPRESS=LZW",overwrite=TRUE) 

# Mask to study area by reclassifying all values to 1 and multiplying 
# (use any raster as template. Here, we used 'recl_thres' from above)
# save for use later
  rcl_df <- c(0, 100, 1)
  rcl_mtrx <- matrix(rcl_df,ncol=3,byrow = TRUE)
  all_1s <- reclassify(recl_thresh,rcl_mtrx)
  all_1s <- reclassify(all_1s, rcl=cbind(0,1)) #don't forget the 0's!
  writeRaster(all_1s,paste0(int.dir,"all_1s_mask.tif"),
              options="COMPRESS=LZW",overwrite=TRUE)
  MMU_sum2 <- prod(MMU_sum, all_1s, na.rm=TRUE)

# Save intermediate raster
  writeRaster(MMU_sum2,paste0(int.dir,"S123_MMUs_sum_all.tif"),
              options="COMPRESS=LZW",overwrite=TRUE) 
```

## 5.3 Final multi-state SDM output: Suitable sites of minimum area size

Since suitable sites of all three behavioural states combined should be of a total MMU size of 562 cells, we reuse the __*clump()*__ function to extract them. This is the final raster output for the multi-state SDM. It will only have suitable sites of this size, while all others are NULL (NA).

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

# Reclassify all values of 0 to NULL (NA) from MMU_sum
  MMU_sum_recl <- reclassify(MMU_sum, rcl=cbind(0,NA),
                             filename = tempfile(fileext = ".tif"),
                             options   = "COMPRESS=LZW")

# save intermediate raster
  writeRaster(MMU_sum_recl,paste0(int.dir,"S123_MMU_sum_recl.tif"),
              options="COMPRESS=LZW",overwrite=TRUE) 

# Add total size of MMU across all 3 states
  mmu_val <- sum(mmu_table$MMU[1:3])

# Detect contiguous groups of suitable pixels and extract by a certain cell count
 
  # clump class raster
    clmp_mmu <- clump(MMU_sum_recl,filename = tempfile(fileext=".tif"),
                             options="COMPRESS=LZW")
  
  # calculate frequency of each clump/patch
    clmp_freq <- as.data.frame(freq(clmp_mmu))
  
  # store unique clump IDs where cell count is < x (MMU)
    remove_ID <- clmp_freq$value[which(clmp_freq$count < mmu_val)]  
  
  # assign NA to all clumps whose ID's have frequency x (unwanted areas deleted)
    MMU_sum_recl[clmp_mmu %in% remove_ID] <- NA

# Reclassify NA to 0
  MMU_sum_recl <- reclassify(MMU_sum_recl, rcl=cbind(NA,0))

# Mask to study area by reclassifying all values to 1 and multiplying 
  MMU_sum_recl <- overlay(MMU_sum_recl, all_1s, fun=function(a,b){return(a*b)})
                                   
# Save and clear from workspace to save memory 
  writeRaster(MMU_sum_recl,paste0(lay.dir,"S123_MMUs_35_females.tif"),
                        options="COMPRESS=LZW",overwrite=TRUE)               
```

Plot.

```{r,warning=FALSE, cache=TRUE, cache.comments=FALSE, results='hide'}
# reload
  MMU_sum_recl <- get.ras(paste0("S123_MMUs_35_females.tif"), ras_dir=lay.dir)

# Use "cut()" to make 7 classifications (warning: this turns 0 into NA)
  MMU_s123_cut <- cut(MMU_sum_recl, c(0,1,10,11,100,101,110,111))

# Reclassify all values to 0 to include 0s in the calculation
  MMU_s123_0s <- calc(MMU_sum_recl, fun=function(x){x[x>=0] <-0;return(x)})
  MMU_s123_NAs <- calc(MMU_sum_recl, fun=function(x){x[x==NA] <-0;return(x)})

# Unify rasters to create classified layer that includes values of 0
  MMU_s123_stack <- stack(MMU_s123_cut,MMU_s123_0s)
  MMU_s123_cut2 <- max(MMU_s123_stack,na.rm=TRUE)

# spplot
  sum_MMU_plot <- spplot(MMU_s123_cut2,
                         main=c("Suitable sites for minimum of 35 females (S1-3)"),
                         width = 1,maxpixels=50000000,
                         cuts=7,
                         col.regions=c(state_cols),
                         colorkey=list(space="bottom",
                                       tick.number=7,
                                       breaks=list(0,1,2,3,4,5,6,7),
                                       labels=list(labels=c("0","1","10","11","100",
                                                            "101","110","111",""),
                                                   cex = 1)),
                        par.settings = list(panel.background=list(col="black"))
                         )

# Save
  png(width=1800,height=1800,
      paste0(fig.dir,"Multi-state_SDM_MMU_35_females.png"),res=300)
  print(sum_MMU_plot)
  dev.off()
```

Plot here. 

```{r, fig.width=9, fig.height=9}
# View graphic
  sum_MMU_plot
```

As shown here, the sites that were smaller than the MMUs were deleted. 

## 5.4 Conversion to polygons and points

### 5.4.1 Conversion to polygons

Next, we read the overall suitability raster layer again (*`S123_MMUs_35_females.tif`*) and reclassify it, where MMU clusters are valued to 1, and all other pixels are reclassified to NA (NULL).

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

# Read raster again
  MMU_sum_recl <- raster(paste0(lay.dir,"S123_MMUs_35_females.tif"))

# Reclassify, where MMU clusters are valued at 1, and all other pixels are NA
# Reclassify all values of 0 to NULL (NA) from MMU_sum
  MMU_s123_NA <- reclassify(MMU_sum_recl, rcl=cbind(0,NA),
                             filename = tempfile(fileext = ".tif"),
                             options   = "COMPRESS=LZW")
  # save intermediate raster
  writeRaster(MMU_s123_NA,paste0(int.dir,"MMU_S123.tif"),
              options="COMPRESS=LZW",overwrite=TRUE) 

  MMU_s123_NA <- calc(MMU_s123_NA, fun=function(x){x[x>0] <-1;return(x)})
  # save intermediate raster
  writeRaster(MMU_s123_NA,paste0(int.dir,"MMU_S123_all_1s.tif"),
              options="COMPRESS=LZW",overwrite=TRUE) 
```

We create a polygon that will be used to extract values from all clusters, making sure to first remove any polygons smaller than the minimum area size for all three states (562 pixels; 351,250 m\textsuperscript{2}). 

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

# Convert raster to polygon and convert from multipart to singlepart (currently, features = 1)
  MMU_poly <- rasterToPolygons(MMU_s123_NA, na.rm=TRUE, digits = 12, dissolve=TRUE)
  MMU_poly <- disaggregate(MMU_poly)

# Project polygons to NZDG2000
  proj4string(MMU_poly) <- CRS("+init=epsg:2193")

# Calculate the area of each individual polygon and round to nearest integer
  MMU_poly$area <- round(area(MMU_poly))
  
# Extract polygons >= 351250 m^2 in size.
  MMU_s123_poly_FINAL <- subset(MMU_poly, area >= 351250)

# Add a column assigning a unique ID number to each polygon 
# (later used for polygon centroids)
  MMU_s123_poly_FINAL$id <- 1:nrow(MMU_s123_poly_FINAL)

```

For the NZSL, suitable sites are areas of a minimum size that contain a location for at least S1, followed by S2 and/or S3. Hence, we prioritise values of 1, 11 and 111, and reclassify a new raster, where values 1, 11 and 111 are converted to 1, and 0 is NA.

```{r,warning=FALSE,cache=TRUE,cache.comments=FALSE}
# reclassify
  MMU_s123_1s <- calc(MMU_sum_recl, fun=function(x){x[x<=0] <-NA;return(x)})
  MMU_s123_1s <- calc(MMU_s123_1s, fun=function(x){x[x==11] <-1;return(x)})
  MMU_s123_1s <- calc(MMU_s123_1s, fun=function(x){x[x==111] <-1;return(x)})
  
# save intermediate raster
  writeRaster(MMU_s123_1s,paste0(int.dir,"MMU_S123_1_11_111.tif"),
              options="COMPRESS=LZW",overwrite=TRUE) 
```

Next, we add a column to the polygons, where the minimum value is extracted from the MMU rasters. We then extract all polygons with a minimum value of "1". This is the final polygon used that highlights suitable sites for a minimum of 35 individuals across all three behavioural states.

```{r,warning=FALSE,cache=TRUE,cache.comments=FALSE}
# Add a column where minimum MMU raster value within each polygon is extracted
  MMU_s123_poly_FINAL$min_val <- as.numeric(extract(MMU_s123_1s, MMU_s123_poly_FINAL,
                                                    fun=min, na.rm=TRUE))

# Make a subset
  MMU_poly_1s_NZ <- subset(MMU_s123_poly_FINAL, min_val == 1)

# Save
  temp_dir <- file.path(paste0("G:/R/NZSL_MSSDM_MCA_2019_NZ/layers"))
  writeOGR(obj=MMU_poly_1s_NZ, dsn=temp_dir, layer="MMU_s123_poly_NZ",
           driver="ESRI Shapefile", overwrite_layer = TRUE)  
```

Also save the intermediate layer of all polygons, `MMU_poly` (unfiltered by size), in case they're needed. 

```{r,eval=FALSE}
# Save intermediate polygon as shapefile
  temp_dir <- file.path(paste0("G:/R/NZSL_MSSDM_MCA_2019_NZ/data/intermediate"))
  writeOGR(obj=MMU_poly, dsn=temp_dir, layer="MMU_S123_poly_UNFILTERED",
           driver="ESRI Shapefile",overwrite_layer=TRUE)
```

Plot. 

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

# Upload New Zealand coastline polygon
  New_Zealand <-readShapePoly(paste0(dat.dir,"Region.shp"))

# Save image
  png(width=1600,height=1800,
      paste0(fig.dir,"Multi-state_SDM_MMU_35_females_polygons.png"),res=300)
  {plot(New_Zealand,
        col="gray96",bg="#8DCBE4",
        axes=TRUE,
        xlim=c(1000000,2005000))
  plot(MMU_poly_1s_NZ,col="#E65518",cex=.5, border="#E65518", add=TRUE)
  title("Suitable Sites for a minimum of 35 females")}
  dev.off()
  
```

Plot here.

```{r polyplot,warning=FALSE,cache=TRUE,cache.comments=FALSE,fig.height=9,fig.width=9}
# plot here
  {plot(New_Zealand,
        col="gray96",bg="#8DCBE4",
        axes=TRUE, 
        xlim=c(1000000,2005000))
  plot(MMU_poly_1s_NZ,col="#E65518",cex=.5, border="#E65518", add=TRUE)
  title("Suitable Sites for minimum of 35 females (multi-state SDM)")}
```

### 5.4.2 Conversion to Points

We convert the extracted polygons from above into polygon centroids.

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

# Polygon centroids by "id"
  MMU_points_NZ = gCentroid(MMU_poly_1s, byid=TRUE)

# Save shapefile
  shapefile(MMU_points_NZ, paste0(lay.dir,'MMU_S123_point_NZ.shp'),overwrite=TRUE)
  
```

Plot.

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

# Upload New Zealand coastline polygon
  New_Zealand <-readShapePoly(paste0(dat.dir,"Region.shp"))

# Save image
  png(width=1600,height=1800,
      paste0(fig.dir,"Multi-state_SDM_MMU_35_females_points.png"),res=300)
  {plot(New_Zealand,
        col="gray96",bg="#8DCBE4",
        axes=TRUE, 
        xlim=c(1000000,2005000))
  points(MMU_points_NZ,col="#E65518",pch=20,cex=.4)
  title("Suitable Sites for a minimum of 35 females")}
  dev.off()
  
```

Plot here.

```{r,warning=FALSE,cache=TRUE,cache.comments=FALSE,fig.height=9,fig.width=9}
# plot here
  {plot(New_Zealand,
        col="gray96",bg="#8DCBE4",
        axes=TRUE, 
        xlim=c(1000000,2005000))
  points(MMU_points_NZ,col="#E65518",pch=20,cex=.4)
  title("Suitable Sites for minimum of 35 females (multi-state SDM)")}
  
```


# 6. Save workspace

```{r}
# save workspace to load later if needed
  save.image("G:/R/mssdm_NZ.RData")
```

This concludes the multi-state SDM.

Next, we conducted CV, MESS, MoD and limiting factor evaluations. **See Appendix S3**.


# 7. References

Augé, A. A., Chilvers, B. L., Moore, A., Mathieu, R., & Robertson, B. C. (2009). Aggregation and dispersion of female New Zealand sea lions at the Sandy Bay breeding colony, Auckland Islands: how unusual is their spatial behaviour? Behaviour, 146(9), 1287–1311. https://doi.org/10.1163/15683909X427687

Frans, V.F., Augé, A.A., Edelhoff, H., Erasmi, S., Balkenhol, N. & Engler, J.O. (2018) Quantifying apart what belongs together: A multi-state species distribution modelling framework for species using distinct habitats. Methods in Ecology and Evolution, 9 (1), 98-108. https://doi.org/10.1111/2041-210X.12847

Phillips, S.J. & Dudik, M. (2008) Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography, 31 (2), 161-175. https://doi.org/10.1111/j.0906-7590.2008.5203.x
