---
title: "GDP_plots"
output: html_document
date: "2025-10-27"
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = "C:/Users/ignasicortesar/OneDrive - Delft University of Technology/Paper3_Private_Adaptation_CGE")
```

```{r, message=FALSE}
library(tidyverse)
library(readxl)
library(ggplot2)
library(gamstransfer)

```

## Setup Components

```{r}
# define function to quickly read in specific variables from the gdx

get_gdx_records <- function(gdx_list, var_name) {
  if (!var_name %in% names(gdx_list)) {
    stop(paste("Variable", var_name, "not found in GDX data."))
  }
  gdx_list[[var_name]]$records
}

```


```{r}

# for mapping purposes

library(eurostat)
library(sf)
library(RColorBrewer)
library(giscoR)

# add Croatia

geodata <- get_eurostat_geospatial(output_class = "sf", resolution = "1", nuts_level = 2, year = 2010, cache = TRUE, update_cache = TRUE)
geodata_1 <- get_eurostat_geospatial(output_class = "sf", resolution = "1", nuts_level = 1, year = 2013, cache = TRUE, update_cache = TRUE)
geodata_1 <- subset(geodata_1, id == "HR0")
geodata_1$id[1] <- "HR00"
geodata_1$NUTS_ID[1] <- "HR00"
geodata_1$FID[1] <- "HR00"
geodata_1$geo[1] <- "HR00"
geodata_1$LEVL_CODE[1] <- 2
geodata_final <- rbind(geodata, geodata_1)
geodata <- geodata_final
```


# Read in all data for 3C scenarios

```{r}

# call from gdx files for all relevant scenarios (not included due to size)

# 3C with adaptation (Ad_Growth disabled)
results_3C_Ad1_r1_gdx <- readGDX("model/EU-EMS-FIRMADAPT/Results/SCGE/FIRM_ADAPT_NUTS2/SimulationResults_Dyn_3C_BCR_Priv_1_Ad_1_AdGrowth_0_AdSpike_0_RegDiff_1_Optimal_0SMALL.gdx")

results_3C_Ad1_r0_gdx <- readGDX("model/EU-EMS-FIRMADAPT/Results/SCGE/FIRM_ADAPT_NUTS2/SimulationResults_Dyn_3C_BCR_Priv_1_Ad_1_AdGrowth_0_AdSpike_0_RegDiff_0_Optimal_0SMALL.gdx")

# 3C without adaptation
results_3C_Ad0_r1_gdx <- readGDX("model/EU-EMS-FIRMADAPT/Results/SCGE/FIRM_ADAPT_NUTS2/SimulationResults_Dyn_3C_BCR_Priv_1_Ad_0_AdGrowth_0_AdSpike_0_RegDiff_1_Optimal_0SMALL.gdx")

results_3C_Ad0_r0_gdx <- readGDX("model/EU-EMS-FIRMADAPT/Results/SCGE/FIRM_ADAPT_NUTS2/SimulationResults_Dyn_3C_BCR_Priv_1_Ad_0_AdGrowth_0_AdSpike_0_RegDiff_0_Optimal_0SMALL.gdx")


# 3C with adaptation (Ad_Growth enabled)
results_3C_Ad1_r1_g1_gdx <- readGDX("model/EU-EMS-FIRMADAPT/Results/SCGE/FIRM_ADAPT_NUTS2/SimulationResults_Dyn_3C_BCR_Priv_1_Ad_1_AdGrowth_1_AdSpike_0_RegDiff_1_Optimal_0SMALL.gdx")

results_3C_Ad1_r0_g1_gdx <- readGDX("model/EU-EMS-FIRMADAPT/Results/SCGE/FIRM_ADAPT_NUTS2/SimulationResults_Dyn_3C_BCR_Priv_1_Ad_1_AdGrowth_1_AdSpike_0_RegDiff_0_Optimal_0SMALL.gdx")

# no climate change
results_nodamages <- readGDX("model/EU-EMS-FIRMADAPT/Results/SCGE/FIRM_ADAPT_NUTS2/SimulationResults_Dyn_3C_NODAMAGES.gdx")

```


# No Adaptation and No Damages

```{r}

GDP_3C_Ad0_r0 <- get_gdx_records(results_3C_Ad0_r0_gdx, "GDP_dyn") %>% 
  rename(GDP_Ad0_r0 = "value")
GDP_3C_Ad0_r1 <- get_gdx_records(results_3C_Ad0_r1_gdx, "GDP_dyn") %>% 
  rename(GDP_Ad0_r1 = "value")
GDP_baseline <- get_gdx_records(results_nodamages, "GDP_dyn") %>% 
  rename(GDP_baseline = "value")

```

# Constant Adaptation

```{r}

GDP_3C_Ad1_r1 <- get_gdx_records(results_3C_Ad1_r1_gdx, "GDP_dyn") %>% 
  rename(GDP_Ad1_r1 = "value")
GDP_3C_Ad1_r0 <- get_gdx_records(results_3C_Ad1_r0_gdx, "GDP_dyn") %>% 
  rename(GDP_Ad1_r0 = "value")


```

# Growing Adaptation

```{r}

GDP_3C_Ad1_r1_g1 <- get_gdx_records(results_3C_Ad1_r1_g1_gdx, "GDP_dyn") %>% 
  rename(GDP_Ad1_r1_g1 = "value")
GDP_3C_Ad1_r0_g1 <- get_gdx_records(results_3C_Ad1_r0_g1_gdx, "GDP_dyn") %>% 
  rename(GDP_Ad1_r0_g1 = "value")


```

```{r}
cluster_classification <- readxl::read_excel("output-processor/cluster_definition.xlsx", sheet = "clustering_output", range = "A1:B272")

cluster_classification <- cluster_classification %>% 
  rename(geo = "cnt")
```


#

```{r}
GDP_3C_allscenarios <- merge(GDP_3C_Ad0_r0, GDP_3C_Ad0_r1, by = c("cnt", "t"))
GDP_3C_allscenarios <- merge(GDP_3C_allscenarios, GDP_3C_Ad1_r0, by = c("cnt", "t"))
GDP_3C_allscenarios <- merge(GDP_3C_allscenarios, GDP_3C_Ad1_r1, by = c("cnt", "t"))
GDP_3C_allscenarios <- merge(GDP_3C_allscenarios, GDP_3C_Ad1_r0_g1, by = c("cnt", "t"))
GDP_3C_allscenarios <- merge(GDP_3C_allscenarios, GDP_3C_Ad1_r1_g1, by = c("cnt", "t"))

# isolate data for 2100

GDP_3C_allscenarios_2100 <- GDP_3C_allscenarios %>% 
  filter(t == 2100) %>% 
  transform(GDP_3C_Ad1_r0_rel = 100 * (GDP_Ad1_r0-GDP_Ad0_r0)/GDP_Ad0_r0) %>%
  transform(GDP_3C_Ad1_r1_rel = 100 * (GDP_Ad1_r1-GDP_Ad0_r1)/GDP_Ad0_r1) %>%
  transform(GDP_3C_Ad1_r0_g1_rel = 100 * (GDP_Ad1_r0_g1-GDP_Ad0_r0)/GDP_Ad0_r0) %>%
  transform(GDP_3C_Ad1_r1_g1_rel = 100 * (GDP_Ad1_r1_g1-GDP_Ad0_r1)/GDP_Ad0_r1) %>%
  select(-GDP_Ad0_r0, -GDP_Ad0_r1, -GDP_Ad1_r0, -GDP_Ad1_r1, -GDP_Ad1_r0_g1, -GDP_Ad1_r1_g1, -t) %>%
  rename("A1 + R0 + G0" = "GDP_3C_Ad1_r0_rel") %>% 
  rename("A1 + R1 + G0" = "GDP_3C_Ad1_r1_rel") %>% 
  rename("A1 + R0 + G1" = "GDP_3C_Ad1_r0_g1_rel") %>% 
  rename("A1 + R1 + G1" = "GDP_3C_Ad1_r1_g1_rel") %>% 
  rename(geo = "cnt")

# compute relative changes data for all years

GDP_3C_allscenarios <- GDP_3C_allscenarios %>% 
  transform(GDP_3C_Ad1_r0_rel = 100 * (GDP_Ad1_r0-GDP_Ad0_r0)/GDP_Ad0_r0) %>%
  transform(GDP_3C_Ad1_r1_rel = 100 * (GDP_Ad1_r1-GDP_Ad0_r1)/GDP_Ad0_r1) %>%
  transform(GDP_3C_Ad1_r0_g1_rel = 100 * (GDP_Ad1_r0_g1-GDP_Ad0_r0)/GDP_Ad0_r0) %>%
  transform(GDP_3C_Ad1_r1_g1_rel = 100 * (GDP_Ad1_r1_g1-GDP_Ad0_r1)/GDP_Ad0_r1) %>%
  select(-GDP_Ad0_r0, -GDP_Ad0_r1, -GDP_Ad1_r0, -GDP_Ad1_r1, -GDP_Ad1_r0_g1, -GDP_Ad1_r1_g1) %>%
  rename("A1 + R0 + G0" = "GDP_3C_Ad1_r0_rel") %>% 
  rename("A1 + R1 + G0" = "GDP_3C_Ad1_r1_rel") %>% 
  rename("A1 + R0 + G1" = "GDP_3C_Ad1_r0_g1_rel") %>% 
  rename("A1 + R1 + G1" = "GDP_3C_Ad1_r1_g1_rel") %>% 
  rename(geo = "cnt")

GDP_3C_allscenarios <- GDP_3C_allscenarios %>% 
  merge(cluster_classification, by = c("geo")) %>% 
  pivot_longer("A1 + R0 + G0":"A1 + R1 + G1", names_to = "Scenario", values_to = "GDP_Change")

```


# Map out GDP changes

```{r}

library(ggplot2)
library(plotly)
library(viridis)

# with regionalised change
GDP_3C_r1_rel <- merge(GDP_3C_Ad1_r1, GDP_3C_Ad0_r1, by = c("cnt", "t")) %>% 
  transform(GDP_change = 100 * (GDP_Ad1_r1-GDP_Ad0_r1)/GDP_Ad0_r1) %>% 
  rename(geo = "cnt") %>% 
  rename(Year = "t")

# without regionalised change
GDP_3C_r0_rel <- merge(GDP_3C_Ad1_r0, GDP_3C_Ad0_r0, by = c("cnt", "t")) %>% 
  transform(GDP_change = 100 * (GDP_Ad1_r0-GDP_Ad0_r0)/GDP_Ad0_r0) %>% 
  rename(geo = "cnt") %>% 
  rename(Year = "t")

# create manual ranges to visualise in map
manual_ranges_GDP <- c(-5, -1, -0.5, 0, 0.5, 1, 2, 5)

# create categorical variables related to the ranges
GDP_3C_r1_rel <- merge(GDP_3C_r1_rel, geodata, by = c("geo")) %>% 
  mutate(cat = cut_to_classes(GDP_change, n = 7, manual = TRUE, manual_breaks = manual_ranges_GDP, decimals = 1))

GDP_3C_r1_rel$cat2<-factor(GDP_3C_r1_rel$cat,
                        levels = c("-5 ~< -1","-1 ~< -0.5", "-0.5 ~< 0", "0 ~< 0.5", "0.5 ~< 1", "1 ~< 2", "2 ~< 5"),
                        labels = c("Less than -1%","-1% to -0.5%", "-0.5% to 0%","0% to 0.5%","0.5% to 1%","1% to 2%","2% to 5%"))

# create categorical variables related to the ranges
GDP_3C_r0_rel <- merge(GDP_3C_r0_rel, geodata, by = c("geo")) %>% 
  mutate(cat = cut_to_classes(GDP_change, n = 7, manual = TRUE, manual_breaks = manual_ranges_GDP, decimals = 1))

GDP_3C_r0_rel$cat2<-factor(GDP_3C_r0_rel$cat,
                        levels = c("-5 ~< -1","-1 ~< -0.5", "-0.5 ~< 0", "0 ~< 0.5", "0.5 ~< 1", "1 ~< 2", "2 ~< 5"),
                        labels = c("Less than -1%","-1% to -0.5%", "-0.5% to 0%","0% to 0.5%","0.5% to 1%","1% to 2%","2% to 5%"))


```

# Plot maps for both scenarios (A1 + R0 and A1 + R1)

```{r}

#Map for 2100 (r0)

GDP_3C_r0_rel_2100 <- GDP_3C_r0_rel %>% 
  filter(Year == 2100)

GDP_3C_r0_rel_2100$GDP_change <- as.numeric(as.character(GDP_3C_r0_rel_2100$GDP_change))

## Joining, by = "geo"
GDP_change_Ad1_r0_2100_map = 
  ggplot(data = GDP_3C_r0_rel_2100) +
  geom_sf(aes(geometry = geometry, fill = GDP_change), color = "dim grey", size = 0.1) +
  scale_fill_gradientn(  colours = RColorBrewer::brewer.pal(11, "RdYlGn"),
  breaks = seq(-2.5, 2.5, by = 1),
  limits = c(-2.5, 2.5),
  guide = guide_colorbar(nbin = 100,
  barheight = unit(10, "cm"),   # increase height
  barwidth  = unit(0.8, "cm")),  # optional
  name = "GDP change in 2100 (%)") +
#  guides(fill = guide_legend(reverse = TRUE, title = "GDP change in 2100")) +
  labs(
    title = "GDP percentage (%) change in 2100 from current adaptation (No Regionalisation)",
    caption = ""
  ) +
  theme_light() +
  theme(legend.position = c(.85, .5),
        panel.border = element_blank(),
        panel.grid.major = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank()) +
  coord_sf(xlim = c(-12, 44), ylim = c(35, 70))

#Map for 2100 (r1)

GDP_3C_r1_rel_2100 <- GDP_3C_r1_rel %>% 
  filter(Year == 2100)

GDP_3C_r1_rel_2100$GDP_change <- as.numeric(as.character(GDP_3C_r1_rel_2100$GDP_change))

## Joining, by = "geo"
GDP_change_Ad1_r1_2100_map = 
  ggplot(data = GDP_3C_r1_rel_2100) +
  geom_sf(aes(geometry = geometry, fill = GDP_change), color = "dim grey", size = 0.1) +
  scale_fill_gradientn(  colours = RColorBrewer::brewer.pal(11, "RdYlGn"),
  breaks = seq(-2.5, 2.5, by = 1),
  limits = c(-2.5, 2.5),
  oob = scales::squish,
  guide = guide_colorbar(nbin = 100,
  barheight = unit(10, "cm"),   # increase height
  barwidth  = unit(0.8, "cm")),  # optional
  name = "GDP change in 2100 (%)") +
#  guides(fill = guide_legend(reverse = TRUE, title = "GDP change in 2100")) +
  labs(
    title = "GDP percentage (%) change in 2100 from current adaptation (Regionalisation)",
    caption = ""
  ) +
  theme_light() +
  theme(legend.position = c(.85, .5),
        panel.border = element_blank(),
        panel.grid.major = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank()) +
  coord_sf(xlim = c(-12, 44), ylim = c(35, 70))

# save underlying data in the figures (no shapefile geometry)

writexl::write_xlsx(GDP_3C_r0_rel_2100, path = "output-analysis/figures/Figure2_data.xlsx")
writexl::write_xlsx(GDP_3C_r1_rel_2100, path = "output-analysis/figures/Figure3_data.xlsx")

```

# Cleveland dot plot

```{r}

GDP_3C_allscenarios_2060_graph <- GDP_3C_allscenarios_2060 %>% 
  merge(cluster_classification, by = c("geo")) %>% 
  pivot_longer("A1 + R0 + G0":"A1 + R1 + G1", names_to = "Scenario", values_to = "GDP_Change")

GDP_3C_allscenarios_2100_graph <- GDP_3C_allscenarios_2100 %>% 
  merge(cluster_classification, by = c("geo")) %>% 
  pivot_longer("A1 + R0 + G0":"A1 + R1 + G1", names_to = "Scenario", values_to = "GDP_Change")


scenario_comparison <- ggplot(GDP_3C_allscenarios_2100_graph, aes(x = GDP_Change, y = geo, color = Scenario)) +
  geom_point(size = 1) +
  facet_wrap(~ Cluster, nrow = 2, ncol = 2, scales = "free_y") +
  labs(
    title = "Cleveland Dot Plot by Cluster",
    x = "Value",
    y = "Region"
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(size = 12, face = "bold")
  )

```


```{r}

# Get the unique clusters (sorted for consistency)
clusters <- sort(unique(GDP_3C_allscenarios_2100_graph$Cluster))


# Loop through clusters and save each plot as p_1, p_2, ...
for (i in seq_along(clusters)) {
  
  cl <- clusters[i]
  
  # Subset and force clean factor levels (alphabetical)
  df_sub <- GDP_3C_allscenarios_2100_graph %>%
    filter(Cluster == cl) %>%
    mutate(
      geo = as.character(geo)                # drop any factor metadata
    )
  
  # Create factor with explicit sorted levels and drop unused levels
  df_sub$geo <- factor(df_sub$geo, levels = rev(sort(unique(df_sub$geo))))
  df_sub <- droplevels(df_sub)
  df_sub <- df_sub %>% 
    filter(Scenario == "A1 + R0 + G0" | Scenario == "A1 + R0 + G1")
  
  p <- ggplot(df_sub, aes(x = GDP_Change, y = geo, color = Scenario)) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "grey60") +
    geom_point(size = 2.5) +
    scale_y_discrete(limits = sort(unique(df_sub$geo))) + 
    labs(
#      title = paste("GDP changes from adaptation by scenario in 2100 - Cluster ", cl),
      x = "GDP Change in 2100",
      y = "NUTS 2 Region"
    ) +
    theme_minimal()+
    theme(
         axis.text.y = element_text(size = 7)
    )

  
  # Create object name like "p_1", "p_2", etc.
  obj_name <- paste0("p_2100_", i)
  
  # Save the plot to the environment
  assign(obj_name, p)
}

```

```{r}

# generate figure 5 comparing GDP changes in the average region per cluster across scenarios

# compute average changes in GDP in 2100 across scenarios

allscenarios_summary <- GDP_3C_allscenarios %>%
  filter(t != 2015 & t !=2020) %>% 
  group_by(Cluster, Scenario, t) %>% 
  summarise(Avg_GDP_Change = mean(GDP_Change)) %>% 
  ungroup()


Figure5_scenariocomparison <- ggplot(allscenarios_summary, 
       aes(x = t, 
           y = Avg_GDP_Change, 
           color = Scenario,
           group = Scenario)) +   # ensures lines connect correctly
  
  geom_line(size = 1) +
  
  facet_wrap(~ Cluster, nrow = 2, ncol = 2, scales = "free_y") +
  
  labs(
    #title = "Avg GDP Change Over Time Across Scenarios",
    x = "Time (t)",
    y = "Average GDP Change (in %)",
    color = "Scenario"
  ) +
  
  theme_minimal() +
  theme(
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "bottom"
  )

writexl::write_xlsx(allscenarios_summary, path = "output-analysis/figures/Figure5_data.xlsx")


```


```{r}

# generate figure 4 portraying the relative change in GDP across scenarios for the 30 most affected regions


allscenarios_2100_summary <- GDP_3C_allscenarios_2100_graph %>% 
   group_by(Cluster, Scenario) %>% 
   summarise(Avg_GDP_Change = mean(GDP_Change)) %>% 
   ungroup()

top_30 <- readxl::read_excel("output-analysis/top_losses.xlsx", sheet = "Sheet1", range = "A1:A30")

GDP_3C_allscenarios_2100_graph <- GDP_3C_allscenarios_2100_graph %>% 
  merge(top_30, by = c("geo"))


GDP_3C_allscenarios_2100_graph$geo <- factor(GDP_3C_allscenarios_2100_graph$geo, levels = rev(sort(unique(GDP_3C_allscenarios_2100_graph$geo))))

Figure4_dotplot <- ggplot(GDP_3C_allscenarios_2100_graph, aes(x = GDP_Change, y = geo, color = Scenario)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey60") +
  geom_line(aes(group = geo), color = "grey70") +
  geom_point(size = 4) +
  scale_y_discrete(limits = sort(unique(GDP_3C_allscenarios_2100_graph$geo))) + 
  labs(
#      title = paste("GDP changes from adaptation by scenario in 2100 - Cluster ", cl),
    x = "GDP Change (in %) in 2100",
    y = "NUTS 2 Region"
    ) +
  theme_minimal()+
  theme(
      axis.text.y = element_text(size = 12)
    )

writexl::write_xlsx(GDP_3C_allscenarios_2100_graph, path = "output-analysis/figures/Figure4_data.xlsx")


```


