# =============================================================================
# Title:
# - Global Moran’s I (terra) with Multiple Normalized Neighborhood Configurations
#
# Author:
# - Oliver Gutiérrez-Hernández, Department of Geography, University of Málaga (UMA)
#
# Environment:
# - R version 4.5.1 (2025-10-22)
# - Package: terra 1.8-80
#
# Objective:
# - Evaluate spatial dependency (Moran’s I) in a raster dataset 
#   using different normalized spatial neighborhood sizes 
#   (3×3, 5×5, 7×7, 9×9, 11×11) consistent with terra defaults.
#
# Reference:
# - Moran, P. A. P. (1950). "Notes on Continuous Stochastic Phenomena". 
#   Biometrika, 37(1–2), 17–23. https://doi.org/10.1093/biomet/37.1-2.17
# =============================================================================


# --- Package check and loading ------------------------------------------------
install_if_missing <- function(package) {
  if (!requireNamespace(package, quietly = TRUE)) {
    install.packages(package, dependencies = TRUE)
  } else {
    cat(paste("Package", package, "is already installed.\n"))
  }
}

install_if_missing("terra")
library(terra)


# --- Load raster data ---------------------------------------------------------
raster_path <- "C:/data/p-values.tif"   # ⚠️ Replace with your actual file path
raster_pvalues <- rast(raster_path)

cat("Total number of pixels:", ncell(raster_pvalues), "\n")

# Validate raster content
p_values <- values(raster_pvalues)
if (is.null(p_values)) stop("❌ The raster contains no valid values or is improperly loaded.")

# Count valid and invalid pixels
num_valid_pixels <- sum(!is.na(p_values))
num_invalid_pixels <- sum(is.na(p_values))
cat("Valid pixels:", num_valid_pixels, "\n")
cat("Invalid pixels:", num_invalid_pixels, "\n")

# Replace NA values with mean of valid pixels
if (num_valid_pixels > 0) {
  p_values[is.na(p_values)] <- mean(p_values, na.rm = TRUE)
  values(raster_pvalues) <- p_values
} else {
  stop("❌ The raster contains no valid pixels to process.")
}


# --- Global Moran's I for multiple neighborhoods ------------------------------
window_sizes <- c(3, 5, 7, 9, 11)

results <- data.frame(
  Window = character(),
  Moran_I = numeric(),
  P_value = numeric(),
  stringsAsFactors = FALSE
)

for (w_size in window_sizes) {
  cat("\n--- Computing Moran's I with", w_size, "×", w_size, "window ---\n")
  
  # Create normalized Queen’s case kernel (center = 0)
  w_matrix <- matrix(1, nrow = w_size, ncol = w_size)
  center <- ceiling(w_size / 2)
  w_matrix[center, center] <- 0
  w_matrix <- w_matrix / sum(w_matrix)
  
  # Safety check: skip if raster too small or constant
  if (ncell(raster_pvalues) < w_size^2 || var(values(raster_pvalues), na.rm = TRUE) == 0) {
    cat("⚠️ Skipped:", w_size, "×", w_size, "- insufficient variability or raster too small.\n")
    results <- rbind(results, data.frame(
      Window = paste0(w_size, "×", w_size),
      Moran_I = NA,
      P_value = NA
    ))
    next
  }
  
  # Compute Moran’s I (wrapped in tryCatch to avoid stop on error)
  I_value <- tryCatch({
    autocor(raster_pvalues, method = "moran", global = TRUE, w = w_matrix)
  }, error = function(e) {
    cat("⚠️ Error computing Moran’s I for window", w_size, ":", e$message, "\n")
    NA
  })
  
  # Approximate p-value if result is numeric
  p_value <- if (is.numeric(I_value)) 2 * (1 - pnorm(abs(I_value))) else NA
  
  cat("Moran’s I:", round(I_value, 4), " | Approx. p-value:", round(p_value, 4), "\n")
  
  results <- rbind(results, data.frame(
    Window = paste0(w_size, "×", w_size),
    Moran_I = round(I_value, 4),
    P_value = round(p_value, 4)
  ))
}


# --- Summary ------------------------------------------------------------------
cat("\n====================================\n")
cat(" Summary of Global Moran’s I (terra, normalized kernels) \n")
cat("====================================\n")
print(results)
cat("====================================\n")

# --- Optional: Save results ---------------------------------------------------
output_csv <- "C:/data/moran_results_terra_normalized.csv"
write.csv(results, output_csv, row.names = FALSE)
cat("\n✅ Results saved to:", output_csv, "\n")
