# Script to create figures from a manuscript
# Manuscript title: Breeding changes water use of winter wheat across Europe
# Author: Dominik Behrend, University of Bonn
# Date: 2026-01-12
#
# Description:
# This script reads the and reproduces Figure 5 from the main manuscript and
# Figures S11 from the supplementary materials.
# It does not reproduce the crop model simulations; it only reproduces data processing and plotting.
# Besides the data set it also requires a shape file for Europe, the crop mask and the 
# climate classifications according to Beck et al. 2018.
# R version 4.5.1 was used for the original visualizations.
# Packages used: sf (version 1.0-21), tidyverse (version 2.0.0), terra (version 1.8-70)
#                data.table (version 1.17.8), lubridate (version 1.9.4)
#
# Inputs:
# - simulation results of the modern cultivar tommi:
#   data/winter_wheat_simulation_outputs_europe_tommi_1991-2020_v1.0.csv.gz
# - simulation results of the historic cultivar S. Dickkopf:
#   data/winter_wheat_simulation_outputs_europe_dickkopf_1991-2020_v1.0.csv
# - shape file country level for the world:
#   Made with Natural Earth. Free vector and raster map data @ naturalearthdata.com.
#   data/shapefile2/ne_10m_admin_0_countries.shp
# - winter wheat crop map for the domain, which was generated using the data by
#   Baumert et al. 2025:  https://doi.org/10.5281/zenodo.14409497
#   data/wheat_cropmask_europe_v1.0.tif
# - the köppen geiger climate classes from Beck et al. 2018:
#   https://doi.org/10.1038/sdata.2018.214
#   Dataset used was the classification according to the present climate:
#   Beck_KG_V1_present_0p0083.tif 
# 
# Outputs:
# - original figures were outputted using the RStudio Export function from the plots section in the
#   RStudio Interface

#- Loading necessary packages and data -####
library(sf)
library(tidyverse)
library(terra)
library(data.table)
library(lubridate)

project_tommi <- fread("F:/20250109_CultivarUpscalingDataPublish/winter_wheat_simulation_outputs_europe_tommi_1991-2020_v1.0.csv.gz", sep = ",")
project_dickkopf <- fread("F:/20250109_CultivarUpscalingDataPublish/winter_wheat_simulation_outputs_europe_dickkopf_1991-2020_v1.0.csv.gz", sep = ",")
shpf <- st_read("F:/20250109_CultivarUpscalingDataPublish/data/shapefile2/ne_10m_admin_0_countries.shp")
cropmap <- rast("F:/20250109_CultivarUpscalingDataPublish/data/wheat_cropmask_europe_v1.0.tif")
kg_raster3 <- rast("F:/20250109_CultivarUpscalingDataPublish/data/Beck_KG_V1_present_0p0083.tif")

#- Preparing the data by adding the climate classes ####
project_dickkopf[, year := as.integer(sub(".*\\.(\\d{4})$", "\\1", CURRENT.DATE))]
project_tommi[, year := as.integer(sub(".*\\.(\\d{4})$", "\\1", CURRENT.DATE))]

project_dickkopf[, WUE := (AbGrBM*10)/ActTran]
project_tommi[, WUE := (AbGrBM*10)/ActTran]

violin_df <- project_dickkopf %>% 
  inner_join(project_tommi,
             by = c("RunningFile_ID", "WGS84_lon", "WGS84_lat", "year"),
             suffix = c("_dick", "_tommi"))
names(violin_df)[14:15] <- c("LON", "LAT")

# first defining different climatic zones based on the Köppen-Geiger classification
kg_raster3 <- rast("D:/PhD/Projects/2_Meine/1_EuroCordexSimulations/KG_Classification/Beck_KG_V1_present_0p0083.tif")
kg_raster3 <- crop(kg_raster3, ext(-11.00348, 30.25406, 35.99601, 67.71814))

pts_df <- violin_df[,13:15]
pts <- vect(pts_df, geom=c("LON","LAT"), crs="EPSG:4326")

kg_raster3 <- project(kg_raster3, cropmap, method = "mode")
violin_df$KG_zone <- terra::extract(kg_raster3, pts)[,2]

# checking how many climate zones there are and how relevant they aare
kg_labels <- data.frame(
  KG_zone = c(0,      4,    5,      7,     8,    9,     14,     15,   18,    25,    26,    27),
  KG_name = c("NaN", "BWh", "BWk", "BSk", "Csa", "Csb", "Cfa", "Cfb", "Dsb", "Dfa", "Dfb", "Dfc")
)

#- Plotting climate zones (Figure S11) -####
shpf <- st_read("F:/working_folder/world/world/borders.shp")
borders <- vect(shpf)                     
borders <- crop(borders, ext(kg_raster3))

kg_small <- terra::aggregate(kg_raster3, fact = 3, fun = "modal")

# 4) Raster -> data.frame for ggplot
df <- as.data.frame(kg_small, xy = TRUE, na.rm = FALSE)
names(df)[3] <- "KG_zone"
df <- df %>%
  left_join(kg_labels, by = "KG_zone") %>%
  filter(!is.na(KG_name), KG_name != "NaN")

# 1) Provide color list 
col_map <- tribble(
  ~KG_name, ~label,                                                 ~hex,
  "BWh",    "BWh - Dry, Arid desert, Hot",                          "#FF0000",
  "BWk",    "BWk - Dry, Arid desert, Cold",                         "#FF9695",
  "BSk",    "BSk - Dry, Semi-arid steppe, Cold",                    "#FFDB63",
  "Csa",    "Csa - Temperate, Dry summer, Hot summer",              "#FFFF00",
  "Csb",    "Csb - Temperate, Dry summer, Warm summer",             "#C6C700",
  "Cfa",    "Cfa - Temperate, No dry season, Hot summer",           "#C6FF4E",
  "Cfb",    "Cfb - Temperate, No dry season, Warm summer",          "#66FF33",
  "Dsb",    "Dsb - Continental, Dry summer, Warm summer",           "#C600C7",
  "Dfa",    "Dfa - Continental, No dry season, Hot summer",         "#00FFFF",
  "Dfb",    "Dfb - Continental, No dry season, Warm summer",        "#38C7FF",
  "Dfc",    "Dfc - Continental, No dry season, Cold summer",        "#007E7D")

# clean hex codes if needed
col_map <- col_map %>% mutate(hex = toupper(trimws(gsub(";$", "", hex))))

df <- df %>% mutate(KG_name = as.character(KG_name))
col_map <- col_map %>% mutate(KG_name = as.character(KG_name),
                              hex = toupper(trimws(gsub(";$", "", hex))))

# attach labels by name 
df$label <- col_map$label[ match(df$KG_name, col_map$KG_name) ]

# if any are NA, quickly see which names didn’t match
setdiff(unique(df$KG_name), col_map$KG_name)  # <- should be character(0)

# factor in your desired legend order
df$label <- factor(df$label, levels = col_map$label)
pal <- setNames(col_map$hex, col_map$label)

ggplot() +
  geom_raster(data = df, aes(x = x, y = y, fill = label)) +
  geom_sf(data = shpf, fill = NA, color = "black", linewidth = 0.25) +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  scale_fill_manual(
    values = setNames(col_map$hex, col_map$label),
    name   = "Köppen–Geiger Zone") +
  theme_bw(base_size = 12) +
  theme(
    panel.grid     = element_blank(),
    axis.title     = element_blank(),
    legend.title   = element_text(size = 20, face = "bold"),
    legend.text    = element_text(size = 18),
    legend.key.size= unit(1.2, "cm"),    # <-- size of the coloured boxes
    legend.spacing.y = unit(1.2, "cm")   # <-- vertical spacing between items
  )

zone_counts <- violin_df %>%
  group_by(KG_zone) %>%
  summarize(
    n_unique_locations = n_distinct(RunningFile_ID),
    n_total_years = n(),  # total rows (locations × years)
    .groups = "drop"
  ) %>%
  arrange(desc(n_unique_locations))
zone_counts <- left_join(zone_counts, kg_labels, by = "KG_zone")

ggplot(zone_counts, aes(x = KG_name, y = n_unique_locations)) +
  geom_col(fill = "steelblue") +
  labs(x = "Köppen-Geiger Climate Zone Code", y = "Number of Unique Grid Locations",
       title = "Distribution of Simulated Locations by Climate Zone") +
  theme_minimal(base_size = 14)

# removing points from zones that are not relevant (less than 1000 individual points)
violin_df <- left_join(violin_df, kg_labels, by = "KG_zone")
bad_ids <- c(0, 5, 18, 4)
violin_df <- violin_df[!KG_zone %in% bad_ids]

#- Plotting Violin plots (Figure 5) -####
# changing into an even longer format based on cultivar
violin_tran <- violin_df[, c("LON","LAT", "ActTran_tommi","ActTran_dick", "KG_name")]
violin_tran_plot <- violin_tran %>% 
  pivot_longer(cols=c(ActTran_dick, ActTran_tommi), 
               names_to="cultivar", values_to="transpiration")

# adding significant differences
stat_test_results <- violin_tran %>%
  group_by(KG_name) %>%
  summarize(
    n = n(),
    mean_diff = round(((mean(ActTran_tommi) - mean(ActTran_dick))/ mean(ActTran_dick) * 100), 1),
    mean_tommi = round((mean(ActTran_tommi)),1),
    mean_dick = round((mean(ActTran_dick)),1),
    sd_tommi = round((sd(ActTran_tommi)),1),
    sd_dick = round((sd(ActTran_dick)),1),
    p_value = tryCatch(
      t.test(ActTran_tommi, ActTran_dick, paired = TRUE)$p.value,
      error = function(e) NA_real_
    ),
    .groups = "drop"
  ) %>%
  mutate(sig_label = case_when(
    p_value < 0.001 ~ "***",
    p_value < 0.01  ~ "**",
    p_value < 0.05  ~ "*",
    TRUE            ~ ""
  ))

stat_test_results$mean_diff2 <- paste0(stat_test_results$mean_diff, "%")
# plotting differences
ggplot(violin_tran_plot, aes(x=KG_name, y=transpiration, fill=cultivar)) +
  geom_violin(trim=FALSE, scale="width") +
  geom_boxplot(width=0.15, position=position_dodge(0.9), alpha=0.7) +
  scale_fill_manual(values=c("ActTran_dick"="#e31a1c", "ActTran_tommi"="#1f78b4"),
                    labels=c("Historic Cultivar", "Modern Cultivar")) +
  labs(x="Climate Zone (Köppen-Geiger)",
       y="Transpiration sum (mm season⁻¹)",
       fill="") + 
  geom_text(data = stat_test_results,
            aes(x = factor(KG_name), y = 750, label = sig_label),
            inherit.aes = FALSE,
            size = 8) + 
  geom_text(data = stat_test_results,
            aes(x = factor(KG_name), y = 700, label = mean_diff2),
            inherit.aes = FALSE,
            size = 8) + 
  theme_bw(base_size=14) +
  theme(axis.text.x = element_text(angle=45, hjust=1, size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 22),
        legend.position = c(0.85, 0.90))