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 METHODOLOGICAL GRAPHICS (EXTRACTION & PLOTTING)
# ==============================================================================
# Purpose: This script extracts climatic variables at occurrence points and 
# generates VERA proof-of-concept graphics (Gaussian vs. Quadratic logic).
# Outputs: Lc_with_Env.csv and high-resolution PNG graphics.
# ==============================================================================

library(terra)
library(ggplot2)
library(dplyr)
library(gridExtra)

# --- 1. GLOBAL DIRECTORY SETUP ---
# Ensure paths are correct for your local C: drive structure
base_path <- "C:/DATA"
path_vars <- file.path(base_path, "Variables_Lc")
path_occ  <- 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) - L. canadensis
# ==============================================================================
message(">>> PART 1: Starting L. canadensis processing...")

# 1.1. Data Loading & Extraction
vars_lc <- c("minTempWarmest", "PETColdestQuarter", "PETDriestQuarter", 
             "PETseasonality", "wc2.1_2.5m_bio_3", "wc2.1_2.5m_bio_8")

layers_lc <- list.files(path_vars, pattern = "\\.tif$", full.names = TRUE)

# Safety check: Stop if no tiffs found
if(length(layers_lc) == 0) stop("ERROR: No .tif files found in: ", path_vars)

stack_lc <- rast(layers_lc)[[vars_lc]]

# Load occurrences and extract climatic values
occ_file <- file.path(path_occ, "Lc.csv")
if(!file.exists(occ_file)) stop("ERROR: Occurrence file not found: ", occ_file)

occ_lc <- read.csv(occ_file)
env_lc <- terra::extract(stack_lc, occ_lc[, c("X_WGS84", "Y_WGS84")])

# Clean data: Drop ID and remove NAs
df_lc <- cbind(occ_lc, env_lc[, -1]) 
df_lc <- na.omit(df_lc)

# Export extracted data
write.csv(df_lc, file.path(path_out, "Lc_with_Env.csv"), row.names = FALSE)
message(">>> Lc_with_Env.csv saved successfully to: ", path_out)

# 1.2. VERA Plotting Function (Gaussian vs. Quadratic Logic)
dict_names <- list(
  "minTempWarmest"    = "Min Temp of Warmest Month (MTWM)",
  "PETColdestQuarter" = "PET Coldest Quarter (PETCoQ)",
  "PETDriestQuarter"  = "PET Driest Quarter (PETDrQ)",
  "PETseasonality"    = "PET Seasonality (PETs)",
  "wc2.1_2.5m_bio_3"  = "Isothermality (Bio3)",
  "wc2.1_2.5m_bio_8"  = "Mean Temp of Wettest Quarter (Bio8)"
)

plot_vera_logic <- function(df, var_col, formal_name) {
  mu    <- mean(df[[var_col]], na.rm = TRUE)
  sigma <- sd(df[[var_col]], na.rm = TRUE)
  
  # Generate theoretical range for curves
  x_range <- seq(min(df[[var_col]]) - (1.5 * sigma), 
                 max(df[[var_col]]) + (1.5 * sigma), length.out = 300)
  
  theo_df <- data.frame(
    x = x_range,
    suitability = dnorm(x_range, mean = mu, sd = sigma),
    restriction = ((x_range - mu) / sigma)^2
  )
  
  # Colors: Suitability (Blue), Restriction (Red), Optimum (Green)
  color_suit <- "#2c7bb6" 
  color_restr <- "#d7191c" 
  color_opt   <- "#1a9641" 
  
  # Panel A: Gaussian Suitability
  p_gauss <- ggplot() +
    geom_line(data = theo_df, aes(x = x, y = suitability), color = color_suit, linewidth = 1.2) +
    geom_vline(xintercept = mu, linetype = "dashed", color = color_opt, linewidth = 1) +
    labs(title = "A) Gaussian Suitability Profile", 
         subtitle = paste("Climatic Optimum =", round(mu, 2)), 
         x = formal_name, y = "Probability Density") +
    theme_minimal()
  
  # Panel B: VERA Quadratic Restriction
  p_vera <- ggplot() +
    geom_line(data = theo_df, aes(x = x, y = restriction), color = color_restr, linewidth = 1.2) +
    geom_point(data = df, aes(x = !!sym(var_col), y = ((!!sym(var_col) - mu) / sigma)^2), 
               alpha = 0.2, color = color_suit, size = 1.2) +
    geom_vline(xintercept = mu, linetype = "dashed", color = color_opt, linewidth = 1) +
    labs(title = "B) VERA Quadratic Restriction Intensity", 
         subtitle = expression(Z^2 == ((Value - Mean) / SD)^2), 
         x = formal_name, y = "Restriction Intensity (Stress)") +
    theme_minimal()
  
  return(arrangeGrob(p_gauss, p_vera, ncol = 2))
}

# 1.3. Execution Loop for Graphics
message(">>> Generating VERA methodological proof graphics...")

for(var in names(dict_names)) {
  png_name <- paste0("VERA_Proof_Lc_", var, ".png")
  message(">>> Processing: ", var)
  
  g <- plot_vera_logic(df_lc, var, dict_names[[var]])
  
  ggsave(filename = file.path(path_out, png_name), plot = g, 
         width = 12, height = 6, dpi = 300, units = "in")
}

message(">>> All L. canadensis VERA graphics generated successfully!")




# ==============================================================================
# PART 2: BOBCAT (Lynx rufus) - L. rufus (FIXED VERSION)
# ==============================================================================
message(">>> PART 2: Starting L. rufus processing...")

# --- 2.1. Directory Setup ---
path_vars_lr <- file.path(base_path, "DATA/Variables_Lr")

# --- 2.2. Data Loading & Extraction ---
# ÖNEMLİ: Dosya isimlerini klasördekiyle birebir eşle (2.5m olduğuna emin ol)
vars_lr <- c("PETColdestQuarter", "minTempWarmest", "PETWarmestQuarter", 
             "wc2.1_2.5m_bio_19", "wc2.1_2.5m_bio_7", "wc2.1_2.5m_bio_13")

layers_lr <- list.files(path_vars_lr, pattern = "\\.tif$", full.names = TRUE)

if(length(layers_lr) == 0) stop("ERROR: No .tif files found in: ", path_vars_lr)

stack_lr <- rast(layers_lr)

# Hata Önleyici: Sadece mevcut olan katmanları seç
available_vars <- intersect(vars_lr, names(stack_lr))
if(length(available_vars) < length(vars_lr)){
  message("WARNING: Missing layers! Found: ", paste(available_vars, collapse=", "))
}
stack_lr <- stack_lr[[available_vars]]

# Load occurrences
occ_file_lr <- file.path(path_occ, "Lr.csv")
occ_lr <- read.csv(occ_file_lr)
env_lr <- terra::extract(stack_lr, occ_lr[, c("X_WGS84", "Y_WGS84")])

# Data Cleaning
df_lr <- cbind(occ_lr, env_lr[, -1]) 
df_lr <- na.omit(df_lr)

# DEBUG: Sütun isimlerini kontrol et (Hata alırsan buraya bak)
message(">>> Extracted columns for Lr: ", paste(colnames(df_lr), collapse=", "))

# --- 2.3. Dictionary (TYPO FIXED: 2.1_2.1 -> 2.1_2.5m) ---
dict_lr <- list(
  "PETColdestQuarter" = "PET Coldest Quarter (PETCoQ)",
  "minTempWarmest"    = "Min Temp of Warmest Month (MTWM)",
  "PETWarmestQuarter" = "PET Warmest Quarter (PETWaQ)",
  "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)"
)

# --- 2.4. Graphics Generation Loop ---
message(">>> Generating VERA methodological proof graphics for L. rufus...")

for(var in names(dict_lr)) {
  # Ekstra Kontrol: Değişken df içinde var mı?
  if(!var %in% colnames(df_lr)){
    message("!!! SKIPPING: Variable '", var, "' not found in dataframe.")
    next
  }
  
  png_name <- paste0("VERA_Proof_Lr_", var, ".png")
  message(">>> Processing: ", var)
  
  # Part 1'deki plot_vera_logic fonksiyonunu kullanır
  g <- plot_vera_logic(df_lr, var, dict_lr[[var]])
  
  ggsave(filename = file.path(path_out, png_name), plot = g, 
         width = 12, height = 6, dpi = 300, units = "in")
}

message(">>> Process for L. rufus complete.")
