R-4.5.2 
SOFTWARE REQUIREMENTS
These scripts were developed and tested in R version 4.2.0+. 
The following libraries are required:
- terra (Spatial data handling)
- biomod2 (Species Distribution Modeling)
- humboldt (Niche overlap and divergence tests)
- dplyr, tidyr (Data manipulation)

# ==============================================================================
# VERA & MESS PROTOCOL: Spatial Diagnostics (MLV, MESS, and MoD)
# ==============================================================================
# This script calculates:
# 1. VERA Most Limiting Variable (MLV) based on niche-centroid distance.
# 2. Multivariate Environmental Similarity Surface (MESS).
# 3. Most Dissimilar Variable (MoD) representing the primary extrapolation driver.
# Outputs are saved directly to the target directory as .tif files.
# ==============================================================================

library(terra)
library(predicts) # Required for MESS and MoD analysis

# ------------------------------------------------------------------------------
# PART 1: CANADA LYNX (Lynx canadensis)
# ------------------------------------------------------------------------------
message(">>> PART 1: Starting spatial diagnostics for L. canadensis...")

# --- 1. Directory Setup ---
path_vars_lc <- "C:/DATA/Variables_Lc"
path_occ_lc  <- "C:/DATA/Occurrences/Lc.csv"
path_out     <- "C:/OUTPUTS"

if(!dir.exists(path_out)) dir.create(path_out, recursive = TRUE)

# Core Climatic Predictors for L. canadensis
vars_lc <- c("minTempWarmest", "PETColdestQuarter", "PETDriestQuarter", 
             "PETseasonality", "wc2.1_2.5m_bio_3", "wc2.1_2.5m_bio_8")

# --- 2. Data Loading ---
layers_lc <- list.files(path_vars_lc, pattern = ".tif$", full.names = TRUE)
stack_lc  <- rast(layers_lc)[[vars_lc]]

occ_lc <- read.csv(path_occ_lc)
occ_vals_lc <- terra::extract(stack_lc, occ_lc[, c("X_WGS84", "Y_WGS84")])
occ_vals_lc <- na.omit(occ_vals_lc[, -1])

# --- 3. VERA: Z-Score (Penalty) Calculation ---
means_lc <- colMeans(occ_vals_lc)
sds_lc   <- apply(occ_vals_lc, 2, sd)

stack_z_lc <- stack_lc
for(i in 1:length(vars_lc)) {
  stack_z_lc[[i]] <- ((stack_lc[[i]] - means_lc[i]) / sds_lc[i])^2
}
names(stack_z_lc) <- vars_lc

# --- 4. VERA: MLV Identification & Export ---
vera_mlv_lc <- which.max(stack_z_lc)
levels(vera_mlv_lc) <- data.frame(ID = 1:length(vars_lc), Variable = vars_lc)

writeRaster(vera_mlv_lc, 
            filename = file.path(path_out, "Lc_MLV_VERA.tif"), 
            overwrite = TRUE, datatype = "INT1U", gdal = c("COMPRESS=LZW", "TFW=YES"))
# --- 5. MESS & MoD Analysis ---
message(">>> Calculating MESS and MoD for L. canadensis...")
# Compute full MESS (includes individual variable surfaces and the composite index)
mess_full_lc <- predicts::mess(x = stack_lc, v = occ_vals_lc, full = TRUE)

# 1. Extract Global MESS Index (The final 'mess' layer)
mess_index_lc <- mess_full_lc[["mess"]]

# 2. Calculate MoD (Most Dissimilar Variable)
# Identifies the variable (1 to 6) with the minimum MESS value for each pixel
mod_index_lc <- which.min(mess_full_lc[[1:length(vars_lc)]])

# Assign variable names as levels to the MoD raster for categorical mapping
levels(mod_index_lc) <- data.frame(ID = 1:length(vars_lc), Variable = vars_lc)

# --- 6. Export Results ---
# Export MESS Index map with LZW compression
writeRaster(mess_index_lc, 
            filename = file.path(path_out, "Lc_MESS_Index.tif"), 
            overwrite = TRUE, gdal = c("COMPRESS=LZW", "TFW=YES"))

# Export MoD Variable map as an 8-bit unsigned integer raster
writeRaster(mod_index_lc, 
            filename = file.path(path_out, "Lc_MoD_MESS.tif"), 
            overwrite = TRUE, datatype = "INT1U", gdal = c("COMPRESS=LZW", "TFW=YES"))

message(">>> Lynx canadensis MLV, MESS, and MoD TIFs successfully generated!")



# ------------------------------------------------------------------------------
# PART 2: BOBCAT (Lynx rufus)
# ------------------------------------------------------------------------------
message(">>> PART 2: Starting spatial diagnostics for L. rufus...")

# --- 1. Directory Setup ---
path_vars_lr <- "C:/DATA/Variables_Lr"
path_occ_lr  <- "C:/DATA/Occurrences/Lr.csv"

# Core Climatic Predictors for L. rufus
vars_lr <- c("PETColdestQuarter", "minTempWarmest", "PETWarmestQuarter", 
             "wc2.1_2.5m_bio_19", "wc2.1_2.5m_bio_7", "wc2.1_2.5m_bio_13")

# --- 2. Data Loading ---
layers_lr <- list.files(path_vars_lr, pattern = ".tif$", full.names = TRUE)
stack_lr  <- rast(layers_lr)[[vars_lr]]

occ_lr <- read.csv(path_occ_lr)
occ_vals_lr <- terra::extract(stack_lr, occ_lr[, c("X_WGS84", "Y_WGS84")])
occ_vals_lr <- na.omit(occ_vals_lr[, -1])

# --- 3. VERA: Z-Score (Penalty) Calculation ---
means_lr <- colMeans(occ_vals_lr)
sds_lr   <- apply(occ_vals_lr, 2, sd)

stack_z_lr <- stack_lr
for(i in 1:length(vars_lr)) {
  stack_z_lr[[i]] <- ((stack_lr[[i]] - means_lr[i]) / sds_lr[i])^2
}
names(stack_z_lr) <- vars_lr

# --- 4. VERA: MLV Identification & Export ---
vera_mlv_lr <- which.max(stack_z_lr)
levels(vera_mlv_lr) <- data.frame(ID = 1:length(vars_lr), Variable = vars_lr)

writeRaster(vera_mlv_lr, 
            filename = file.path(path_out, "Lr_MLV_VERA.tif"), 
            overwrite = TRUE, datatype = "INT1U", gdal = c("COMPRESS=LZW", "TFW=YES"))

# --- 5. MESS & MoD Analysis ---
message(">>> Calculating MESS and MoD for L. rufus...")
# Compute full MESS (includes individual variable surfaces and the composite index)
mess_full_lr <- predicts::mess(x = stack_lr, v = occ_vals_lr, full = TRUE)

# 1. Extract Global MESS Index (The final 'mess' layer)
mess_index_lr <- mess_full_lr[["mess"]]

# 2. Calculate MoD (Most Dissimilar Variable)
# Identifies the variable index (1 to 6) with the minimum MESS value for each pixel
mod_index_lr <- which.min(mess_full_lr[[1:length(vars_lr)]])

# Assign variable names as levels to the MoD raster for categorical mapping
levels(mod_index_lr) <- data.frame(ID = 1:length(vars_lr), Variable = vars_lr)

# --- 6. Export Results ---
# Export MESS Index map with LZW compression and TFW world file
writeRaster(mess_index_lr, 
            filename = file.path(path_out, "Lr_MESS_Index.tif"), 
            overwrite = TRUE, gdal = c("COMPRESS=LZW", "TFW=YES"))

# Export MoD Variable map as an 8-bit unsigned integer raster
writeRaster(mod_index_lr, 
            filename = file.path(path_out, "Lr_MoD_MESS.tif"), 
            overwrite = TRUE, datatype = "INT1U", gdal = c("COMPRESS=LZW", "TFW=YES"))

message(">>> Lynx rufus MLV, MESS, and MoD TIFs successfully generated!")
