# 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 Figures S2, S7 and S8 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), lme4 (version 1.1-37),
#                viridis (0.6.5)
#
# 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
# - Seasonal (from sowing to harvest) mean weather variables from Era5 reanalysis
#   https://doi.org/10.24381/cds.adbb2d47
#   data/winter_wheat_simulation_inputs_seasonal_weather_1990-2020_v1.0.csv
# - Yearly mean weather variables from Era5 reanalysis
#   https://doi.org/10.24381/cds.adbb2d47
#   data/winter_wheat_simulation_inputs_yearly_weather_1990-2020_v1.0.csv
# - Soil profile mean variables from the SoilGrids2.0 1000 m dataset
#   https://doi.org/10.5194/soil-7-217-2021
#   data/winter_wheat_simulation_inputs_europe_soilprofile_mean_1990-2020_v1.0.csv.gz
# - Average yearly CO2 concentration measured at Mauna  Loa, Hawaii
#   Can be downloaded here: https://doi.org/10.15138/9N0H-ZH07
#   Data was not included here, due to unclear licensing.
# - 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)
library(lme4)
library(viridis)

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

season_weather <- fread("F:/20250109_CultivarUpscalingDataPublish/data/winter_wheat_simulation_inputs_europe_seasonal_weather_1990-2020_v1.0.csv.gz")
yearly_means <- fread("F:/20250109_CultivarUpscalingDataPublish/data/winter_wheat_simulation_inputs_europe_yearly_weather_1990-2020_v1.0.csv.gz")
profile_means <- fread("F:/20250109_CultivarUpscalingDataPublish/data/winter_wheat_simulation_inputs_europe_soil_profile_mean_v1.0.csv.gz")

c02 <- fread("D:/simplace/workspace/simplace_run/simulation/dominik/EU_Cordex_Simulations/data/202504418_C02Content.csv")

#- processing the data -####
names(c02) <- c("year", "C02")
input_df <- yearly_means %>% 
  left_join(c02, by = c("year"))

mean_input <- input_df %>% 
  group_by(year) %>%
  summarize(
    mean_C02 = mean(C02),
    mean_radiation = mean(RDD_total),
    mean_PRECT = mean(PRECT_total),
    mean_Temp = mean(Temp_mean),
    mean_VP = mean(VP_mean)
  )


mean30 <- yearly_means[
  , .(
    Temp_mean30   = mean(Temp_mean,   na.rm=TRUE),
    RDD_total30   = mean(RDD_total,   na.rm=TRUE),
    VP_mean30     = mean(VP_mean,     na.rm=TRUE),
    PRECT_total30 = mean(PRECT_total, na.rm=TRUE)
  ),
  by = .(WGS84_lon, WGS84_lat)
]

mean30$RDD_total30 <- mean30$RDD_total30*1e-9
names(mean30)[1:2] <- c("LON", "LAT")

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

#- plotting figure S7 -####
mean_input$mean_radiation <- as.numeric(mean_input$mean_radiation)*0.000000001
mean_input$mean_C02 <- as.numeric(mean_input$mean_C02)
plot_input <- pivot_longer(mean_input, cols = 2:6, names_to = "variable", values_to = "climate")

plot_input$variable[plot_input$variable == "mean_C02"] <- "CO2 Concentration (ppm)"
plot_input$variable[plot_input$variable == "mean_radiation"] <- "Radiation (GJ m⁻²)"
plot_input$variable[plot_input$variable == "mean_PRECT"] <- "Precipitation (mm m⁻²)"
plot_input$variable[plot_input$variable == "mean_Temp"] <- "Temperature (°C)"
plot_input$variable[plot_input$variable == "mean_VP"] <- "Vapor Pressure (kPa)"

ggplot(plot_input, aes(x = year, y = climate)) +
  # Raw cultivar time‐series
  geom_line(size = 1.8) +
  facet_wrap(~variable, scales = "free_y") +
  labs(x = "Year", y = "Annual Mean Climate Variable") +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(size = 18),
        strip.text = element_text(size = 18, face = "bold"),
        axis.text.y = element_text(hjust=1, size = 18),
        axis.title = element_text(size = 20, face = "bold"), 
        legend.text = element_text(size = 18), 
        legend.title = element_text(size = 20, face = "bold"))

#- plotting figure S2 -####
ggplot(mean30, aes(x = LON, y = LAT, fill = Temp_mean30)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  #facet_wrap(~ variable, ncol = 2, scales = "free") +
  scale_fill_viridis(option = "D", name   = "Mean Temp (°C)\n 1990-2020") +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 14) +
  theme(strip.text = element_text(face = "bold"),
        axis.text  = element_text(size = 8),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))

ggplot(mean30, aes(x = LON, y = LAT, fill = RDD_total30)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  #facet_wrap(~ variable, ncol = 2, scales = "free") +
  scale_fill_viridis(option = "D", name   = "Mean Radiation \n(GJ m⁻² Year⁻¹)\n1990-2020") +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 14) +
  theme(strip.text = element_text(face = "bold"),
        axis.text  = element_text(size = 8),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))

ggplot(mean30, aes(x = LON, y = LAT, fill = PRECT_total30)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  #facet_wrap(~ variable, ncol = 2, scales = "free") +
  scale_fill_viridis(option = "D", name   = "Mean Precipitation \n(mm m⁻² Year⁻¹)\n1990-2020") +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 14) +
  theme(strip.text = element_text(face = "bold"),
        axis.text  = element_text(size = 8),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))


#- Plotting Figure S8 -####
names(yearly_means)[7:8] <- 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 <- yearly_means[,7:8]
pts <- vect(pts_df, geom=c("LON","LAT"), crs="EPSG:4326")

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

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

yearly_means <- yearly_means %>% group_by(KG_name, year) %>%
  summarise(Temperature = mean(Temp_mean),
            Radiation = mean(RDD_total),
            VP = mean(VP_mean),
            Precipitation = mean(PRECT_total))

# plotting mean environmental variables re-scaled
yearly_means$Radiation <- as.numeric(yearly_means$Radiation)*0.000000001
yearly_means$Temperature <- as.numeric(yearly_means$Temperature)
plot_input2 <- pivot_longer(yearly_means, cols = 3:6, names_to = "variable", values_to = "climate")

# precipitation
Csa <- filter(plot_input2, KG_name == "Cfa" & variable == "Precipitation")

# calculate linear model
lm_prec  <- lm(climate  ~ year, data = Csa)
summary(lm_prec)

# create annotation strings
eqn_climate <- with(summary(lm_prec),
                    paste0("y = ", round(coef(lm_prec)[2],2), "x + ", round(coef(lm_prec)[1],1),
                           "\np = ", format.pval(coef(summary(lm_prec))[2,4], digits=2)))

ggplot(Csa, aes(x = year, y = climate)) +
  # Raw cultivar time‐series
  geom_line(size = 2) +
  facet_wrap(KG_name~variable, scales = "free_y", nrow = 1) +
  labs(x = "Year", y = "Mean Annual Precipitation (mm)") +
  geom_smooth(data = Csa, aes(y = climate),
              method  = "lm", se = FALSE, linetype = "dashed", size = 1.2) +
  annotate("text", x = 1991, y = max(Csa$climate) - 8,
           label = eqn_climate, hjust = 0, size = 7) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        strip.text = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 18), 
        legend.title = element_text(size = 20, face = "bold"))

# radiation
Csa <- filter(plot_input2, KG_name == "Cfa" & variable == "Radiation")

# calculate linear model
lm_prec  <- lm(climate  ~ year, data = Csa)
summary(lm_prec)

# create annotation strings
eqn_climate <- with(summary(lm_prec),
                    paste0("y = ", round(coef(lm_prec)[2],2), "x + ", round(coef(lm_prec)[1],1),
                           "\np = ", format.pval(coef(summary(lm_prec))[2,4], digits=2)))

ggplot(Csa, aes(x = year, y = climate)) +
  # Raw cultivar time‐series
  geom_line(size = 2) +
  facet_wrap(KG_name~variable, scales = "free_y", nrow = 1) +
  labs(x = "Year", y = "Mean Annual Radiation Sum (GJ m⁻² Year⁻¹)") +
  geom_smooth(data = Csa, aes(y = climate),
              method  = "lm", se = FALSE, linetype = "dashed", size = 1.2) +
  annotate("text", x = 1991, y = max(Csa$climate)-0.05,
           label = eqn_climate, hjust = 0, size = 7) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        strip.text = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 18), 
        legend.title = element_text(size = 20, face = "bold"))

# temperature
Csa <- filter(plot_input2, KG_name == "Cfa" & variable == "Temperature")

# calculate linear model
lm_prec  <- lm(climate  ~ year, data = Csa)
summary(lm_prec)

# create annotation strings
eqn_climate <- with(summary(lm_prec),
                    paste0("y = ", round(coef(lm_prec)[2],2), "x + ", round(coef(lm_prec)[1],1),
                           "\np = ", format.pval(coef(summary(lm_prec))[2,4], digits=2)))

ggplot(Csa, aes(x = year, y = climate)) +
  # Raw cultivar time‐series
  geom_line(size = 2) +
  facet_wrap(KG_name~variable, scales = "free_y", nrow = 1) +
  labs(x = "Year", y = "Mean Annual Temperature (°C)") +
  geom_smooth(data = Csa, aes(y = climate),
              method  = "lm", se = FALSE, linetype = "dashed", size = 1.2) +
  annotate("text", x = 1991, y = max(Csa$climate)-0.05,
           label = eqn_climate, hjust = 0, size = 7) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        strip.text = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 18), 
        legend.title = element_text(size = 20, face = "bold"))

# VP
Csa <- filter(plot_input2, KG_name == "Cfa" & variable == "VP")

# calculate linear model
lm_prec  <- lm(climate  ~ year, data = Csa)
summary(lm_prec)

# create annotation strings
eqn_climate <- with(summary(lm_prec),
                    paste0("y = ", round(coef(lm_prec)[2],2), "x + ", round(coef(lm_prec)[1],1),
                           "\np = ", format.pval(coef(summary(lm_prec))[2,4], digits=2)))
Csa$variable[Csa$variable == "VP"] <- "Vapor Pressure"

ggplot(Csa, aes(x = year, y = climate)) +
  # Raw cultivar time‐series
  geom_line(size = 2) +
  facet_wrap(KG_name~variable, scales = "free_y", nrow = 1) +
  labs(x = "Year", y = "Mean Vapor Pressure (kPa)") +
  geom_smooth(data = Csa, aes(y = climate),
              method  = "lm", se = FALSE, linetype = "dashed", size = 1.2) +
  annotate("text", x = 1991, y = max(Csa$climate)-0.01,
           label = eqn_climate, hjust = 0, size = 7) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        strip.text = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 18), 
        legend.title = element_text(size = 20, face = "bold"))


# precipitation
Csa <- filter(plot_input2, KG_name == "Dfb" & variable == "Precipitation")

# calculate linear model
lm_prec  <- lm(climate  ~ year, data = Csa)
summary(lm_prec)

# create annotation strings
eqn_climate <- with(summary(lm_prec),
                    paste0("y = ", round(coef(lm_prec)[2],2), "x + ", round(coef(lm_prec)[1],1),
                           "\np = ", format.pval(coef(summary(lm_prec))[2,4], digits=2)))

ggplot(Csa, aes(x = year, y = climate)) +
  # Raw cultivar time‐series
  geom_line(size = 2) +
  facet_wrap(KG_name~variable, scales = "free_y", nrow = 1) +
  labs(x = "Year", y = "Mean Annual Precipitation (mm)") +
  geom_smooth(data = Csa, aes(y = climate),
              method  = "lm", se = FALSE, linetype = "dashed", size = 1.2) +
  annotate("text", x = 1991, y = max(Csa$climate) - 8,
           label = eqn_climate, hjust = 0, size = 7) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        strip.text = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 18), 
        legend.title = element_text(size = 20, face = "bold"))

# radiation
Csa <- filter(plot_input2, KG_name == "Dfb" & variable == "Radiation")

# calculate linear model
lm_prec  <- lm(climate  ~ year, data = Csa)
summary(lm_prec)

# create annotation strings
eqn_climate <- with(summary(lm_prec),
                    paste0("y = ", round(coef(lm_prec)[2],2), "x + ", round(coef(lm_prec)[1],1),
                           "\np = ", format.pval(coef(summary(lm_prec))[2,4], digits=2)))

ggplot(Csa, aes(x = year, y = climate)) +
  # Raw cultivar time‐series
  geom_line(size = 2) +
  facet_wrap(KG_name~variable, scales = "free_y", nrow = 1) +
  labs(x = "Year", y = "Mean Annual Radiation Sum (GJ m⁻² Year⁻¹)") +
  geom_smooth(data = Csa, aes(y = climate),
              method  = "lm", se = FALSE, linetype = "dashed", size = 1.2) +
  annotate("text", x = 1991, y = max(Csa$climate)-0.05,
           label = eqn_climate, hjust = 0, size = 7) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        strip.text = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 18), 
        legend.title = element_text(size = 20, face = "bold"))

# temperature
Csa <- filter(plot_input2, KG_name == "Dfb" & variable == "Temperature")

# calculate linear model
lm_prec  <- lm(climate  ~ year, data = Csa)
summary(lm_prec)

# create annotation strings
eqn_climate <- with(summary(lm_prec),
                    paste0("y = ", round(coef(lm_prec)[2],2), "x + ", round(coef(lm_prec)[1],1),
                           "\np = ", format.pval(coef(summary(lm_prec))[2,4], digits=2)))

ggplot(Csa, aes(x = year, y = climate)) +
  # Raw cultivar time‐series
  geom_line(size = 2) +
  facet_wrap(KG_name~variable, scales = "free_y", nrow = 1) +
  labs(x = "Year", y = "Mean Annual Temperature (°C)") +
  geom_smooth(data = Csa, aes(y = climate),
              method  = "lm", se = FALSE, linetype = "dashed", size = 1.2) +
  annotate("text", x = 1991, y = max(Csa$climate)-0.05,
           label = eqn_climate, hjust = 0, size = 7) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        strip.text = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 18), 
        legend.title = element_text(size = 20, face = "bold"))

# VP
Csa <- filter(plot_input2, KG_name == "Dfb" & variable == "VP")

# calculate linear model
lm_prec  <- lm(climate  ~ year, data = Csa)
summary(lm_prec)

# create annotation strings
eqn_climate <- with(summary(lm_prec),
                    paste0("y = ", round(coef(lm_prec)[2],2), "x + ", round(coef(lm_prec)[1],1),
                           "\np = ", format.pval(coef(summary(lm_prec))[2,4], digits=2)))

ggplot(Csa, aes(x = year, y = climate)) +
  # Raw cultivar time‐series
  geom_line(size = 2) +
  facet_wrap(KG_name~variable, scales = "free_y", nrow = 1) +
  labs(x = "Year", y = "Mean Vapor Pressure (kPa)") +
  geom_smooth(data = Csa, aes(y = climate),
              method  = "lm", se = FALSE, linetype = "dashed", size = 1.2) +
  annotate("text", x = 1991, y = max(Csa$climate)-0.01,
           label = eqn_climate, hjust = 0, size = 7) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(hjust=1, size = 20),
        axis.title = element_text(size = 22, face = "bold"), 
        strip.text = element_text(size = 22, face = "bold"), 
        legend.text = element_text(size = 18), 
        legend.title = element_text(size = 20, face = "bold"))


#- processing the data for figure S1 -####
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 <- 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")
)

# 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]

# 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")

csa <- filter(violin_tran_plot, KG_name == "Csa")
#csa <- filter(csa, transpiration<= 90)
#csa <- filter(csa, cultivar == "ActTran_tommi")

csa_locations <- csa |>
  distinct(LON, LAT) |>
  mutate(location = paste(LON, LAT, sep = "_"))

project_locations_tommi <- project_tommi |>
  mutate(location = paste(WGS84_lon, WGS84_lat, sep = "_"))
project_locations_dickkopf <- project_dickkopf |>
  mutate(location = paste(WGS84_lon, WGS84_lat, sep = "_"))

# Step 2: Filter project_tommi for rows not present in csa
project_filtered_tommi <- project_locations_tommi |>
  filter(location %in% csa_locations$location)
project_filtered_dickkopf <- project_locations_dickkopf |>
  filter(location %in% csa_locations$location)

filtered_dickkopf <- project_filtered_dickkopf[, .(RunningFile_ID, ActTran, PotTran, year, AbGrBM, cultivar = "Historic Cultivar")]
filtered_tommi <- project_filtered_tommi[, .(RunningFile_ID, ActTran, PotTran, year, AbGrBM, cultivar = "Modern Cultivar")]

combined <- rbindlist(list(filtered_dickkopf, filtered_tommi))

# adding weather and soil
combined <- left_join(combined, season_weather)
combined <- left_join(combined, profile_means)

#- Plotting Figure S1 -####
combined[
  , AWC_class := cut(
    nfc_profile,
    breaks = quantile(nfc_profile, probs = c(0, 0.25, 0.75, 1), na.rm = TRUE),
    labels = c("AWC 0.083-0.106", "AWC 0.106-0.116", "AWC 0.116-0.208"),
    include.lowest = TRUE
  )
]

ggplot(combined, aes(x = PRECT_total, y = ActTran)) +
  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 = "Precipitation sum (mm Season⁻¹)", y = "Transpiration sum (mm Season⁻¹)") +
  #annotate("text", x = -8, y = 100, label = eqn, hjust = 0, vjust = 1,
  #         size = 8, color = "black") +
  facet_grid(AWC_class ~ cultivar) + 
  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.70))