################################################################################
# Process data & create figures ################################################
################################################################################

# load required R packages

library(cowplot)
library(dplyr)
library(ggplot2)
library(grid)
library(patchwork)
library(RColorBrewer)
library(reshape2)
library(sf)
library(trend)
library(tidyr)

# load required input datasets

home         = "/.../.../" # define the path to your home directory
ampts_values = read.csv(paste0(home,"AMPTs_values.csv"))
ampts_dates  = read.csv(paste0(home,"AMPTs_dates.csv"))
aic          = read.csv(paste0(home,"AIC.csv"))
metadata     = read.csv(paste0(home,"Metadata.csv"))

# GWL-REA classification is not redistributed with this repository due to its data policy,
# please download the dataset from: https://opendata.dwd.de/climate_environment/CDC/event_catalogues/europe/weather_types/GWL-REA/v1.6/historical/
# GWL_cat      = read.csv(paste0(home,"GWL-REA_v1.6_ERA5_GLP_catalogue.csv"),sep=";", header=FALSE)
# GWL_des      = read.csv(paste0(home,"GWL-REA_v1.6_ERA5_GLP_description.csv"),sep=";", header=TRUE)

# add spatial data 
# classification of natural regions in Germany of Meynem and Schmithausen (1959)
# is not with this repository due to its data policy

WGS84          = ("+proj=longlat +datum=WGS84")
LAEA           = ("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")

################################################################################
# Process AMPTs values #########################################################
################################################################################

# for each station (CODE) and each pcp duration (DUR*), check whether at least 
# 30 yrs of valid data are available. Missing values are coded as -99.9 and are 
# treated as NA. If a duration has fewer than 30 valid yrs at a station, all 
# values of this duration for that station are set to NA. Otherwise, the data 
# are kept unchanged

dur_cols               <- grep("^DUR", names(ampts_values), value = TRUE)
ampts_values[dur_cols] <- lapply(ampts_values[dur_cols], \(x) replace(x, x==-99.9, NA_real_))

for (code in unique(ampts_values$CODE)){
  idx <- ampts_values$CODE == code
  for (dur in dur_cols){
    n_valid <- sum(!is.na(ampts_values[idx,dur]))
    if(n_valid < 30){
      ampts_values[idx,dur] <- NA_real_
      cat(
        "Station:", code,
        "| Duration:", dur,
        "| Valid yrs:", n_valid, "\n"
      )
    }
  }
}

################################################################################
# Figure 1 #####################################################################
################################################################################

stations_per_year <- ampts_values %>%
  group_by(YYYY) %>%
  summarise(
    across(
      all_of(c("DUR5","DUR420","DUR1440")),
      ~ sum(!is.na(.) & . != -99.9)
    ),
    .groups = "drop"
  )

stations_per_year <- melt(
  stations_per_year,
  id = "YYYY",
  variable.name = "DURATION",
  value.name    = "N_STATIONS"
)

stations_per_year <- stations_per_year %>%
  select(YYYY, DURATION, N_STATIONS) %>%
  pivot_wider(
    names_from = DURATION,
    values_from = N_STATIONS
  ) %>%
  arrange(YYYY)

Figure1 <- ggplot(stations_per_year, aes(x = YYYY)) +
  # Area from 0 to 5-minute data
  geom_ribbon(aes(ymin = 0,ymax = DUR5),fill = "#95C4C0",alpha = 0.5) +
  # Area between 5-minute and subdaily data
  geom_ribbon(aes(ymin = DUR5,ymax = DUR420),fill = "#307970",alpha = 0.5) +
  # Area between subdaily and daily data
  geom_ribbon(aes(ymin = DUR420,ymax = DUR1440),fill = "#B78543",alpha = 0.5) +
  # Original lines
  geom_line(aes(y = DUR5, color = "DUR5"),linewidth = 0.5) +
  geom_line(aes(y = DUR420, color = "DUR420"),linewidth = 0.5) +
  geom_line(aes(y = DUR1440, color = "DUR1440"),linewidth = 0.5) +
  # Line colors and legend labels
  scale_color_manual(values = c("DUR5" = "#95C4C0","DUR420" = "#307970","DUR1440" = "#B78543"),
                     breaks = c("DUR5", "DUR420","DUR1440"),
                     labels = c("5-minute data","subdaily data","daily data")) +
  labs(title = "",x = "\nYears",y = "Number of stations\n",color = "") +
  scale_x_continuous(limits = c(1895, 2025),
                     breaks = seq(1900, 2020, by = 20),
                     expand = c(0, 0) ) +
  coord_cartesian(ylim = c(0, 1000)) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.title = element_text(size = 8),
    legend.text = element_text(size = 8),
    axis.title.x = element_text(size = 8),
    axis.title.y = element_text(size = 8),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8)
  )

plot(Figure1)
ggsave(paste0(home,"Figure1.png"),Figure1,width = 14,height = 14,units = "cm",dpi = 300)

################################################################################
# Figure 2 #####################################################################
################################################################################

ampts_values_long <- melt(ampts_values,id=c("YYYY","CODE"), variable="DURATION",value.name = "AMPTs")

ampts_max_st_dur  <- ampts_values_long %>%
  group_by(CODE, DURATION) %>%
  summarise(
    AMPTs_max = if (all(is.na(AMPTs))) NA_real_ else max(AMPTs, na.rm = TRUE),
    .groups = "drop"
  )

ampts_max_summary <- ampts_max_st_dur %>%
  group_by(DURATION) %>%
  summarise(
    Mean   = mean(AMPTs_max,na.rm=TRUE),
    Median = median(AMPTs_max, na.rm=TRUE),
    P01    = quantile(AMPTs_max, 0.01, na.rm=TRUE),
    P99    = quantile(AMPTs_max, 0.99, na.rm=TRUE), ,.groups = "drop"
  )

ampts_max_summary <- subset(ampts_max_summary,
                            ampts_max_summary$DURATION %in% c("DUR5","DUR10",
                                                              "DUR15","DUR30",
                                                              "DUR60","DUR120",
                                                              "DUR360","DUR720",
                                                              "DUR1440"))
ampts_max_summary$DURATION <- as.numeric(gsub("DUR","",ampts_max_summary$DURATION))


# create df with max recordable intensities for various measurement systems for 
# selected durations

recordable_int_durs <- c(1, 5, 10, 15, 30, 60, 120, 360, 720, 1440)
recordable_int      <- data.frame(
  DURATION        = recordable_int_durs,
  RainRec_daily   = 2.67 * recordable_int_durs,
  RainRec_weekly  = 0.38 * recordable_int_durs,
  RainRec_monthly = 3.33 * recordable_int_durs,
  Pluvio_early    = 10 * recordable_int_durs,
  Pluvio_modern   = 20 * recordable_int_durs
)

ampts_max_summary <- left_join(recordable_int,ampts_max_summary, by="DURATION")

draw_key_ampts <- function(data, params, size) {
  grobTree(
    rectGrob(width = unit(0.9, "npc"), height = unit(0.65, "npc"),
             gp = gpar(fill = "#3876C8", col = NA, alpha = 0.3)),
    segmentsGrob(x0 = unit(0.05, "npc"), x1 = unit(0.95, "npc"),
                 y0 = unit(0.5, "npc"), y1 = unit(0.5, "npc"),
                 gp = gpar(col = "#145CE0", lwd = 1.5)),
    pointsGrob(x = unit(0.5, "npc"), y = unit(0.5, "npc"),
               pch = 19, size = unit(2.5, "mm"),
               gp = gpar(col = "#145CE0"))
  )
}

Figure2 <- ggplot(ampts_max_summary, aes(x = DURATION)) +
  geom_line(data = recordable_int, aes(x = DURATION, y = RainRec_weekly, color = "weekly rain recorder"), linewidth = 0.8, inherit.aes = FALSE) +
  geom_point(data = recordable_int, aes(x = DURATION, y = RainRec_weekly, color = "weekly rain recorder"), size = 2, shape = 19, inherit.aes = FALSE) +
  geom_line(data = recordable_int, aes(x = DURATION, y = RainRec_daily, color = "daily rain recorder"), linetype = 2, linewidth = 0.8, inherit.aes = FALSE) +
  geom_point(data = recordable_int, aes(x = DURATION, y = RainRec_daily, color = "daily rain recorder"), size = 2, shape = 19, inherit.aes = FALSE) +
  geom_line(data = recordable_int, aes(x = DURATION, y = Pluvio_early, color = "early pluviometer"), linewidth = 0.8, inherit.aes = FALSE) +
  geom_point(data = recordable_int, aes(x = DURATION, y = Pluvio_early, color = "early pluviometer"), size = 2, shape = 19, inherit.aes = FALSE) +
  geom_line(data = recordable_int, aes(x = DURATION, y = Pluvio_modern, color = "modern pluviometer"), linetype = 2, linewidth = 0.8, inherit.aes = FALSE) +
  geom_point(data = recordable_int, aes(x = DURATION, y = Pluvio_modern, color = "modern pluviometer"), size = 2, shape = 19, inherit.aes = FALSE) +
  geom_ribbon(aes(ymin = P01, ymax = P99), fill = "#3876C8", alpha = 0.3, na.rm = TRUE, show.legend = FALSE) +
  geom_line(aes(y = Median, color = "recorded AMPTs"), linewidth = 0.8, na.rm = TRUE, key_glyph = draw_key_ampts) +
  geom_point(aes(y = Median), color = "#145CE0", size = 2, shape = 19, na.rm = TRUE, show.legend = FALSE) +
  scale_color_manual(
    values = c("early pluviometer" = "#307970", "modern pluviometer" = "#307970",
               "daily rain recorder" ="#B78543", "weekly rain recorder" = "#B78543",
               "recorded AMPTs" = "#145CE0"),
    breaks = c("daily rain recorder","weekly rain recorder",
               "early pluviometer","modern pluviometer","recorded AMPTs")) +
  scale_x_log10(breaks = recordable_int$DURATION) +
  scale_y_log10() +
  labs(x = "\nDuration [min]", y = "Precipitation [mm]\n", color = "") +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.box.just = "center",
    legend.title = element_blank(),
    legend.spacing.x = unit(0.5, "cm"),
    legend.spacing.y = unit(0.2, "cm"),
    legend.margin = margin(0, 0, 0, 0),
    legend.key.width = unit(1.5, "cm"),
    legend.text = element_text(size = 8),
    axis.title = element_text(size = 8),
    axis.text = element_text(size = 8)
  ) +
  guides(
    color = guide_legend(
      nrow = 3,
      byrow = TRUE
    )
  )

plot(Figure2)
ggsave(paste0(home,"Figure2.png"),Figure2,width = 14,height = 16,units = "cm",dpi = 300)

################################################################################
# Figure 3 #####################################################################
################################################################################

# case study for station Erfurt-Weimar "KO00240" with a known date of the system
# replacement on 1 January 1993

one_stat_code  <- "KO00240" # Erfurt-Weimar with change date:01/01/1993
one_stat_data  <- subset(ampts_values,ampts_values$CODE == one_stat_code)
one_stat_data  <- subset(one_stat_data,one_stat_data$YYYY>=1959)
one_stat_chyr  <- 1993

legend_colors <- c("recorded AMPTs"="black",
                   "Sen's slope: 1959-2020"="#145CE0",
                   "Sen's slope: before replacement"="#B78543",
                   "Sen's slope: after replacement"="#307970",
                   "Median: before replacement"="#B78543",
                   "Median: after replacement"="#307970")

legend_linetypes <- c("recorded AMPTs"="solid",
                      "Sen's slope: 1959-2020"="solid",
                      "Sen's slope: before replacement"="solid",
                      "Sen's slope: after replacement"="solid",
                      "Median: before replacement"="dashed",
                      "Median: after replacement"="dashed")

legend_order <- c("recorded AMPTs",
                  "Sen's slope: 1959-2020",
                  "Sen's slope: before replacement",
                  "Sen's slope: after replacement",
                  "Median: before replacement",
                  "Median: after replacement")

get_sen_results <- function(data){
  data <- data %>% arrange(YYYY) %>% filter(!is.na(YYYY),!is.na(AMPT))
  if(nrow(data)<3) return(list(slope=NA_real_,intercept=NA_real_))
  ss <- trend::sens.slope(data$AMPT)
  slope <- unname(ss$estimates)
  intercept <- median(data$AMPT-slope*data$YYYY,na.rm=TRUE)
  list(slope=slope,intercept=intercept)
}

make_station_plot <- function(data,duration_col,plot_title,change_year=1993){
  plot_data   <- data %>% 
    filter(YYYY>=1951,YYYY<=2020,!is.na(.data[[duration_col]])) %>%
    transmute(YYYY,AMPT=.data[[duration_col]])
  
  data_before <- plot_data %>%
    filter(YYYY<change_year); data_after <- plot_data %>%
    filter(YYYY>=change_year)
  
  med_before  <- if(all(is.na(data_before$AMPT))) NA_real_ else median(data_before$AMPT,na.rm=TRUE) 
  med_after   <- if(all(is.na(data_after$AMPT))) NA_real_ else median(data_after$AMPT,na.rm=TRUE)
  sen_all     <- get_sen_results(plot_data)
  sen_before  <- get_sen_results(data_before)
  sen_after   <- get_sen_results(data_after)
  slope_label <- sprintf("Sen's slope [mm yr⁻¹]\n1959-2020: %+.2f\n1951–%s: %+.2f\n%s–2020: %+.2f",
                         sen_all$slope,change_year-1,sen_before$slope,change_year,sen_after$slope)
  trend_data  <- data.frame(x=c(min(plot_data$YYYY),min(data_before$YYYY),
                                min(data_after$YYYY)),
                            xend=c(max(plot_data$YYYY),max(data_before$YYYY),max(data_after$YYYY)),
                            y=c(sen_all$intercept+sen_all$slope*min(plot_data$YYYY),
                                sen_before$intercept+sen_before$slope*min(data_before$YYYY),
                                sen_after$intercept+sen_after$slope*min(data_after$YYYY)),
                            yend=c(sen_all$intercept+sen_all$slope*max(plot_data$YYYY),
                                   sen_before$intercept+sen_before$slope*max(data_before$YYYY),
                                   sen_after$intercept+sen_after$slope*max(data_after$YYYY)),
                            legend_label=c("Sen's slope: 1959-2020","Sen's slope: before replacement","Sen's slope: after replacement"))
  
  median_data <- data.frame(x=c(min(data_before$YYYY),min(data_after$YYYY)),
                            xend=c(max(data_before$YYYY),max(data_after$YYYY)),
                            y=c(med_before,med_after),yend=c(med_before,med_after),
                            legend_label=c("Median: before replacement","Median: after replacement"))
  
  ggplot(plot_data,aes(x=YYYY,y=AMPT))+
    geom_line(aes(color="recorded AMPTs",linetype="recorded AMPTs"),linewidth=0.5)+
    geom_segment(data=trend_data,aes(x=x,
                                     xend=xend,
                                     y=y,
                                     yend=yend,
                                     color=legend_label,
                                     linetype=legend_label),linewidth=0.8,inherit.aes=FALSE)+
    geom_segment(data=median_data,aes(x=x,
                                      xend=xend,
                                      y=y,
                                      yend=yend,
                                      color=legend_label,linetype=legend_label),
                 linewidth=0.8,inherit.aes=FALSE)+
    geom_vline(xintercept=change_year,color="grey40",linewidth=0.5)+
    annotate("text",x=change_year+2,y=max(plot_data$AMPT,na.rm=TRUE),
             label=paste0("01/01/",change_year),
             color="grey40",angle=90,hjust=1,size=3.5)+
    annotate("text",x=1948,y=Inf,label=slope_label,hjust=0,vjust=1.15,size=3.2,lineheight=1.1,color="black")+
    scale_color_manual(values=legend_colors,breaks=legend_order,name=NULL)+
    scale_linetype_manual(values=legend_linetypes,breaks=legend_order,name=NULL)+
    scale_y_continuous(breaks=function(x) seq(min(x),max(x),length.out=5),
                       labels=function(x) sprintf("%.0f",x),expand=expansion(mult=c(0.05,0.20)))+
    scale_x_continuous(limits=c(1945,2025),breaks=seq(1960,2020,by=20),expand=c(0,0))+
    labs(title=plot_title,y="AMPTs [mm]\n",x="\nYears")+
    guides(color=guide_legend(nrow=1,byrow=TRUE,override.aes=list(linewidth=0.9)),
           linetype=guide_legend(nrow=1,byrow=TRUE,override.aes=list(linewidth=0.9)))+
    theme_bw()+
    theme(plot.title=element_text(face="bold",size=10),
          legend.title=element_blank(),
          legend.text=element_text(size=8),
          axis.title.x=element_text(size=8),
          axis.title.y=element_text(size=8),
          axis.text.x=element_text(size=8),
          axis.text.y=element_text(size=8),
          legend.position="bottom")
}

p_5min    <- make_station_plot(one_stat_data,"DUR5","Duration: 5 minutes",one_stat_chyr)
p_60min   <- make_station_plot(one_stat_data,"DUR60","Duration: 1 h",one_stat_chyr)
p_1440min <- make_station_plot(one_stat_data,"DUR1440","Duration: 1 day",one_stat_chyr)

Figure3 <- (p_5min|p_60min|p_1440min)+plot_layout(guides="collect")&theme(legend.position="bottom",legend.justification="center",legend.box="horizontal",legend.margin=margin(t=4,r=0,b=0,l=0))
plot(Figure3)
ggsave(filename=paste0(home,"Figure3.png"),plot=Figure3,width=36,height=17,units="cm",dpi=300)

# compare medians before & after

med_before_5min    <- median(one_stat_data$DUR5[one_stat_data$YYYY<one_stat_chyr],na.rm=TRUE)
med_after_5min     <- median(one_stat_data$DUR5[one_stat_data$YYYY>=one_stat_chyr ],na.rm=TRUE)

med_before_60min   <- median(one_stat_data$DUR60[one_stat_data$YYYY<one_stat_chyr],na.rm=TRUE)
med_after_60min    <- median(one_stat_data$DUR60[one_stat_data$YYYY>=one_stat_chyr],na.rm=TRUE)

med_before_1440min <- median(one_stat_data$DUR1440[one_stat_data$YYYY<one_stat_chyr],na.rm=TRUE)
med_after_1440min  <- median(one_stat_data$DUR1440[one_stat_data$YYYY>=one_stat_chyr],na.rm=TRUE)

med_change_5min    <- round(100*(med_after_5min-med_before_5min)/med_before_5min,1)
med_change_60min   <- round(100*(med_after_60min-med_before_60min)/med_before_60min,1)
med_change_1440min <- round(100*(med_after_1440min-med_before_1440min)/med_before_1440min,1)

################################################################################
# Table 1 ######################################################################
################################################################################
# 
# colnames(GWL_cat) <- c("Date","GWL_number","GWL_abbreviation")
# GWL_cat$YYYY      <- substr(GWL_cat$Date,1,4)
# GWL_cat$MM        <- substr(GWL_cat$Date,6,7)
# GWL_cat$DD        <- substr(GWL_cat$Date,9,10)
# GWL_cat           <- subset(GWL_cat,YYYY<=2020) 
# 
# GWL_fq            <- prop.table(table(GWL_cat$GWL_abbreviation))*100
# GWL_fq            <- as.data.frame(GWL_fq)
# colnames(GWL_fq)  <- c("GWL_abbreviation","Frequency")
# GWL_fq$Frequency  <- round(GWL_fq$Frequency,1)
# 
# missing <- round(100 - sum(GWL_fq$Frequency), 10)
# x       <- as.integer(round(missing / 0.1))
# top_idx <- order(GWL_fq$Frequency, decreasing = TRUE)[1:x]
# 
# GWL_fq$Frequency[top_idx] <- GWL_fq$Frequency[top_idx] + 0.1
# Table1                    <- left_join(GWL_des,GWL_fq, by="GWL_abbreviation")
# write.csv(Table1,paste0(home,"GWL-REA_FQ.csv"),row.names = FALSE)

################################################################################
# Figure 4 #####################################################################
################################################################################

ampts_med_st_dur <- ampts_values_long %>%
  group_by(CODE, DURATION) %>%
  summarise(
    AMPTs_median = median(AMPTs, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  pivot_wider(
    names_from = DURATION,
    values_from = AMPTs_median
  )

ampts_med_st_dur <- left_join(ampts_med_st_dur, metadata, by = "CODE")
ampts_med_st_dur <- st_as_sf(ampts_med_st_dur, coords = c("LON","LAT"), crs = WGS84)
ampts_med_st_dur <- st_transform(ampts_med_st_dur, crs = LAEA)

# colorbar solution 1:
# define breaks,labels & colors

BreaksDefined_med_DUR5     <- round(c(seq(3,8,length.out=10),Inf),1)
BreaksDefined_med_DUR15    <- round(c(seq(7,16,length.out=10),Inf),1)
BreaksDefined_med_DUR30    <- round(c(seq(9,21,length.out=10),Inf),1)
BreaksDefined_med_DUR60    <- round(c(seq(13,24,length.out=10),Inf),1)
BreaksDefined_med_DUR180   <- round(c(seq(15,34,length.out=10),Inf),1)
BreaksDefined_med_DUR420   <- round(c(seq(15,45,length.out=10),Inf),1)
BreaksDefined_med_DUR1440  <- round(c(seq(20,70,length.out=10),Inf),1)
BreaksDefined_med_DUR4320  <- round(c(seq(30,90,length.out=10),Inf),1)
BreaksDefined_med_DUR10080 <- round(c(seq(40,120,length.out=10),Inf),1)

LabelsDefined_med_DUR5     <- c(format(BreaksDefined_med_DUR5[1:(length(BreaksDefined_med_DUR5) - 1)], nsmall = 1))
LabelsDefined_med_DUR15    <- c(format(BreaksDefined_med_DUR15[1:(length(BreaksDefined_med_DUR15) - 1)], nsmall = 1))
LabelsDefined_med_DUR30    <- c(format(BreaksDefined_med_DUR30[1:(length(BreaksDefined_med_DUR30) - 1)], nsmall = 1))
LabelsDefined_med_DUR60    <- c(format(BreaksDefined_med_DUR60[1:(length(BreaksDefined_med_DUR60) - 1)], nsmall = 1))
LabelsDefined_med_DUR180   <- c(format(BreaksDefined_med_DUR180[1:(length(BreaksDefined_med_DUR180) - 1)], nsmall = 1))
LabelsDefined_med_DUR420   <- c(format(BreaksDefined_med_DUR420[1:(length(BreaksDefined_med_DUR420) - 1)], nsmall = 1))
LabelsDefined_med_DUR1440  <- c(format(BreaksDefined_med_DUR1440[1:(length(BreaksDefined_med_DUR1440) - 1)], nsmall = 1))
LabelsDefined_med_DUR4320  <- c(format(BreaksDefined_med_DUR4320[1:(length(BreaksDefined_med_DUR4320) - 1)], nsmall = 1))
LabelsDefined_med_DUR10080 <- c(format(BreaksDefined_med_DUR10080[1:(length(BreaksDefined_med_DUR10080) - 1)], nsmall = 1))

med_colors   <- c("#FFF9C4","#D8ECBC","#B1E0B4","#91CBB0","#77AEB0","#5D91B0","#3876C8","#145CE0","#124BCF","#334393","#543C58")
dur_selected <- c("DUR5","DUR15","DUR30","DUR60","DUR180","DUR420","DUR1440","DUR4320","DUR10080")
dur_labels   <- c(
  DUR5     = "5 minutes",
  DUR15    = "15 minutes",
  DUR30    = "30 minutes",
  DUR60    = "1 hour",
  DUR180   = "3 hours",
  DUR420   = "7 hours",
  DUR720   = "12 hours",
  DUR1440  = "1 day",
  DUR4320  = "3 days",
  DUR10080 = "7 days"
)

for (i in 1:length(dur_selected)){
  
  plt_col       <- dur_selected[i]
  plt_input     <- ampts_med_st_dur
  plt_input     <- subset(plt_input,plt_input[[plt_col]] !=0)
  
  plt_breaks    <- get(paste0("BreaksDefined_med_",dur_selected[i]))
  plt_labels    <- get(paste0("LabelsDefined_med_",dur_selected[i]))
  plt_input$cut <- cut(plt_input[[plt_col]],
                      breaks = plt_breaks, 
                      include.lowest = TRUE,
                      right = FALSE, 
                      labels = plt_labels)
  
  # colorbar solution 2:
  # # define breaks & labels & colorsbased on mix/max values
  #  plt_breaks    <- round(seq(floor(min(plt_input[[plt_col]],na.rm=TRUE)),
  #                             ceiling(max(plt_input[[plt_col]],na.rm=TRUE)),
  #                             length.out = 10),1)
  #  plt_labels    <- c(format(plt_breaks[1:(length(plt_breaks) - 1)], nsmall = 1))
  #  plt_input$cut <- cut(plt_input[[plt_col]],
  #                       breaks = plt_breaks, 
  #                       include.lowest = TRUE,
  #                       right = FALSE, 
  #                       labels = plt_labels)
   
  plt <- ggplot() +
    geom_sf(data = GER_regions, fill = NA, colour = "black", linewidth = 0.35) +
    geom_sf(data = plt_input, aes(fill = cut),size=3, pch = 21, show.legend = TRUE) +
    scale_fill_manual(values = med_colors, 
                       labels = {
                       labs <- plt_labels
                       labs[length(labs)] <- paste0("≥ ", labs[length(labs)])
                       labs
                     },
                     drop = FALSE, na.value = "red", na.translate = TRUE,
                     guide = guide_legend(reverse = TRUE)) +
    theme_bw() +
    theme(legend.position = 'right') +
    theme(axis.title = element_blank(), 
          axis.text = element_text(size = 8),
          legend.text = element_text(size = 10),
          legend.title = element_text(size = 10),
          #legend.text = element_text(hjust),
          #legend.key.size = unit(1.2, "cm"),  
          #plot.margin = unit(c(0.25, 0, 0.25, 0), "cm"),
          plot.title = element_text(color = "black", size = 10, face = "bold"),
          plot.subtitle = element_text(color = "black", size = 10, face = "bold")) +
    labs(fill="[mm]") +
    ggtitle(paste0("Duration: ", dur_labels[dur_selected[i]]))+
    geom_sf_label(data=GER_regions,aes(label=ID),size=3, alpha=0.8)
  #plot(plt)
  
  assign(paste0("Figure4_",i), plt)
  
}

# combine part 1-9 to final Figure 4

Figure4 <- plot_grid(Figure4_1,Figure4_2,Figure4_3,Figure4_4,Figure4_5,Figure4_6,
                              Figure4_7,Figure4_8,Figure4_9,
                              ncol=3,nrow=3)+
  theme(plot.background = element_rect(fill="white",colour = "white"))

#plot(Figure4)
ggsave(paste0(home,"Figure4.png"),Figure4,width = 36,height = 36,units = "cm",dpi = 300)

################################################################################
# Figure 5 #####################################################################
################################################################################

ampts_cv_st_dur <- ampts_values_long %>%
  group_by(CODE, DURATION) %>%
  summarise(
    AMPTs_cv = {
      x <- AMPTs[!is.na(AMPTs)]
      if (length(x) < 2 || mean(x) == 0) {
        NA_real_
      } else {
        100 * sd(x) / mean(x)
      }
    },
    .groups = "drop"
  ) %>%
  pivot_wider(
    names_from = DURATION,
    values_from = AMPTs_cv
  )

ampts_cv_st_dur <- left_join(ampts_cv_st_dur, metadata, by = "CODE")
ampts_cv_st_dur <- st_as_sf(ampts_cv_st_dur, coords = c("LON","LAT"), crs = WGS84)
ampts_cv_st_dur <- st_transform(ampts_cv_st_dur, crs = LAEA)

# define legend properties for plotting (breaks,labels & colors) - median

BreaksDefined_cv <- c(10,20,30,40,50,Inf)
LabelsDefined_cv <- c(format(BreaksDefined_cv[1:(length(BreaksDefined_cv) - 1)], nsmall = 1))

cv_colors   <- c("#E3D7C8","#DDC8AC","#CFA56D","#825822","#543005")
dur_selected <- c("DUR5","DUR15","DUR30","DUR60","DUR180","DUR420","DUR1440","DUR4320","DUR10080")

for (i in 1:length(dur_selected)){
  
  plt_col       <- dur_selected[i]
  plt_input     <- ampts_cv_st_dur
  plt_input     <- subset(plt_input,plt_input[[plt_col]] !=0)

  plt_breaks    <- BreaksDefined_cv
  plt_labels    <- LabelsDefined_cv
  plt_input$cut <- cut(plt_input[[plt_col]],
                       breaks = plt_breaks, 
                       include.lowest = TRUE,
                       right = FALSE, 
                       labels = plt_labels)
  
  plt <- ggplot() +
    geom_sf(data = GER_regions, fill = NA, colour = "black", linewidth = 0.35) +
    geom_sf(data = plt_input, aes(fill = cut),size=3, pch = 21, show.legend = TRUE) +
    scale_fill_manual(values = cv_colors, 
                      labels = {
                        labs <- plt_labels
                        labs[length(labs)] <- paste0("≥ ", labs[length(labs)])
                        labs
                      },
                      drop = FALSE, na.value = "red", na.translate = TRUE,
                      guide = guide_legend(reverse = TRUE)) +
    theme_bw() +
    theme(legend.position = 'right') +
    theme(axis.title = element_blank(), 
          axis.text = element_text(size = 8),
          legend.text = element_text(size = 10),
          legend.title = element_text(size = 10),
          #legend.text = element_text(hjust),
          #legend.key.size = unit(1.2, "cm"),  
          #plot.margin = unit(c(0.25, 0, 0.25, 0), "cm"),
          plot.title = element_text(color = "black", size = 10, face = "bold"),
          plot.subtitle = element_text(color = "black", size = 10, face = "bold")) +
    labs(fill="[%]") +
    ggtitle(paste0("Duration: ", dur_labels[dur_selected[i]]))+
    geom_sf_label(data=GER_regions,aes(label=ID),size=3, alpha=0.8)
  #plot(plt)
  
  assign(paste0("Figure5_",i), plt)
  
}

# combine part 1-9 to final Figure 5

Figure5 <- plot_grid(Figure5_1,Figure5_2,Figure5_3,Figure5_4,Figure5_5,Figure5_6,
                     Figure5_7,Figure5_8,Figure5_9,ncol=3,nrow=3)+
  theme(plot.background = element_rect(fill="white",colour = "white"))

#plot(Figure5)
ggsave(paste0(home,"Figure5.png"),Figure5,width = 36,height = 36,units = "cm",dpi = 300)

################################################################################
# Process AMPTs dates ##########################################################
################################################################################

# for each station (CODE) and each pcp duration (DUR*), check whether at least 
# 30 yrs of valid data are available. Missing values are coded as -99.9 and are 
# treated as NA. If a duration has fewer than 30 valid yrs at a station, all 
# values of this duration for that station are set to NA. Otherwise, the data 
# are kept unchanged

dur_cols               <- grep("^DUR", names(ampts_dates), value = TRUE)
ampts_dates[dur_cols]  <- lapply(ampts_dates[dur_cols], \(x) replace(x, x==-99, NA_real_))

for (code in unique(ampts_dates$CODE)){
  idx <- ampts_dates$CODE == code
  for (dur in dur_cols){
    n_valid <- sum(!is.na(ampts_dates[idx,dur]))
    if(n_valid < 30){
      ampts_dates[idx,dur] <- NA_real_
      cat(
        "Station:", code,
        "| Duration:", dur,
        "| Valid yrs:", n_valid, "\n"
      )
    }
  }
}

################################################################################
# Figure 6 #####################################################################
################################################################################

ampts_dates_long      <- melt(ampts_dates,id = c("YYYY", "CODE"),variable.name = "DURATION", value.name = "Date")
ampts_dates_long$YYYY <- substr(ampts_dates_long$Date, 1, 4)
ampts_dates_long$MM   <- substr(ampts_dates_long$Date, 5, 6)
ampts_dates_long$DD   <- substr(ampts_dates_long$Date, 7, 8)
ampts_dates_long      <- ampts_dates_long[, c("CODE","DURATION","YYYY","MM","DD")]

ampts_dates_gwl       <- left_join(ampts_dates_long,GWL_cat,by = c("YYYY", "MM", "DD"))

ampts_dates_gwl <- ampts_dates_gwl %>%
  filter(
    !is.na(GWL_abbreviation),
    !is.na(DURATION),
    !is.na(YYYY),
    !is.na(MM),
    !is.na(DD)
  ) %>%
  mutate(
    DURATION = as.character(DURATION)
  )

GWL_fq_SB_durations <- ampts_dates_gwl %>%
  count(
    DURATION,
    GWL_abbreviation,
    name = "Count"
  ) %>%
  group_by(DURATION) %>%
  mutate(
    Frequency = 100 * Count / sum(Count)
  ) %>%
  ungroup()

# CP-based frequencies

GWL_fq_DB_durations <- ampts_dates_gwl %>%
  distinct(
    YYYY,
    MM,
    DD,
    DURATION,
    GWL_abbreviation
  ) %>%
  count(
    DURATION,
    GWL_abbreviation,
    name = "Count"
  ) %>%
  group_by(DURATION) %>%
  mutate(
    Frequency = 100 * Count / sum(Count)
  ) %>%
  ungroup()

GWL_fq_SB_grouped <- GWL_fq_SB_durations %>%
  mutate(
    Duration_group = case_when(
      DURATION %in% c("DUR5", "DUR15", "DUR30", "5", "15", "30") ~
        "short-duration",
      
      DURATION %in% c("DUR60", "DUR180", "DUR420", "60", "180", "420") ~
        "medium-duration",
      
      DURATION %in% c("DUR1440", "DUR4320", "DUR10080",
                      "1440", "4320", "10080") ~
        "long-duration",
      
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(Duration_group)) %>%
  group_by(
    Duration_group,
    GWL_abbreviation
  ) %>%
  summarise(
    Count = sum(Count),
    .groups = "drop"
  ) %>%
  group_by(Duration_group) %>%
  mutate(
    Frequency = 100 * Count / sum(Count)
  ) %>%
  ungroup() %>%
  mutate(
    Statistic = "SB-CPFQ"
  )

GWL_fq_DB_grouped <- GWL_fq_DB_durations %>%
  mutate(
    Duration_group = case_when(
      DURATION %in% c("DUR5", "DUR15", "DUR30", "5", "15", "30") ~
        "short-duration",
      
      DURATION %in% c("DUR60", "DUR180", "DUR420", "60", "180", "420") ~
        "medium-duration",
      
      DURATION %in% c("DUR1440", "DUR4320", "DUR10080",
                      "1440", "4320", "10080") ~
        "long-duration",
      
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(Duration_group)) %>%
  group_by(
    Duration_group,
    GWL_abbreviation
  ) %>%
  summarise(
    Count = sum(Count),
    .groups = "drop"
  ) %>%
  group_by(Duration_group) %>%
  mutate(
    Frequency = 100 * Count / sum(Count)
  ) %>%
  ungroup() %>%
  mutate(
    Statistic = "DB-CPFQ"
  )

# normalised CP-based frequencies

GWL_overall_frequency <- Table1 %>%
  select(
    GWL_abbreviation,
    Overall_Frequency = Frequency
  )

GWL_fq_NSB_grouped <- GWL_fq_SB_grouped %>%
  left_join(GWL_overall_frequency, by = "GWL_abbreviation") %>%
  mutate(
    NormalisedFrequency = Frequency / Overall_Frequency,
    Statistic = "NSB-CPFQ"
  )

GWL_fq_NDB_grouped <- GWL_fq_DB_grouped %>%
  left_join(GWL_overall_frequency, by = "GWL_abbreviation") %>%
  mutate(
    NormalisedFrequency = Frequency / Overall_Frequency,
    Statistic = "NDB-CPFQ"
  )

# plotting data

GWL_fq_NSB_plot <- GWL_fq_NSB_grouped %>%
  select(
    Duration_group,
    GWL_abbreviation,
    Count,
    Frequency = NormalisedFrequency,
    Statistic
  )

GWL_fq_NDB_plot <- GWL_fq_NDB_grouped %>%
  select(
    Duration_group,
    GWL_abbreviation,
    Count,
    Frequency = NormalisedFrequency,
    Statistic
  )


# Combine all four statistics
GWL_plot_data <- bind_rows(
  GWL_fq_SB_grouped,
  GWL_fq_DB_grouped,
  GWL_fq_NSB_plot,
  GWL_fq_NDB_plot
) %>%
  complete(
    Statistic,
    Duration_group,
    GWL_abbreviation,
    fill = list(
      Count = 0,
      Frequency = 0
    )
  )


GWL_plot_data$Duration_group <- factor(
  GWL_plot_data$Duration_group,
  levels = c(
    "short-duration",
    "medium-duration",
    "long-duration"
  )
)

GWL_plot_data$Statistic <- factor(
  GWL_plot_data$Statistic,
  levels = c(
    "SB-CPFQ",
    "NSB-CPFQ",
    "DB-CPFQ",
    "NDB-CPFQ"
  ),
  labels = c(
    "SB-CPFQ",
    "NSB-CPFQ",
    "DB-CPFQ",
    "NDB-CPFQ"
  )
)

gwl_order <- GWL_des %>%
  arrange(GWL_number) %>%
  pull(GWL_abbreviation) %>%
  rev()

GWL_plot_data$GWL_abbreviation <- factor(
  GWL_plot_data$GWL_abbreviation,
  levels = gwl_order
)

GWL_plot_data$Statistic <- factor(
  GWL_plot_data$Statistic,
  levels = c(
    "SB-CPFQ",
    "NSB-CPFQ",
    "DB-CPFQ",
    "NDB-CPFQ"
  )
)

GWL_plot_data_SB <- GWL_plot_data %>%
  filter(Statistic %in% c("SB-CPFQ", "NSB-CPFQ")) %>%
  droplevels()

Figure6 <- ggplot(GWL_plot_data_SB,aes(x = GWL_abbreviation,y = Frequency,fill = Duration_group)) +
  geom_col(position = position_dodge(width = 0.8),width = 0.7,colour = "black", linewidth = 0.25) +
  geom_text(aes(label = ifelse(Frequency == 0,"",sprintf("%.1f", Frequency))),
            position = position_dodge(width = 0.8),
            hjust = -0.15,
            size = 2.6) +
  scale_fill_manual(values = c("short-duration"  = "#145CE0",
                               "medium-duration" = "#307970",
                               "long-duration"   = "#B78543")) +
  coord_flip() +
  facet_wrap(~Statistic,nrow = 1,scales = "free_x",
             labeller = as_labeller(c("SB-CPFQ"  = "SB-CP[FQ]",
                                      "NSB-CPFQ" = "NSB-CP[FQ]"),
                                    label_parsed)) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.17))) +
  labs(x = NULL,y = NULL,fill = NULL) +
  geom_hline(data = data.frame(Statistic = factor("NSB-CPFQ",
                                                  levels = levels(GWL_plot_data_SB$Statistic))),
             aes(yintercept = 1),
             linetype = "dashed",
             colour = "black",
             linewidth = 0.5,
             inherit.aes = FALSE)+
  theme_bw() +
  theme(legend.position = "bottom",
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        strip.text = element_text(face = "bold"),
        panel.spacing = grid::unit(1, "cm"))

print(Figure6)

ggsave(paste0(home, "Figure6.png"), Figure6,width = 24,height = 36,units = "cm",dpi = 300)

GWL_plot_data_DB <- GWL_plot_data %>%
  filter(Statistic %in% c("DB-CPFQ", "NDB-CPFQ")) %>%
  droplevels()

Figure_App1 <- ggplot(GWL_plot_data_DB,aes(x = GWL_abbreviation,y = Frequency,fill = Duration_group)) +
  geom_col(position = position_dodge(width = 0.8),width = 0.7,colour = "black", linewidth = 0.25) +
  geom_text(aes(label = ifelse(Frequency == 0,"",sprintf("%.1f", Frequency))),
            position = position_dodge(width = 0.8),
            hjust = -0.15,
            size = 2.6) +
  scale_fill_manual(values = c("short-duration"  = "#145CE0",
                               "medium-duration" = "#307970",
                               "long-duration"   = "#B78543")) +
  coord_flip() +
  facet_wrap(~Statistic,nrow = 1,scales = "free_x",
             labeller = as_labeller(c("DB-CPFQ"  = "DB-CP[FQ]",
                                      "NDB-CPFQ" = "NDB-CP[FQ]"),
                                    label_parsed)) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.17))) +
  labs(x = NULL,y = NULL,fill = NULL) +
  geom_hline(data = data.frame(Statistic = factor("NDB-CPFQ",
                                                  levels = levels(GWL_plot_data_DB$Statistic))),
             aes(yintercept = 1),
             linetype = "dashed",
             colour = "black",
             linewidth = 0.5,
             inherit.aes = FALSE)+
  theme_bw() +
  theme(legend.position = "bottom",
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        strip.text = element_text(face = "bold"),
        panel.spacing = grid::unit(1, "cm"))

print(Figure_App1)
ggsave(paste0(home, "Figure_App1.png"), Figure_App1,width = 24,height = 36, units = "cm",dpi = 300)

################################################################################
# Trend analysis: Mann-Kendall & Sen's slope + FDR #############################
################################################################################

ta_results           <- as.data.frame(matrix(ncol = 11, nrow = 0))
colnames(ta_results) <- c("DURATION", "CODE", "SYEAR", "EYEAR", "PVALUE",
                          "PVALUE_FDR", "SIG_RAW", "SIG_FDR",
                          "TAU", "SENSLOPE", "SENSLOPE_PCT")

k_min <- 30
c_min <- 0.9

for (i in 1:length(dur_cols)) {
  
  ta_dur                 <- dur_cols[i]
  ta_input_dur           <- ampts_values[, c("CODE", "YYYY", ta_dur)]
  names(ta_input_dur)[3] <- "AMPTs"
  
  for (j in 1:length(unique(ampts_values$CODE))) {
    
    ta_code           <- unique(ampts_values$CODE)[j]
    ta_input_dur_code <- subset(ta_input_dur, ta_input_dur$CODE == ta_code)
    ta_input_dur_code <- ta_input_dur_code[!is.na(ta_input_dur_code$AMPTs), ]
    
    if (nrow(ta_input_dur_code) > 0) {
      
      ta_input_dur_code <- ta_input_dur_code[order(ta_input_dur_code$YYYY), ]
      ta_first_year_avl <- min(ta_input_dur_code$YYYY)
      ta_last_year_avl  <- max(ta_input_dur_code$YYYY)
      
      if ((ta_last_year_avl - ta_first_year_avl + 1) > k_min) {
        
        ta_looping_years <- seq(ta_first_year_avl, ta_last_year_avl - (k_min - 1))
        
        for (k in 1:length(ta_looping_years)) {
          
          ta_syear                <- ta_looping_years[k]
          ta_input_dur_code_syear <- subset(ta_input_dur_code, ta_input_dur_code$YYYY >= ta_syear)
          
          ta_noyrs_expected <- ta_last_year_avl - ta_syear + 1
          ta_noyrs_real     <- length(unique(ta_input_dur_code_syear$YYYY))
          
          if (1 - (ta_noyrs_expected - ta_noyrs_real) / ta_noyrs_expected >= c_min) {
            
            MannKendall <- mk.test(ta_input_dur_code_syear$AMPTs)
            SenSlope    <- sens.slope(ta_input_dur_code_syear$AMPTs)
            
            mk_pvalue       <- as.numeric(MannKendall$p.value)
            mk_tau          <- as.numeric(MannKendall$estimates["tau"])
            mk_senslope     <- as.numeric(SenSlope$estimates)
            mk_senslope_pct <- (mk_senslope / median(ta_input_dur_code_syear$AMPTs)) * 100
            
            ta_results[nrow(ta_results) + 1, ]        <- NA
            ta_results$DURATION[nrow(ta_results)]     <- ta_dur
            ta_results$CODE[nrow(ta_results)]         <- ta_code
            ta_results$SYEAR[nrow(ta_results)]        <- ta_syear
            ta_results$EYEAR[nrow(ta_results)]        <- ta_last_year_avl
            ta_results$PVALUE[nrow(ta_results)]       <- mk_pvalue
            ta_results$TAU[nrow(ta_results)]          <- mk_tau
            ta_results$SENSLOPE[nrow(ta_results)]     <- mk_senslope
            ta_results$SENSLOPE_PCT[nrow(ta_results)] <- mk_senslope_pct
          }
        }
      }
    }
  }
  
  cat("Duration:", gsub("DUR", "", ta_dur), "min - ready\n")
}

ta_results$PVALUE       <- as.numeric(ta_results$PVALUE)
ta_results$TAU          <- as.numeric(ta_results$TAU)
ta_results$SENSLOPE     <- as.numeric(ta_results$SENSLOPE)
ta_results$SENSLOPE_PCT <- as.numeric(ta_results$SENSLOPE_PCT)
ta_results$SYEAR        <- as.numeric(ta_results$SYEAR)
ta_results$EYEAR        <- as.numeric(ta_results$EYEAR)

# Apply FDR within each family of simultaneous tests:
# one family = one DURATION + one exact analysis period (SYEAR–EYEAR)

ta_results <- ta_results %>%
  group_by(DURATION, SYEAR, EYEAR) %>%
  mutate(
    PVALUE_FDR = p.adjust(PVALUE, method = "BH"),
    
    SIG_RAW = case_when(
      is.na(PVALUE) ~ NA,
      TRUE          ~ PVALUE <= 0.05
    ),
    
    SIG_FDR = case_when(
      is.na(PVALUE_FDR) ~ NA,
      TRUE              ~ PVALUE_FDR <= 0.05
    )
  ) %>%
  ungroup() %>%
  mutate(
    TDIR_RAW = case_when(
      is.na(SENSLOPE_PCT) | is.na(SIG_RAW) ~ NA_character_,
      abs(SENSLOPE_PCT) < 0.01             ~ "no_trend",
      SENSLOPE_PCT <= -0.01 & SIG_RAW      ~ "negative_sign",
      SENSLOPE_PCT <= -0.01 & !SIG_RAW     ~ "negative",
      SENSLOPE_PCT >=  0.01 & SIG_RAW      ~ "positive_sign",
      SENSLOPE_PCT >=  0.01 & !SIG_RAW     ~ "positive",
      TRUE                                 ~ NA_character_
    ),
    
    TDIR_FDR = case_when(
      is.na(SENSLOPE_PCT) | is.na(SIG_FDR) ~ NA_character_,
      abs(SENSLOPE_PCT) < 0.01             ~ "no_trend",
      SENSLOPE_PCT <= -0.01 & SIG_FDR      ~ "negative_sign",
      SENSLOPE_PCT <= -0.01 & !SIG_FDR     ~ "negative",
      SENSLOPE_PCT >=  0.01 & SIG_FDR      ~ "positive_sign",
      SENSLOPE_PCT >=  0.01 & !SIG_FDR     ~ "positive",
      TRUE                                 ~ NA_character_
    )
  )

################################################################################
# Figure 7 #####################################################################
################################################################################

ta_results$NOYEARS <- ta_results$EYEAR - ta_results$SYEAR + 1

ta_results_longtpd <- ta_results %>%
  group_by(DURATION, CODE) %>%
  slice_max(
    order_by = NOYEARS,
    n = 1,
    with_ties = FALSE
  ) %>%
  ungroup() %>%
  filter(
    DURATION %in% c("DUR60","DUR180","DUR420","DUR1440","DUR4320","DUR10080")) %>%
  arrange(DURATION, CODE, SYEAR)


make_summary <- function(data, tdir_col) {
  
  out <- data %>%
    filter(!is.na(.data[[tdir_col]])) %>%
    group_by(
      DURATION,
      TDIR = .data[[tdir_col]]
    ) %>%
    summarise(
      COUNT = n(),
      .groups = "drop_last"
    ) %>%
    mutate(
      FREQ = COUNT / sum(COUNT) * 100
    ) %>%
    ungroup()
  
  out$DURATION <- gsub("DUR", "", out$DURATION)
  
  out$DURATION <- factor(
    out$DURATION,
    levels = c("60", "180", "420", "1440", "4320", "10080"),
    labels = c("1 h", "3 h", "7 h", "1 day", "3 days", "7 days")
  )
  
  out$TDIR <- factor(
    out$TDIR,
    levels = c(
      "negative_sign",
      "negative",
      "no_trend",
      "positive",
      "positive_sign"
    )
  )
  
  return(out)
}


raw_df <- make_summary(
  ta_results_longtpd,
  "TDIR_RAW"
)

fdr_df <- make_summary(
  ta_results_longtpd,
  "TDIR_FDR"
)

plot_trends <- function(df, title_txt) {
  
  ggplot(df,aes(x = DURATION, y = FREQ,fill = TDIR)) +
    geom_col(position = "stack",colour = "black",linewidth = 0.3) +
    scale_fill_manual(
      values = c(
        "negative_sign" = "#543005",
        "negative"      = "#D6B991",
        "no_trend"      = "grey90",
        "positive"      = "#95C4C0",
        "positive_sign" = "#003C30"
      ),
      labels = c(
        "negative_sign" = "negative significant",
        "negative"      = "negative non-significant",
        "no_trend"      = "no trend",
        "positive"      = "positive non-significant",
        "positive_sign" = "positive significant"
      ),
      drop = FALSE
    ) +
    scale_y_continuous(
      limits = c(0, 100),
      breaks = seq(0, 100, 20),
      labels = scales::number_format(accuracy = 0.1),
      expand = expansion(mult = c(0, 0.02))
    ) +
    labs(
      x = "\nDuration\n",
      y = "Fraction of stations [%]\n",
      fill = NULL,
      title = title_txt
    ) +
    guides(
      fill = guide_legend(
        nrow = 1,
        byrow = TRUE
      )
    ) +
    theme_bw() +
    theme(
      legend.position = "bottom",
      legend.title = element_text(size = 8),
      legend.text = element_text(size = 8),
      axis.title = element_text(size = 10),
      axis.text = element_text(size = 10),
      plot.title = element_text(
        size = 10,
        face = "bold"
      )
    )
}

Figure7_raw <- plot_trends(raw_df,"Uncorrected")
Figure7_fdr <- plot_trends(fdr_df,"FDR-corrected")

Figure7 <- Figure7_raw + Figure7_fdr +
  patchwork::plot_layout(
    ncol = 2,
    guides = "collect"
  ) &
  theme(
    legend.position = "bottom"
  )
plot(Figure7)

ggsave(paste0(home, "Figure7.png"), Figure7, width = 32,  height = 16,dpi = 300,units = "cm")

################################################################################
# Figure 8 #####################################################################
################################################################################

# retain only station that have ampts available till 2020

ta_results_till2020         <- subset(ta_results,ta_results$EYEAR == 2020)

ta_results_till2020_summary <- ta_results_till2020 %>%
  group_by(DURATION,SYEAR,TDIR_FDR) %>%
  summarise(COUNT = n(), .groups = "drop_last")%>%
  mutate(FREQ = COUNT / sum(COUNT)*100) %>%
  ungroup()

durs_selected_min1h         <- c("DUR60","DUR180","DUR420","DUR1440","DUR4320","DUR10080")
ta_results_till2020_summary <- subset(ta_results_till2020_summary,
                                      ta_results_till2020_summary$DURATION %in% durs_selected_min1h)

for (i in 1:length(durs_selected_min1h)){
  
  plt_dur      <- durs_selected_min1h[i]
  plt_input    <- subset(ta_results_till2020_summary,ta_results_till2020_summary$DURATION == plt_dur)
  
  plt_input <- plt_input %>%
    group_by(SYEAR) %>%
    filter(sum(COUNT) >= 10) %>%
    ungroup()
  
  plt_input_fq <- dcast(plt_input, SYEAR ~ TDIR_FDR, value.var = "FREQ") 
  plt_input_ct <- dcast(plt_input, SYEAR ~ TDIR_FDR, value.var = "COUNT") 
  plt_input_fq[is.na(plt_input_fq)] <- 0
  
  # add mising categories
  
  trend_categories <- c(
    "negative_sign",
    "negative",
    "no_trend",
    "positive",
    "positive_sign"
  )
  
  missing_categories <- setdiff(
    trend_categories,
    names(plt_input_fq)
  )
  
  plt_input_fq[missing_categories] <- 0
  plt_input_fq$Fq_sum              <- rowSums(plt_input_fq[,-1])
  plt_input_fq                     <- melt(plt_input_fq[,1:6], id.vars = "SYEAR") 
  plt_input_fq$variable            <- factor(plt_input_fq$variable,
                                              levels = c("negative_sign","negative",
                                                         "no_trend","positive","positive_sign"))
  
                                                         
  missing_categories <- setdiff(
    trend_categories,
    names(plt_input_ct)
  )                                                       
  plt_input_ct[missing_categories] <- 0                                              
  plt_input_ct$Ct_sum              <- rowSums(plt_input_ct[,-1],na.rm=TRUE)
  
  
  plt <- ggplot() +
    geom_area(data=plt_input_fq,aes(x = SYEAR, y = value, fill = variable), color="black") +
    geom_line(data=plt_input_ct,aes(x = SYEAR, y = Ct_sum/10), color="black", linetype = "dashed",size=0.5) +  # scaled
    geom_point(data=plt_input_ct,aes(x = SYEAR, y = Ct_sum/10), fill="white",alpha=0.5,color="black",pch=21, size=4) +  # scaled
    scale_y_continuous(
      name = "Fraction of stations [%]",
      sec.axis = sec_axis(~.*10, name = "Number of stations")   # invert scaling
    ) +
    scale_x_continuous(limits = c(1900, 1992), expand = c(0, 0)) +
    theme_bw() +
    scale_fill_manual(
      values = c(
        "negative_sign" = "#543005",
        "negative"      = "#D6B991",
        "no_trend"      = "grey90",
        "positive"      = "#95C4C0",
        "positive_sign" = "#003C30"
      ),
      labels = c(
        "negative_sign" = "negative significant",
        "negative"      = "negative non-significant",
        "no_trend"      = "no trend",
        "positive"      = "positive non-significant",
        "positive_sign" = "positive significant"
      ),
      drop = FALSE
    ) +
    labs(x="\nStart Year", fill="") +
    ggtitle(paste("Duration:", dur_labels[plt_dur])) +
    theme(plot.title = element_text(face = "bold", size = 10)) +
    coord_cartesian(clip = "off")
  #plot(plt)
  assign(paste0("Figure8_",i), plt)

  
}

# combine part 1-9 to final Figure 7

TDIR_fill <- scale_fill_manual(
  name   = NULL,
  values = c(
    "negative_sign" = "#543005",
    "negative"      = "#D6B991",
    "no_trend"      = "grey90",
    "positive"      = "#95C4C0",
    "positive_sign" = "#003C30"
  ),
  labels = c(
    "negative_sign" = "negative significant",
    "negative"      = "negative non-significant",
    "no_trend"      = "no trend",
    "positive"      = "positive non-significant",
    "positive_sign" = "positive significant"
  ),
  drop   = FALSE
)

Figure8 <- (
   Figure8_1 + Figure8_2 + Figure8_3 +
    Figure8_4 + Figure8_5 + Figure8_6) +
  plot_layout(ncol = 3, guides = "collect") &
  TDIR_fill &
  theme(
    legend.position   = "bottom",
    legend.box        = "vertical",
    legend.spacing.y  = unit(-1, "mm"),
    legend.spacing.x  = unit(3, "mm"),
    legend.key.width  = unit(7, "mm"),
    plot.margin       = margin(3, 3, 3, 3, unit = "mm")
  )
#plot(Figure8)
ggsave(paste0(home, "Figure8.png"),Figure8,
  width = 36, height = 24, dpi = 300, units = "cm"
)

################################################################################
# Figure 9 #####################################################################
################################################################################

# process aic 

aic <- melt(aic, id="CODE", variable.name = "DURATION", value.name = "AIC")

# define legend properties for plotting (breaks,labels & colors) - median

BreaksDefined_SensSlope      <- c(-Inf,-2.0,-1.5, -1.0,-0.5,-0.1,-0.01,0.01,0.1, 0.5, 1.0, 1.5, 2.0,  Inf)
LabelsDefined_SensSlope      <- c(paste("≤", format(BreaksDefined_SensSlope[2],nsmall =1)),
                                  format(BreaksDefined_SensSlope[3:(length(BreaksDefined_SensSlope) - 1)], nsmall = 1),
                                  paste(">", format(BreaksDefined_SensSlope[length(BreaksDefined_SensSlope-1)-1],nsmall=1)))
ColorsDefined_SensSlope      <- c("#543005","#754C19","#96692E","#B78543","#D6B991","grey90",
                                  "grey90", "#95C4C0","#6AB0AA","#489790","#307970","#185A50","#003C30")

ta_results$SENSLOPE_PCT_cut  <- cut(ta_results$SENSLOPE_PCT,
                                    breaks = BreaksDefined_SensSlope, 
                                    include.lowest = TRUE,
                                    right = FALSE, 
                                    labels = LabelsDefined_SensSlope)                          

# define starting years & investigation period

plt_SYEAR_all <- c(1951,1971,1991)
plt_EYEAR     <- 2020

for (i in 1:length(plt_SYEAR_all)){
  
  plt_SYEAR        <- plt_SYEAR_all[i]
  plt_input_syear  <- subset(ta_results,ta_results$SYEAR==plt_SYEAR & ta_results$EYEAR==plt_EYEAR) 
  
  for (j in 1:length(durs_selected_min1h)){
    
    plt_dur             <- durs_selected_min1h[j]
    plt_input_syear_dur <- subset(plt_input_syear,plt_input_syear$DURATION==plt_dur)
    plt_input_syear_dur <- left_join(plt_input_syear_dur, aic, by=c("CODE","DURATION"))
    plt_input_syear_dur <- left_join(plt_input_syear_dur, metadata, by="CODE")
    plt_input_syear_dur <- st_as_sf(plt_input_syear_dur, coords = c("LON","LAT"), crs = WGS84)
    plt_input_syear_dur <- st_transform(plt_input_syear_dur, crs = LAEA)
    
    plt_input_notsign_notjump <- subset(plt_input_syear_dur,
                                        plt_input_syear_dur$PVALUE_FDR>0.05 & (plt_input_syear_dur$AIC!="Jump"| is.na(plt_input_syear_dur$AIC)))
    plt_input_notsign_jump    <- subset(plt_input_syear_dur,
                                        plt_input_syear_dur$PVALUE_FDR>0.05 & plt_input_syear_dur$AIC=="Jump")
    plt_input_sign_notjump    <- subset(plt_input_syear_dur,
                                        plt_input_syear_dur$PVALUE_FDR<=0.05 & (plt_input_syear_dur$AIC!="Jump"| is.na(plt_input_syear_dur$AIC)))
    plt_input_sign_jump       <- subset(plt_input_syear_dur,
                                        plt_input_syear_dur$PVALUE_FDR<=0.05 & plt_input_syear_dur$AIC=="Jump")
    
    plt <- ggplot()+
      geom_sf(data = GER_regions, fill = NA, colour = "black", linewidth = 0.35) +
      # layer: not-sign & jump
      geom_sf(data = plt_input_notsign_jump, aes(fill = SENSLOPE_PCT_cut), colour="white", size = 2.5, pch = 24, show.legend = TRUE)+
      # layer: not-sign & not-jump
      geom_sf(data = plt_input_notsign_notjump, aes(fill = SENSLOPE_PCT_cut), colour="white", size = 3, pch = 21, show.legend = FALSE)+
      # layer: sign & jump
      geom_sf(data = plt_input_sign_jump, aes(fill = SENSLOPE_PCT_cut), colour="black", size = 2.5, pch = 24, show.legend = FALSE)+
      # layer: sign & not-jump
      geom_sf(data = plt_input_sign_notjump, aes(fill = SENSLOPE_PCT_cut), colour="black", size = 3, pch = 21, show.legend = FALSE)+
      # further settings
      scale_fill_manual(values = ColorsDefined_SensSlope, 
                        drop = FALSE, na.value = "red", na.translate = TRUE,
                        guide = guide_legend(reverse = TRUE)) +
      theme_bw() +
      theme(legend.position = 'bottom') +
      theme(axis.title = element_blank(), 
            axis.text = element_text(size = 8),
            legend.text = element_text(size = 10),
            legend.title = element_text(size = 10),
            plot.margin = unit(c(0.25, 0, 0.25, 0), "cm"),
            plot.title = element_text(color = "black", size = 10, face = "bold"),
            plot.subtitle = element_text(color = "black", size = 10, face = "bold")) +
      labs(fill=NULL) +
      ggtitle(paste("Duration:", dur_labels[plt_dur])) +
      geom_sf_label(data=GER_regions,aes(label=ID),size=3, alpha=0.8)
    #plot(plt)
    assign(paste0("Figure9_",j), plt)
    
  }
  
  Sen_levels <- LabelsDefined_SensSlope  # same order as your cut labels
  SenSlope_fill <- scale_fill_manual(
    name   = NULL,
    values = ColorsDefined_SensSlope,
    limits = LabelsDefined_SensSlope,
    breaks = LabelsDefined_SensSlope,
    drop   = FALSE,
    na.translate = TRUE,
    na.value = "red",
    guide = guide_legend(
      nrow = 1, byrow = TRUE, reverse = FALSE,
      override.aes = list(shape = 21, size = 4, colour = "black")
    )
  )
  
  SenSlope_shape <- scale_shape_manual(
    name   = NULL,
    values = c("no jump" = 21, "jump" = 24),  # 21 = circle, 24 = triangle
    breaks = c("no jump","jump"),
    drop   = FALSE,
    na.translate = TRUE,
    na.value = "red",
    guide = guide_legend(
      nrow = 1, byrow = TRUE, reverse = FALSE,
      override.aes = list(size = 4, colour = "black")
    )
  ) 
  
  Figure9 <- (
    Figure9_1 + Figure9_2 + Figure9_3 +
      Figure9_4 + Figure9_5 + Figure9_6 ) +
    plot_layout(ncol = 3, guides = "collect") &
    SenSlope_fill &
    theme(
      legend.position  = "bottom",
      legend.direction = "horizontal",  # <-- horizontal legend
      legend.box       = "horizontal",  # <-- fix typo ("vertikal")
      legend.spacing.y = unit(0, "mm"),
      legend.spacing.x = unit(3, "mm"),
      legend.key.width = unit(7, "mm"),
      plot.margin      = margin(3, 3, 3, 3, unit = "mm")
    )
  
  if(plt_SYEAR==1951){
    ggsave(Figure9, 
           file=paste0(home, "Figure9.png"),
           width = 32 , height = 24,  dpi = 300, units = "cm")
  } else if(plt_SYEAR==1971){
    ggsave(Figure9, 
           file=paste0(home, "Figure_App2.png"),
           width = 32 , height = 24,  dpi = 300, units = "cm")
  }else if(plt_SYEAR==1991){
    ggsave(Figure9, 
           file=paste0(home, "Figure_App3.png"),
           width = 32 , height = 24,  dpi = 300, units = "cm")
  }

}