01_Wallace.rmd

---
title: "Modeling in Wallace"
author: "Bethany A. Johnson"
output:
  html_notebook: default
date: "19 Jan 2023"
---

Bethany A. Johnson  
Here is the session code from modeling *H. chapmani* on Nov 1, 2022.  
Updated 19 Jan 2023. 

```{r, eval=FALSE}
# Install Wallace from GitHub
install.packages("remotes")
remotes::install_github("wallaceEcoMod/wallace")
# Load Wallace
library(wallace)
# Run Wallace
run_wallace()
```

# Modeling in Wallace

Wallace Session 2022-11-01
================

Please find below the R code history from your *Wallace* v1.9.4-5
session.

You can reproduce your session results by running this R Markdown file
in RStudio.

Each code block is called a “chunk”, and you can run them either
one-by-one or all at once by choosing an option in the “Run” menu at the
top-right corner of the “Source” pane in RStudio.

For more detailed information see <http://rmarkdown.rstudio.com>).

### Package installation

Wallace uses the following R packages that must be installed and loaded
before starting.

```{r}
library(spocc)
library(spThin)
library(dismo)
library(rgeos)
library(ENMeval)
library(wallace)
```

My analyses are below.

------------------------------------------------------------------------

## Analysis for *Handleyomys chapmani* (Hc)

User CSV path with occurrence data. If the CSV file is not in the
current workspace, change to the correct file path
(e.g.“/Users/darwin/Documents/occs/”).

INPUTS  
**occs_full.csv** : This is the full dataset of occurrence points for *H. chapmani* (116 total) obtained from ILCB & LG  
**occs_processed.csv** : this dataset has been processed (41 total); removed two outliers and applied a spatial thin of 10km  
**CHELSA_bio10_xx** : Chelsa v1.2 bioclimatic variables 05, 06, 13, & 14 (these were not included due to file size)  

```{r}
# NOTE: provide the folder path of the .csv file
occs_path <- "/Users/bajohnson135/Desktop/github/Handleyomys/inputs/"
occs_path <- file.path(occs_path, "occs_processed.csv")
# get a list of species occurrence data
userOccs_Hc <- occs_userOccs(
  txtPath = occs_path, 
  txtName = "occs_processed.csv", 
  txtSep = ",", 
  txtDec = ".")
occs_Hc <- userOccs_Hc$Handleyomys_chapmani$cleaned
```

### Obtain environmental data

Using user-specified variables.

```{r}
## Specify the directory with the environmental variables
dir_envs_Hc <- "/Users/bajohnson135/Desktop/github/Handleyomys/inputs"
envs_path <- file.path(dir_envs_Hc, c('CHELSA_bio10_05.tif', 'CHELSA_bio10_06.tif', 'CHELSA_bio10_13.tif', 'CHELSA_bio10_14.tif'))
# Create environmental object 
envs_Hc <- envs_userEnvs(
  rasPath = envs_path,
  rasName = c('CHELSA_bio10_05.tif', 'CHELSA_bio10_06.tif', 'CHELSA_bio10_13.tif', 'CHELSA_bio10_14.tif'),
  doBrick = FALSE)
occs_xy_Hc <- occs_Hc[c('longitude', 'latitude')]
occs_vals_Hc <- as.data.frame(raster::extract(envs_Hc, occs_xy_Hc, cellnumbers = TRUE))
# Remove duplicated same cell values
occs_Hc <- occs_Hc[!duplicated(occs_vals_Hc[, 1]), ]
occs_vals_Hc <- occs_vals_Hc[!duplicated(occs_vals_Hc[, 1]), -1]
# remove occurrence records with NA environmental values
occs_Hc <- occs_Hc[!(rowSums(is.na(occs_vals_Hc)) >= 1), ]
# also remove variable value rows with NA environmental values
occs_vals_Hc <- na.omit(occs_vals_Hc)
# add columns for env variable values for each occurrence record
occs_Hc <- cbind(occs_Hc, occs_vals_Hc)
```

### Process environmental data

Sampling of 131680 background points and corresponding environmental
data using a “point buffers” method with a 0.5 degree buffer.

```{r}
# Generate background extent 
bgExt_Hc <- penvs_bgExtent(
  occs = occs_Hc,
  bgSel = "point buffers",
  bgBuf = 0.5)
```

```{r}
# Mask environmental data to provided extent
bgMask_Hc <- penvs_bgMask(
  occs = occs_Hc,
  envs = envs_Hc,
  bgExt = bgExt_Hc)
```

```{r}
# Sample background points from the provided area
bgSample_Hc <- penvs_bgSample(
  occs = occs_Hc,
  bgMask =  bgMask_Hc,
  bgPtsNum = 131680)
```

```{r}
# Extract values of environmental layers for each background point
bgEnvsVals_Hc <- as.data.frame(raster::extract(bgMask_Hc,  bgSample_Hc))
##Add extracted values to background points table
bgEnvsVals_Hc <- cbind(scientific_name = paste0("bg_", "Handleyomys chapmani"), bgSample_Hc,
                            occID = NA, year = NA, institution_code = NA, country = NA,
                            state_province = NA, locality = NA, elevation = NA,
                            record_type = NA, bgEnvsVals_Hc)
```

### Partition occurrence data

Partition occurrences and background points for model training and
validation using “spatial block”, a spatial partition method with an
aggregation factor of 2.

```{r}
# R code to get partitioned data
groups_Hc <- part_partitionOccs(
  occs = occs_Hc ,
  bg =  bgSample_Hc, 
  method = "block",
  bgMask = bgMask_Hc,
  aggFact = 2) 
```

### Build and Evaluate Niche Model

Generating a species distribution model using the maxnet algorithm as
implemented in ENMeval V2.0 (with clamping = TRUE). For tuning using L,
LQ, H, LQH feature classes and regularization multipliers in the 1, 4
range increasing by 0.5. Not using any categorical predictor variables.

```{r}
# Run maxent model for the selected species
model_Hc <- model_maxent(
  occs = occs_Hc,
  bg = bgEnvsVals_Hc,
  user.grp = groups_Hc, 
  bgMsk = bgMask_Hc,
  rms = c(1, 5), 
  rmsStep =  0.5,
  fcs = c('L', 'LQ', 'H', 'LQH'),
  clampSel = TRUE,
  algMaxent = "maxnet",
  parallel = FALSE,
  numCores = 3)
```

H2 was selected as the optimal model.

### Visualize

Generate a map of the maxnet generated model with with a “mtp” threshold
rule of 0.236196536870172.

```{r}
# Select current model and obtain raster prediction
m_Hc <- model_Hc@models[["fc.H_rm.2"]] 
predSel_Hc <- predictMaxnet(m_Hc, bgMask_Hc,
                                          type = "cloglog", 
                                          clamp = TRUE)
# extract the suitability values for all occurrences
occs_xy_Hc <- occs_Hc[c('longitude', 'latitude')]
# determine the threshold based on the current prediction
occPredVals_Hc <- raster::extract(predSel_Hc, occs_xy_Hc)
# Define probability of quantile based on selected threshold
thresProb_Hc <- switch("mtp", 
                              "mtp" = 0, "p10" = 0.1, "qtp" = 0)
# Define threshold value
thres_Hc <- stats::quantile(occPredVals_Hc, probs = thresProb_Hc)
# Applied selected threshold
predSel_Hc <- predSel_Hc > thres_Hc
```

```{r}
# Get values of prediction
mapPredVals_Hc <- getRasterVals(predSel_Hc, "cloglog")
# Define colors and legend  
rasCols <- c("#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c")
legendPal <- colorNumeric(rev(rasCols), mapPredVals_Hc, na.color = 'transparent')
rasPal <- c('gray', 'blue')
```

```{r}
# Generate map
m <- leaflet() %>% addProviderTiles(providers$Esri.WorldTopoMap) 
m  %>%
  leaflet::addLegend("bottomright", colors = c('gray', 'blue'),
            title = "Thresholded Suitability<br>(Training)",
            labels = c("predicted absence", "predicted presence"),
            opacity = 1, layerId = "train") %>% 
  #add occurrence data
  addCircleMarkers(data = occs_Hc, lat = ~latitude, lng = ~longitude,
                   radius = 5, color = 'red', fill = TRUE, fillColor = "red",
                   fillOpacity = 0.2, weight = 2, popup = ~pop) %>% 
  ##Add model prediction
  addRasterImage(predSel_Hc, colors = rasPal, opacity = 0.7,
                 group = 'vis', layerId = 'mapPred', method = "ngb") %>%
 ##add background polygons
  addPolygons(data = bgExt_Hc, fill = FALSE,
              weight = 4, color = "blue", group = 'proj')
```


```{r}
sdm_mtp <- predSel_Hc
plot(sdm_mtp)
```

```{r}
# set NA values to 0, keep values above 0
sdm_mtp[is.na(sdm_mtp)] <- 0
sdm_mtp <- sdm_mtp > 0
plot(sdm_mtp)
```

```{r}
writeRaster(sdm_mtp, "/Users/bajohnson135/Desktop/Manuscript/H_chapmani/sdm_mtp.tif")
```

OUTPUTS  

1. **SDM.tif**: the continuous prediction from H2  
2. **SDM_mtp.tif**: the minimum training presence threshold prediction
3. **Handleyomys_chapmani_bgShp.dbf**, **Handleyomys_chapmani_bgShp.shp**, **Handleyomys_chapmani_bgShp.shx**: the shapefile of the background extent background extent  
4. **Handleyomys_chapmani_maxnet_evalTbl.csv**: Maxent evaluation table for candidate models