Intro

First, we go to the Wallace web-page (https://wallaceecomod.github.io/) that is used for species distribution and niche Modeling. It’s open source!

Required Libraries

library(shiny)
library(leaflet)
library(raster)
library(sp)
library(sf)
library(ENMeval)
library(wallace)
library(spocc)
library(spThin)
library(dismo)

Analysis for Morpho sulkowskyi

Occurrence Data

queryDb_Ms <- occs_queryDb(
  spNames = "Morpho sulkowskyi", 
  occDb = "gbif", 
  occNum = 300,
  RmUncertain = FALSE)
occs_Ms <- queryDb_Ms$Morpho_sulkowskyi$cleaned

Environmental Data

envs_Ms <- envs_worldclim(
  bcRes = 2.5, 
  bcSel = c('bio02', 'bio05', 'bio08', 'bio09', 'bio11', 'bio18'), 
  mapCntr = c(-75.356, 4.748), 
  doBrick = TRUE)

occs_xy_Ms <- occs_Ms[c('longitude', 'latitude')]
occs_vals_Ms <- as.data.frame(raster::extract(envs_Ms, occs_xy_Ms, cellnumbers = TRUE))

occs_Ms <- occs_Ms[!duplicated(occs_vals_Ms[, 1]), ]
occs_vals_Ms <- occs_vals_Ms[!duplicated(occs_vals_Ms[, 1]), -1]

occs_Ms <- occs_Ms[!(rowSums(is.na(occs_vals_Ms)) >= 1), ]
occs_vals_Ms <- na.omit(occs_vals_Ms)

occs_Ms <- cbind(occs_Ms, occs_vals_Ms)

Process Background

bgExt_Ms <- penvs_bgExtent(
  occs = occs_Ms,
  bgSel = "point buffers",
  bgBuf = 2)

bgMask_Ms <- penvs_bgMask(
  occs = occs_Ms,
  envs = envs_Ms,
  bgExt = bgExt_Ms)

bgSample_Ms <- penvs_bgSample(
  occs = occs_Ms,
  bgMask =  bgMask_Ms,
  bgPtsNum = 10000)

bgEnvsVals_Ms <- as.data.frame(raster::extract(bgMask_Ms,  bgSample_Ms))

bgEnvsVals_Ms <- cbind(scientific_name = "bg_Morpho sulkowskyi", bgSample_Ms,
                       occID = NA, year = NA, institution_code = NA, country = NA,
                       state_province = NA, locality = NA, elevation = NA,
                       record_type = NA, bgEnvsVals_Ms)

Partition Data

groups_Ms <- part_partitionOccs(
  occs = occs_Ms ,
  bg =  bgSample_Ms, 
  method = "cb2",
  bgMask = bgMask_Ms,
  aggFact = 3)

Modeling

model_Ms <- model_maxent(
  occs = occs_Ms,
  bg = bgEnvsVals_Ms,
  user.grp = groups_Ms, 
  bgMsk = bgMask_Ms,
  rms = c(1, 2), 
  rmsStep = 1,
  fcs = c('L', 'LQ'),
  clampSel = TRUE,
  algMaxent = "maxnet",
  parallel = FALSE,
  numCores = 2)
##   |                                                                              |                                                                      |   0%  |                                                                              |==================                                                    |  25%  |                                                                              |===================================                                   |  50%  |                                                                              |====================================================                  |  75%  |                                                                              |======================================================================| 100%

Prediction Map

m_Ms <- model_Ms@models[["fc.L_rm.1"]] 
predSel_Ms <- predictMaxnet(m_Ms, bgMask_Ms, type = "cloglog", clamp = TRUE)

occPredVals_Ms <- raster::extract(predSel_Ms, occs_xy_Ms)
thres_Ms <- quantile(occPredVals_Ms, probs = 0.1)
predSel_Ms_thres <- predSel_Ms > thres_Ms

leaflet() %>% 
  addProviderTiles(providers$Esri.WorldTopoMap) %>%
  addRasterImage(predSel_Ms_thres, colors = c("gray", "blue"), opacity = 0.7) %>%
  addCircleMarkers(data = occs_Ms, lat = ~latitude, lng = ~longitude,
                   radius = 3, color = "red", fill = TRUE, fillOpacity = 0.6)

Response Curves

n <- mxNonzeroCoefs(model_Ms@models[["fc.L_rm.1"]], "maxnet")

for (i in n) {
  maxnet::response.plot(model_Ms@models[["fc.L_rm.1"]], v = i, type = "cloglog")
}