# 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:
# Figures reproduced include: Figures 4, 7 from the main manuscript and 
# Figures S2 and S4 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), 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
# - 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(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")

#- 1. Calculating 30 year means for each location accross Europe -####
# there are 12 location-year combinations with negative yields/biomass
# we decided to discard them for the analysis
project_dickkopf <- project_dickkopf %>% filter(!(Yield < 0 | AbGrBM < 0))
project_tommi <- project_tommi %>% filter(!(Yield < 0 | AbGrBM < 0))

setnames(project_dickkopf, old = c("WGS84_lon","WGS84_lat"), new = c("LON","LAT"))
setnames(project_tommi, old = c("WGS84_lon","WGS84_lat"), new = c("LON","LAT"))

# adding climate zones, to check differences in the trend between zones
# first: adapting extend
kg_raster3 <- crop(kg_raster3, ext(-11.00348, 30.25406, 35.99601, 67.71814))

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

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

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

kg_raster3 <- project(kg_raster3, cropmap, method = "mode")
project_dickkopf$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)
project_dickkopf <- left_join(project_dickkopf, kg_labels, by = "KG_zone")
project_tommi <- left_join(project_tommi, kg_labels, by = "KG_zone")

# calculating means
tommi_mean <- project_tommi[, .(mean_AGBM = mean(AbGrBM[AbGrBM >= 0], na.rm = TRUE),
                                mean_Yield = mean(Yield[Yield >= 0], na.rm = TRUE),
                                mean_E = mean(ActEvap[ActEvap >= 0], na.rm = TRUE),
                                mean_T = mean(ActTran[ActTran >= 0], na.rm = TRUE),
                                mean_LAI_max = mean(LAIMax[LAIMax >= 0], na.rm = TRUE),
                                mean_Outflow = mean(SubsOutfl[SubsOutfl >= 0], na.rm = TRUE),
                                mean_LAI_anth = mean(LAIAnt[LAIAnt >= 0], na.rm = TRUE),
                                mean_PotTran = mean(PotTran[PotTran >= 0], na.rm = TRUE),
                                mean_RootHydrCond = mean(RootHydrCond[RootHydrCond >= 0], na.rm = TRUE),
                                mean_PotET0 = mean(PotET0[PotET0 >= 0], na.rm = TRUE), 
                                mean_WUE = (mean(AbGrBM[AbGrBM >= 0], na.rm = TRUE)*10)/mean(ActTran[ActTran >= 0], na.rm = TRUE)),
                            by = c("RunningFile_ID", "LON", "LAT")]

dickkopf_mean <- project_dickkopf[, .(mean_AGBM = mean(AbGrBM[AbGrBM >= 0], na.rm = TRUE),
                                      mean_Yield = mean(Yield[Yield >= 0], na.rm = TRUE),
                                      mean_E = mean(ActEvap[ActEvap >= 0], na.rm = TRUE),
                                      mean_T = mean(ActTran[ActTran >= 0], na.rm = TRUE),
                                      mean_LAI_max = mean(LAIMax[LAIMax >= 0], na.rm = TRUE),
                                      mean_Outflow = mean(SubsOutfl[SubsOutfl >= 0], na.rm = TRUE),
                                      mean_LAI_anth = mean(LAIAnt[LAIAnt >= 0], na.rm = TRUE),
                                      mean_PotTran = mean(PotTran[PotTran >= 0], na.rm = TRUE),
                                      mean_RootHydrCond = mean(RootHydrCond[RootHydrCond >= 0], na.rm = TRUE), 
                                      mean_PotET0 = mean(PotET0[PotET0 >= 0], na.rm = TRUE), 
                                      mean_WUE = (mean(AbGrBM[AbGrBM >= 0], na.rm = TRUE)*10)/mean(ActTran[ActTran >= 0], na.rm = TRUE)),
                                  by = c("RunningFile_ID", "LON", "LAT")]

# calculating cultivar differences
diff_df <- dickkopf_mean %>% 
  inner_join(tommi_mean,
             by = c("RunningFile_ID", "LON", "LAT"),
             suffix = c("_dick", "_tommi")) %>%
  mutate(pctdiff_E = (mean_E_tommi - mean_E_dick) / mean_E_dick   * 100,
         pctdiff_T = (mean_T_tommi - mean_T_dick) / mean_T_dick   * 100,
         pctdiff_LAI_max = (mean_LAI_max_tommi - mean_LAI_max_dick) / mean_LAI_max_dick * 100,
         pctdiff_LAI_anth = (mean_LAI_anth_tommi - mean_LAI_anth_dick) / mean_LAI_anth_dick * 100,
         pctdiff_Yield = (mean_Yield_tommi - mean_Yield_dick) / mean_Yield_dick * 100,
         pctdiff_Outflow = (mean_Outflow_tommi - mean_Outflow_dick) / mean_Outflow_dick * 100,
         pctdiff_PotTran = (mean_PotTran_tommi - mean_PotTran_dick) / mean_PotTran_dick * 100,
         pctdiff_RootHydrCond = (mean_RootHydrCond_tommi - mean_RootHydrCond_dick) / mean_RootHydrCond_dick * 100,
         pctdiff_PotET0 = (mean_PotET0_tommi - mean_PotET0_dick) / mean_PotET0_dick * 100,
         pctdiff_AGBM = (mean_AGBM_tommi - mean_AGBM_dick) / mean_AGBM_dick * 100,
         pctdiff_WUE = (mean_WUE_tommi - mean_WUE_dick) / mean_WUE_dick * 100)
names(diff_df)[2:3] <- c("LON", "LAT")

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

diff_df$mean_ET_dick <- diff_df$mean_E_dick + diff_df$mean_T_dick
diff_df$mean_ET_tommi <- diff_df$mean_E_tommi + diff_df$mean_T_tommi
diff_df$pctdiff_ET <- ((diff_df$mean_ET_tommi - diff_df$mean_ET_dick) / diff_df$mean_ET_dick * 100)

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

#- 2. Plotting Results -####
# Plotting Transpiration ####
plot_t <- diff_df[, c("LON","LAT", "pctdiff_T", "mean_T_tommi","mean_T_dick")]
names(plot_t)[3:5] <- c("Cultivar Difference", "Modern Cultivar", "Historic Cultivar")
plot_t <- pivot_longer(plot_t, cols = 3:5, values_to = "values", names_to = "names")
max_value <- max(max(diff_df$mean_T_tommi), max(diff_df$mean_T_dick))

plot <- filter(plot_t, names == "Modern Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean transpiration\nsum 1990-2020\n(mm season⁻¹)",option = "D", limits = c(0, max_value)) +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 18) +
  theme(strip.text = element_text(face = "bold", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
round(mean(plot$values),1) # 266.5498
round(sd(plot$values),1) # 41.24926

#[1] 266.7
#[1] 41.1

plot <- filter(plot_t, names == "Historic Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean transpiration\nsum 1990-2020\n(mm season⁻¹)",option = "D", limits = c(0, max_value)) +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 18) +
  theme(strip.text = element_text(face = "bold", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
round(mean(plot$values),1) # 320.1697
round(sd(plot$values),1) # 32.4
sd(plot$values)/abs(mean(plot$values)) * 100 # 13.5

plot <- filter(plot_t, names == "Cultivar Difference")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "∆Transpiration (%)",option = "D", limits = c(-100, 0)) +
  #scale_fill_gradient2(name = "∆Transpiration (%)", low = "#0571b0", mid = "#f7f7f7", high = "#ca0020",
  #                     midpoint = 0, limits = c(-100,100)) +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 18) +
  theme(strip.text = element_text(face = "bold", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
round(mean(plot$values),2) # -17
round(sd(plot$values),1) # 8.4
sd(plot$values)/abs(mean(plot$values)) * 100 # 49.70216

# Plotting leaf area index ####
plot_LAI <- diff_df[, c("LON","LAT", "pctdiff_LAI_max", "mean_LAI_max_tommi","mean_LAI_max_dick")]
names(plot_LAI)[3:5] <- c("Cultivar Difference", "Modern Cultivar", "Historic Cultivar")
plot_LAI <- pivot_longer(plot_LAI, cols = 3:5, values_to = "values", names_to = "names")
max_value <- max(max(diff_df$mean_LAI_max_tommi), max(diff_df$mean_LAI_max_dick))

plot <- filter(plot_LAI, names == "Modern Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean Yearly\nMaximum LAI\n1990-2020",option = "D", limits = c(0, max_value)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # 3.8
sd(plot$values) # 0.5
sd(plot$values)/abs(mean(plot$values)) * 100 # 13.5

plot <- filter(plot_LAI, names == "Historic Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean Yearly\nMaximum LAI\n1990-2020",option = "D", limits = c(0, max_value)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # 5.93
sd(plot$values) # 0.37
sd(plot$values)/abs(mean(plot$values)) * 100 # 6.27

plot <- filter(plot_LAI, names == "Cultivar Difference")
plot$names[plot$names == "Cultivar Difference"] <- "Mean Maximum LAI"
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  #scale_fill_viridis(name = "∆Maximum LAI (%)",option = "D", limits = c(-100, 0)) +
  scale_fill_gradient2(name = "∆Maximum LAI (%)", low = "#0571b0", mid = "#f7f7f7", high = "#ca0020",
                       midpoint = 0, limits = c(-100,25)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # -36.22
sd(plot$values) # 6.68
sd(plot$values)/abs(mean(plot$values)) * 100 # 18.45

# Plotting above ground biomass #### 
plot_ABG <- diff_df[, c("LON","LAT", "pctdiff_AGBM", "mean_AGBM_tommi","mean_AGBM_dick")]
names(plot_ABG)[3:5] <- c("Cultivar Difference", "Modern Cultivar", "Historic Cultivar")
plot_ABG <- pivot_longer(plot_ABG, cols = 3:5, values_to = "values", names_to = "names")
max_value <- max(max(diff_df$mean_AGBM_tommi), max(diff_df$mean_AGBM_dick))

plot <- filter(plot_ABG, names == "Modern Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean Yearly\nBiomass\n1990-2020 t ha⁻¹",option = "D", limits = c(0, max_value)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # 17.71
sd(plot$values) # 3.04
sd(plot$values)/abs(mean(plot$values)) * 100 # 17.18

plot <- filter(plot_ABG, names == "Historic Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean Yearly\nBiomass\n1990-2020 t ha⁻¹",option = "D", limits = c(0, max_value)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # 18.13
sd(plot$values) #  2.02
sd(plot$values)/abs(mean(plot$values)) * 100 # 11.12

plot <- filter(plot_ABG, names == "Cultivar Difference")
plot$names[plot$names == "Cultivar Difference"] <- "Mean Above Ground Biomass"
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  #scale_fill_viridis(name = "∆Biomass (%)",option = "D", limits = c(-100, 26)) +
  scale_fill_gradient2(name = "∆Biomass (%)", low = "#0571b0", mid = "#f7f7f7", high = "#ca0020",
                       midpoint = 0, limits = c(-100,25)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # -2.84
sd(plot$values) #  10.98
sd(plot$values)/abs(mean(plot$values)) * 100 # 387.21

# Plotting Root hydraulic conductance ####
plot_ABG <- diff_df[, c("LON","LAT", "pctdiff_RootHydrCond", "mean_RootHydrCond_tommi","mean_RootHydrCond_dick")]
names(plot_ABG)[3:5] <- c("Cultivar Difference", "Modern Cultivar", "Historic Cultivar")
plot_ABG <- pivot_longer(plot_ABG, cols = 3:5, values_to = "values", names_to = "names")
max_value <- max(max(diff_df$mean_RootHydrCond_tommi), max(diff_df$mean_RootHydrCond_dick))

plot <- filter(plot_ABG, names == "Modern Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean Yearly\nMaximum Krs\n1990-2020)" ,option = "D", limits = c(0, max_value*0.014)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # 0.0002289236
sd(plot$values) #  0.0001252863


plot <- filter(plot_ABG, names == "Historic Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean Yearly\nMaximum Krs\n1990-2020)",option = "D", limits = c(0, max_value*0.085)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # 0.001272098
sd(plot$values) #  0.0006223623


plot <- filter(plot_ABG, names == "Cultivar Difference")
plot$names[plot$names == "Cultivar Difference"] <- "Mean Root Hydraulic Conductance"
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  #scale_fill_viridis(name = "∆Krs (%)",option = "D", limits = c(-100, -0)) +
  scale_fill_gradient2(name = "∆Krs (%)", low = "#0571b0", mid = "#f7f7f7", high = "#ca0020",
                       midpoint = 0, limits = c(-100,25)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))

mean(plot$values) # -82.46
sd(plot$values) #  2.65

# Plotting Potential ETranspiration ####
plot_ABG <- diff_df[, c("LON","LAT", "pctdiff_ET", "mean_ET_tommi","mean_ET_dick")]
names(plot_ABG)[3:5] <- c("Cultivar Difference", "Modern Cultivar", "Historic Cultivar")
plot_ABG <- pivot_longer(plot_ABG, cols = 3:5, values_to = "values", names_to = "names")
max_value <- max(max(diff_df$mean_ET_tommi), max(diff_df$mean_ET_dick))

plot <- filter(plot_ABG, names == "Modern Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean seasonal ET\n1990-2020 (mm)" ,option = "D", limits = c(0, max_value)) +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 18) +
  theme(strip.text = element_text(face = "bold", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # 332.7
sd(plot$values) # 42.6
sd(plot$values)/abs(mean(plot$values)) * 100 # 12.8

plot <- filter(plot_ABG, names == "Historic Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean seasonal ET\n1990-2020 (mm)",option = "D", limits = c(0, max_value)) +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 18) +
  theme(strip.text = element_text(face = "bold", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # 378
sd(plot$values) # 34.7
sd(plot$values)/abs(mean(plot$values)) * 100 # 9.1

plot <- filter(plot_ABG, names == "Cultivar Difference")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  #scale_fill_viridis(name = "∆Potential ET (%)",option = "D", limits = c(-100, -0)) +
  scale_fill_gradient2(name = "∆ET (%)", low = "#0571b0", mid = "#f7f7f7", high = "#ca0020",
                       midpoint = 0, limits = c(-100, 25)) +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 18) +
  theme(strip.text = element_text(face = "bold", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # -12.5
sd(plot$values) # 6.9
sd(plot$values)/abs(mean(plot$values)) * 100 # 55

# Plotting Outflow ####
plot_ABG <- diff_df[, c("LON","LAT", "pctdiff_Outflow", "mean_Outflow_tommi","mean_Outflow_dick")]
names(plot_ABG)[3:5] <- c("Cultivar Difference", "Modern Cultivar", "Historic Cultivar")
plot_ABG <- pivot_longer(plot_ABG, cols = 3:5, values_to = "values", names_to = "names")
max_value <- max(max(diff_df$mean_PotET0_tommi), max(diff_df$mean_PotET0_dick))

plot <- filter(plot_ABG, names == "Modern Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean Yearly\nSubsurface Outflow\n1990-2020 (mm)" ,option = "D", limits = c(0, max_value)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))

plot <- filter(plot_ABG, names == "Historic Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean Yearly\nSubsurface Outflow\n1990-2020 (mm)",option = "D", limits = c(0, max_value)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))

plot <- filter(plot_ABG, names == "Cultivar Difference")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  #scale_fill_viridis(name = "∆Subsurface Outflow (%)",option = "D", limits = c(-25, 25)) +
  scale_fill_gradient2(
    name = "∆Subsurface Outflow (%)",
    low      = "blue",
    mid      = "white",
    high     = "red",
    midpoint = 0,
    limits   = c(-25, 25)
  ) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))

# Plotting WUE ####
plot_WUE <- diff_df[, c("LON","LAT", "pctdiff_WUE", "mean_WUE_tommi","mean_WUE_dick")]
names(plot_WUE)[3:5] <- c("Cultivar Difference", "Modern Cultivar", "Historic Cultivar")
plot_WUE <- pivot_longer(plot_WUE, cols = 3:5, values_to = "values", names_to = "names")
max_value <- max(max(diff_df$mean_WUE_tommi), max(diff_df$mean_WUE_dick))

plot <- filter(plot_WUE, names == "Modern Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean WUE 1990-\n2020 (kg mm⁻¹)",option = "D", limits = c(0, max_value)) +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 18) +
  theme(strip.text = element_text(face = "bold", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),        
        legend.text = element_text(size=16),  
        legend.title = element_text(size=18,face = "bold"),                  
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # 0.67
sd(plot$values) # 0.1
sd(plot$values)/abs(mean(plot$values)) * 100 # 15.52

plot <- filter(plot_WUE, names == "Historic Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean WUE 1990-\n2020 (kg mm⁻¹)",option = "D", limits = c(0, max_value)) +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 18) +
  theme(strip.text = element_text(face = "bold", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),        
        legend.text = element_text(size=16),  
        legend.title = element_text(size=18,face = "bold"),                  
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) #  0.57
sd(plot$values) # 0.09
sd(plot$values)/abs(mean(plot$values)) * 100 # 15.11

plot <- filter(plot_WUE, names == "Cultivar Difference")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  #scale_fill_viridis(name = "∆WUE (%)",option = "D", limits = c(0, 50)) +
  scale_fill_gradient2(name = "∆WUE (%)", low = "#0571b0", mid = "#f7f7f7", high = "#ca0020",
                       midpoint = 0, limits = c(-100,100)) +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 18) +
  theme(strip.text = element_text(face = "bold", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),        
        legend.text = element_text(size=16),  
        legend.title = element_text(size=18,face = "bold"),                  
        legend.background  = element_rect(fill = alpha("white", 0.99)))
mean(plot$values) # 16.71
sd(plot$values) # 3.77
sd(plot$values)/abs(mean(plot$values)) * 100 # 22.54

plot <- filter(plot_WUE, names == "Cultivar Difference")
plot$names[plot$names == "Cultivar Difference"] <- "Mean Water Use Efficiency"
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  #scale_fill_viridis(name = "∆WUE (%)",option = "D", limits = c(0, 50)) +
  scale_fill_gradient2(name = "∆WUE (%)", low = "#0571b0", mid = "#f7f7f7", high = "#ca0020",
                       midpoint = 0, limits = c(-100,100)) +
  coord_sf(xlim = c(-11.00348,30.25406), ylim = c(35.99601,67.71814), expand = FALSE) +
  labs(x = "", y = "") +
  theme_bw(base_size = 18) +
  theme(strip.text = element_text(face = "bold", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),        
        legend.text = element_text(size=16),  
        legend.title = element_text(size=18,face = "bold"),                  
        legend.background  = element_rect(fill = alpha("white", 0.99)))

# Plotting Actual ET ####
plot_ABG <- diff_df[, c("LON","LAT", "pctdiff_PotET0", "mean_PotET0_tommi","mean_PotET0_dick")]
names(plot_ABG)[3:5] <- c("Cultivar Difference", "Modern Cultivar", "Historic Cultivar")
plot_ABG <- pivot_longer(plot_ABG, cols = 3:5, values_to = "values", names_to = "names")
max_value <- max(max(diff_df$mean_PotET0_tommi), max(diff_df$mean_PotET0_dick))

plot <- filter(plot_ABG, names == "Modern Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean Yearly\nPotential ET\n1990-2020 (mm)" ,option = "D", limits = c(0, max_value)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))

plot <- filter(plot_ABG, names == "Historic Cultivar")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  scale_fill_viridis(name = "Mean Yearly\nPotential ET\n1990-2020 (mm)",option = "D", limits = c(0, max_value)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))

plot <- filter(plot_ABG, names == "Cultivar Difference")
ggplot(plot, aes(x = LON, y = LAT, fill = values)) +
  geom_sf(data = shpf, fill = "grey90", color = "black", size = 0.3) +
  geom_raster(interpolate = FALSE) +           # or use geom_tile()
  facet_wrap(~ names, ncol = 3) +
  #scale_fill_viridis(name = "∆Potential ET (%)",option = "D", limits = c(-100, -0)) +
  scale_fill_gradient2(name = "∆Potential ET (%)", low = "#0571b0", mid = "#f7f7f7", high = "#ca0020",
                       midpoint = 0, limits = c(-100, 25)) +
  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", size = 18),
        axis.text  = element_text(size = 6),
        legend.position    = c(0.01, 0.99),        
        legend.justification = c(0, 1),            
        legend.background  = element_rect(fill = alpha("white", 0.99)))