Wallace Session 2023-05-01
================

Please find below the R code history from your *Wallace* v2.0.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)
```

The *Wallace* session code .Rmd file is composed of a chain of module
functions that are internal to *Wallace*. Each of these functions
corresponds to a single module that the user ran during the session. To
see the internal code for these module functions, click on the links in
the .Rmd file. Users are encouraged to write custom code in the .Rmd
directly to modify their analysis, and even modify the module function
code to further customize. To see the source code for any module
function, just type its name into the R console and press Return.

```{r}
# example:
# just type the function name and press Return to see its source code
# paste this code into a new script to edit it
occs_queryDb
```

Your analyses are below.

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

## Analysis for *Vespa velutina* (Vv)

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

```{r}
# NOTE: provide the folder path of the .csv file
occs_path <- ""
occs_path <- file.path(occs_path, "Vespa_velutina_FINAL.txt")
# get a list of species occurrence data
userOccs_Vv <- occs_userOccs(
  txtPath = occs_path, 
  txtName = "Vespa_velutina_FINAL.txt", 
  txtSep = ";", 
  txtDec = ",")
occs_Vv <- userOccs_Vv$Vespa_velutina$cleaned
```

### Obtain environmental data

Using user-specified variables.

```{r}
## Specify the directory with the environmental variables
dir_envs_Vv <- ""
envs_path <- file.path(dir_envs_Vv, c('Bio2_AMDR.tif', 'Bio3_Isothermality.tif', 'Bio4_Temp_seasonality.tif', 'Bio10_MTWQ.tif', 'Bio11_MTCQ.tif', 'Bio15_PPT_seasonality.tif', 'Bio18_PPTWQ.tif', 'Bio19_PPTCQ.tif', 'LC.tif'))
# Create environmental object 
envs_Vv <- envs_userEnvs(
  rasPath = envs_path,
  rasName = c('Bio2_AMDR.tif', 'Bio3_Isothermality.tif', 'Bio4_Temp_seasonality.tif', 'Bio10_MTWQ.tif', 'Bio11_MTCQ.tif', 'Bio15_PPT_seasonality.tif', 'Bio18_PPTWQ.tif', 'Bio19_PPTCQ.tif', 'LC.tif'),
  doBrick = TRUE)
occs_xy_Vv <- occs_Vv[c('longitude', 'latitude')]
occs_vals_Vv <- as.data.frame(raster::extract(envs_Vv, occs_xy_Vv, cellnumbers = TRUE))
# Remove duplicated same cell values
occs_Vv <- occs_Vv[!duplicated(occs_vals_Vv[, 1]), ]
occs_vals_Vv <- occs_vals_Vv[!duplicated(occs_vals_Vv[, 1]), -1]
# remove occurrence records with NA environmental values
occs_Vv <- occs_Vv[!(rowSums(is.na(occs_vals_Vv)) >= 1), ]
# also remove variable value rows with NA environmental values
occs_vals_Vv <- na.omit(occs_vals_Vv)
# add columns for env variable values for each occurrence record
occs_Vv <- cbind(occs_Vv, occs_vals_Vv)
```

### Process Occurrence Data

Thinning the occurrences to 5 km

```{r}
# Thin occurrences 
occs_Vv <- poccs_thinOccs(
  occs = occs_Vv, 
  thinDist = 5)
```

### Process environmental data

Sampling of 50000 background points and corresponding environmental data
using a user provided background extent with a 0 degree buffer.

```{r}
# Load the user provided shapefile or csv file with the desired extent.
  ##User must input the path to shapefile or csv file and the file name 
# Define path
bgPath_Vv <- ""
bgExt_Vv <- penvs_userBgExtent(
  bgShp_path = paste0(bgPath_Vv, "Study_area", ".shp"),
  bgShp_name = paste0("Study_area", c(".shp", ".shx", ".dbf")),
  userBgBuf = 0,
  occs = occs_Vv)
# Mask environmental data to provided extent
bgMask_Vv <- penvs_bgMask(
  occs = occs_Vv,
  envs = envs_Vv,
  bgExt = bgExt_Vv)
# Sample background points from the provided area
bgSample_Vv <- penvs_bgSample(
  occs = occs_Vv,
  bgMask =  bgMask_Vv,
  bgPtsNum = 50000)
# Extract values of environmental layers for each background point
bgEnvsVals_Vv <- as.data.frame(raster::extract(bgMask_Vv,  bgSample_Vv))
##Add extracted values to background points table
bgEnvsVals_Vv <- cbind(scientific_name = paste0("bg_", "Vespa velutina"), bgSample_Vv,
                            occID = NA, year = NA, institution_code = NA, country = NA,
                            state_province = NA, locality = NA, elevation = NA,
                            record_type = NA, bgEnvsVals_Vv)
```

### 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_Vv <- part_partitionOccs(
  occs = occs_Vv ,
  bg =  bgSample_Vv, 
  method = "block",
  bgMask = bgMask_Vv,
  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 LQ,
H, LQH, LQHP feature classes and regularization multipliers in the 0.5,
2.5 range increasing by 0.5. Using a total of 1 categorical predictor
variables.

```{r}
# Run maxent model for the selected species
model_Vv <- model_maxent(
  occs = occs_Vv,
  bg = bgEnvsVals_Vv,
  user.grp = groups_Vv, 
  bgMsk = bgMask_Vv,
  rms = c(0.5, 2.5), 
  rmsStep =  0.5,
  fcs = c('LQ', 'H', 'LQH', 'LQHP'),
  clampSel = TRUE,
  algMaxent = "maxnet",
  catEnvs = "LC",
  parallel = FALSE,
  numCores = 31)
```

### Visualize

Generate a map of the maxnet generated model with no threshold

```{r}
# Select current model and obtain raster prediction
m_Vv <- model_Vv@models[["fc.LQHP_rm.1"]]
predSel_Vv <- predictMaxnet(m_Vv, bgMask_Vv,
                                          type = "cloglog", 
                                          clamp = TRUE)
#Get values of prediction
mapPredVals_Vv <- getRasterVals(predSel_Vv, "cloglog")
#Define colors and legend  
rasCols <- c("#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c")
legendPal <- colorNumeric(rev(rasCols), mapPredVals_Vv, na.color = 'transparent')
rasPal <- colorNumeric(rasCols, mapPredVals_Vv, na.color = 'transparent')
#Generate map
m <- leaflet() %>% addProviderTiles(providers$Esri.WorldTopoMap) 
m  %>%
  leaflet::addLegend("bottomright", pal = legendPal,
            title = "Predicted Suitability<br>(Training)",
            values = mapPredVals_Vv, layerId = "train",
            labFormat = reverseLabel(2, reverse_order = TRUE)) %>% 
  #add occurrence data
  addCircleMarkers(data = occs_Vv, lat = ~latitude, lng = ~longitude,
                   radius = 5, color = 'red', fill = TRUE, fillColor = "red",
                   fillOpacity = 0.2, weight = 2, popup = ~pop) %>% 
  ##Add model prediction
  addRasterImage(predSel_Vv, colors = rasPal, opacity = 0.7,
                 group = 'vis', layerId = 'mapPred', method = "ngb") %>%
 ##add background polygons
  addPolygons(data = bgExt_Vv,fill = FALSE,
              weight = 4, color = "blue", group = 'proj')
```

### Visualize

Generate a Maxent evaluation plot using “auc.val” as evaluation
statistic.

```{r}
# Generate an evaluation plot
maxentEvalPlot_Vv <- ENMeval::evalplot.stats(
   model_Vv,
   "auc.val",
   "rm",
   "fc")
#plot
maxentEvalPlot_Vv
```

```{r}
# First must generate the transfer area according to the drawn polygon in the GUI
xfer_draw_Vv <-xfer_draw(
  polyXfXY = matrix(c(19.087939, 11.351981, 1.857851, -2.098036, -9.746085, -12.558963, -12.558963, -4.647188, -1.746204, 4.583216, 9.989595, 15.527838, 17.90137, 18.956274, 19.087939, 50.260377, 54.647619, 54.494771, 53.799841, 49.209525, 41.923059, 36.866806, 35.554949, 35.982824, 37.983538, 38.191066, 38.707306, 39.219844, 48.505992, 50.260377),ncol=2,byrow=FALSE),
  polyXfID = 81398,
  drawXfBuf = 0)
# Create object of transfer variables
xferAreaEnvs_Vv <- envs_Vv
# Generate a transfer of the model to the desired area
xfer_area_Vv <- xfer_area(
  evalOut = model_Vv,
  curModel = "fc.LQHP_rm.1",
  envs = xferAreaEnvs_Vv , 
  outputType = "cloglog",
  alg = "maxnet",
  clamp = TRUE,
  xfExt = xfer_draw_Vv) 
# store the cropped transfer variables
xferExt_Vv <- xfer_area_Vv$xferExt
 #map result
 
 ###Make map of transfer
bb_Vv <-  bgExt_Vv@bbox
bbZoom <- polyZoom(bb_Vv[1, 1], bb_Vv[2, 1], bb_Vv[1, 2], 
                   bb_Vv[2, 2], fraction = 0.05)
mapXferVals_Vv <- getRasterVals(xfer_area_Vv$xferArea,"cloglog")
rasCols_Vv <- c("#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c")
# if no threshold specified
legendPal <- colorNumeric(rev(rasCols_Vv), mapXferVals_Vv, na.color = 'transparent')
rasPal_Vv <- colorNumeric(rasCols_Vv, mapXferVals_Vv, na.color = 'transparent')
m <- leaflet() %>% addProviderTiles(providers$Esri.WorldTopoMap) 
m %>%
  fitBounds(bbZoom[1], bbZoom[2], bbZoom[3], bbZoom[4]) %>%
  leaflet::addLegend("bottomright", pal = legendPal,
            title = "Predicted Suitability<br>(Transferred)",
            values = mapXferVals_Vv, layerId = 'xfer',
            labFormat = reverseLabel(2, reverse_order = TRUE)) %>%
# map model prediction raster and transfer polygon
  clearMarkers() %>% clearShapes() %>% removeImage('xferRas') %>%
  addRasterImage(xfer_area_Vv$xferArea, colors = rasPal_Vv, opacity = 0.7,
                 layerId = 'xferRas', group = 'xfer', method = "ngb") %>%
 ##add transfer polygon (user drawn area)
  addPolygons(data = xfer_draw_Vv, fill = FALSE,
              weight = 4, color = "blue", group = 'xfer')
```