library(sf)
library(ggplot2)
library(dplyr)
library(ggpubr)
library(tidyr)  # Load tidyr for drop_na()
library(rlang)  # Load rlang for sym

# Set the working directory
setwd("D:/Whiplash_Paper_version_2/Future_Climate_Results_with_all_models/Shape_file")

# Read the shapefile
ssp585 <- st_read("ssp245_with_risks_8.shp")
print(head(ssp585))

# Specify the columns to keep
columns_to_keep <- c('min', 'max', 'RISK_RATNG', 'EAL_RATNG', 'SOVI_RATNG', 'RESL_RATNG', 'DRGT_RISKR', 'RFLD_EALR','DRGT_EALR','RFLD_RISKR')

columns <- c('RISK_RATNG', 'RISK_RATNG', 'EAL_RATNG', 'SOVI_RATNG', 'RESL_RATNG', 'DRGT_RISKR', 'RFLD_EALR','DRGT_EALR','RFLD_RISKR')

# Define the palette for the boxplots
palette <- c(
  "Low" = "blue",
  "High" = "red"
)

# Define the order for the categories
order <- c("Low", "High")

# Define the comparisons for p-value calculation
my_comparisons <- list(
  c("Low", "High")
)

for (col in columns) {
  print(col)
  ssp585_subset <- ssp585 %>%
    select(all_of(c(col, 'max'))) %>%
    filter(!is.na(!!sym(col)) & !!sym(col) != "" & !(!!sym(col) %in% c("none", "No Rating", "No Expected Annual Losses")))
  
  # Check if there is data left after filtering
  if (nrow(ssp585_subset) == 0) {
    print(paste("No data left for column:", col))
    next
  }
  
  # Relabel the categories
  ssp585_subset[[col]] <- recode(ssp585_subset[[col]],
                                 "Very Low" = "Low",
                                 "Relatively Low" = "Low",
                                 "Relatively High" = "High",
                                 "Very High" = "High")
  
  # Filter to only include "Low" and "High"
  ssp585_subset <- ssp585_subset %>%
    filter(!!sym(col) %in% c("Low", "High"))
  
  # Convert the column to factor and keep only levels present in the data
  levels_present <- unique(ssp585_subset[[col]])
  levels_to_keep <- intersect(order, levels_present)
  ssp585_subset[[col]] <- factor(ssp585_subset[[col]], levels = levels_to_keep)
  print(levels(ssp585_subset[[col]]))
  
  # Check if there are any valid levels left for plotting
  if (length(levels_to_keep) < 2) {
    print(paste("Not enough levels for column:", col))
    next
  }
  
  p <- ggplot(ssp585_subset, aes_string(x = col, y = "max")) +
    geom_boxplot(aes_string(color = col), width = 0.2, fill = "white", outlier.shape = NA) +  # Change color of box lines
    #geom_point(position = position_dodge(width = 0.1), size = 1.0, color = "black") +  # Set point color to black
    scale_color_manual(values = palette) +
    stat_compare_means(comparisons = my_comparisons, label = "p.signif", step.increase = 0.1, label.y = c(1.5, 0.2, 0.3, 0.4)) +  # Decrease spacing for p-value labels
    theme_minimal() +
    scale_y_continuous(limits = c(0, 2.5), breaks = seq(0, 2.5, by = 1)) +  # Set y-axis limits and ticks
    labs(title = " ",
         x = " ",
         y = "Surplus") +
    theme(
      legend.position = "none",  # Remove legend
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.background = element_rect(fill = "white", color = NA),  # Set panel background to white
      plot.background = element_rect(fill = "white", color = NA),  # Set plot background to white
      axis.line = element_line(color = "black", linewidth = 0.5),
      panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
      axis.ticks = element_line(color = "black", linewidth = 0.5),
      axis.ticks.length = unit(0.3, "cm"),  # Increase the length of the tick lines
      #axis.text.x = element_blank(),  # Remove x tick labels
      #axis.text.y = element_blank()   # Remove y tick labels
    )
  print(p)
  
  ggsave(paste0("D:/Whiplash_Paper_version_2/Future_Climate_Results_with_all_models/Vul_2/", col, "_surplus_ssp245", ".png"), 
         plot = p + 
           theme_minimal() + 
           theme(plot.margin = margin(0, 0, 0, 0), 
                 axis.title.y = element_text(size = 12, color = "black"),
                 axis.text = element_text(size = 10, color = "black"),
                 axis.line = element_line(color = "black"),
                 axis.ticks = element_line(color = "black"),
                 panel.grid = element_blank(),
                 panel.border = element_blank(),
                 panel.background = element_rect(fill = "white", color = NA),  # Set panel background to white
                 plot.background = element_rect(fill = "white", color = NA),  # Set plot background to white
                 legend.position = "none") + 
           scale_y_continuous(limits = c(-2, 2), breaks = seq(-2, 2, by = 1), expand = c(0, 0)),
         dpi = 400, 
         width = 1.2, 
         height = 1.5, 
         units = "in")
}