Gap analysis workflow for multiple species

# ── Load libraries ─────────────────────────────────────────────────────────────
library(GapAnalysis)
library(dplyr)

library(tidyr)
library(stringr)
library(terra)

library(leaflet)
library(openxlsx)
library(ggplot2)
library(ggtext)
library(cowplot)
library(htmltools)
library(htmlwidgets)
library(kableExtra)

# ── Set outputs folder ────────────────────────────────────────────────────────
# NOTE: This path must match the path set in the knit: field in the YAML header
# All files will be saved here. Update both locations if changing the path, e.g.:
#outputs_folder <- "C:/Users/you/Desktop/GapAnalysis/outputs"
outputs_folder <- "outputs"
if (!dir.exists(outputs_folder)) dir.create(outputs_folder)

# ── Source functions from GapAnalysis repo (alternative method) ───────────────
#files <- list.files("R", pattern = "\\.[rR]$", full.names = TRUE)
#for (f in files) source(f)

# ── Load data from GapAnalysis repo (alternative method) ─────────────────────
#load("data/CucurbitaData.rda")      # Occurrence data
#load("data/CucurbitaRasts.rda")     # Raster list
#load("data/ProtectedAreas.rda")     # Protected areas raster
#load("data/ecoregions.rda")         # Ecoregions data

# ── Load your own data ────────────────────────────────────────────────────────
# NOTE: To use your own occurrence data and SDM rasters, remove the hashes
# below and update the file paths.
#occData <- read.csv("path/to/your_occurrences.csv") # Local occurrence data
#sdms <- list(                                  # Local rasters read in as list
#   rast("path/to/species1_sdm.tif"),
#   rast("path/to/species2_sdm.tif"))
# NOTE: To download protected areas and ecoregions data locally, run getDatasets(). 
# This requires an active internet connection.
#getDatasets()          # Local download of protected areas and ecoregions data

# ── Load example Cucurbita data (GapAnalysis library) ────────────────────────
# NOTE: To use your own data, hash out the two lines below.
occData <- CucurbitaData                  
sdms    <- terra::unwrap(CucurbitaRasts)

# ── Prepare spatial layers ────────────────────────────────────────────────────
ecos     <- terra::vect(ecoregions)
proAreas <- terra::unwrap(ProtectedAreas)

# Standardize CRS across all spatial layers to prevent mismatch warnings
target_crs <- "+proj=longlat +datum=WGS84"
ecos     <- terra::project(ecos, target_crs)
proAreas <- terra::project(proAreas, target_crs)

# ── Prepare species list and SDM rasters ──────────────────────────────────────
taxa    <- unique(occData$species)
sdmList <- lapply(seq_along(taxa), function(i) {
  r <- sdms[[i]]
  if (terra::crs(r, proj = TRUE) != terra::crs(ecos, proj = TRUE)) {
    r <- terra::project(r, terra::crs(ecos))
  }
  r <- terra::subst(r, 0, NA)
  r[!(is.na(r) | r == 1)] <- NA
  r
})
names(sdmList) <- taxa

# Confirm species names
message("\nConfirm species for gap analysis:\n", paste(taxa, collapse = "\n"))

results_list <- list()

for (i in seq_along(taxa)) {
  taxon          <- taxa[i]
  occurrenceData <- occData[occData$species == taxon, ]
  sdm            <- sdms[[i]]

  if (terra::crs(sdm, proj = TRUE) != terra::crs(ecos, proj = TRUE)) {
    sdm <- terra::project(sdm, terra::crs(ecos))
  }
  sdm <- terra::subst(sdm, 0, NA)
  sdm[!(is.na(sdm) | sdm == 1)] <- NA

  # 1. Run SRSex first (on all records)
  srsex <- SRSex(taxon = taxon, occurrenceData = occurrenceData)

  # 2. Check inputs
  occurrences     <- checkOccurrences(csv = occurrenceData, taxon = taxon)
  if (!inherits(sdm, "SpatRaster")) stop("sdm is not a SpatRaster")
  sdm_checked     <- checksdm(sdm)
  proArea_cropped <- terra::crop(proAreas, terra::ext(sdm_checked))
  proArea_checked <- checkProtectedAreas(protectedArea = proArea_cropped,
                                         sdm           = sdm_checked)

  # Skip species if SDM raster has no values
  if (all(is.na(terra::values(sdm_checked)))) {
    message("Skipping ", taxon, ": SDM raster has no values.")
    next
  }

  # 3. Generate G buffers
  gBuffer <- generateGBuffers(
    taxon          = taxon,
    occurrenceData = occurrences$data,
    bufferDistM    = 50000
  )

  # 4. Ex-situ analysis
  grsex <- GRSex(taxon = taxon, sdm = sdm_checked, gBuffer = gBuffer)
  ersex <- ERSex(
    taxon          = taxon,
    sdm            = sdm_checked,
    occurrenceData = occurrences$data,
    gBuffer        = gBuffer,
    ecoregions     = ecos,
    idColumn       = "ECO_NAME"
  )
  fcsex <- FCSex(taxon = taxon, srsex = srsex, grsex = grsex, ersex = ersex)

  # 5. In-situ analysis
  srsin <- SRSin(
    taxon          = taxon,
    sdm            = sdm_checked,
    occurrenceData = occurrences$data,
    protectedAreas = proArea_checked
  )
  grsin <- GRSin(taxon = taxon, sdm = sdm_checked, protectedAreas = proArea_checked)
  ersin <- ERSin(
    taxon          = taxon,
    sdm            = sdm_checked,
    occurrenceData = occurrences$data,
    protectedAreas = proArea_checked,
    ecoregions     = ecos,
    idColumn       = "ECO_NAME"
  )
  fcsin   <- FCSin(taxon = taxon, srsin = srsin, grsin = grsin, ersin = ersin)
  fcsmean <- FCSc_mean(taxon = taxon, fcsin = fcsin, fcsex = fcsex)

  # 6. Store results
  results_list[[taxon]] <- list(
    srsex           = srsex,
    occurrences     = occurrences,
    sdm_checked     = sdm_checked,
    proArea_checked = proArea_checked,
    eco_checked     = ecoregions,
    grsex           = grsex,
    ersex           = ersex,
    fcsex           = fcsex,
    srsin           = srsin,
    grsin           = grsin,
    ersin           = ersin,
    fcsin           = fcsin,
    fcsmean         = fcsmean
  )
}

# ── Ex-situ metrics ───────────────────────────────────────────────────────────
exsitu_table <- do.call(rbind, lapply(names(results_list), function(taxon) {
  r <- results_list[[taxon]]
  data.frame(
    Species          = taxon,
    SRS_exsitu       = r$srsex$`SRS exsitu`,
    GRS_exsitu       = r$fcsex$`GRS exsitu`,
    ERS_exsitu       = r$fcsex$`ERS exsitu`,
    FCS_exsitu       = r$fcsex$`FCS exsitu`,
    # NOTE: "FCS existu score" is a known typo in the GapAnalysis package.
    # Update to "FCS exsitu score" when fixed upstream.
    FCS_exsitu_score = r$fcsex$`FCS existu score`
  )
}))

# ── In-situ metrics ───────────────────────────────────────────────────────────
insitu_table <- do.call(rbind, lapply(names(results_list), function(taxon) {
  r <- results_list[[taxon]]
  data.frame(
    Species          = taxon,
    SRS_insitu       = r$fcsin$`SRS insitu`,
    GRS_insitu       = r$fcsin$`GRS insitu`,
    ERS_insitu       = r$fcsin$`ERS insitu`,
    FCS_insitu       = r$fcsin$`FCS insitu`,
    FCS_insitu_score = r$fcsin$`FCS insitu score`
  )
}))

# ── Combined metrics ──────────────────────────────────────────────────────────
combined_table <- do.call(rbind, lapply(names(results_list), function(taxon) {
  r <- results_list[[taxon]]
  data.frame(
    Species         = taxon,
    FCSc_min        = r$fcsmean$FCSc_min,
    FCSc_max        = r$fcsmean$FCSc_max,
    FCSc_mean       = r$fcsmean$FCSc_mean,
    FCSc_min_class  = r$fcsmean$FCSc_min_class,
    FCSc_max_class  = r$fcsmean$FCSc_max_class,
    FCSc_mean_class = r$fcsmean$FCSc_mean_class
  )
}))

knitr::kable(exsitu_table,   digits = 2, caption = "<span style='text-align:left; display:block;'>Ex-situ gap analysis metrics</span>",  format = "html") %>% kableExtra::kable_styling()

knitr::kable(insitu_table,   digits = 2, caption = "<span style='text-align:left; display:block;'>In-situ gap analysis metrics</span>",   format = "html") %>% kableExtra::kable_styling()

knitr::kable(combined_table, digits = 2, caption = "<span style='text-align:left; display:block;'>Combined gap analysis metrics</span>",  format = "html") %>% kableExtra::kable_styling()