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)

# ==============================================================================
# SPATIAL EXTRACTION PROTOCOL: Occurrences & Environmental Predictors
# ==============================================================================
# Purpose: Extract climatic values from raster surfaces (TIFF) to species 
# occurrence points (CSV) and save the consolidated datasets.
# ==============================================================================

library(terra)

# --- 1. GLOBAL DIRECTORY SETUP ---
# Defining paths for raw occurrences and environmental variables
base_path <- "C:/DATA"
path_occ  <- file.path(base_path, "Occurrences")
path_lc_vars <- file.path(base_path, "Variables_Lc")
path_lr_vars <- file.path(base_path, "Variables_Lr")

# --- 2. DATA EXTRACTION: CANADA LYNX (Lc) ---
message(">>> Starting environmental extraction for L. canadensis...")

# Define target predictors for Lc
vars_lc <- c("minTempWarmest", "PETColdestQuarter", "PETDriestQuarter", 
             "PETseasonality", "wc2.1_2.5m_bio_3", "wc2.1_2.5m_bio_8")

# Load Lc rasters and occurrences
stack_lc <- rast(list.files(path_lc_vars, pattern = "\\.tif$", full.names = TRUE))[[vars_lc]]
occ_lc   <- read.csv(file.path(path_occ, "Lc.csv"))

# Perform spatial extraction
env_values_lc <- terra::extract(stack_lc, occ_lc[, c("X_WGS84", "Y_WGS84")])

# Combine coordinates with extracted values and remove NA records (sampling gaps)
df_lc_final <- cbind(occ_lc, env_values_lc[, -1])
df_lc_final <- na.omit(df_lc_final)

# Save the consolidated Lc dataset
write.csv(df_lc_final, file.path(path_occ, "Lc_with_Env.csv"), row.names = FALSE)
message(">>> Lc_with_Env.csv saved successfully.")


# --- 3. DATA EXTRACTION: BOBCAT (Lr) ---
message(">>> Starting environmental extraction for L. rufus...")

# Define target predictors for Lr
vars_lr <- c("PETColdestQuarter", "minTempWarmest", "PETWarmestQuarter", 
             "wc2.1_2.5m_bio_19", "wc2.1_2.5m_bio_7", "wc2.1_2.5m_bio_13")

# Load Lr rasters and occurrences
stack_lr <- rast(list.files(path_lr_vars, pattern = "\\.tif$", full.names = TRUE))[[vars_lr]]
occ_lr   <- read.csv(file.path(path_occ, "Lr.csv"))

# Perform spatial extraction
env_values_lr <- terra::extract(stack_lr, occ_lr[, c("X_WGS84", "Y_WGS84")])

# Combine coordinates with extracted values and remove NA records
df_lr_final <- cbind(occ_lr, env_values_lr[, -1])
df_lr_final <- na.omit(df_lr_final)

# Save the consolidated Lr dataset
write.csv(df_lr_final, file.path(path_occ, "Lr_with_Env.csv"), row.names = FALSE)
message(">>> Lr_with_Env.csv saved successfully.")

message(">>> ALL PROCESSES COMPLETE. Datasets are ready for analysis.")





# ==============================================================================
# STANDARDIZED NICHE DENSITY PROFILES (RIDGE PLOTS)
# ==============================================================================
# This script reads the extracted climatic values, standardizes them (Z-scores), 
# and generates stacked density profiles for both L. canadensis and L. rufus.
# ==============================================================================

library(ggplot2)
library(dplyr)
library(tidyr)
library(ggridges)

# --- DIRECTORY SETUP ---
base_path <- "C:/DATA"
path_in   <- file.path(base_path, "Occurrences")
path_out  <- file.path(base_path, "OUTPUTS")

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

# ==============================================================================
# PART 1: CANADA LYNX (Lynx canadensis)
# ==============================================================================
message(">>> Generating Niche Density Profile for L. canadensis...")

df_lc <- read.csv(file.path(path_in, "Lc_with_Env.csv"))

vars_lc <- c("minTempWarmest", "PETColdestQuarter", "PETDriestQuarter", 
             "PETseasonality", "wc2.1_2.5m_bio_3", "wc2.1_2.5m_bio_8")

df_scaled_lc <- as.data.frame(scale(df_lc[, vars_lc]))

df_long_lc <- df_scaled_lc %>%
  pivot_longer(cols = everything(), names_to = "Variable", values_to = "Z_Score") %>%
  mutate(Variable = factor(Variable, levels = rev(vars_lc)))

labels_lc <- c(
  "minTempWarmest"    = "Min Temp of Warmest Month",
  "PETColdestQuarter" = "PET Coldest Quarter",
  "PETDriestQuarter"  = "PET Driest Quarter",
  "PETseasonality"    = "PET Seasonality",
  "wc2.1_2.5m_bio_3"  = "Isothermality (Bio3)",
  "wc2.1_2.5m_bio_8"  = "Mean Temp of Wettest Quarter"
)

# Standardized color palette
colors_lc <- c(
  "minTempWarmest"    = "#e41a1c",
  "PETColdestQuarter" = "#377eb8",
  "PETDriestQuarter"  = "#4daf4a",
  "PETseasonality"    = "#984ea3",
  "wc2.1_2.5m_bio_3"  = "#ff7f00",
  "wc2.1_2.5m_bio_8"  = "#a65628"
)

p_lc <- ggplot(df_long_lc, aes(x = Z_Score, y = Variable, fill = Variable)) +
  geom_density_ridges(alpha = 0.85, color = "white", linewidth = 0.8, scale = 1.3) +
  scale_fill_manual(values = colors_lc) +
  scale_y_discrete(labels = labels_lc) +
  coord_cartesian(xlim = c(-3.5, 3.5)) +
  labs(
    title = "Standardized Niche Density Profile",
    subtitle = "Lynx canadensis",
    x = "Standardized Variable Values (Z-scores)",
    y = NULL
  ) +
  theme_minimal(base_family = "sans") +
  theme(
    legend.position = "none",
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5, color = "grey30", face = "italic"),
    axis.text.y = element_text(face = "bold", size = 11, color = "black"),
    axis.text.x = element_text(size = 10, color = "black"),
    axis.title.x = element_text(face = "bold", size = 12, margin = margin(t = 10)),
    panel.grid.major.x = element_line(color = "grey80", linetype = "dashed"),
    panel.grid.minor.x = element_blank()
  )

ggsave(file.path(path_out, "Lc_Niche_Density_Profile.png"), plot = p_lc, 
       width = 10, height = 7, dpi = 300, bg = "white")

# ==============================================================================
# PART 2: BOBCAT (Lynx rufus)
# ==============================================================================
message(">>> Generating Niche Density Profile for L. rufus...")

df_lr <- read.csv(file.path(path_in, "Lr_with_Env.csv"))

vars_lr <- c("PETColdestQuarter", "minTempWarmest", "PETWarmestQuarter", 
             "wc2.1_2.5m_bio_19", "wc2.1_2.5m_bio_7", "wc2.1_2.5m_bio_13")

df_scaled_lr <- as.data.frame(scale(df_lr[, vars_lr]))

df_long_lr <- df_scaled_lr %>%
  pivot_longer(cols = everything(), names_to = "Variable", values_to = "Z_Score") %>%
  mutate(Variable = factor(Variable, levels = rev(vars_lr)))

labels_lr <- c(
  "PETColdestQuarter" = "PET Coldest Quarter",
  "minTempWarmest"    = "Min Temp of Warmest Month",
  "PETWarmestQuarter" = "PET Warmest Quarter",
  "wc2.1_2.5m_bio_19" = "Precip of Coldest Quarter (Bio19)",
  "wc2.1_2.5m_bio_7"  = "Temperature Annual Range (Bio7)",
  "wc2.1_2.5m_bio_13" = "Precip of Wettest Month (Bio13)"
)

colors_lr <- c(
  "PETColdestQuarter" = "#377eb8",
  "minTempWarmest"    = "#e41a1c",
  "PETWarmestQuarter" = "#ff7f00",
  "wc2.1_2.5m_bio_19" = "#4daf4a",
  "wc2.1_2.5m_bio_7"  = "#984ea3",
  "wc2.1_2.5m_bio_13" = "#a65628"
)

p_lr <- ggplot(df_long_lr, aes(x = Z_Score, y = Variable, fill = Variable)) +
  geom_density_ridges(alpha = 0.85, color = "white", linewidth = 0.8, scale = 1.3) +
  scale_fill_manual(values = colors_lr) +
  scale_y_discrete(labels = labels_lr) +
  coord_cartesian(xlim = c(-3.5, 3.5)) +
  labs(
    title = "Standardized Niche Density Profile",
    subtitle = "Lynx rufus",
    x = "Standardized Variable Values (Z-scores)",
    y = NULL
  ) +
  theme_minimal(base_family = "sans") +
  theme(
    legend.position = "none",
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5, color = "grey30", face = "italic"),
    axis.text.y = element_text(face = "bold", size = 11, color = "black"),
    axis.text.x = element_text(size = 10, color = "black"),
    axis.title.x = element_text(face = "bold", size = 12, margin = margin(t = 10)),
    panel.grid.major.x = element_line(color = "grey80", linetype = "dashed"),
    panel.grid.minor.x = element_blank()
  )

ggsave(file.path(path_out, "Lr_Niche_Density_Profile.png"), plot = p_lr, 
       width = 10, height = 7, dpi = 300, bg = "white")

message(">>> Niche Density Profiles successfully generated and saved to OUTPUTS.")

