# 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 8 from the main manuscript and
# Figures S5 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")

#- Processing the Data -####
project_dickkopf[, year := as.integer(sub(".*\\.(\\d{4})$", "\\1", CURRENT.DATE))]
project_tommi[, year := as.integer(sub(".*\\.(\\d{4})$", "\\1", CURRENT.DATE))]

# there are 12 location-year combinations with negative yields/biomass, not sure how this happened, removing them for now
project_dickkopf <- project_dickkopf %>% filter(!(Yield < 0 | AbGrBM < 0))
project_tommi <- project_tommi %>% filter(!(Yield < 0 | AbGrBM < 0)) 

physio_df <- project_dickkopf %>% 
  inner_join(project_tommi,
             by = c("RunningFile_ID", "WGS84_lon", "WGS84_lat", "year"),
             suffix = c("_dick", "_tommi")) %>% 
  mutate(pctdiff_T = (ActTran_tommi - ActTran_dick) / ActTran_dick   * 100,
         pctdiff_LAI_max = (LAIMax_tommi - LAIMax_dick) / LAIMax_dick * 100,
         pctdiff_LAI_anth = (LAIAnt_tommi - LAIAnt_dick) / LAIAnt_dick * 100,
         pctdiff_Yield = (Yield_tommi - Yield_dick) / Yield_dick * 100,
         pctdiff_Outflow = (SubsOutfl_tommi - SubsOutfl_dick) / SubsOutfl_dick * 100,
         pctdiff_PotTran = (PotTran_tommi - PotTran_dick) / PotTran_dick * 100,
         pctdiff_RootHydrCond = (RootHydrCond_tommi - RootHydrCond_dick) / RootHydrCond_dick * 100,
         pctdiff_PotET0 = (PotET0_tommi - PotET0_dick) / PotET0_dick * 100,
         pctdiff_AGBM = (AbGrBM_tommi - AbGrBM_dick) / AbGrBM_dick * 100)

names(physio_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 <- physio_df[,13:15]
pts <- vect(pts_df, geom=c("LON","LAT"), crs="EPSG:4326")

kg_raster3 <- project(kg_raster3, cropmap, method = "mode")
physio_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")
)

zone_counts <- physio_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")


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

# changing dates to idate and then into doy
physio_df[, AntDate_tommi := as.IDate(AntDate_tommi, format = "%d.%m.%Y")]
physio_df[, AntDate_dick := as.IDate(AntDate_dick, format = "%d.%m.%Y")]

physio_df$doy_tommi <- yday(physio_df$AntDate_tommi)
physio_df$doy_dick <- yday(physio_df$AntDate_dick)

physio_df$pctdiff_AnthDoy <- (physio_df$doy_tommi - physio_df$doy_dick) / physio_df$doy_dick * 100

shpf$LAT <- shpf$LABEL_Y 
shpf$LON <- shpf$LABEL_X 

#- Plotting for the whole domain (Figure 8)-#####
# Plotting Phenology ####
physio_phenology <- physio_df[, c("pctdiff_T","pctdiff_AnthDoy", "KG_name")]

# calculating statistics
mod <- lm(pctdiff_T ~ pctdiff_AnthDoy, data = physio_phenology)
s <- summary(mod)

eqn <- paste0("y = ", round(coef(mod)[2], 3), " x + ", round(coef(mod)[1], 3),
              "\np = ", format.pval(coef(s)[2,4], digits=2), 
              "\nR-squared: ", round(s$r.squared, 3))
physio_phenology$physio <- "Phenology"

ggplot(physio_phenology, aes(x = pctdiff_AnthDoy, y = pctdiff_T)) +
  stat_binhex(bins = 80) +
  scale_fill_viridis_c(name = "Count\n(log10)", trans = "log10") +
  geom_smooth(method = "lm", color = "black", se = FALSE, linewidth = 2) +
  labs(x = "Anthesis DOY difference (%)", y = "Transpiration difference (%)") +
  annotate("text", x = -8, y = 100, label = eqn, hjust = 0, vjust = 1,
           size = 8, color = "black") +
  facet_grid(~physio) + 
  theme_bw()+ 
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 22, face = "bold"),
        strip.text = element_text(size = 22, face = "bold"),
        legend.position = c(0.90, 0.80))

# Plotting LAI ####
physio_lai <- physio_df[, c("pctdiff_T","pctdiff_LAI_max", "KG_name")]

# calculating statistics
mod <- lm(pctdiff_T ~ pctdiff_LAI_max, data = physio_lai)
s <- summary(mod)

eqn <- paste0("y = ", round(coef(mod)[2], 3), " x + ", round(coef(mod)[1], 3),
              "\np = ", format.pval(coef(s)[2,4], digits=2), 
              "\nR-squared: ", round(s$r.squared, 3))
physio_lai$physio <- "Maximum LAI"

ggplot(physio_lai, aes(x = pctdiff_LAI_max, y = pctdiff_T)) +
  stat_binhex(bins = 80) +
  scale_fill_viridis_c(name = "Count\n(log10)", trans = "log10") +
  geom_smooth(method = "lm", color = "black", se = FALSE, linewidth = 2) +
  labs(x = "Maximum LAI difference (%)", y = "Transpiration difference (%)") +
  annotate("text", x = -100, y = 100, label = eqn, hjust = 0, vjust = 1,
           size = 8, color = "black") +
  facet_grid(~physio) + 
  theme_bw()+ 
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 22, face = "bold"),
        strip.text = element_text(size = 22, face = "bold"),
        legend.position = c(0.90, 0.80))

# Plotting Root hydraulic conductance ####
physio_krs <- physio_df[, c("pctdiff_T","pctdiff_RootHydrCond", "KG_name")]

# calculating statistics
mod <- lm(pctdiff_T ~ pctdiff_RootHydrCond, data = physio_krs)
s <- summary(mod)

eqn <- paste0("y = ", round(coef(mod)[2], 3), " x + ", round(coef(mod)[1], 3),
              "\np = ", format.pval(coef(s)[2,4], digits=2), 
              "\nR-squared: ", round(s$r.squared, 3))
physio_krs$physio <- "Maximum Krs"

ggplot(physio_krs, aes(x = pctdiff_RootHydrCond, y = pctdiff_T)) +
  stat_binhex(bins = 80) +
  scale_fill_viridis_c(name = "Count\n(log10)", trans = "log10") +
  geom_smooth(method = "lm", color = "black", se = FALSE, linewidth = 2) +
  labs(x = "Maximum Krs difference (%)", y = "Transpiration difference (%)") +
  annotate("text", x = -100, y = 100, label = eqn, hjust = 0, vjust = 1,
           size = 8, color = "black") +
  facet_grid(~physio) + 
  theme_bw()+ 
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 22, face = "bold"),
        strip.text = element_text(size = 22, face = "bold"),
        legend.position = c(0.90, 0.80))

# Plotting AGBM ####
physio_agbm <- physio_df[, c("pctdiff_T","pctdiff_AGBM", "KG_name")]

# calculating statistics
mod <- lm(pctdiff_T ~ pctdiff_AGBM, data = physio_agbm)
s <- summary(mod)

eqn <- paste0("y = ", round(coef(mod)[2], 3), " x + ", round(coef(mod)[1], 3),
              "\np = ", format.pval(coef(s)[2,4], digits=2), 
              "\nR-squared: ", round(s$r.squared, 3))
physio_agbm$physio <- "Above Ground Biomass"

ggplot(physio_agbm, aes(x = pctdiff_AGBM, y = pctdiff_T)) +
  stat_binhex(bins = 80) +
  scale_fill_viridis_c(name = "Count\n(log10)", trans = "log10") +
  geom_smooth(method = "lm", color = "black", se = FALSE, linewidth = 2) +
  labs(x = "Above Ground Biomass difference (%)", y = "Transpiration difference (%)") +
  annotate("text", x = -100, y = 100, label = eqn, hjust = 0, vjust = 1,
           size = 6.5, color = "black") +
  facet_grid(~physio) + 
  theme_bw()+ 
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 22, face = "bold"),
        strip.text = element_text(size = 22, face = "bold"),
        legend.position = c(0.90, 0.80))

#- Plotting accross different climate zones (Figure S5) -####
# Plotting Phenology Csa ####
physio_df2 <- physio_df %>% filter(KG_name == "Csa")
physio_phenology <- physio_df2[, c("pctdiff_T","pctdiff_AnthDoy", "KG_name")]

# calculating statistics
mod <- lm(pctdiff_T ~ pctdiff_AnthDoy, data = physio_phenology)
s <- summary(mod)

eqn <- paste0("y = ", round(coef(mod)[2], 3), " x + ", round(coef(mod)[1], 3),
              "\np = ", format.pval(coef(s)[2,4], digits=2), 
              "\nR-squared: ", round(s$r.squared, 3))
physio_phenology$physio <- "Phenology"

ggplot(physio_phenology, aes(x = pctdiff_AnthDoy, y = pctdiff_T)) +
  stat_binhex(bins = 80) +
  scale_fill_viridis_c(name = "Count\n(log10)", trans = "log10") +
  geom_smooth(method = "lm", color = "black", se = FALSE, linewidth = 2) +
  labs(x = "Anthesis DOY difference (%)", y = "Transpiration difference (%)") +
  annotate("text", x = -8, y = 100, label = eqn, hjust = 0, vjust = 1,
           size = 8, color = "black") +
  facet_grid(KG_name ~ physio) + 
  theme_bw()+ 
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 22, face = "bold"),
        strip.text = element_text(size = 22, face = "bold"),
        legend.position = c(0.90, 0.80))

# Plotting LAI Csa ####
physio_lai <- physio_df2[, c("pctdiff_T","pctdiff_LAI_max", "KG_name")]

# calculating statistics
mod <- lm(pctdiff_T ~ pctdiff_LAI_max, data = physio_lai)
s <- summary(mod)

eqn <- paste0("y = ", round(coef(mod)[2], 3), " x + ", round(coef(mod)[1], 3),
              "\np = ", format.pval(coef(s)[2,4], digits=2), 
              "\nR-squared: ", round(s$r.squared, 3))
physio_lai$physio <- "Maximum LAI"

ggplot(physio_lai, aes(x = pctdiff_LAI_max, y = pctdiff_T)) +
  stat_binhex(bins = 80) +
  scale_fill_viridis_c(name = "Count\n(log10)", trans = "log10") +
  geom_smooth(method = "lm", color = "black", se = FALSE, linewidth = 2) +
  labs(x = "Maximum LAI difference (%)", y = "Transpiration difference (%)") +
  annotate("text", x = -100, y = 100, label = eqn, hjust = 0, vjust = 1,
           size = 8, color = "black") +
  facet_grid(~physio) + 
  theme_bw()+ 
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 22, face = "bold"),
        strip.text = element_text(size = 22, face = "bold"),
        legend.position = c(0.90, 0.80))

# Plotting Root hydraulic conductance Csa ####
physio_krs <- physio_df2[, c("pctdiff_T","pctdiff_RootHydrCond", "KG_name")]

# calculating statistics
mod <- lm(pctdiff_T ~ pctdiff_RootHydrCond, data = physio_krs)
s <- summary(mod)

eqn <- paste0("y = ", round(coef(mod)[2], 3), " x + ", round(coef(mod)[1], 3),
              "\np = ", format.pval(coef(s)[2,4], digits=2), 
              "\nR-squared: ", round(s$r.squared, 3))
physio_krs$physio <- "Maximum Krs"

ggplot(physio_krs, aes(x = pctdiff_RootHydrCond, y = pctdiff_T)) +
  stat_binhex(bins = 80) +
  scale_fill_viridis_c(name = "Count\n(log10)", trans = "log10") +
  geom_smooth(method = "lm", color = "black", se = FALSE, linewidth = 2) +
  labs(x = "Maximum Krs difference (%)", y = "Transpiration difference (%)") +
  annotate("text", x = -100, y = 100, label = eqn, hjust = 0, vjust = 1,
           size = 8, color = "black") +
  facet_grid(~physio) + 
  theme_bw()+ 
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 22, face = "bold"),
        strip.text = element_text(size = 22, face = "bold"),
        legend.position = c(0.90, 0.80))

# Plotting Phenology Dfc ####
physio_df2 <- physio_df %>% filter(KG_name == "Dfc")
physio_phenology <- physio_df2[, c("pctdiff_T","pctdiff_AnthDoy", "KG_name")]

# calculating statistics
mod <- lm(pctdiff_T ~ pctdiff_AnthDoy, data = physio_phenology)
s <- summary(mod)

eqn <- paste0("y = ", round(coef(mod)[2], 3), " x + ", round(coef(mod)[1], 3),
              "\np = ", format.pval(coef(s)[2,4], digits=2), 
              "\nR-squared: ", round(s$r.squared, 3))
physio_phenology$physio <- "Phenology"

ggplot(physio_phenology, aes(x = pctdiff_AnthDoy, y = pctdiff_T)) +
  stat_binhex(bins = 80) +
  scale_fill_viridis_c(name = "Count\n(log10)", trans = "log10") +
  geom_smooth(method = "lm", color = "black", se = FALSE, linewidth = 2) +
  labs(x = "Anthesis DOY difference (%)", y = "Transpiration difference (%)") +
  annotate("text", x = -8, y = 100, label = eqn, hjust = 0, vjust = 1,
           size = 8, color = "black") +
  facet_grid(KG_name ~ physio) + 
  theme_bw()+ 
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 22, face = "bold"),
        strip.text = element_text(size = 22, face = "bold"),
        legend.position = c(0.90, 0.80))

# Plotting LAI Dfc ####
physio_lai <- physio_df2[, c("pctdiff_T","pctdiff_LAI_max", "KG_name")]

# calculating statistics
mod <- lm(pctdiff_T ~ pctdiff_LAI_max, data = physio_lai)
s <- summary(mod)

eqn <- paste0("y = ", round(coef(mod)[2], 3), " x + ", round(coef(mod)[1], 3),
              "\np = ", format.pval(coef(s)[2,4], digits=2), 
              "\nR-squared: ", round(s$r.squared, 3))
physio_lai$physio <- "Maximum LAI"

ggplot(physio_lai, aes(x = pctdiff_LAI_max, y = pctdiff_T)) +
  stat_binhex(bins = 80) +
  scale_fill_viridis_c(name = "Count\n(log10)", trans = "log10") +
  geom_smooth(method = "lm", color = "black", se = FALSE, linewidth = 2) +
  labs(x = "Maximum LAI difference (%)", y = "Transpiration difference (%)") +
  annotate("text", x = -100, y = 100, label = eqn, hjust = 0, vjust = 1,
           size = 8, color = "black") +
  facet_grid(~physio) + 
  theme_bw()+ 
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 22, face = "bold"),
        strip.text = element_text(size = 22, face = "bold"),
        legend.position = c(0.90, 0.80))

# Plotting Root hydraulic conductance Dfc ####
physio_krs <- physio_df2[, c("pctdiff_T","pctdiff_RootHydrCond", "KG_name")]

# calculating statistics
mod <- lm(pctdiff_T ~ pctdiff_RootHydrCond, data = physio_krs)
s <- summary(mod)

eqn <- paste0("y = ", round(coef(mod)[2], 3), " x + ", round(coef(mod)[1], 3),
              "\np = ", format.pval(coef(s)[2,4], digits=2), 
              "\nR-squared: ", round(s$r.squared, 3))
physio_krs$physio <- "Maximum Krs"

ggplot(physio_krs, aes(x = pctdiff_RootHydrCond, y = pctdiff_T)) +
  stat_binhex(bins = 80) +
  scale_fill_viridis_c(name = "Count\n(log10)", trans = "log10") +
  geom_smooth(method = "lm", color = "black", se = FALSE, linewidth = 2) +
  labs(x = "Maximum Krs difference (%)", y = "Transpiration difference (%)") +
  annotate("text", x = -100, y = 100, label = eqn, hjust = 0, vjust = 1,
           size = 8, color = "black") +
  facet_grid(~physio) + 
  theme_bw()+ 
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 22, face = "bold"),
        strip.text = element_text(size = 22, face = "bold"),
        legend.position = c(0.90, 0.80))