Maxent_SDM.Rmd

---
title: "Maxent SDM using BioClim and Land Use data"
  
author: "Pablo Deschepper"
date: "`r Sys.Date()`"
output: html_document

params:
  Species:
    label: "Species:"
    value: "Vipera latastei"
    input: text
  Country:
    label: "2 letter country code (optional):"
    value: ""
    input: text
  Continent:
    label: "Continent (optional):"
    value: all
    input: select
    choices: [all, africa, antarctica, asia, europe, north_america, south_america, oceania]
  Database:
    label: "Database:"
    value: gbif
    input: select
    choices: [gbif, inaturalist, observation]
  Limit:
    label: "Max number of records:"
    value: "1000"
    input: numeric
  Workdir: 
    label: "Select your working directory:"
    value: "D:/temp/SDM_Herpetofauna/Maxent_HighRes_SSP245SSP585/"
    input: text
  Ncores:
    label: "Number of computing cores:"
    value: 10
    input: numeric
  Range:
    label: "Apply an overlay of the species' range to the predictions:"
    value: "YES"
    input: select
    choices: ["YES","NO"]
  ExportRaster:
    label: "Export raster objects for predictions:"
    value: "NO"
    input: select
    choices: ["NO", "YES"]
  Climate.model:
    label: "Select the CMIP6 climate model:"
    value: MPI-ESM1-2-HR
    input: select
    choices: [ACCESS-CM2, ACCESS-ESM1-5, AWI-CM-1-1-MR, BCC-CSM2-MR, CanESM5, CanESM5-CanOE, CMCC-ESM2, CNRM-CM6-1,CNRM-CM6-1-HR,CNRM-ESM2-1,EC-Earth3-Veg, EC-Earth3-Veg-LR, FIO-ESM-2-0, GFDL-ESM4, GISS-E2-1-G, GISS-E2-1-H, HadGEM3-GC31-LL, INM-CM4-8,INM-CM5-0,IPSL-CM6A-LR,MIROC-ES2L, MIROC6, MPI-ESM1-2-HR, MPI-ESM1-2-LR, MRI-ESM2-0, UKESM1-0-LL]
  SSP:
    label: "Select the SSP code:"
    value: 585
    input: select
    choices: [245, 585]
  Plot.xmin:
    label: "The numeric longitude minimum for SDM plots (optional):"
    value: ""
    input: numeric
  Plot.xmax:
    label: "The numeric longitude maximum for SDM plots (optional):"
    value: ""
    input: numeric
  Plot.ymin:
    label: "The numeric latitude minimum for SDM plots (optional):"
    value: ""
    input: numeric
  Plot.ymax:
    label: "The numeric latitude minimum for SDM plots (optional):"
    value: ""
    input: numeric
  

---

```{r setup, include = TRUE, eval = TRUE, echo = FALSE}
library(knitr)
knitr::opts_chunk$set(warning = FALSE, message = TRUE)
knitr::opts_knit$set(root.dir = params$Workdir)
logo_url <- "https://github.com/thesnakeguy/SDM_BioClim-and-LandUse/blob/main/Logo_SDM_resize.png?raw=true"
include_graphics(logo_url)
```
<br> 
<br> 
*This Markdown script was developed to model species distributions and ecological niches under both current and future climate conditions. The species distribution models (SDMs) for present-day conditions are trained using BioClim data, land use data, or a combination of both. For future projections, the models use BioClim data based on climate projections alone or combined with contemporary land use data.*<br>
*The modeled distributions reflect the species’ ecological niche, although several factors—such as dispersal barriers, geological history, and biotic interactions—may create differences between the realized and modeled niches.*<br>
*The script allows flexible model training, enabling users to specify geographic boundaries (country or continent) for occurrence data. Users can visualize the predicted species distribution or ecological niche within customized coordinate ranges, making this tool adaptable to various research needs. The script automatically removes strongly correlated (r > 0.75) predictors within the range of the species occurrences.*<br> 
<br> 
<br> 

- **Species = *`r params$Species`* **
- Using records from following database: `r params$Database`
- The number of records is limited to `r params$Limit` records
- Climate model for future predictions: `r params$Climate.model`
- Shared Socioeconomic Pathway (SSP): `r params$SSP`  

The working directory = *`r params$Workdir`*<br>


note: The script will download all BioClim and Land Use data upon first use of the script, translating into a longer processing times. The script will run faster after first use.

```{r libraries, include=FALSE}
library(dismo)
library(raster)
library(rgbif)
library(sf)
library(rmaxent)
library(rJava)
library(rasterVis)
library(viridis)
library(geodata)
library(ggplot2)
library(stringr)
library(rnaturalearth)
library(caret)
library(dichromat)
library(dplyr)
library(kableExtra)
library(rangeBuilder)
library(ENMeval)
library(cowplot)
```

```{r Get species records, include=FALSE, echo=FALSE}
##################### PARAMETERS ##################################

# set georeference obscurity in meters for fetching gbif records
georefobsc <- 500
georefobsc_gbif <- paste0("0,",georefobsc)
# minimum records needed for processing
minrecs <- 100

###################################################################

#### Search for species occurrences in GBIF ####


ContinentParam <- NULL
if (params$Continent == "all") {
  ContinentParam <- ""} else {ContinentParam <- params$Continent}

taxonKey <- name_backbone(params$Species)$usageKey

occurrences <- rgbif::occ_data(taxonKey = taxonKey,
                              limit = params$Limit,
                              hasCoordinate = TRUE,
                              hasGeospatialIssue = FALSE,
                              coordinateUncertaintyInMeters = georefobsc_gbif,
                              country = params$Country,
                              continent = ContinentParam)


occurrences <- data.frame(occurrences$data[, c("decimalLongitude", "decimalLatitude", "occurrenceID", "country")])

if (params$Database == "inaturalist") {
  occurrences <- occurrences[str_detect(occurrences$occurrenceID, "inaturalist"),] %>% na.omit()
  recs <- nrow(occurrences)
  print(paste("Using",recs,"inaturalist records"))
  } else if (params$Database == "observation") {
    occurrences <- occurrences[str_detect(occurrences$occurrenceID, "observation.org"),] %>% na.omit()
    recs <- nrow(occurrences)
    print(paste("Using",recs,"observation.org records"))
  } else {
    recs <- nrow(occurrences)
    print(paste("Using",recs,"gbif records"))
  }
  
if (recs < minrecs) {
  stop("Fewer than ", minrecs," species records: Process aborted.")
}

# Create a SpatialPointsDataFrame object for the occurrence data
spdf <- SpatialPointsDataFrame(coords = occurrences[, c("decimalLongitude", "decimalLatitude")], data = occurrences)
# Bounding box of distribution
bbox <- st_bbox(spdf) 
points <- data.frame(occurrences[, c("decimalLongitude", "decimalLatitude")])
# Create an extent object for the study area using bounding box coordinates
# By adding 2, we add 1 degree at each side
study_area_ext <- extent(bbox)
# Create a raster object for the study area
study_area_ras <- raster::raster(study_area_ext)
# Set the resolution of the raster
res(study_area_ras) <- 2.5

```

```{r Estimating a species range map, include=FALSE, echo=FALSE, message=FALSE}

# Create range map
if (params$Range == "YES"){
  # Get alphahull for the species range
  range_polygon <- getDynamicAlphaHull(occurrences[,1:2], buff = 20000)
  range_polygon <- range_polygon[[1]]
  # Convert sf polygon to SpatVector
  range_spatvector <- terra::vect(range_polygon)
  # Step 1: Create an empty raster spanning the extent of the polygon
  raster_template <- rast(terra::ext(range_spatvector), resolution = 0.01, crs = "EPSG:4326")
  # Step 2: Rasterize the polygon
  range <- rasterize(range_spatvector, raster_template, field = 1)
}

```

## Species records
This plot shows the georeferenced species records that will be used to train the Species Distribution Models (SDM) using MaxEnt. It is good to maximize the number of records.
When the option of the species' range map is toggled, a range map is plotted by an alpha hull polygon.

**`r recs` records** of `r params$Species` are used for processing (georeference obscurity < `r georefobsc`m).

```{r Records plot, include=TRUE, echo = FALSE, dev='png'}

all_countries <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
gg <- ggplot() + geom_point(data = points, aes(x=decimalLongitude, y = decimalLatitude), col="blue", size = 0.4) +
  geom_sf(data = all_countries, fill = "transparent", color="black") +
  coord_sf(xlim = c(bbox$xmin-2, bbox$xmax+2), ylim = c(bbox$ymin-2, bbox$ymax+2), expand = TRUE) + theme_bw()

if (params$Range == "YES") {
  range_df <- as.data.frame(range, xy = TRUE, na.rm = TRUE)
  gg <- ggplot() + geom_point(data = points, aes(x=decimalLongitude, y = decimalLatitude), col="blue", size = 0.4) +
    geom_sf(data = all_countries, fill = "transparent", color="black") +
    coord_sf(xlim = c(bbox$xmin-2, bbox$xmax+2), ylim = c(bbox$ymin-2, bbox$ymax+2), expand = TRUE) + 
    geom_raster(data = range_df, aes(x = x, y = y), fill = "limegreen", alpha = 0.5) + theme_bw() 
} 
plot(gg)

```

```{r Load Land Use and BioClimatic data , include=FALSE}

############## FIRST TIME USE ONLY ####################

# UNDER DEVELOPMENT --> in a future update the process of downloading, resampling and merging predictive raster layers will be automated
# USE THIS CODE TO DOWNLOAD BIOCLIM AND LAND COVER DATA. CHANGE THE SSP TO THE DESIRED VALUE. SET IT TO 245 and 585 AND RUN THE SCRIPT TO 
# DOWNLAD THE RESPECTIVE SSP RASTER FILES. ONLY 245 and 585 ARE CURRENTLY SUPPORTED AS PARAMETERS FOR MODELLING USING THE .RMD SCRIPT.


## LAND USE ###
### load land use data (if res is not and option, it is 0.5 minutes) ####
# elevation <- geodata::elevation_global(path = "Elevation", res = 2.5)
# trees <- geodata::landcover(var = "trees", path = "Trees")
# grassland <- geodata::landcover(var = "grassland", path = "Grassland")
# shrubs <- geodata::landcover(var = "shrubs", path = "Shrubs")
# built <- geodata::landcover(var = "built", path = "Built")
# bare <- geodata::landcover(var = "bare", path = "Bare")
# wetland <- geodata::landcover(var = "wetland", path = "Wetland")
# water <- geodata::landcover(var = "water", path = "Water")
# names(elevation) <- c("elevation") # this name was awkwardly long
# 
# # Create a list of the raster objects
# LU_list <- list(elevation, trees, grassland, shrubs, built, bare, wetland, water)
# # After checking resolutions and extent, resample where needed
# sapply(LU_list, res)
# sapply(LU_list, dim)
# LU_list[[1]] <- terra::resample(LU_list[[1]], LU_list[[2]], threads = TRUE)
# # Half the resolution to avoid overkill and long computing times
# LU_list_aggregated <- lapply(LU_list, function(raster) {
#   aggregate(raster, fact = 2, fun = "mean")
# })
# 
# LandUse <- rast(LU_list_aggregated)
# 
# # save
# terra::writeRaster(x = LandUse, filename = paste0(getwd(),"/LandUse.tif"))
# 
# 
# ### BioClim ###
# if (!dir.exists("worldclim_data")) {
#   # Download WorldClim data if the directory does not exist
#   dir.create("worldclim_data")
#   worldclim_data <- geodata::worldclim_global(var = 'bio', res = 2.5, download = TRUE, path = "worldclim_data")
#   print("BioClim data has been downloaded and saved to the working directory.")
# } else {
#   print("BioClim data already exists and was loaded from the working directory.")
#   rast_list <- c(list.files("worldclim_data/climate/wc2.1_2.5m/", full.names = TRUE)) %>%
#     .[grepl("bio_([1-9]|1[0-9])", .)]
#   worldclim_data <- terra::rast(rast_list)
#   }
# # rename bio variables
# names(worldclim_data) <- gsub("wc2.1_2.5m_", "", names(worldclim_data))
# # Resample to match resololution with landuse
# worldclim_data_resampled <- terra::resample(worldclim_data, LandUse, threads = TRUE)
# 
# # save
# terra::writeRaster(x = worldclim_data_resampled, filename = paste0(getwd(),"/worldclim.tif"), overwrite = TRUE)
# 
# #Download BioClim data for two future time periods
# fut <- geodata::cmip6_world(model=params$Climate.model, ssp=params$SSP, time='2041-2060', var='bioc', download=F, res=2.5, path='worldclim_data')
# fut2 <- geodata::cmip6_world(model=params$Climate.model, ssp=params$SSP, time='2061-2080', var='bioc', download=F, res=2.5, path='worldclim_data')
# names(fut) <- str_replace(names(fut), "wc2.1_2.5m_bioc_MPI-ESM1-2-HR_ssp245_2041-2060_(\\d+)", "bio_\\1")
# names(fut2) <- str_replace(names(fut2), "wc2.1_2.5m_bioc_MPI-ESM1-2-HR_ssp245_2061-2080_(\\d+)", "bio_\\1")
# 
# fut <- terra::resample(fut, LandUse, threads = TRUE)
# fut2 <- terra::resample(fut2, LandUse, threads = TRUE)
# 
# terra::writeRaster(x = fut, filename = paste0(getwd(),"/worldclim_fut.tif"), overwrite = TRUE)
# terra::writeRaster(x = fut2, filename = paste0(getwd(),"/worldclim_fut2.tif"), overwrite = TRUE)

```


```{r Combine BioClim and Land use data into predictors, include=FALSE}
############## FIRST TIME USE ONLY ####################
### Make predictors dataset ####
# worldclim_data_aligned <- worldclim_data_resampled
# predictors <- c(worldclim_data_aligned, LandUse)
# fut_predictors <- c(fut, LandUse)
# fut2_predictors <- c(fut2, LandUse)
# 
# terra::writeRaster(x = predictors, filename = paste0(getwd(),"/All_predictors.tif"), overwrite = TRUE)
# terra::writeRaster(x = fut_predictors, filename = paste0(getwd(),"/fut_LUpredictors.tif"), overwrite = TRUE)
# terra::writeRaster(x = fut2_predictors, filename = paste0(getwd(),"/fut2_LUpredictors.tif"), overwrite = TRUE)

```


```{r Read all SpatRaster object into predictor sets, include=TRUE, echo=FALSE}
landuse_rasterbrick <- rast(paste0(params$Workdir,"LandUse.tif"))
worldclim_data_aligned <- rast(paste0(params$Workdir,"worldclim.tif"))
predictors <- rast(paste0(params$Workdir,"All_predictors.tif"))

# Check the SSP scenario specified by the user
if (params$SSP == "245") {
  fut_predictors <- rast(paste0(params$Workdir, "fut_LUpredictors_SSP245.tif"))
  fut2_predictors <- rast(paste0(params$Workdir, "fut2_LUpredictors_SSP245.tif"))
  fut <- rast(paste0(params$Workdir, "worldclim_fut_SSP245.tif"))
  fut2 <- rast(paste0(params$Workdir, "worldclim_fut2_SSP245.tif"))
} else if (params$SSP == "585") {
  fut_predictors <- rast(paste0(params$Workdir, "fut_LUpredictors_SSP585.tif"))
  fut2_predictors <- rast(paste0(params$Workdir, "fut2_LUpredictors_SSP585.tif"))
  fut <- rast(paste0(params$Workdir, "worldclim_fut_SSP585.tif"))
  fut2 <- rast(paste0(params$Workdir, "worldclim_fut2_SSP585.tif"))
} else {
  stop("Invalid SSP scenario specified. Please use '245' or '585'.")
}

```

```{r Clean the Full, BioClim and Land Use dataset by removing correlation, include = TRUE, echo = FALSE}

# Clean dataset for correlation > 0.75
datasets <- list(predictors = predictors, 
                 worldclim_data_aligned = worldclim_data_aligned, 
                 landuse_rasterbrick = landuse_rasterbrick)
cleaned_datasets <- list()
# Loop through each dataset and remove highly correlated variables within the range of the species -> there could be correlation within the region, but not worldwide!
for (dataset_name in names(datasets)) {
  dataset_df <- crop(datasets[[dataset_name]], study_area_ext) %>% as.data.frame(na.rm = TRUE)
  cor_mat <- cor(dataset_df, method = 'spearman', use = 'complete.obs')
  highly_correlated <- findCorrelation(cor_mat, cutoff = 0.75)
  # Only remove layers if there are correlated layers to remove
  if (length(highly_correlated) > 0) {
    cleaned_dataset <- datasets[[dataset_name]][[ -highly_correlated ]]
  } else {
    cleaned_dataset <- datasets[[dataset_name]]
  }
  cleaned_datasets[[dataset_name]] <- cleaned_dataset
  # print(paste("Cleaned layers in", dataset_name, ":"))
  # print(names(cleaned_dataset))
}

# Access the cleaned datasets
predictors_cleaned <- cleaned_datasets$predictors 
worldclim_data_aligned_cleaned <- cleaned_datasets$worldclim_data_aligned 
landuse_rasterbrick_cleaned <- cleaned_datasets$landuse_rasterbrick 
# Also subset (clean) the future BioClim data and combine to full predictors set
fut_bioclim_cleaned <- fut[[names(worldclim_data_aligned_cleaned)]] 
fut2_bioclim_cleaned <- fut2[[names(worldclim_data_aligned_cleaned)]] 
fut_predictors_cleaned <- fut_predictors[[names(predictors_cleaned)]] 
fut2_predictors_cleaned <- fut2_predictors[[names(predictors_cleaned)]]

```


```{r Run the maxent models after tuning with ENMevaluate, include=FALSE, echo=FALSE}
#######################

# Tuning with ENMeval

#######################


occs <- data.frame(spdf$decimalLongitude, spdf$decimalLatitude)
# List of environmental datasets
env_datasets <- list(
  predictors = terra::crop(predictors_cleaned, study_area_ras),
  worldclim = terra::crop(worldclim_data_aligned_cleaned, study_area_ras),
  landuse = terra::crop(landuse_rasterbrick_cleaned, study_area_ras)
)

# Initialize lists to store evaluation results and best models
eval_results <- list()
best_models <- list()
best_params <- list()

# Loop through each environmental dataset
for (env_name in names(env_datasets)) {
  # Evaluate the model
  eval_results[[env_name]] <- ENMevaluate(
    occ = occs,
    env = env_datasets[[env_name]], 
    partitions = "block", 
    parallel = TRUE, 
    algorithm = "maxnet",
    numCores = params$Ncores, 
    tune.args = list(fc = c("L", "Q", "LQ"), rm = seq(0.5, 2, 0.5))
  )
  
  # Get the best model index based on AICc
  bestmod <- which(eval_results[[env_name]]@results$AICc == min(eval_results[[env_name]]@results$AICc))
  
  # Store the best model and parameters
  best_params[[env_name]] <- eval_results[[env_name]]@results[bestmod,]
}


# Store results
params_predictors <- best_params$predictors
params_worldclim <- best_params$worldclim
params_landuse <- best_params$landuse

opt_params <- list(
  predictors = list(fc = params_predictors$fc, rm = params_predictors$rm),
  worldclim = list(fc = params_worldclim$fc, rm = params_worldclim$rm),
  landuse = list(fc = params_landuse$fc, rm = params_landuse$rm)
)

##################################################

# Run dismo::maxent with tuned parameters

#################################################


#create test and train data #
# Now let us divide this file into 5 parts. 75% of the data will be used to train the model (to create the mode), and 25% will be used to test if it is a good model.
group <- kfold(points,5)
pres_train <- points[group!=1,]
pres_test <- points[group==1,]

# Function to create maxent arguments based on optimized feature classes and regularization
fc_to_args <- function(fc, rm) {
  args <- c("linear=false", "quadratic=false", "product=false", "threshold=false", "hinge=false")
  
  if (grepl("L", fc)) args[1] <- "linear=true"
  if (grepl("Q", fc)) args[2] <- "quadratic=true"
  # Add other feature classes if needed, e.g., "P" for product, "T" for threshold, "H" for hinge
  
  # Add regularization parameter
  args <- c(args, paste0("beta_lqp=", rm))
  
  return(args)
}

# Applying this for each model setup using dismo::maxent with optimal parameters
model_full <- maxent(
  raster::brick(crop(predictors_cleaned, study_area_ras)),
  pres_train,
  removeDuplicates = TRUE,
  replicates = 10,
  args = fc_to_args(opt_params$predictors$fc, opt_params$predictors$rm)
)

model_bioclim <- maxent(
  raster::brick(crop(worldclim_data_aligned_cleaned, study_area_ras)),
  pres_train,
  removeDuplicates = TRUE,
  replicates = 10,
  args = fc_to_args(opt_params$worldclim$fc, opt_params$worldclim$rm)
)

model_landuse <- maxent(
  raster::brick(crop(landuse_rasterbrick_cleaned, study_area_ras)),
  pres_train,
  removeDuplicates = TRUE,
  replicates = 10,
  args = fc_to_args(opt_params$landuse$fc, opt_params$landuse$rm)
)

################################################################################################################

# create pseudo-absences to evaluate the model using AUC
predictorset <- list(raster::brick(crop(predictors_cleaned, study_area_ras)), raster::brick(crop(worldclim_data_aligned_cleaned, study_area_ras)), raster::brick(crop(landuse_rasterbrick_cleaned, study_area_ras)))
models <- c(model_full, model_bioclim, model_landuse)
plot_titles <- c("Full model", "BioClim model", "Land Use model")

evalplots <- list()
e.metrics <- list()
AUC <- c()
for (i in c(1:length(models))) {
  backg = randomPoints(predictorset[[i]], n=10000, ext = study_area_ext, extf=1.25)
  colnames(backg) <- c("lon","lat")
  group=kfold(backg, 5)
  backg_train <- backg[group!=1,]
  backg_test <- backg[group==1,]
  e = evaluate(pres_test, backg_test, models[[i]], predictorset[[i]])
  e.metrics[[i]] <- e
  # Extract TPR and FPR for ROC curve
  TPR <- e@TPR  # True Positive Rate
  FPR <- e@FPR  # False Positive Rate
  AUC <- round(e@auc, 3)
  # Store the ROC plot as a ggplot object in evalplots
  evalplots[[i]] <- ggplot() +
    geom_line(aes(x = FPR, y = TPR), color = "blue") +
    xlab("False Positive Rate") + 
    ylab("True Positive Rate") +
    labs(title = plot_titles[i],
         subtitle = paste("AUC =",AUC)) +
    theme_minimal()
}

```

# MaxEnt SDM: the models {.tabset}
## Predictor contributions for every model

 - *Full model*: This model includes all 19 BioClim variables all land use variables as predictors.
 - *BioClim model*: This model exclusively uses the BioClim variables to train the model.
 - *Land Use model*: This model exclusively uses the land use variables to train the model.<br>
 
 All models are employed with 10-fold cross validation. The ENMevaluate function from the ENMeval package was used to tune parameters: tune.args = list(fc = c("L", "Q", "LQ"), rm = seq(0.5, 2, 0.5)). <br>
 
 
**BioClimatic variables**<br>
BIO1 = Annual Mean Temperature<br>
BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))<br>
BIO3 = Isothermality (BIO2/BIO7) (×100)<br>
BIO4 = Temperature Seasonality (standard deviation ×100)<br>
BIO5 = Max Temperature of Warmest Month<br>
BIO6 = Min Temperature of Coldest Month<br>
BIO7 = Temperature Annual Range (BIO5-BIO6)<br>
BIO8 = Mean Temperature of Wettest Quarter<br>
BIO9 = Mean Temperature of Driest Quarter<br>
BIO10 = Mean Temperature of Warmest Quarter<br>
BIO11 = Mean Temperature of Coldest Quarter<br>
BIO12 = Annual Precipitation<br>
BIO13 = Precipitation of Wettest Month<br>
BIO14 = Precipitation of Driest Month<br>
BIO15 = Precipitation Seasonality (Coefficient of Variation)<br>
BIO16 = Precipitation of Wettest Quarter<br>
BIO17 = Precipitation of Driest Quarter<br>
BIO18 = Precipitation of Warmest Quarter<br>
BIO19 = Precipitation of Coldest Quarter

**Land use variables**<br>
Values are the fraction of a landcover class in each cell. The values are derived from the ESA WorldCover data set at 0.3-seconds resolution.<br>
Landcover variables:<br>
  - Trees<br>
  - Grassland<br> 
  - Shrubs<br>
  - Water<br>
  - Wetland<br>
  - Bare<br>
  - Built<br>

more info on: https://github.com/rspatial/geodata

```{r Plotting predictor contributions, include=TRUE, echo=FALSE, dev='jpeg'}

par(mfrow = c(3, 1), mar = c(2, 10, 1.5, 1))
terra::plot(model_full, main = "Full model") 
terra::plot(model_bioclim, main = "BioClim model")
terra::plot(model_landuse, main = "Land Use model")
par(mfrow = c(1,1))

```

Predictors within every model have been purged for correlation if r > 0.75. When pairwise correlation is detected above this threshold, one of the variables within the pair is omitted before further processing.


### Partial response curve

For each predictor the fitted species-environment relationship is plotted along the entire gradient while keeping the other predictors at their mean value. Partial response plots are oversimplifications and therefore hard to interpret as they do not represent the full species-environment relationship.

#### Full model

```{r Plotting predictor response curve: full model, include=TRUE, echo=FALSE, dev='jpeg'}

# Function to plot response curves for a given model
plot_response_curves <- function(model) {
  # Set up the plotting layout
  par(mfrow = c(3, ceiling(length(names(model@presence)) / 3)), mar = c(4, 4, 2, 1))
  
  # Loop through predictors and plot response curves with titles
  for (var in names(model@presence)) {
    dismo::response(model, var = var, range = 'p', col = "grey", lwd = 4, expand = 0)
    # title(main = var)
  }
  # Reset plotting layout back to default
  par(mfrow = c(1, 1))
}
plot_response_curves(model_full)

```

#### BioClim model

```{r Plotting predictor response curve: BioClim model, include=TRUE, echo=FALSE, dev='jpeg'}
plot_response_curves(model_bioclim)
```

#### Land Use model

```{r Plotting predictor response curve: Land Use model, include=TRUE, echo=FALSE, dev='jpeg'}
plot_response_curves(model_landuse)
```


## Model validation

In order to validate the models and make a comparison between the three models trained on different predictor sets, 75% of the data is kept to train the MaxEnt (Maximum Entropy) models, while the remaining 25% is kept aside for model validation.
The ROC curve is drawn by calculating the true positive rate (TPR) and false positive rate (FPR) at every possible threshold. A perfect model has a TPR of 1.0 and a FPR of 0.0. The larger the Area Under the Curve (AUC), the better a model performs.

```{r Plotting model accuracy, include=TRUE, echo=FALSE, dev='jpeg'}
cowplot::plot_grid(plotlist = evalplots, ncol = 3)
```

Metrics for model evaluation:

```{r Print model evaluation metrics, include=TRUE, echo=FALSE}

names(e.metrics) <- c("Full model","BioClim model","Land use model")
print(e.metrics[1:length(models)])

```


```{r Maxent predictions, include=TRUE, echo=FALSE}

#### Create a map of the predicted distributions in the present and future ####
# Extent of predictions
bbox.plot <- c(as.numeric(params$Plot.xmin), as.numeric(params$Plot.xmax), 
                as.numeric(params$Plot.ymin), as.numeric(params$Plot.ymax))
ext.pred <- ""
if (length(bbox.plot) != 4 || any(is.na(bbox.plot))) {
  ext.pred <- study_area_ext + 2
} else {
  ext.pred <- raster::extent(bbox.plot)
}

# PRESENT SDM for all predictors, bioclim only, and landuse only
present_full <- terra::predict(object = terra::crop(predictors_cleaned, ext.pred), model = model_full, cores = params$Ncores, na.rm=TRUE)
present_bioclim <- terra::predict(object = terra::crop(worldclim_data_aligned_cleaned, ext.pred), model = model_bioclim, cores = params$Ncores, na.rm=TRUE)
present_landuse <- terra::predict(object = terra::crop(landuse_rasterbrick_cleaned, ext.pred), model = model_landuse, cores = params$Ncores, na.rm=TRUE)

# FUTURE SDM for future BioClim predictors, future BioClim combined with contemporary land use
fut_bioclim_pred <- terra::predict(object = terra::crop(fut_bioclim_cleaned, ext.pred), model = model_bioclim, cores = params$Ncores, na.rm=TRUE)
fut2_bioclim_pred <- terra::predict(object = terra::crop(fut2_bioclim_cleaned, ext.pred), model = model_bioclim, cores = params$Ncores, na.rm=TRUE)
fut_full_pred <- terra::predict(object = terra::crop(fut_predictors_cleaned, ext.pred), model = model_full, cores = params$Ncores, na.rm=TRUE)
fut2_full_pred <- terra::predict(object = terra::crop(fut2_predictors_cleaned, ext.pred), model = model_full, cores = params$Ncores, na.rm=TRUE)

```


# MaxEnt SDM: plots {.tabset}

## Present

### SDM using BioClim and land use data

```{r cropping and resampling the species range map to map the predictions raster, include=TRUE, echo=FALSE}

# crop to extent and resample the range map
if (params$Range == "YES") {
range_raster_resampled <- resample(range, present_full, method = "ngb")

present_full <- terra::mask(present_full, range_raster_resampled)
present_bioclim <- terra::mask(present_bioclim, range_raster_resampled)
present_landuse <- terra::mask(present_landuse, range_raster_resampled)

fut_full_pred <- terra::mask(fut_full_pred, range_raster_resampled)
fut2_full_pred <- terra::mask(fut2_full_pred, range_raster_resampled)
fut_bioclim_pred <- terra::mask(fut_bioclim_pred, range_raster_resampled)
fut2_bioclim_pred <- terra::mask(fut2_bioclim_pred, range_raster_resampled)
}

```


```{r plotting the present 1/3, include=TRUE, echo=FALSE, dev='jpeg'}
# Get elevation data (this code here is to add a elevation layer as the bottom layer if wanted)
# slope <- terra::terrain(elevation_rc, "slope", unit="radians")
# aspect <- terra::terrain(elevation_rc, "aspect", unit="radians")
# hill <- terra::shade(slope, aspect, 45, 270)
# all_countries <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")

# plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4),
#      xlim = c(ext.pred@xmin, ext.pred@xmax),
#      ylim = c(ext.pred@ymin, ext.pred@ymax))
raster::plot(present_full, col = terrain.colors(25, alpha = 1, rev = TRUE))
title(bquote('MaxEnt full model ' * italic(.(params$Species))))
# shapefile of countries
plot(all_countries, add = TRUE, col = "transparent", border = "black")
# add occurrences with correct projection
plot(spdf, add = TRUE, pch = 16, cex = 0.3, col = rgb(red = 0, green = 0, blue = 1, alpha = 0.3))
```

### SDM using BioClim data only


```{r plotting the present 2/3, include=TRUE, echo=FALSE, dev='jpeg'}
# plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4),
#      xlim = c(ext.pred@xmin, ext.pred@xmax),
#      ylim = c(ext.pred@ymin, ext.pred@ymax))
raster::plot(present_bioclim, col = terrain.colors(25, alpha = 1, rev = TRUE))
title(bquote('MaxEnt BioClim model ' * italic(.(params$Species))))
# shapefile of countries
plot(all_countries, add = TRUE, col = "transparent", border = "black")
# add occurrences with correct projection
plot(spdf, add = TRUE, pch = 16, cex = 0.3, col = rgb(red = 0, green = 0, blue = 1, alpha = 0.3))
```

### SDM using land use data only


```{r plotting the present 3/3, include=TRUE, echo=FALSE, dev='jpeg'}
# plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4),
#      xlim = c(ext.pred@xmin, ext.pred@xmax),
#      ylim = c(ext.pred@ymin, ext.pred@ymax))
raster::plot(present_landuse, col = terrain.colors(25, alpha = 1, rev = TRUE))
title(bquote('MaxEnt landuse model ' * italic(.(params$Species))))
# shapefile of countries
plot(all_countries, add = TRUE, col = "transparent", border = "black")
# add occurrences with correct projection
plot(spdf, add = TRUE, pch = 16, cex = 0.3, col = rgb(red = 0, green = 0, blue = 1, alpha = 0.3))

```

## Future

### SDM using future BioClim data and contemporary land use data combined


```{r plotting the future 1/2, include=TRUE, echo=FALSE, dev='jpeg'}
# plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4),
#      xlim = c(ext.pred@xmin, ext.pred@xmax),
#      ylim = c(ext.pred@ymin, ext.pred@ymax))
raster::plot(fut_full_pred, col = terrain.colors(25, alpha = 1, rev = TRUE))
title(bquote('MaxEnt full model ' * italic(.(params$Species)) * ' : 2041-2060'))
plot(all_countries, add = TRUE, col = "transparent", border = "black")
plot(spdf, add = TRUE, pch = 16, cex = 0.3, col = rgb(red = 0, green = 0, blue = 1, alpha = 0.3))

# plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4),
#     xlim = c(ext.pred@xmin, ext.pred@xmax),
#      ylim = c(ext.pred@ymin, ext.pred@ymax))
raster::plot(fut2_full_pred, col = terrain.colors(25, alpha = 1, rev = TRUE))
title(bquote('MaxEnt full model ' * italic(.(params$Species)) * ' : 2061-2080'))
plot(all_countries, add = TRUE, col = "transparent", border = "black")
plot(spdf, add = TRUE, pch = 16, cex = 0.3, col = rgb(red = 0, green = 0, blue = 1, alpha = 0.3))
```


### SDM using future BioClim data


```{r plotting the future 2/2, include=TRUE, echo=FALSE, dev='jpeg'}
# plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4),
#      xlim = c(ext.pred@xmin, ext.pred@xmax),
#      ylim = c(ext.pred@ymin, ext.pred@ymax))
raster::plot(fut_bioclim_pred, col = terrain.colors(25, alpha = 1, rev = TRUE))
title(bquote('MaxEnt BioClim model ' * italic(.(params$Species)) * ' : 2041-2060'))
plot(all_countries, add = TRUE, col = "transparent", border = "black")
plot(spdf, add = TRUE, pch = 16, cex = 0.3, col = rgb(red = 0, green = 0, blue = 1, alpha = 0.3))

# plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4),
#     xlim = c(ext.pred@xmin, ext.pred@xmax),
#      ylim = c(ext.pred@ymin, ext.pred@ymax))
raster::plot(fut2_bioclim_pred, col = terrain.colors(25, alpha = 1, rev = TRUE))
title(bquote('MaxEnt BioClim model ' * italic(.(params$Species)) * ' : 2061-2080'))
plot(all_countries, add = TRUE, col = "transparent", border = "black")
plot(spdf, add = TRUE, pch = 16, cex = 0.3, col = rgb(red = 0, green = 0, blue = 1, alpha = 0.3))
```

## SDM change

### Change in climatic niche

Change in climatic niche availability using BioClim predictors between the current climate and future climate. Increase in climatic niche size is indicated in blue, reduction is indicated in red. 
Additionally, a trend plot is shown that indicates the change in climatic suitability. For every grid cell of 5x5km that has at least one species record, the average value of change is calculated. Density plots give context to the variability across sites. Averaging values across grid cells mitigates biased results towards extensively visited regions.


```{r plotting the change, include=TRUE, echo=FALSE, dev='jpeg'}
# Calculate change: negative = loss, positive = gain
change <- fut_bioclim_pred-present_bioclim
change2 <- fut2_bioclim_pred-present_bioclim
colrange <- colorRampPalette(c("darkred", "white", "darkblue"))

# plot change raster object 
terra::plot(change, col = colrange(255), range = c(-1, 1))
title(bquote('Response to climate change ' * italic(.(params$Species)) * ' : 2041-2060'))
plot(all_countries, add = TRUE, col = "transparent", border = "black")
# plot change2 raster object
raster::plot(change2, col = colrange(255), range = c(-1, 1))
title(bquote('Response to climage change ' * italic(.(params$Species)) * ' : 2061-2080'))
plot(all_countries, add = TRUE, col = "transparent", border = "black")


# Desired cell size in kilometers
desired_resolution_km <- 5
# Calculate aggregation factor
original_res <- terra::res(change) # Resolution in degrees (x, y)
original_res_km <- original_res * 111.32 # Convert to kilometers (approximation)
agg_factor <- desired_resolution_km / original_res_km
# Aggregate rasters to 5x5 km resolution
agg_change <- terra::aggregate(change, fact = agg_factor, fun = mean, na.rm = TRUE)
agg_change2 <- terra::aggregate(change2, fact = agg_factor, fun = mean, na.rm = TRUE)
# Mask aggregated rasters to cells with occurrences
occurrences <- terra::vect(spdf) # Spatial points for occurrences
occ_mask_change <- terra::mask(agg_change, occurrences)
occ_mask_change2 <- terra::mask(agg_change2, occurrences)
# Calculate mean change values directly from the masked layers
mean_change <- terra::global(occ_mask_change, fun = mean, na.rm = TRUE)[[1]]
mean_change2 <- terra::global(occ_mask_change2, fun = mean, na.rm = TRUE)[[1]]
# Create data frame for plotting
trend_data <- data.frame(
  Time = factor(c("Present", "2041-2060", "2061-2080"), levels = c("Present", "2041-2060", "2061-2080")),
  MeanChange = c(0, mean_change, mean_change2)
)

# Trend plot
trend_plot <- ggplot(data = trend_data, aes(x = Time, y = MeanChange, group = 1)) +
  geom_line(size = 2, color = "lightgrey") +
  geom_point(size = 4, color = "tomato") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_text(
    data = trend_data, # Add labels only to future time points
    aes(label = round(MeanChange, 2)), # Show rounded values
    nudge_y = 0.05, # Adjust position vertically
    size = 4,
    color = ifelse(trend_data$MeanChange == 0, "white", "black")
  ) +
  labs(
    x = "Time Period",
    y = "Mean Change in Bioclimatic Suitability",
    title = bquote(italic(.(params$Species)))
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14)
  )

# Density plots
# Future 1
change_values <- na.omit(values(occ_mask_change))
data <- data.frame(change_values)
density_plot_2041 <- ggplot(data, aes(x = lyr1)) +
  geom_density(
    data = subset(data, lyr1 < 0), 
    aes(x = lyr1), 
    fill = "tomato", alpha = 0.6
  ) +
  geom_density(
    data = subset(data, lyr1 >= 0), 
    aes(x = lyr1), 
    fill = "skyblue", alpha = 0.6
  ) +
  geom_point(
    aes(x = lyr1, y = 0), 
    position = position_jitter(height = 0), 
    color = "darkblue",
    alpha = 0.6
  ) +
  labs(
    x = "Change in Suitability",
    y = "Density",
    title = "Density of Suitability Changes (2041-2060)"
  ) +
  theme_minimal(base_size = 8)


# Future 2
change_values <- na.omit(values(occ_mask_change2))
data <- data.frame(change_values)
density_plot_2061 <- ggplot(data, aes(x = lyr1)) +
  geom_density(
    data = subset(data, lyr1 < 0), 
    aes(x = lyr1), 
    fill = "tomato", alpha = 0.6
  ) +
  geom_density(
    data = subset(data, lyr1 >= 0), 
    aes(x = lyr1), 
    fill = "skyblue", alpha = 0.6
  ) +
  geom_point(
    aes(x = lyr1, y = 0), 
    position = position_jitter(height = 0), 
    color = "darkblue",
    alpha = 0.6
  ) +
  labs(
    x = "Change in Suitability",
    y = "Density",
    title = "Density of Suitability Changes (2061-2080)"
  ) +
  theme_minimal(base_size = 8)

# Display plots
species_trend <- plot_grid(trend_plot, plot_grid(density_plot_2041, density_plot_2061, ncol = 1, nrow = 2), ncol = 2, rel_widths = c(1, 1))
plot(species_trend)

```


## Limiting variables

The variable most responsible for unfavorable conditions in respect of the species' climatic niche in the future (2041-2060) is plotted for every grid cell within the estimated species' range. 
<br>
<br>

```{r plotting the limiting factors, include=TRUE, echo=FALSE, dev='jpeg'}
lim <- rmaxent::limiting(raster::brick(terra::crop(fut_bioclim_cleaned, ext.pred)), model_bioclim)
rasterVis::levelplot(lim, col.regions = rainbow, main = "Limiting variables under future climatic conditions (2041-2060)") 

# Limiting predictors
lim2 <- rmaxent::limiting(raster::brick(terra::crop(fut2_bioclim_cleaned, ext.pred)), model_bioclim)
rasterVis::levelplot(lim2, col.regions = rainbow, main = "Limiting variables under future climatic conditions (2061-2080)")

```

This table shows the relative contribution of every BioClim predictor to the species' climatic niche contraction where the strength of range reducion is > -0.1. (present vs 2041-2060). The relative contribution is estimated by: <br>
1) Extracting the limiting predictor variables that overlap with a value lower than -0.1 for range contraction.<br>
2) Overlaying the species' spatial occurrence points with the extracted raster image from step 1 and count the relative
   contribution of every predictor variable.
<br>

```{r plotting a table for limiting factors, include=TRUE, echo=FALSE}
# Estimate which limiting variables contribute to the species' climatic niche contraction
Limiting_Predictors <- function(ChangeRaster, LimitRaster, SpatPoints, NegChange_limit = -0.1) {
  # Identify areas of negative change
  negative_change <- ChangeRaster < NegChange_limit
  negative_change <- raster::raster(negative_change)
  # Mask the LimitRaster to show only negative change areas
  lim_masked <- LimitRaster
  lim_masked[!negative_change] <- NA
  # Extract raster values at each spatial point
  spdf_values <- extract(lim_masked, SpatPoints)
  # Create a summary data frame with counts of each raster value
  summary_df <- table(spdf_values) %>% as.data.frame()
  # Retrieve the attribute table with predictor names
  attribute_table <- lim_masked@data@attributes[[1]]
  # Map raster values to predictors
  summary_df$spdf_values <- setNames(attribute_table$predictor[match(summary_df$spdf_values, attribute_table$ID)], 
                                     summary_df$spdf_values)
  # Rename columns for clarity
  names(summary_df) <- c("Predictor", "Freq")
  # Add relative percentage column and order by frequency
  summary_df <- summary_df %>%
    mutate(Perc = round((as.numeric(Freq) / sum(as.numeric(Freq))) * 100, 2)) %>%
    arrange(desc(Freq))
  # Output the summary table with kable for a clear display
  kable(summary_df, caption = "Relative Contribution of Predictors for Negative Change Areas") %>%
  kable_paper("hover", full_width = F)
}
Limiting_Predictors(ChangeRaster = change, LimitRaster = lim, SpatPoints = spdf, NegChange_limit = -0.1)

```

```{r exporting raster objects if ExportRaster is TRUE, include=TRUE, echo=FALSE}
# Create directory
if (!dir.exists(file.path(params$Workdir, "Raster_objects")) & params$ExportRaster == "YES") {  dir.create(file.path(params$Workdir,"Raster_objects"), recursive = TRUE)}

# Export relevant raster objects
if (params$ExportRaster == "YES") {
  setwd(dir = file.path(params$Workdir,"Raster_objects"))
    # Present: full model (consistently has the highest AUC)
    terra::writeRaster(x = present_full, filename = paste0(params$Species, "_Present_FullModel.tif"), overwrite = TRUE)
    # Future 1: BioClim
    terra::writeRaster(x = fut_bioclim_pred, filename = paste0(params$Species, "_Future_BioClim.tif"), overwrite = TRUE)
    # Future 2: BioClim
    terra::writeRaster(x = fut2_bioclim_pred, filename = paste0(params$Species, "_Future2_BioClim.tif"), overwrite = TRUE)
    # Change 1
    terra::writeRaster(x = change, filename = paste0(params$Species, "_Change.tif"), overwrite = TRUE)
    # Change 2
    terra::writeRaster(x = change2, filename = paste0(params$Species, "_Change2.tif"), overwrite = TRUE)
  setwd(dir = params$Workdir)
}
```
<br>
<br>
<br>
<br>