analysis.r

################################################################################
## Analysis for Paper "Climate shift", Fillon, Linsenmeier, Wagner #############
### Handling of large datasets are on another file (Columbia research grid) ####
################################################################################
dev.off()
rm(list=ls())

#Step 0: load libraries, define paths
#Step 1: load data compute in pre_analysis and KG areas
#Step 2: prepare some function for data processing
#Step 3: plots

#Load some tools
libraries=c("renv","xtable","stringr","ggplot2","viridis","hrbrthemes","openxlsx","patchwork","ncdf4","ggpubr","lwgeom","gam","raster","caret","sys","sf","terra","stats","exactextractr","spaMM","metR","rgdal","dplyr","tidyr","base","geifer","data.table","hrbrthemes","pROC","DescTools","rgeos","maptools")
lapply(libraries, require, character.only = TRUE)

#Define path
setwd("")

library(extrafont)
font_import()
loadfonts(device = "pdf")

mydirection <- getwd()  

#Define some elements
shape_file="KOPPEN5" #if other Koppen region, all climate outputs from huge data need to be changed
list_SSP=c("ssp126","ssp370","ssp585") #ssp that we wil study
bin= 2 # bin for visualization
sspname="ssp585" # for Figures specialized on one SSP

#if to run subsample data subsample_run =1
#if to create subsample data from large data subsample_creation=1
#illustrative only
subsample_run=0
subsample_creation=0

#if to run with other damage functions (robustness)
#ie damage function in GDP levels rather than GDP growth
robustness=0

#create subsample

#create subsample
if (subsample_creation==1){
file_list <- list.files(file.path(mydirection, "climate_outputs_pretreated/"), full.names = TRUE)
for (file_path in file_list) {
  file_name <- basename(file_path)
  data <- readRDS(file_path) 
  data <- as.data.frame(data)
  data$value=round(data$value,0)
  if ("model" %in% colnames(data)){
  data_aggregate <- data %>%
    group_by(value, region, type, model, SSP) %>%
    summarise(frequency=sum(frequency,na.rm=T))
    }
  else{
  data_aggregate <- data %>%
    group_by(value, region, type, SSP) %>%
    summarise(frequency=sum(frequency,na.rm=T))
    }
  output_path <- file.path(file.path(mydirection, "climate_outputs_pretreated_subsample/"), file_name)
  saveRDS(data_aggregate, file = sub("\\.csv$", ".rds", output_path))}}

################################################################################
## STEP 1 - LOAD GENERAL DATASETS AND SHAPEFILES  ##############################
################################################################################
list_models=c("GFDL","IPSL","MPI","MRI","UK")

###Take data computed in preanalysis
#These are files for historical, synthetic general, synthetic control, KOPPEN level, all climate models
for (model_name in list_models){
  dataset_SSP=data.frame()
  for (sspname in list_SSP){
    dataset= readRDS(file.path(mydirection,paste0("climate_outputs_pretreated/KOPPEN",model_name,"_KOPPEN5_picontrol_synthetic_",sspname,".rds")))
    if (subsample_run==1){
      dataset= readRDS(file.path(mydirection,paste0("climate_outputs_pretreated_subsample/KOPPEN",model_name,"_KOPPEN5_picontrol_synthetic_",sspname,".rds")))
      dataset = dataset[!is.na(dataset$region),]}
    dataset=as.data.frame(dataset)
    dataset$SSP="Control"
    dataset$SSP[dataset$type==paste0("sm_",sspname)]=substr(dataset$type[dataset$type==paste0("sm_",sspname)],4,nchar(dataset$type[dataset$type==paste0("sm_",sspname)]))
    dataset$SSP[dataset$type==paste0("sg_",sspname)]=substr(dataset$type[dataset$type==paste0("sg_",sspname)],4,nchar(dataset$type[dataset$type==paste0("sg_",sspname)]))
    dataset$type[dataset$type==paste0("sg_",sspname)]="sg"
    dataset$type[dataset$type==paste0("sm_",sspname)]="sm"
    dataset=dataset[!is.na(dataset$frequency),]
    dataset_SSP=rbind(dataset_SSP,dataset)}
  dataset_SSP <- dataset_SSP %>%
    group_by(region,type,value,SSP) %>%
    summarise(frequency=sum(frequency,na.rm=T))
  dataset_SSP$frequency[dataset_SSP$type=="picontrol"]=dataset_SSP$frequency[dataset_SSP$type=="picontrol"]/length(list_SSP)
  assign(paste0("hist_",model_name),dataset_SSP)}

#These are files for projections for SSP, KOPPEN level, all climate models
for (model_name in list_models){
  dataset_SSP=data.frame()
  for (sspname in list_SSP){
    dataset<-data.frame()
    for (SSP in list_SSP){
    dataset= readRDS(file.path(mydirection,paste0("/climate_outputs_pretreated/",model_name,"_KOPPEN_KOPPEN5_",sspname,".rds")))
    if (subsample_run==1){dataset= readRDS(file.path(mydirection,paste0("/climate_outputs_pretreated_subsample/",model_name,"_KOPPEN_KOPPEN5_",sspname,".rds")))}
    dataset=as.data.frame(dataset)
    dataset_SSP=rbind(dataset_SSP,dataset)}
assign(model_name,dataset_SSP)}}

#Take Koppen Boundaries 5 KG regions
boundaries <- rast(file.path(mydirection,"/data_additional/koppen_geiger_0p5.tif"))
name_regions_complete=c(1:5)
name_regions=c(1:5)
for (k in c(2,3)){boundaries <- subst(boundaries, k, 1)}
for (k in c(5,6,7)){boundaries <- subst(boundaries, k, 4)}
for (k in c(9,10,11,12,13,14,15,16)){boundaries <- subst(boundaries, k, 8)}
for (k in c(18,19,20,21,22,23,24,25,26,27,28)){boundaries <- subst(boundaries, k, 17)}
boundaries <- subst(boundaries, 30, 29)
boundaries = as.polygons(boundaries)
boundaries=st_as_sf(boundaries)

################################################################################
###################### STEP 2  - PREPARE SOME FUNCTIONS  #######################
################################################################################
#Function1_processing2 is to process data to clean distributions 
function1_processing<-function(dataset, bin, location, name){
  #transform weight to frequency
  dataset <- dataset %>%
    group_by(region,SSP) %>%
    mutate(frequency=365*(frequency)/sum(frequency))

  #location's choice and round for the bins
  dataset=dataset[dataset$region==location,]
  dataset$value=round(dataset$value, bin)

  #Frequency of a daily temp per location*bin 
  dataset <- dataset %>%
    group_by(value,SSP) %>%
    mutate(frequency=sum(frequency)) %>%
    distinct(value, .keep_all=TRUE)
  dataset$model=name
  return(dataset)}

process_dataset <- function(dataset, bin, location, modelname, type_label, type_label2) {
  dataset_filtered <- dataset %>%
    filter(type == type_label) %>%
    function1_processing(bin = bin, location = location, name = modelname) %>%
    mutate(type = type_label2)
  return(dataset_filtered)}

# Main function to handle the synthetic, projection, and picontrol datasets for a given location, model, and SSP
function1 <- function(dataset1, dataset2, location, modelname, bin) {
  # Process each dataset type using the helper function
  dataset1_picontrol <- process_dataset(dataset1, bin, location, modelname, "picontrol","Control")
  dataset1_sg <- process_dataset(dataset1, bin, location, modelname, "sg","Synth. General")
  dataset1_sm <- process_dataset(dataset1, bin, location, modelname, "sm","Synth. Model")
  dataset2_processed <- process_dataset(dataset2, bin, location, modelname, "projections", "Proj.")
  dataset_combined <- bind_rows(dataset1_picontrol, dataset1_sg, dataset1_sm, dataset2_processed)
  return(as.data.frame(dataset_combined))}

#load and bin synthetic, projection, control, for all models all RCP one location one bin
function2 <- function(location, bin){
  plot1=function1(hist_GFDL,GFDL,location,"GFDL",bin) 
  plot2=function1(hist_MPI,MPI,location,"MPI",bin) 
  plot3=function1(hist_IPSL,IPSL,location,"IPSL",bin) 
  plot4=function1(hist_MRI,MRI,location,"MRI",bin) 
  plot5=function1(hist_UK,UK,location,"UK",bin) 
  value=rbind(plot1,plot2,plot3,plot4,plot5)
  return(value)}

#This function plots the shift index of interest for all models for a location
processing_plot2 <- function(dataset,location,specification){
  location_name=name_regions_complete[location]
  plot<-ggplot(plot1, aes(x = value, y = frequency, color=model)) +
    geom_line(size = 1.5) +
    scale_color_viridis(discrete = TRUE) +
    labs(colour="Model")+
    theme_ipsum() +
    theme(axis.title.x = element_blank(), axis.title.y = element_blank(),legend.position = "none")
    return(plot)}

#prepare data for plot climate shift with winsoring
plot2_prepare <- function(dataset, sspname, location, bin, xmin, xmax, ymin, ymax) {
  plot_averagessp <- dataset[dataset$SSP == sspname, ]
  plot_averagehist <- dataset[dataset$SSP == "Control", ]
  dataset <- rbind(plot_averagessp, plot_averagehist)

  dataset <- dataset %>%
    group_by(value, type) %>%
    summarize(frequency = sum(frequency, na.rm = TRUE)) %>%
    group_by(type) %>%
    mutate(total_frequency = sum(frequency, na.rm = TRUE))

  dataset$frequency <- 365 * dataset$frequency / dataset$total_frequency
  dataset <- dataset[dataset$type != "Synth. General", ]
  dataset$type[dataset$type == "Synth. Model"] <- "Synthetic"
  dataset$type[dataset$type == "Proj."] <- "Projections"

  winsorize_weighted <- function(df) {
    df <- df[order(df$value), ]
    df$cumfreq <- cumsum(df$frequency) / sum(df$frequency)
    lower_cut <- min(df$value[df$cumfreq >= 0.01])
    upper_cut <- max(df$value[df$cumfreq <= 0.99])
    df <- df %>%
      filter(value >= lower_cut, value <= upper_cut)
    return(df)}

  dataset <- dataset %>%
    group_by(type) %>%
    group_modify(~winsorize_weighted(.x)) %>%
    ungroup()

  plot1 <- ggplot(dataset, aes(x = value, y = frequency, color = type)) +
    geom_line(size = 1.5) +
    scale_color_viridis(discrete = TRUE) +
    theme_ipsum() +
    theme(
      axis.text.x = element_text(size = 18),
      axis.text.y = element_text(size = 18),
       axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      legend.position = "none")

  if (location == 1) {
    plot1 <- ggplot(dataset, aes(x = value, y = frequency, color = type)) +
      geom_line(size = 1.5) +
      labs(x = "Daily mean temperature (1°C bin)", colour = "Landscape ") +
      scale_color_viridis(discrete = TRUE) +
      theme_ipsum() +
      theme(
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        axis.title.x = element_text(size = 22, face = "bold"),
        axis.title.y = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size = 22),
        legend.title = element_text(size = 22, face = "bold"))}

  if (location == 2) {
    plot1 <- ggplot(dataset, aes(x = value, y = frequency, color = type)) +
      geom_line(size = 1.5) +
      scale_color_viridis(discrete = TRUE) +
      ggtitle("Distribution of daily mean temperatures") +
      theme_ipsum() +
      theme(
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.title = element_text(size = 24, face = "bold", hjust = 1),
        legend.position = "none")}
  return(plot1)}

#Function to prepare graph - climate shift
function3 <- function(location, sspname, bin){
  dataset=function2(location,bin)
  plot_ssp=dataset[dataset$SSP==sspname & dataset$type=="Proj.",]
  plot_synthetic_model=dataset[dataset$SSP==sspname & dataset$type=="Synth. Model",]
  plot_synthetic_general=dataset[dataset$SSP==sspname & dataset$type=="Synth. General",]  
  colnames(plot_ssp)=c("region","type_proj", "value","SSP","frequency_proj","model")
  dataset=merge(plot_ssp, plot_synthetic_model,by=c("model","value","region","SSP"), all=TRUE)
  dataset$frequency_proj[is.na(dataset$frequency_proj)] <- 0
  dataset$frequency[is.na(dataset$frequency)] <- 0
  dataset$frequency_diff=dataset$frequency_proj - dataset$frequency
  dataset$frequency_diff=dataset$frequency_diff*3.65 #nb days
  return(dataset)}

#plot climate shift (Figure 2 middle)
plot2_prepare2 <- function(dataset, sspname, location, bin, xmin, xmax, ymin, ymax) {
  winsorize_weighted <- function(df) {
    df <- df[order(df$value), ]
    weights <- abs(df$frequency_diff)
    weights[is.na(weights)] <- 0
    if (sum(weights) == 0) return(df)  
    df$cumw <- cumsum(weights) / sum(weights)
    lower_cut <- min(df$value[df$cumw >= 0.01])
    upper_cut <- max(df$value[df$cumw <= 0.99])
    df <- df %>% filter(value >= lower_cut, value <= upper_cut)
    return(df)}

  dataset <- dataset %>%
    group_by(model) %>%
    group_modify(~winsorize_weighted(.x)) %>%
    ungroup()

  dataset_binned <- dataset %>%
    mutate(temp_bin = floor(value)) %>%
    group_by(model, temp_bin) %>%
    summarise(freq_diff_mean = mean(frequency_diff, na.rm = TRUE), .groups = "drop")

  plot1 <- ggplot(dataset_binned, aes(x = temp_bin, y = freq_diff_mean, color = model)) +
    geom_line(size = 1.5) +
    scale_color_brewer(palette = "Dark2") +
    theme_ipsum() +
    theme(
      axis.text.x = element_text(size = 18),
      axis.text.y = element_text(size = 18),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      legend.position = "none")

  if (location == 1) {
    plot1 <- ggplot(dataset_binned, aes(x = temp_bin, y = freq_diff_mean, color = model)) +
      geom_line(size = 1.5) +
      labs(x = "Daily mean temperature (1°C bin)", colour = "ESM") +
      scale_color_brewer(palette = "Dark2") +
      theme_ipsum() +
      theme(
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        axis.title.x = element_text(size = 22, face = "bold"),
        axis.title.y = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size = 22),
        legend.title = element_text(size = 22, face = "bold"))}

  if (location == 2) {
    plot1 <- ggplot(dataset_binned, aes(x = temp_bin, y = freq_diff_mean, color = model)) +
    geom_line(size = 1.5) +
    scale_color_brewer(palette = "Dark2") +
    theme_ipsum() +
    ggtitle("Difference in future change in distribution") +
    theme(
      axis.text.x = element_text(size = 18),
      axis.text.y = element_text(size = 18),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      plot.title = element_text(size = 24, face = "bold", hjust = 1),
      legend.position = "none")}
  return(plot1)}

plot2_prepare3 <- function(dataset_global, bin, xmin_plot, xmax_plot, xmin, xmax, ymin = NULL, ymax = NULL) {
  dataset_global <- dataset_global[dataset_global$Temperature >= xmin, ]
  dataset <- dataset_global[dataset_global$Temperature <= xmax, ]
  ymin = min(dataset_global$Min)*(1+0.1)
  ymax = max(dataset_global$Max)*(1+0.1)
  dataset <- dataset %>%
    pivot_longer(cols = c(Central, Min, Max),
                 names_to = "estimate",
                 values_to = "value") %>%
    dplyr::select(location, Temperature, estimate, value)

  plot1 <- ggplot(dataset, aes(x = Temperature, y = value, linetype = estimate)) +
    geom_hline(yintercept = 0, linetype = "solid", color = "gray30", size = 0.5) +  # ligne y=0
    geom_line(aes(colour = estimate), size = 1) +
    xlim(xmin_plot, xmax_plot) +
    ylim(ymin, ymax) +
    scale_color_manual(values = c("Central" = "black", "Min" = "gray50", "Max" = "gray50"),
                       labels = c("Central" = "Mean", "Min" = "Lower CI", "Max" = "Upper CI")) +
    scale_linetype_manual(values = c("Central" = "solid", "Min" = "dashed", "Max" = "dashed"),
                          labels = c("Central" = "Mean", "Min" = "Lower CI", "Max" = "Upper CI")) +
    labs(
      x = "Daily mean temperature (°C)",
      y = "Relative growth rate",
      colour = "Estimate",
      linetype = "Estimate",
      title = "Growth rate relative to baseline [20°C–22°C]"
    ) +
    theme_classic(base_size = 16) +
    theme(
      plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
      axis.title.y = element_text(size = 16, face = "bold"),
      axis.title.x = element_text(size = 16, face = "bold"),
      axis.text = element_text(size = 14),
      legend.position = "bottom",
      axis.line.x = element_blank(),
      legend.title = element_text(size = 14, face = "bold"),
      legend.text = element_text(size = 13),
      axis.line.y = element_line(size = 0.3))
  return(plot1)}

#Global damage function
function_damage <- function(xmin_value, xmax_value,location_name,zone_name){
  if (zone_name=="Global"){dataset=damage[damage$type=="Global",]}

  #smoothing
  lo <- loess(coefficient ~ Temperature, data=dataset, weights=1/SE, span=0.75, degree=2)
  lo_min <- loess(min~Temperature, data=dataset, weights=1/SE, span=0.75, degree=2)
  lo_max <- loess(max~Temperature, data=dataset, weights=1/SE, span=0.75, degree=2)

  #from 2°C bin to appropriate bin
  temp_seq= seq(floor(min(dataset$Temperature)), ceiling(max(dataset$Temperature)),by = 1/10^(bin))
                                        
  #predict
  dataset <- data.frame(
    Temperature = temp_seq,
    Central = predict(lo, newdata = data.frame(Temperature = temp_seq)),
    Min = predict(lo_min, newdata = data.frame(Temperature = temp_seq)),
    Max = predict(lo_max, newdata = data.frame(Temperature = temp_seq))
  )

  #fill with missing values
  #this allows to know how many dots we have to add below
  diff_min=xmin_value-min(dataset$Temperature)
  diff_max=xmax_value-max(dataset$Temperature)

  coeff_lowerbound=do.call("rbind", replicate(abs(diff_min), dataset[dataset$Temperature==min(dataset$Temperature),], simplify = FALSE))
  coeff_lowerbound=coeff_lowerbound[,c("Min","Max","Temperature","Central")]
  coeff_lowerbound$Temperature=seq(from=xmin_value, to =min(dataset$Temperature)-1, by=1)
  coeff_higherbound=do.call("rbind", replicate(abs(diff_max), dataset[dataset$Temperature==max(dataset$Temperature),], simplify = FALSE))
  coeff_higherbound=coeff_higherbound[,c("Min","Max","Temperature","Central")]
  coeff_higherbound$Temperature=seq(from=max(dataset$Temperature)+1,to=xmax_value, by=1)
  dataset=rbind(dataset, coeff_higherbound)
  dataset=rbind(dataset, coeff_lowerbound)
  dataset=dataset[dataset$Temperature>=xmin_value,]
  dataset=dataset[dataset$Temperature<=xmax_value,]
  dataset$zone=zone_name
  if (zone_name=="Regional"){ dataset$location=location_name}
  return(dataset)}

################################################################################
###################### STEP 3 - RUN THE FUNCTION FOR PLOTS #####################
################################################################################
#in /damage_frompython/, a is for tropical, b is for arid, c temperate, d continental, e polar
#in the rest, 1 is tropical, 2 is arid, 3 is temperate, 4 is continental, 5 is polar

if (robustness==0){
global_damage <- read.csv(file.path(mydirection,"/damage_frompython/coefficients_model_global_01.csv"))
colnames(global_damage)=c("bin","coefficient","SE","P")
global_damage=global_damage[,c("bin","coefficient","SE")]}

if (robustness==1){
global_damage <- read.csv(file.path(mydirection,"/damage_frompython/coefficients_model_global_01fd.csv"))
colnames(global_damage)=c("bin","coefficient","SE","P")
global_damage=global_damage[,c("bin","coefficient","SE")]}

global_damage <- global_damage %>%
  filter(bin != "tp_total") %>%
  mutate(type = "Global")

bins <- read.xlsx(file.path(mydirection,"/data_additional/bins.xlsx"), skipEmptyRows = FALSE, colNames = TRUE)
colnames(bins)=c("bin","Lowerbound","Upperbound")

#last and first bin are winsorized
bins$Lowerbound[bins$Lowerbound=="- infinite"]=-41
bins$Upperbound[bins$Upperbound=="+ infinite"]=41

damage <- bind_rows(
  global_damage %>%
    merge(bins, by = "bin") %>%
    rename(Temperature = Lowerbound) %>%
    mutate(min = coefficient - 2* SE, max = coefficient + 2* SE, location=0))
damage$Temperature <- as.numeric(damage$Temperature)

plot1_data=function2(1,bin)
plot2_data=function2(2,bin)
plot3_data=function2(3,bin)
plot4_data=function2(4,bin)
plot5_data=function2(5,bin)

#plot the distribution of changes
# Calculate xmin, xmax, ymin, ymax across all plots
plot_values <- list(plot1_data$value, plot2_data$value, plot3_data$value, plot4_data$value, plot5_data$value)
plot_frequencies <- list(plot1_data$frequency, plot2_data$frequency, plot3_data$frequency, plot4_data$frequency, plot5_data$frequency)
xmin <- min(unlist(plot_values), na.rm = TRUE) - 0.01
xmax <- max(unlist(plot_values), na.rm = TRUE) + 0.01
ymin <- min(unlist(plot_frequencies), na.rm = TRUE) - 0.01
ymax <- max(unlist(plot_frequencies), na.rm = TRUE) + 0.01

# Compute rounded min/max for each dataset
xmin <- min(unlist(plot_values)) - 0.01
xmax <- max(unlist(plot_values)) + 0.01
ymin <- min(unlist(plot_frequencies)) - 0.01
ymax <- max(unlist(plot_frequencies)) + 0.01
xmin_values <- sapply(plot_values, function(x) round(min(x), bin))
xmax_values <- sapply(plot_values, function(x) round(max(x), bin))
xmin1=round(min(plot1_data$value),bin)
xmin2=round(min(plot2_data$value),bin)
xmin3=round(min(plot3_data$value),bin)
xmin4=round(min(plot4_data$value),bin)
xmin5=round(min(plot5_data$value),bin)
xmax1=round(max(plot1_data$value),bin)
xmax2=round(max(plot2_data$value),bin)
xmax3=round(max(plot3_data$value),bin)
xmax4=round(max(plot4_data$value),bin)
xmax5=round(max(plot5_data$value),bin)

#smooth damage functions over right interval
#do it absolute (same global damage functions for all regions)
#smoothing of coefficients with polynomials
#for global : on whole interval possible
#then add coefficents for bins that are in data but not in damage function interval
smooth_list <- lapply(1:5, function(loc) {
  df <- function_damage(xmin, xmax, loc, "Global")
  df$location <- loc  # ajouter explicitement la colonne location
  return(df)
})
smooth_global <- do.call(rbind, smooth_list)

#alternative xmin and xmax for visualization
xmin=-20
xmax=40 
 
plot1=plot2_prepare(plot1_data,sspname,1,bin,xmin,xmax,ymin,ymax)
plot2=plot2_prepare(plot2_data,sspname,2,bin,xmin,xmax,ymin,ymax)
plot3=plot2_prepare(plot3_data,sspname,3,bin,xmin,xmax,ymin,ymax)
plot4=plot2_prepare(plot4_data,sspname,4,bin,xmin,xmax,ymin,ymax)
plot5=plot2_prepare(plot5_data,sspname,5,bin,xmin,xmax,ymin,ymax)

plot6=function3(1,sspname,bin)
plot7=function3(2,sspname,bin)
plot8=function3(3,sspname,bin)
plot9=function3(4,sspname,bin)
plot10=function3(5,sspname,bin)

ymin= min(sign(plot6$frequency_diff) * log10(1 + abs(plot6$frequency_diff)))
ymax= max(sign(plot6$frequency_diff) * log10(1 + abs(plot6$frequency_diff)))

plot6=plot2_prepare2(plot6,sspname,1,bin,xmin,xmax,ymin,ymax)
plot7=plot2_prepare2(plot7,sspname,2,bin,xmin,xmax,ymin,ymax)
plot8=plot2_prepare2(plot8,sspname,3,bin,xmin,xmax,ymin,ymax)
plot9=plot2_prepare2(plot9,sspname,4,bin,xmin,xmax,ymin,ymax)
plot10=plot2_prepare2(plot10,sspname,5,bin,xmin,xmax,ymin,ymax)

#then plot
y_min_damage=min(smooth_global$Min)
y_max_damage=max(smooth_global$Max)

plot11=plot2_prepare3(smooth_global[smooth_global$location==1,],bin,-30,40, -30, 40, y_min_damage,y_max_damage)

################## FIGURE CLIMATE SHIFT AND DAMAGE FUNCTIONS ###################################
plotA=plot2 + plot7 
plotB=plot4 + plot9 
plotC=plot5 + plot10 
plotD=plot3 + plot8 
plotE=plot1 + plot6 
p_all<-cowplot::plot_grid(plotA,plotB,plotC,plotD,plotE,nrow=5,labels=c("ARID","CONTINENTAL","POLAR","TEMPERATE","TROPICAL"),label_size = 24,label_x=0.01,rel_heights=c(0.85,0.85,0.85,0.85,1))
p_all
ggsave(file.path(mydirection, paste0("/Figures/plot2_shift_",robustness,".pdf")), plot = p_all, width = 10, height = 10,  scale = 1.6, family = "Helvetica")

plot11
ggsave(file.path(mydirection, paste0("/Figures/damage_function.pdf")), plot = last_plot(), width = 10, height = 6, family = "Helvetica")


################## FIGURE UNCERTAINTY ACROSS DIMENSIONS  ###################################
calculate_damage_diff <- function(damage_data, control_values, type_column, type_hist, type_synth) {
  # Calculate damage differences and merge control values
  damage_data <- merge(damage_data, control_values, by = c("location", "model"))
  # Calculate damage differences
  damage_data$damage_diff_central_synth <- damage_data$damage_central 
  damage_data$damage_diff_min_synth <- damage_data$damage_min 
  damage_data$damage_diff_max_synth <- damage_data$damage_max 
  names(damage_data)[names(damage_data) == "type"] <- type_column
  names(damage_data)[names(damage_data) == "SSP"] <- "SSP.x"
  return(damage_data)}

# Function to calculate damage for both synth models
process_synth_models <- function(damage, type_value, control) {
  synth_values <- damage[damage$type == type_value, ]
  synth_values <- calculate_damage_diff(synth_values, control, "type_synth", "Synth. Model", "Synth. General")
  synth_values_model <- synth_values[, c("location", "model", "type_synth", "SSP.x", "damage_diff_central_synth", "damage_diff_min_synth", "damage_diff_max_synth")]
  return(synth_values_model)}

damage_final <- rbind(plot1_data, plot2_data, plot3_data, plot4_data, plot5_data)
names(damage_final)[names(damage_final) == "value"] <- "Temperature"
names(damage_final)[names(damage_final) == "region"] <- "location"
damage_smooth <- smooth_global
damage_final <- merge(damage_final, damage_smooth, by = c("Temperature", "location"))

#save damage with details for future work
damage_final_save_decomposition <- damage_final

# Summarize damage ie over whole temperature distribution
damage_final <- damage_final %>% 
  group_by(model, type, SSP, location) %>%
  summarize(damage_central = sum(frequency * Central, na.rm = TRUE),
            damage_min = sum(frequency * Min, na.rm = TRUE),
            damage_max = sum(frequency * Max, na.rm = TRUE))

# Save a copy of damage_final for future work
damage_final_save <- damage_final

#quick checll
# Process control values (damage_hist)
control_values <- damage_final[damage_final$type == "Control", ]
names(control_values)[names(control_values) == "type"] <- "type_hist"
names(control_values)[names(control_values) == "damage_central"] <- "damage_central_hist"
names(control_values)[names(control_values) == "damage_min"] <- "damage_min_hist"
names(control_values)[names(control_values) == "damage_max"] <- "damage_max_hist"
damage_final <- damage_final[damage_final$type != "Control", ]

# Merge projections with control values
damage_final <- merge(damage_final, control_values, by = c("location", "model"))

# Calculate differences (ratio of change between damage in SSP vs hist)
damage_final$damage_diff_central <- damage_final$damage_central 
damage_final$damage_diff_min <- damage_final$damage_min 
damage_final$damage_diff_max <- damage_final$damage_max 

# Keep only relevant columns
damage_final <- damage_final[, c("location", "model", "type", "SSP.x", "damage_diff_central", "damage_diff_min", "damage_diff_max","damage_central_hist","damage_min_hist","damage_max_hist")]

# Process synthesized models (Synth. Model and Synth. General)
synth_values_model <- process_synth_models(damage_final_save, "Synth. Model", control_values)
synth_values_general <- process_synth_models(damage_final_save, "Synth. General", control_values)

# Combine synthesized model and general values
synth_values <- rbind(synth_values_model, synth_values_general)

# Process projection values
damage_final <- damage_final[damage_final$type == "Proj.", ]
damage_final <- merge(damage_final, synth_values, by = c("location", "model", "SSP.x"))

#compute DD
damage_final$damage_DD_central <- - 100 * (damage_final$damage_diff_central - damage_final$damage_diff_central_synth) / abs(damage_final$damage_diff_central_synth - damage_final$damage_central_hist)
damage_final$damage_DD_min <- - 100 * (damage_final$damage_diff_min - damage_final$damage_diff_min_synth) / abs(damage_final$damage_diff_min_synth-damage_final$damage_min_hist)
damage_final$damage_DD_max <- - 100 * (damage_final$damage_diff_max - damage_final$damage_diff_max_synth) / abs(damage_final$damage_diff_max_synth - damage_final$damage_max_hist)

damage_final$location[damage_final$location==1]="TROPICAL"
damage_final$location[damage_final$location==2]="ARID"
damage_final$location[damage_final$location==3]="TEMPERATE"
damage_final$location[damage_final$location==4]="CONTINENTAL"
damage_final$location[damage_final$location==5]="POLAR"

#in addition, we compute for damage_final the 
#damage_final_weighted
damage_final$damage_weight=abs(damage_final$damage_diff_central_synth-damage_final$damage_central_hist)

# Keep relevant columns for the final output
damage_final <- damage_final[, c("location", "model", "type", "type_synth", "SSP.x", "damage_DD_central", "damage_DD_min", "damage_DD_max","damage_weight")]

#prepare dataset to plot average vs model specific omitted damages
damage_final2 <- damage_final %>%
  group_by(location, type, type_synth, SSP.x) %>%
  mutate(damage_DD_central= mean(damage_DD_central,na.rm=T),
  damage_DD_min= mean(damage_DD_min,na.rm=T),
  damage_DD_max= mean(damage_DD_max,na.rm=T))

damage_final$model_type="Model"
damage_final2$model_type="Average"
damage_final2$model=NA
damage_final=rbind(damage_final,damage_final2)

library(scales)
signed_log_trans <- trans_new(
  "signed_log",
  transform = function(x) sign(x) * log1p(abs(x)),
  inverse = function(x) sign(x) * (expm1(abs(x))),
  domain = c(-Inf, Inf))

breaks_custom <- c(-100, -10, 1, 10, 100, 4000)

plot1 =ggplot(damage_final[damage_final$type_synth=="Synth. Model",], aes(x = damage_DD_central, y = SSP.x, colour=model_type)) +
  geom_point(size=4) +
  facet_grid(location ~ ., scales = "free_y", space = "free_y") +
  ggtitle('Between Earth System Models') +
  xlab('% (signed log scale) of damage over | under-estimated') +
  ylab('') +
  xlim(-10,4000)+
  labs(colour="ESM")+
  geom_vline(xintercept=1, linetype="dotted") +
  theme_linedraw() +
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank())+
  scale_color_manual(values=c("darkred", "steelblue"))+
  theme(legend.text= element_text(size=16),axis.title.x= element_text(size=14),axis.text.y= element_text(size=14), axis.text.x= element_text(size=14), strip.text = element_text(size = 14), legend.title=element_text(face='bold', size=16),plot.title=element_text(size= 18, face='bold'))+
  scale_x_continuous(
    trans = signed_log_trans,
    breaks = breaks_custom,
    labels = scales::label_number(accuracy = 1)) 

plot1
ggsave(file.path(mydirection, paste0("/Figures/plot3_shift_",robustness,".pdf")), plot = last_plot(), width = 6, height = 6,scale=1.4)

#################  FIGURE AGGREGATE DAMAGE  ###################################
#load economic data 2015
DOSE=read.csv(file.path(mydirection,"/data_additional/DOSE/DOSE_V2.csv"))
DOSE=DOSE[DOSE$year==2015,]
DOSE_shapefiles <- st_read(file.path(mydirection,"/data_additional/DOSE/geometries_all_simplified.shp"))
DOSE_shapefiles=as.data.frame(as.data.frame(DOSE_shapefiles)[,c("GID_1")])
colnames(DOSE_shapefiles)=c("GID_1")
DOSE_shapefiles$DOSE_region=1:nrow(DOSE_shapefiles)
DOSE=DOSE[,c("GID_1","grp_pc_usd","pop")]
DOSE=left_join(DOSE,DOSE_shapefiles,by="GID_1")

#merge with dose koppen weights
a=read.csv(file.path(mydirection, "/climate_weighting/a.csv"))
a<- a%>%
  group_by(DOSE_region,region)%>%
  summarise(weight_DOSElevel=sum(weight_DOSElevel,na.rm=T))

GDP_impact=left_join(DOSE,a,by=c("DOSE_region"))

#gdp not pc
GDP_impact$GDP=GDP_impact$pop*GDP_impact$grp_pc_usd
GDP_impact=GDP_impact[,c("DOSE_region","region","weight_DOSElevel","GDP")]

#weight_DOSElevel is multiplication share of zone in DOSE by share of DOSE in Koppen
GDP_impact <- GDP_impact %>%
  group_by(DOSE_region)%>%
  mutate(weight_ID=sum(weight_DOSElevel,na.rm=T))
GDP_impact$weight_DOSElevel=GDP_impact$weight_DOSElevel/GDP_impact$weight_ID

#total GDP (all DOSE)
GDP_impact$value_total=sum(GDP_impact$GDP,na.rm=T)

#for each dose*koppen, multiply share of zone in dose * share of DOSE in koppen by GDP in this DOSE region
#gives the share of each 
GDP_impact <- GDP_impact %>%
  group_by(DOSE_region,region,value_total) %>%
  summarise(GDP=sum(weight_DOSElevel*GDP,na.rm=T))

#then, give GDP of each Dose*Koppen region as a total of all GDP
GDP_impact$share=GDP_impact$GDP/GDP_impact$value_total

#then apply the damage associated to each DOSE_KOPPEN
GDP_impact$region[GDP_impact$region==1]="TROPICAL"
GDP_impact$region[GDP_impact$region==2]="ARID"
GDP_impact$region[GDP_impact$region==3]="TEMPERATE"
GDP_impact$region[GDP_impact$region==4]="CONTINENTAL"
GDP_impact$region[GDP_impact$region==5]="POLAR"
names(GDP_impact)[names(GDP_impact) == "region"] <- "location"
GDP_impact=GDP_impact[,c("location","share")]

#total share, express share as a %
GDP_impact$share_total=sum(GDP_impact$share, na.rm=T)
GDP_impact$share=GDP_impact$share/GDP_impact$share_total

#total share of each location (Koppen Geiger)
GDP_impact <- GDP_impact %>%
  group_by(location) %>%
  summarise(share=sum(share,na.rm=T))
GDP_impact=GDP_impact[!is.na(GDP_impact$location),]

#now merge with DD omitted damages (growth effect) to have global estimates
damage_final=damage_final[damage_final$model_type=="Model",]
damage=merge(GDP_impact,damage_final[,c("location","damage_DD_central","damage_weight","model","model_type","type_synth","SSP.x")],by="location")

#for damage uncertainty, take SSP method, global damages, synthetic model
#average over all models
damage3=damage[damage$type_synth=="Synth. Model",]

damage3 <- damage3 %>%
  group_by(SSP.x, model) %>%
  summarise(change_final=sum(share*damage_DD_central,na.rm=T), change_final_weighted=sum(share*damage_DD_central*damage_weight)/sum(share*damage_weight))

damage3 <- damage3 %>%
  group_by(SSP.x) %>%    # Group by the variable SSP.x
  summarise(change_final = mean(change_final, na.rm = TRUE), change_final_weighted=mean(change_final_weighted,na.rm=T)) %>%  # Calculate the mean
  ungroup()  # Ungroup if needed, to remove grouping after summarising

mean_global=mean(damage3$change_final_weighted,na.rm=T)

#for Global
#for ESM uncertainty, take SSP method, synthetic model
damage2=damage[damage$type_synth=="Synth. Model",]
damage2 <- damage2 %>%
  group_by(SSP.x,model) %>%
  summarise(change_final=sum(share*damage_DD_central,na.rm=T), change_final_weighted=sum(share*damage_DD_central*damage_weight)/sum(share*damage_weight))
damage2$model_type="Model"
damage2_new <- damage2 %>%
  group_by(SSP.x) %>%    # Group by the variable SSP.x
  summarise(change_final = mean(change_final, na.rm = TRUE), change_final_weighted=mean(change_final_weighted,na.rm=T)) %>%  # Calculate the mean
  ungroup()  # Ungroup if needed, to remove grouping after summarising
damage2_new$model_type="Average"
damage2 <- bind_rows(damage2[,c("change_final","SSP.x","model_type","change_final_weighted")], damage2_new)

plot1 <- ggplot(damage2, aes(x = change_final_weighted, y = SSP.x, colour=model_type)) +
  geom_point(size=4) +
  ggtitle('Global Dose-Response') +
  xlab('% of damage underestimated') +
  ylab('') +
  xlim(-8,50)+
  geom_vline(xintercept = mean_global, colour = "darkred", linetype = "dotted", size = 1) +  # Example for steelblue
  labs(colour="ESM")+
  geom_vline(xintercept=0, linetype="dotted") +
  theme_linedraw() +
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank())+
  scale_color_manual(values=c("darkred", "steelblue"))+
  theme(legend.text= element_text(size=16),axis.title.x= element_text(size=14),axis.text.y= element_text(size=14), axis.text.x= element_text(size=14), strip.text = element_text(size = 14), legend.title=element_text(face='bold', size=16), plot.title=element_blank())

if (robustness==1){
plot1 <- ggplot(damage2, aes(x = change_final_weighted, y = SSP.x, colour=model_type)) +
  geom_point(size=4) +
  ggtitle('Global Dose-Response') +
  xlab('% of damage underestimated') +
  ylab('') +
  xlim(-8,50)+
  geom_vline(xintercept = mean_global, colour = "darkred", linetype = "dotted", size = 1) +  # Example for steelblue
  labs(colour="ESM")+
  geom_vline(xintercept=0, linetype="dotted") +
  theme_linedraw() +
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank())+
  scale_color_manual(values=c("darkred", "steelblue"))+
  theme(legend.text= element_text(size=16),axis.title.x= element_text(size=14),axis.text.y= element_text(size=14), axis.text.x= element_text(size=14), strip.text = element_text(size = 14), legend.title=element_text(face='bold', size=16), plot.title=element_blank())}

plot1
ggsave(file.path(mydirection, paste0("/Figures/plot4_shift_",robustness,".pdf")), plot = last_plot(), width = 8, height = 6, scale=1)

#now decomposition par temperature bin
################## FIGURE IMPACTS PER BIN ###########################
damage_final_save_decomposition$location[damage_final_save_decomposition$location==1]="TROPICAL"
damage_final_save_decomposition$location[damage_final_save_decomposition$location==2]="ARID"
damage_final_save_decomposition$location[damage_final_save_decomposition$location==3]="TEMPERATE"
damage_final_save_decomposition$location[damage_final_save_decomposition$location==4]="CONTINENTAL"
damage_final_save_decomposition$location[damage_final_save_decomposition$location==5]="POLAR"

#share of each area in GDP 2050
damage_final_save_decomposition=merge(damage_final_save_decomposition, GDP_impact, by="location")
damage_final_save_decomposition$type[damage_final_save_decomposition$type=="Proj."]="Projection"
damage_final_save_decomposition$type[damage_final_save_decomposition$type=="Synth. Model"]="Synthetic"

#just compute a table for illustration in the paper
total_comp <- damage_final_save_decomposition%>%
  group_by(zone, location, type, SSP) %>%
  summarise(total=sum(Central*frequency,na.rm=T))

total_comp <- total_comp %>%  
  filter(type %in% c("Projection", "Synthetic", "Control"))

total_wide <- total_comp %>% 
  pivot_wider(
    names_from = type,
    values_from = total)

total_wide <- total_wide %>%
  group_by(location, SSP) %>%
  mutate(
    Control = ifelse(SSP == "Control", Control, NA)
  ) %>%
  ungroup()  %>% 
  fill(Control, .direction = "down")  # Remplir les valeurs NA de Control avec la valeur la plus proche (si applicable) " plutôt remplir par 0 stp

total_wide=total_wide[total_wide$SSP!="Control",]

total_wide <- total_wide %>%
  mutate(
    absolute = -100*( Synthetic-Control),
    diff = - 100*(Projection-Synthetic)/abs(Synthetic-Control)
  )
total_wide=total_wide[,c("location","SSP","diff","absolute")]
total_wide=merge(total_wide, GDP_impact, by="location")
total_wide=as.data.frame(total_wide)

total_wide_clean <- total_wide %>%
  mutate(
    SSP = str_sub(SSP, -3, -1),  # Prend les 3 derniers caractères
    SSP = str_replace(SSP, "^(.)(.)(.*)$", "\\1-\\2.\\3"),
    location = str_to_title(location),  # Met en majuscule uniquement la première lettre
    share_rounded = 100*round(share, 2),
    Location_share = paste0(location, "\n(", share_rounded, "% of GDP)")
  )

# Étape 2 : Transformer le tableau
total_wide_clean$absolute=round(total_wide_clean$absolute,3)
total_wide_clean$absolute=ifelse(abs(total_wide_clean$absolute)<0.001,"|x|<1e-3",total_wide_clean$absolute)
total_wide_clean$value = paste0(round(total_wide_clean$diff,2), " [", total_wide_clean$absolute,"]")
table_final <- total_wide_clean %>%
  select(SSP, Location_share, value) %>%
  pivot_wider(
    names_from = Location_share,
    values_from = value
  )

# Affichage du tableau
latex_table <- xtable(table_final, 
                      caption="DD (in \\%) for different SSP in each Köppen-Geiger zones (with their share of 2015 GDP) with global dose-response functions",
                      digits = c(0, rep(2, ncol(table_final))))  # Arrondi à 2 décimales

# Afficher le code LaTeX dans la console
print(latex_table,
      file = file.path(mydirection, paste0("/Figures/table_recap",robustness,".tex")),
      include.rownames = FALSE,
      floating = FALSE,
      size = "small",
      add.to.row = list(pos = list(0), command = "\\centering "))

#now, I plot for each SSP/KG the contributors of the final effect
decomposition=damage_final_save_decomposition

decomposition <- decomposition %>%
  filter(type %in% c("Projection", "Synthetic", "Control")) %>%
  select(location, Temperature, SSP, type, Central, frequency, share) %>%
  mutate(total = Central * frequency) %>%
  replace_na(list(total = 0)) %>%
    group_by(location, Temperature, SSP, type, share) %>%
  summarise(total = sum(total), .groups = "drop") %>%
    pivot_wider(
    names_from = type,
    values_from = total, #total was already 
    values_fill = list(total = 0))

control_values <- decomposition %>%
  filter(SSP == "Control") %>%
  select(location, Temperature, Control)

decomposition=decomposition[,c("location","Temperature","SSP","Synthetic","Projection","share")]
decomposition <- decomposition %>%
  filter(SSP != "Control") %>%
  left_join(control_values, by = c("location", "Temperature")) %>%
  mutate(Control = replace_na(Control, 0))

decomposition <- decomposition %>%
  mutate(across(c(Control, Projection, Synthetic), ~replace_na(., 0)))

decomposition$Temperature_binned = floor(decomposition$Temperature)

#first, compute per bin per region
#same as up, per region
decomposition <- decomposition %>%
  group_by(location, SSP) %>%
  mutate(
    damage_r_weight = abs(sum(Synthetic,na.rm=T) - sum(Control,na.rm=T)),
    damage_DD_r = -100 * (sum(Projection,na.rm=T) - sum(Synthetic,na.rm=T)) / damage_r_weight
  )

decomposition_reg <- decomposition %>%
  group_by(location, SSP) %>%
  summarise(
    damage_r_weight = abs(sum(Synthetic - Control,na.rm=T)),
    damage_DD_r = -100 * (sum(Projection - Synthetic,na.rm=T)) / damage_r_weight,
    damage_DD_r_abs = -100 * (sum(abs(Projection - Synthetic),na.rm=T)) / damage_r_weight
  )

#per bin
decomposition <- decomposition %>%
  group_by(location, SSP, Temperature_binned) %>%
  mutate(
    damage_b_weight = abs(sum(Synthetic - Control,na.rm=T)),
    damage_DD_b = ifelse(damage_b_weight!=0, -100 * (sum(Projection - Synthetic,na.rm=T)) / damage_b_weight, 0)
  )

#global estimate 
global_by_ssp <- decomposition %>%
  group_by(SSP) %>%
  summarise(
    denominator_t = sum(share * damage_r_weight, na.rm = TRUE),
    damage_t = sum(share * damage_DD_r * damage_r_weight, na.rm = TRUE) / denominator_t,
    .groups = "drop")

#regional weighted (i.e. the share of each region in total weighted)
decomposition_regional <- decomposition %>%
  group_by(SSP, location) %>%
  summarise(
    contribution_absolute = sum(share * damage_DD_r * damage_r_weight, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  left_join(global_by_ssp, by = "SSP") %>%
  mutate(
    damage_t_r = 100 * contribution_absolute / (damage_t * denominator_t))

damage_by_bin <- decomposition %>%
  group_by(SSP, Temperature_binned, location) %>%
  summarise(
    damage_weighted = sum(damage_DD_b * damage_b_weight, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  group_by(SSP, location) %>%
  mutate(
    total_signed = sum(abs(damage_weighted), na.rm = TRUE),
    contribution_signed_norm = 100 * damage_weighted / total_signed
  ) %>%
  ungroup()

#this does not sum up to 100
#because signed 
damage_by_bin_global <- damage_by_bin %>%
  left_join(decomposition_regional, by = c("location","SSP")) %>%
  mutate(
    contribution_signed_norm = contribution_signed_norm*damage_t_r) %>%
  group_by(Temperature_binned,SSP) %>%
  summarise(
    contribution_percent = sum(contribution_signed_norm,na.rm=T)) %>%
  group_by(SSP) %>%
  mutate(
    contribution_total = sum(contribution_percent,na.rm=T)) %>%
  group_by(SSP, Temperature_binned) %>%
  summarise(
    contribution_percent = 100*contribution_percent/contribution_total)


plot_global <- ggplot(damage_by_bin_global,
       aes(x = Temperature_binned, y = contribution_percent,
           color = contribution_percent > 0)) +
  geom_segment(aes(xend = Temperature_binned, y = 0, yend = contribution_percent), size = 1) +
  geom_point(size = 2) +
  facet_grid(rows = vars(SSP), scales = "free_y") +
  scale_color_manual(values = c("TRUE" = "#D73027", "FALSE" = "#4575B4")) +
  labs(
    x = "Daily mean temperature (°C)",
    y = "Contribution of temperature bin (in % of omitted damages)") +
  coord_cartesian(xlim = c(-20, 40), ylim = c(min(damage_by_bin_global$contribution_percent), max(damage_by_bin_global$contribution_percent)))+
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "none",
    axis.title.x = element_text(size=16, face="bold"),
    axis.title.y = element_text(size=16, face="bold"),
    strip.text = element_text(face = "bold", size = 14),
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )
plot_global
ggsave(file.path(mydirection, paste0("/Figures/plot_decomposition_temperature",robustness,"_global.pdf")), plot = last_plot(), width = 12, height = 8, scale=1)

#now regional #this does not sum up to 100
#because signed 
damage_by_bin_regional <- damage_by_bin %>%
  left_join(decomposition_regional, by = c("location","SSP")) %>%
  mutate(
    contribution_signed_norm = contribution_signed_norm*damage_t_r) %>%
  group_by(Temperature_binned,SSP, location) %>%
  summarise(
    contribution_percent = sum(contribution_signed_norm,na.rm=T)) %>%
  group_by(SSP, location) %>%
  mutate(
    contribution_total = sum(contribution_percent,na.rm=T)) %>%
  group_by(SSP, Temperature_binned, location) %>%
  summarise(
    contribution_percent = 100*contribution_percent/contribution_total)

plot_regional <- ggplot(damage_by_bin_regional,
       aes(x = Temperature_binned, y = contribution_percent,
           color = contribution_percent > 0)) +
  geom_segment(aes(xend = Temperature_binned, y = 0, yend = contribution_percent), size = 1) +
  geom_point(size = 2) +
  facet_grid(rows = vars(location), cols = vars(SSP), scales = "free_y") +
  scale_color_manual(values = c("TRUE" = "#D73027", "FALSE" = "#4575B4")) +
  labs(
    x = "Daily mean temperature (°C)",
    y = "Contribution of temperature bin (in % of omitted damages)") +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "none",
    axis.title.x = element_text(size=16, face="bold"),
    axis.title.y = element_text(size=16, face="bold"),
    strip.text = element_text(face = "bold", size = 14),
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )
plot_regional
ggsave(file.path(mydirection, paste0("/Figures/plot_decomposition_temperature",robustness,"_regional.pdf")), plot = last_plot(), width = 12, height = 8, scale=1)

################## FIGURE DISTRIBUTIONAL IMPACTS ###################################
DOSE=read.csv(file.path(mydirection,"/data_additional/DOSE/DOSE_V2.csv"))
DOSE=DOSE[DOSE$year==2015,]
DOSE_shapefiles <- st_read(file.path(mydirection,"/data_additional/DOSE/geometries_all_simplified.shp"))
DOSE_shapefiles=as.data.frame(as.data.frame(DOSE_shapefiles)[,c("GID_1")])
colnames(DOSE_shapefiles)=c("GID_1")
DOSE_shapefiles$DOSE_region=1:nrow(DOSE_shapefiles)
DOSE=DOSE[,c("GID_1","grp_pc_usd","ag_grp_pc_usd","man_grp_pc_usd","serv_grp_pc_usd","pop")]
DOSE=left_join(DOSE,DOSE_shapefiles,by="GID_1")

#merge with dose*koppen weights
a=read.csv(file.path(mydirection, "/climate_weighting/a.csv"))
a<- a%>%
  group_by(DOSE_region,region)%>%
  summarise(weight_DOSElevel=sum(weight_DOSElevel,na.rm=T))

distribution=left_join(DOSE,a,by=c("DOSE_region"))
distribution=distribution[,c("grp_pc_usd","region","weight_DOSElevel","DOSE_region")]

#check if complete
distribution <- distribution %>%
  group_by(DOSE_region)%>%
  mutate(weight_ID=sum(weight_DOSElevel,na.rm=T))
distribution$weight_DOSElevel=distribution$weight_DOSElevel/distribution$weight_ID

#then apply the damage associated to each DOSE_KOPPEN
distribution$region[distribution$region==1]="TROPICAL"
distribution$region[distribution$region==2]="ARID"
distribution$region[distribution$region==3]="TEMPERATE"
distribution$region[distribution$region==4]="CONTINENTAL"
distribution$region[distribution$region==5]="POLAR"
names(distribution)[names(distribution) == "region"] <- "location"

#merge with losses from the data
distribution=merge(distribution,total_wide,by=c("location"))

distribution <-distribution %>%
  group_by(DOSE_region,SSP, diff, absolute)%>%
  summarise(grp_pc_usd=sum(grp_pc_usd*weight_DOSElevel,na.rm=T),relative=sum(diff*absolute*weight_DOSElevel,na.rm=T), diff=sum(diff*weight_DOSElevel,na.rm=T), absolute=sum(absolute*weight_DOSElevel,na.rm=T),weight_ID=sum(weight_DOSElevel,na.rm=T))

#first, with simple approach
distribution <- distribution %>%
  mutate(
    log_gdp = log(grp_pc_usd),
    signed_log_diff = sign(diff) * log(1 + abs(diff)))

ymin= min(distribution$signed_log_diff, na.rm=T)
ymax= max(distribution$signed_log_diff, na.rm=T)

df_plot <-  distribution %>%
  filter(SSP == "ssp126")
plot1g  <- ggplot(df_plot, aes(x = log_gdp, y = signed_log_diff)) +
  geom_point(color = "#4682B4", alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", color = "#B22222", size = 1.2, se = TRUE) +
  stat_cor(method = "pearson", label.x = -5, label.y = 5,
         size = 5, color = "gray20", label.sep = "\n", r.digits = 3,
           aes(label = ..r.label..)) +
  labs(
    subtitle = "SSP1-2.6",
    x = "Log GDP per Capita",
    y = "Signed Log Omitted Damages (%)"
  ) +
  ylim(ymin, ymax)+
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 14, hjust = 0.5,  face = "bold"),
    axis.title = element_text(size = 13),
        axis.title.x = element_blank(),
    axis.text = element_text(size = 11)
  )

df_plot <-  distribution %>%
  filter(SSP == "ssp370")
plot2g  <- ggplot(df_plot, aes(x = log_gdp, y = signed_log_diff)) +
  geom_point(color = "#4682B4", alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", color = "#B22222", size = 1.2, se = TRUE) +
 stat_cor(method = "pearson", label.x = -5, label.y = 5,
         size = 5, color = "gray20", label.sep = "\n", r.digits = 3,
         aes(label = ..r.label..)) +
    labs(
    subtitle = "SSP3-7.0",
    x = "Log GDP per Capita",
    y = "Signed Log Omitted Damages (%)"
  ) +
    ylim(ymin, ymax)+
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 14, hjust = 0.5,  face = "bold"),
    axis.title = element_text(size = 13),
    axis.title.y = element_blank(),
    axis.title.x = element_blank(),
    axis.text = element_text(size = 11)
  )

df_plot <-  distribution %>%
  filter(SSP == "ssp585")
plot3g <- ggplot(df_plot, aes(x = log_gdp, y = signed_log_diff)) +
  geom_point(color = "#4682B4", alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", color = "#B22222", size = 1.2, se = TRUE) +
  stat_cor(method = "pearson", label.x = -5, label.y = 5,
         size = 5, color = "gray20", label.sep = "\n", r.digits = 3,
         aes(label = ..r.label..)) +
  labs(
    subtitle = "SSP5-8.5",
    x = "Log GDP per Capita",
    y = "Signed Log Omitted Damages (%)"
  ) +
    ylim(ymin, ymax)+
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 14, hjust = 0.5,  face = "bold"),
    axis.title = element_text(size = 13),
      axis.title.y = element_blank(),
              axis.title.x = element_blank(),
    axis.text = element_text(size = 11)
  )

plot = (plot1g + plot2g + plot3g)
plot

ggsave(file.path(mydirection, paste0("/Figures/plot_distributional_diff_",robustness,".pdf")), plot = last_plot(), width = 10, height = 4)

#first, with simple approach
distribution <- distribution %>%
  mutate(
    log_gdp = log(grp_pc_usd),
    signed_log_diff = sign(absolute) * log(1 + abs(100*absolute)))

ymin= min(distribution$signed_log_diff, na.rm=T)
ymax= max(distribution$signed_log_diff, na.rm=T)

df_plot <-  distribution %>%
  filter(SSP == "ssp126")
plot1g  <- ggplot(df_plot, aes(x = log_gdp, y = signed_log_diff)) +
  geom_point(color = "#4682B4", alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", color = "#B22222", size = 1.2, se = TRUE) +
   stat_cor(method = "pearson", label.x = -5, label.y = 5,
         size = 5, color = "gray20", label.sep = "\n", r.digits = 3,
           aes(label = ..r.label..)) +
  labs(
    subtitle = "SSP1-2.6",
    x = "Log GDP per Capita",
    y = "Signed Log Omitted Damages (%)"
  ) +
  ylim(ymin, ymax)+
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 14, hjust = 0.5,  face = "bold"),
    axis.title = element_text(size = 13),
        axis.title.x = element_blank(),
    axis.text = element_text(size = 11)
  )

df_plot <-  distribution %>%
  filter(SSP == "ssp370")
plot2g  <- ggplot(df_plot, aes(x = log_gdp, y = signed_log_diff)) +
  geom_point(color = "#4682B4", alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", color = "#B22222", size = 1.2, se = TRUE) +
   stat_cor(method = "pearson", label.x = -5, label.y = 5,
         size = 5, color = "gray20", label.sep = "\n", r.digits = 3,
         aes(label = ..r.label..)) +
    labs(
    subtitle = "SSP3-7.0",
    x = "Log GDP per Capita",
    y = "Signed Log Omitted Damages (%)"
  ) +
    ylim(ymin, ymax)+
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 14, hjust = 0.5,  face = "bold"),
    axis.title = element_text(size = 13),
    axis.title.y = element_blank(),
    axis.title.x = element_blank(),
    axis.text = element_text(size = 11)
  )

df_plot <-  distribution %>%
  filter(SSP == "ssp585")
plot3g <- ggplot(df_plot, aes(x = log_gdp, y = signed_log_diff)) +
  geom_point(color = "#4682B4", alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", color = "#B22222", size = 1.2, se = TRUE) +
   stat_cor(method = "pearson", label.x = -5, label.y = 5,
         size = 5, color = "gray20", label.sep = "\n", r.digits = 3,
         aes(label = ..r.label..)) +
  labs(
    subtitle = "SSP5-8.5",
    x = "Log GDP per Capita",
    y = "Signed Log Omitted Damages (%)"
  ) +
    ylim(ymin, ymax)+
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 14, hjust = 0.5,  face = "bold"),
    axis.title = element_text(size = 13),
      axis.title.y = element_blank(),
              axis.title.x = element_blank(),
    axis.text = element_text(size = 11)
  )

plot = (plot1g + plot2g + plot3g)
plot

ggsave(file.path(mydirection, paste0("/Figures/plot_distributional_absolute_",robustness,".pdf")), plot = last_plot(), width = 10, height = 4)

#distribution of absolute damage on histogram
distribution_long <- distribution %>%
  pivot_longer(cols = c(absolute), names_to = "type", values_to = "value")

plot <- ggplot(distribution_long, aes(x = value, fill = type)) +
  geom_histogram(position = "identity", alpha = 0.6, bins = 30, color = "white") +
  scale_fill_brewer(palette = "Set2") +
  labs(
    x = "Value (%)",
    y = "Count",
    fill = "Metric",
    title = "Distribution of Absolute Difference in Damages between Synthetic and Control climate for all SSP and DOSE region",
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(size = 12, face = "bold"),
    axis.title = element_text(size = 13),
    axis.text = element_text(size = 11),
    plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
    legend.position = "NA"
  )

plot
ggsave(file.path(mydirection, paste0("/Figures/distribution_absolute_",robustness,".pdf")), plot = last_plot(), width = 12, height = 7)

#distribution of relative diff and absolute damage #global
DOSE_shapefiles <- st_read(file.path(mydirection,"/data_additional/DOSE/geometries_all_simplified.shp"))
DOSE_shapefiles=as.data.frame(as.data.frame(DOSE_shapefiles)[,c("GID_1","geometry")])
DOSE_shapefiles$DOSE_region=1:nrow(DOSE_shapefiles)
DOSE_shapefiles=merge(DOSE_shapefiles,distribution[,c("diff","absolute","DOSE_region")],by="DOSE_region")
DOSE_shapefiles=DOSE_shapefiles[,c("geometry","diff","absolute","DOSE_region")]
DOSE_shapefiles <- st_as_sf(DOSE_shapefiles)
DOSE_shapefiles <- st_make_valid(DOSE_shapefiles)

#remove strange polygon (rm in manuel's work)
DOSE_shapefiles$bbox <- lapply(st_geometry(DOSE_shapefiles), st_bbox)
DOSE_shapefiles <- DOSE_shapefiles %>%
  mutate(
    width = sapply(bbox, function(x) x["xmax"] - x["xmin"]),
    height = sapply(bbox, function(x) x["ymax"] - x["ymin"])) 
suspect_regions <- DOSE_shapefiles %>%
  filter(width > quantile(width, 0.999))
remove_ID=suspect_regions[1,]$DOSE_region
DOSE_shapefiles=DOSE_shapefiles[DOSE_shapefiles$DOSE_region!=remove_ID,]

p3 <- ggplot(data = DOSE_shapefiles) +
  geom_sf(aes(fill = absolute)) +    # 'data' est la colonne utilisée pour le remplissage
  scale_fill_gradient2(
    low = "blue",        # Couleur pour les valeurs négatives
    mid = "white",       # Couleur pour zéro
    high = "red",        # Couleur pour les valeurs positives
    midpoint = 0         # Point de transition entre les valeurs négatives et positives
  ) +
  scale_x_continuous(labels = function(x) ifelse(x == 0, "", x)) +  # Supprimer le label 0°
  labs(fill = "Absolute \n impacts") +    # Modifier le titre de la légende
  theme_minimal() +
theme(
  plot.background = element_rect(fill = "white", color = NA),
  panel.background = element_rect(fill = "white", color = NA),
  legend.title = element_text(size = 14)
)
ggsave(file.path(mydirection, paste0("/Figures/plot6_shift_absolute",robustness,"_global.png")), plot = p3, width = 10, height = 7, dpi = 300)

DOSE_shapefiles$diff_log <- sign(DOSE_shapefiles$diff) * log(1 + abs(DOSE_shapefiles$diff))
lim <- max(abs(DOSE_shapefiles$diff_log), na.rm = TRUE)

p4 <-  ggplot(data = DOSE_shapefiles) +
  geom_sf(aes(fill = diff_log)) +
  scale_fill_gradient2(
    low = "blue",
    mid = "white",
    high = "red",
    midpoint = 0,
    limits = c(-lim, lim),  # Assure une symétrie visuelle
    name = "Omitted impacts \n (signed log scale)",
    guide = guide_colorbar(
      title.position = "top",
      title.hjust = 0.5
    ),
    labels = function(x){x}
  ) +
  scale_x_continuous(labels = function(x) ifelse(x == 0, "", x)) +
  theme_minimal() +
theme(
  plot.background = element_rect(fill = "white", color = NA),
  panel.background = element_rect(fill = "white", color = NA),
  legend.title = element_text(size = 14)
)
ggsave(file.path(mydirection, paste0("/Figures/plot6_shift_diff",robustness,"_global.png")), plot = p4, width = 10, height = 7, dpi = 300)