# mapping MAXENT outputs

library(sf)
library(terra)
library(ggplot2)


rm(list=ls())


# setwd
mywd <- "C:/Users/fjascensao/OneDrive - Universidade de Lisboa/Aulas/MIBC/Classes/Class_3" # <--- change to your folder! 
setwd(file.path(mywd, "outputs_fullModel"))

dir() # what is inside the folder

maxentResults <- read.csv("maxentResults.csv")
names(maxentResults)

# We are goint to use the "10th percentile training presence"
# The 10th percentile training presence (P10) is a threshold that defines a habitat's suitability 
# by omitting the lowest 10% of training occurrences. 
# This approach is used to create a more conservative model, 
# reducing the impact of potential errors in input data or unusual, 
# low-quality, or transient records, making it ideal for native habitat mapping. 
treshold_P10 <- maxentResults$X10.percentile.training.presence.Cloglog.threshold # 

# current
myresult.current <- rast("Pelo_ibe.asc") # <--- change to your file!
plot(myresult.current)

myresult.current01 <- myresult.current >= treshold_P10
names(myresult.current01) <- "Current conditions"
plot(myresult.current01)


# FUTURE predictions
HadGEM3L_ssp245 <- rast(file.path(mywd, "future","outputs", "HadGEM3-GC31-LL_ssp245_2050","Pelo_ibe.asc"))
HadGEM3L_ssp585 <- rast(file.path(mywd, "future","outputs", "HadGEM3-GC31-LL_ssp585_2050","Pelo_ibe.asc"))

MIROC6_ssp245 <- rast(file.path(mywd, "future","outputs", "MIROC6_ssp245_2050","Pelo_ibe.asc"))
MIROC6_ssp585 <- rast(file.path(mywd, "future","outputs", "MIROC6_ssp585_2050","Pelo_ibe.asc"))

MPI_ssp245 <- rast(file.path(mywd, "future","outputs", "MPI-ESM1-2-HR_ssp245_2050","Pelo_ibe.asc"))
MPI_ssp585 <- rast(file.path(mywd, "future","outputs", "MPI-ESM1-2-HR_ssp585_2050","Pelo_ibe.asc"))


predictions245 <- c(HadGEM3L_ssp245, MIROC6_ssp245, MPI_ssp245)
plot(predictions245)

predictions585 <- c(HadGEM3L_ssp585, MIROC6_ssp585, MPI_ssp585)
plot(predictions585)

# Average models
predictions245_mean <- mean(predictions245)
names(predictions245_mean) <- "SSP: 245"
plot(predictions245_mean)

predictions585_mean <- mean(predictions585)
names(predictions585_mean) <- "SSP: 585"
plot(predictions585_mean)


predictions245_01 <- predictions245_mean >= treshold_P10
predictions585_01 <- predictions585_mean >= treshold_P10

plot(c(myresult.current01, predictions245_01, predictions585_01))




# ==============================================================================
# GAIN / LOSS / STABLE ANALYSIS
# ==============================================================================

# Reload binary layers fresh (your NAs above alter them)
current_bin <- (rast("Pelo_ibe.asc") >= treshold_P10) * 1
future245_bin <- (predictions245_mean >= treshold_P10) * 1
future585_bin <- (predictions585_mean >= treshold_P10) * 1


# Classification:
#   0 = absent in both (no change, unsuitable)
#   1 = lost (present now, absent future)
#   2 = gained (absent now, present future)
#   3 = stable (present in both)

change_245 <- current_bin * 2 + future245_bin
# reclassify: 0->0 (absent both), 1->2 (gained), 2->1 (lost), 3->3 (stable)
rcl <- matrix(c(0, 0,
                1, 2,
                2, 1,
                3, 3), ncol = 2, byrow = TRUE)
change_245 <- classify(change_245, rcl)
names(change_245) <- "SSP 245"

change_585 <- current_bin * 2 + future585_bin
change_585 <- classify(change_585, rcl)
names(change_585) <- "SSP 585"


plot(change_245)
plot(change_585)


writeRaster(change_245, file.path(mywd, "change_245.tiff"))
writeRaster(change_585, file.path(mywd, "change_585.tiff"))



# ==============================================================================
# GGPLOT VERSION — PUBLICATION-QUALITY FIGURE
# ==============================================================================

library(tidyterra)

# Define colour palette and labels
change_colours <- c("0" = "grey90",       # absent both
                    "1" = "#d73027",       # lost
                    "2" = "#1a9850",       # gained
                    "3" = "#4575b4")       # stable

change_labels <- c("0" = "Absent",
                   "1" = "Lost",
                   "2" = "Gained",
                   "3" = "Stable")

# Convert to factor rasters for ggplot
change_245_f <- as.factor(change_245)
change_585_f <- as.factor(change_585)

# Optional: load country borders for context
# study area
Iberia <- st_read("C:/Users/fjascensao/OneDrive - Universidade de Lisboa/Aulas/MIBC/Classes/Class_1_20260302/Administrative/Iberia.shp")
Iberia <- st_transform(Iberia, 4326)


# --- Panel 1: Current distribution ---
(p_current <- ggplot() +
  geom_spatraster(data = as.factor(current_bin)) +
  geom_sf(data = Iberia, fill = NA, colour = "black", linewidth = 0.3) +
  scale_fill_manual(
    values = c("0" = "grey90", "1" = "#4575b4"),
    labels = c("Absent", "Present"),
    na.translate = FALSE,
    name = NULL
  ) +
  labs(title = "Current distribution") +
  theme_minimal(base_size = 11) +
  theme(
    legend.position = "bottom",
    axis.text = element_text(size = 7),
    plot.title = element_text(face = "bold", hjust = 0.5)
  ) +
  coord_sf(expand = FALSE))


# --- Panel 2: SSP 245 change ---
(p_245 <- ggplot() +
  geom_spatraster(data = change_245_f) +
  geom_sf(data = Iberia, fill = NA, colour = "black", linewidth = 0.3) +
  scale_fill_manual(
    values = change_colours,
    labels = change_labels,
    na.translate = FALSE,
    name = NULL
  ) +
  labs(title = "SSP2-4.5 (2041–2060)") +
  theme_minimal(base_size = 11) +
  theme(
    legend.position = "bottom",
    axis.text = element_text(size = 7),
    plot.title = element_text(face = "bold", hjust = 0.5)
  ) +
  coord_sf(expand = FALSE))


# --- Panel 3: SSP 585 change ---
(p_585 <- ggplot() +
  geom_spatraster(data = change_585_f) +
  geom_sf(data = Iberia, fill = NA, colour = "black", linewidth = 0.3) +
  scale_fill_manual(
    values = change_colours,
    labels = change_labels,
    na.translate = FALSE,
    name = NULL
  ) +
  labs(title = "SSP5-8.5 (2041–2060)") +
  theme_minimal(base_size = 11) +
  theme(
    legend.position = "bottom",
    axis.text = element_text(size = 7),
    plot.title = element_text(face = "bold", hjust = 0.5)
  ) +
  coord_sf(expand = FALSE))


# --- Combine panels ---
library(patchwork)

(p_combined <- p_current + p_245 + p_585 +
  plot_layout(ncol = 3, guides = "collect") +
  plot_annotation(
    title = expression(italic("Pelodytes ibericus") ~ "— Projected range change under climate scenarios"),
    subtitle = "Ensemble mean of 3 GCMs | Threshold: 10th percentile training presence",
    theme = theme(
      plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
      plot.subtitle = element_text(size = 10, hjust = 0.5, colour = "grey40"),
      legend.position = "bottom"
    )
  ))


ggsave(file.path(mywd, "range_change_map.png"),
       p_combined, width = 14, height = 6, dpi = 300)

p_combined


# ==============================================================================
# AREA SUMMARY TABLE
# ==============================================================================

# Cell area (approximate or use expanse() for accurate values)
cell_km2 <- prod(res(current_bin)) * (111.32)^2  # rough approx

calc_areas <- function(change_raster, scenario_name) {
  freq_tbl <- freq(change_raster)
  get_n <- function(val) {
    n <- freq_tbl$count[freq_tbl$value == val]
    if (length(n) == 0) 0 else n
  }
  data.frame(
    scenario = scenario_name,
    stable_km2 = round(get_n(3) * cell_km2, 0),
    gained_km2 = round(get_n(2) * cell_km2, 0),
    lost_km2   = round(get_n(1) * cell_km2, 0)
  )
}


area_summary <- rbind(
  calc_areas(change_245, "SSP2-4.5"),
  calc_areas(change_585, "SSP5-8.5")
)

area_summary$net_change_km2 <- area_summary$gained_km2 - area_summary$lost_km2
print(area_summary)


write.csv(area_summary, file.path(mywd, "range_change_summary.csv"), row.names = FALSE)