################################################################################ # Process annual maximum precipitation totals (ampts) and create figures ####### ################################################################################ # load required R packages library(cowplot) library(dplyr) library(ggplot2) library(patchwork) library(RColorBrewer) library(reshape2) library(sf) library(trend) # load required input datasets home = "/INSERT/YOUR/PATH/HERE" ampts = read.csv(paste0(home,"AMPTs.csv")) aic = read.csv(paste0(home,"AIC.csv")) metadata = read.csv(paste0(home,"Metadata.csv")) # process AMPTs: # 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), value = TRUE) ampts[dur_cols] <- lapply(ampts[dur_cols], \(x) replace(x, x==-99.9, NA_real_)) for (code in unique(ampts$CODE)){ idx <- ampts$CODE == code for (dur in dur_cols){ n_valid <- sum(!is.na(ampts[idx,dur])) if(n_valid < 30){ ampts[idx,dur] <- NA_real_ cat( "Station:", code, "| Duration:", dur, "| Valid yrs:", n_valid, "\n" ) } } } ################################################################################ # Figure 1 ##################################################################### ################################################################################ dur_keep <- c("DUR5", "DUR420", "DUR1440") dur_labels <- c("DUR5" = "5-minute", "DUR420" = "7-hour", "DUR1440" = "1-day" ) stations_per_year <- ampts %>% group_by(YYYY) %>% summarise( across( all_of(dur_keep), ~ sum(!is.na(.) & . != -99.9) ), .groups = "drop" ) stations_per_year <- melt( stations_per_year, id = "YYYY", variable.name = "DURATION", value.name = "N_STATIONS" ) Figure1 <- ggplot(stations_per_year,aes(x = YYYY, y = N_STATIONS, color = DURATION)) + geom_line()+ scale_color_manual( values = c("orange", "red", "darkblue"), labels = c( "DUR5" = "5-minute", "DUR420" = "7-hour", "DUR1440" = "24-hour")) + labs(title = "", x = "\nYears", y = "Number of stations\n", color="Temporal resolution:") + scale_x_continuous(limits = c(1895, 2025), breaks= seq(1900, 2020, by=20), expand = c(0, 0))+ ylim(0,1500)+ theme_bw()+ theme(legend.position = "bottom", legend.title = element_text(size=8), legend.text = element_text(size=8), axis.title.x = element_text(size=10), axis.title.y = element_text(size=10), axis.text.x = element_text(size=10), axis.text.y = element_text(size=10)) #plot(Figure1) ggsave(Figure1, file=paste0(home, "Figure1.png"), width = 16, height = 16, dpi = 300, units = "cm") ################################################################################ # Figure 2 ##################################################################### ################################################################################ years_per_station <- ampts %>% group_by(CODE) %>% summarise( across( all_of(dur_keep), \(x) sum(!is.na(x) & x != -99.9), .names = "N_YEARS_{.col}" ), .groups = "drop" ) # convert "year_per_station" to spatial object & prepare for plotting # NOTE: physical region shp-files are not included due to the licensing restrictions 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") years_per_station <- left_join(years_per_station, metadata, by = "CODE") years_per_station <- st_as_sf(years_per_station, coords = c("LON","LAT"), crs = WGS84) years_per_station <- st_transform(years_per_station, crs = LAEA) # define legend properties for plotting (breaks,labels & colors) - number of years BreaksDefined_YrsNo <- c(30,40,50,60,70,80,90,100,Inf) LabelsDefined_YrsNo <- c(format(BreaksDefined_YrsNo[1:(length(BreaksDefined_YrsNo) - 1)], nsmall = 1)) ColorsDefined_YrsNo <- brewer.pal(length(LabelsDefined_YrsNo),"YlOrBr") for (i in 1:length(dur_keep)){ plt_col <- paste0("N_YEARS_",dur_keep[i]) plt_input <- years_per_station plt_input <- subset(plt_input,plt_input[[plt_col]] !=0) plt_input$cat <- cut(plt_input[[plt_col]], breaks = BreaksDefined_YrsNo, include.lowest = TRUE, right = FALSE, labels = LabelsDefined_YrsNo) plt <- ggplot() + geom_sf(data = plt_input, aes(fill = cat),size=1.7, pch = 21, show.legend = TRUE) + scale_fill_manual(values = ColorsDefined_YrsNo, labels = c("100.0" = "≥ 100.0"), 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, hjust = 0), legend.title = element_text(size = 10), #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("Temporal resolution: ", dur_labels[i])) #plot(plt) assign(paste0("Figure2_",i), plt) } # define legend properties for plotting (breaks,labels & colors) - altitude BreaksDefined_Alt <- c(0,100,250,500,750,1000,1250,1500,1750,Inf) LabelsDefined_Alt <- trimws(format(BreaksDefined_Alt[1:(length(BreaksDefined_Alt) - 1)], nsmall = 1)) ColorsDefined_Alt <- c("#00A600", "#3EBB00", "#8BD000", "#E6E600", "#E8C32E", "#EBB25E", "#EDB48E", "#D8B29A", "#8B5A2B") years_per_station$ALT_cut <- cut(years_per_station$ALT, breaks = BreaksDefined_Alt, include.lowest = TRUE, right = FALSE, labels = LabelsDefined_Alt) Figure2_4 <- ggplot() + geom_sf(data = years_per_station, aes(fill = ALT_cut),size=1.7, pch = 21, show.legend = TRUE) + scale_fill_manual(values = ColorsDefined_Alt, labels = c("1750.0" = "≥ 1750.0"), drop = FALSE, na.value = NA, 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.title = element_text(size = 10), legend.text = element_text(size = 10, hjust = 0), #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("Altitude above sea level [m]") #plot(Figure2_4) # combine part 1-4 to final Figure 2 Figure2 <- plot_grid(Figure2_1,Figure2_2,Figure2_3, Figure2_4, ncol=2, nrow = 2)+ theme(plot.background = element_rect(fill="white", colour = "white")) #plot(Figure2) ggsave(Figure2, file=paste0(home, "Figure2.png"), width = 30 ,height = 24, dpi = 300, units = "cm") ################################################################################ # Figure 3 ##################################################################### ################################################################################ ampts_long <- melt(ampts, id=c("CODE","YYYY"), variable.name = "DURATION", value.name = "AMPTs") ampts_long$DURATION <- factor(ampts_long$DURATION, levels = levels(ampts_long$DURATION), labels = sub("^DUR","", levels(ampts_long$DURATION))) Figure3 <- ggplot(ampts_long, aes(x = DURATION, y = AMPTs) )+ geom_violin(trim=FALSE, fill="grey80",color="black",linewidth=0.2) + geom_boxplot(width=0.2, color="red", linewidth=0.5,outlier.shape = 1)+ stat_summary(fun = median, aes(group = 1), geom = "line", color = "darkgreen", size = 0.5) + stat_summary(fun = median, aes(group = 1), geom = "point", color = "darkgreen", size = 2) + labs(title = "", x = "\nDurations [min]", y = "AMPTs [mm]\n") + scale_y_continuous(labels = scales::number_format(accuracy = 0.1)) + coord_cartesian(ylim = c(1, 500))+ theme_bw()+ theme(legend.position = "bottom", legend.title = element_text(size=10), legend.text = element_text(size=10), axis.title.x = element_text(size=10), axis.title.y = element_text(size=10), axis.text.x = element_text(size=10), axis.text.y = element_text(size=10)) #plot(Figure3) ggsave(Figure3, file=paste0(home, "Figure3.png"), width = 16 ,height = 16, dpi = 300, units = "cm") ################################################################################ # Figure 4 ##################################################################### ################################################################################ stats_per_station_med <- ampts %>% group_by(CODE) %>% summarise( across(all_of(dur_cols), \(x) median(x, na.rm =TRUE), .names = "MEDIAN_{.col}"), .groups = "drop" ) # convert "stats_per_station_med" to spatial object & prepare for plotting stats_per_station_med <- left_join(stats_per_station_med, metadata, by = "CODE") stats_per_station_med <- st_as_sf(stats_per_station_med, coords = c("LON","LAT"), crs = WGS84) stats_per_station_med <- st_transform(stats_per_station_med, crs = LAEA) # define legend properties for plotting (breaks,labels & colors) - median BreaksDefined_med_DUR5 <- c(seq(1,10,by=1),Inf) BreaksDefined_med_DUR30 <- c(seq(2.7,27,by=2.7),Inf) BreaksDefined_med_DUR60 <- c(seq(3.2,32,by=3.2),Inf) BreaksDefined_med_DUR180 <- c(seq(4,40,by=4),Inf) BreaksDefined_med_DUR420 <- c(seq(5,50,by=5),Inf) BreaksDefined_med_DUR720 <- c(seq(6,60,by=6),Inf) BreaksDefined_med_DUR1440 <- c(seq(7,70,by=7),Inf) BreaksDefined_med_DUR4320 <- c(seq(9,90,by=9),Inf) BreaksDefined_med_DUR10080 <- c(seq(12,120,by=12),Inf) LabelsDefined_med_DUR5 <- c(format(BreaksDefined_med_DUR5[1:(length(BreaksDefined_med_DUR5) - 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_DUR720 <- c(format(BreaksDefined_med_DUR720[1:(length(BreaksDefined_med_DUR720) - 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)) ColorsDefined_med <- c("#FFF9C4","#D8ECBC","#B1E0B4","#91CBB0","#77AEB0","#5D91B0","#3876C8","#145CE0","#124BCF","#334393","#543C58") for (i in 1:length(dur_cols)){ plt_col <- paste0("MEDIAN_",dur_cols[i]) plt_input <- stats_per_station_med plt_input <- subset(plt_input,plt_input[[plt_col]] !=0) plt_breaks <- get(paste0("BreaksDefined_med_",dur_cols[i])) plt_labels <- get(paste0("LabelsDefined_med_",dur_cols[i])) plt_input$cut <- cut(plt_input[[plt_col]], breaks = plt_breaks, include.lowest = TRUE, right = FALSE, labels = plt_labels) plt <- ggplot() + geom_sf(data = plt_input, aes(fill = cut),size=2.5, pch = 21, show.legend = TRUE) + scale_fill_manual(values = ColorsDefined_med, 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: ",gsub("DUR","",dur_cols[i])," min")) #plot(plt) assign(paste0("Figure4_",i), plt) } # combine part 1-9 to final Figure 4 Figure4_combined <- 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_combined) ggsave(Figure4_combined, file=paste0(home, "Figure4.png"), width = 36 , height = 36, dpi = 300, units = "cm") ################################################################################ # Figure 5 ##################################################################### ################################################################################ stats_per_station_cv <- ampts %>% group_by(CODE) %>% summarise( across(all_of(dur_cols), \(x) sd(x, na.rm =TRUE)/mean(x, na.rm =TRUE)*100, .names = "CV_{.col}"), .groups = "drop" ) # convert "stats_per_station_cv" to spatial object & prepare for plotting stats_per_station_cv <- left_join(stats_per_station_cv, metadata, by = "CODE") stats_per_station_cv <- st_as_sf(stats_per_station_cv, coords = c("LON","LAT"), crs = WGS84) stats_per_station_cv <- st_transform(stats_per_station_cv, 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)) ColorsDefined_cv <- brewer.pal(length(LabelsDefined_cv),"YlOrBr") for (i in 1:length(dur_cols)){ plt_col <- paste0("CV_",dur_cols[i]) plt_input <- stats_per_station_cv 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 = plt_input, aes(fill = cut),size=2.5, pch = 21, show.legend = TRUE) + scale_fill_manual(values = ColorsDefined_cv, 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: ",gsub("DUR","",dur_cols[i])," min")) #plot(plt) assign(paste0("Figure5_",i), plt) } # combine part 1-9 to final Figure 5 Figure5_combined <- 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_combined) ggsave(Figure5_combined, file=paste0(home, "Figure5.png"), width = 36 , height = 36, dpi = 300, units = "cm") ################################################################################ # Trend analysis: Mann-Kendall & Sen's Slope ################################### ################################################################################ ta_results <- as.data.frame(matrix(ncol=8,nrow=0)) colnames(ta_results) <- c("DURATION","CODE","SYEAR","EYEAR","PVALUE", "TAU","SENSLOPE","SENSLOPE_PCT") k_min <- 30 # predefined minimum number of yrs to perform trend analysis c_min <- 0.9 # predefined minimum ratio of yrs to perform trend analysis # here ampts for minimum 90% of yrs has to be available for (i in 1:length(dur_cols)){ ta_dur <- dur_cols[i] ta_input_dur <- ampts[,c("CODE","YYYY",ta_dur)] names(ta_input_dur)[3] <- "AMPTs" for(j in 1:length(unique(ampts$CODE))){ ta_code <- unique(ampts$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 <- round(MannKendall$p.value,3) mk_tau <- round(MannKendall$estimates["tau"],3) mk_senslope <- round(SenSlope$estimates,5) mk_senslope_pct <- round((mk_senslope/median(ta_input_dur_code_syear$AMPTs))*100,5) 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") } # group results into classes ta_results <- ta_results %>% mutate( TDIR = case_when( abs(SENSLOPE_PCT) < 0.01 ~ "no_trend", SENSLOPE_PCT <= -0.01 & PVALUE <= 0.05 ~ "negative_sign", SENSLOPE_PCT <= -0.01 & PVALUE > 0.05 ~ "negative", SENSLOPE_PCT >= 0.01 & PVALUE <= 0.05 ~ "positive_sign", SENSLOPE_PCT >= 0.01 & PVALUE > 0.05 ~ "positive", TRUE ~ NA_character_ ) ) ################################################################################ # Figure 6 ##################################################################### ################################################################################ # retain only the longest available time series per station and duration 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() %>% arrange(DURATION, CODE, SYEAR) ta_results_longtpd_summary <- ta_results_longtpd %>% group_by(DURATION,TDIR) %>% summarise(COUNT = n(), .groups = "drop_last")%>% mutate(FREQ = COUNT / sum(COUNT)*100) %>% ungroup() ta_results_longtpd_summary$DURATION <-gsub("DUR","",ta_results_longtpd_summary$DURATION) ta_results_longtpd_summary$DURATION <- factor(ta_results_longtpd_summary$DURATION, levels = c("5","30","60","180","420", "720","1440","4320","10080")) ta_results_longtpd_summary$TDIR <- factor(ta_results_longtpd_summary$TDIR, levels = c("negative_sign","negative", "no_trend","positive","positive_sign")) Figure6 <- ggplot(ta_results_longtpd_summary,aes(x = DURATION, y = FREQ, fill = TDIR)) + geom_rect(xmin = 0.5, xmax = 2.5, ymin = -Inf, ymax = Inf,fill = NA, colour = "red", linewidth = 1, linetype= "dashed") + geom_col(position = "stack", colour="black", linewidth=0.3) + theme_bw() + scale_fill_manual( values = c("negative_sign"="#543005", "negative"="#D6B991", "positive"="#95C4C0", "positive_sign"="#003C30"), labels = c("negative_sign"="negative significant", "negative"="negative non-significant", "positive"="positive non-significant", "positive_sign"="positive significant") ) + labs(x = "\nDurations [min]", y = "Fraction of stations [%]\n", fill = "Trend direction and significance:") + scale_y_continuous(labels = scales::number_format(accuracy = 0.1)) + guides(fill = guide_legend(nrow = 2)) + 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(Figure6) ggsave(Figure6, file=paste0(home, "Figure6.png"), width = 16 ,height = 16, dpi = 300, units = "cm") ################################################################################ # Figure 7 ##################################################################### ################################################################################ # 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) %>% summarise(COUNT = n(), .groups = "drop_last")%>% mutate(FREQ = COUNT / sum(COUNT)*100) %>% ungroup() for (i in 1:length(dur_cols)){ plt_dur <- dur_cols[i] plt_input <- subset(ta_results_till2020_summary,ta_results_till2020_summary$DURATION == plt_dur) plt_input_fq <- dcast(plt_input, SYEAR ~ TDIR, value.var = "FREQ") plt_input_ct <- dcast(plt_input, SYEAR ~ TDIR, value.var = "COUNT") plt_input_fq[is.na(plt_input_fq)] <- 0 plt_input_fq$Fq_sum <- rowSums(plt_input_fq[,2:6]) 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")) 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", "positive"="#95C4C0", "positive_sign"="#003C30"), labels = c( negative_sign = "negative significant", negative = "negative non-significant", positive = "positive non-significant", positive_sign = "positive significant" )) + labs(x="\nStart Year", fill="") + ggtitle(paste("Duration:", gsub("DUR","",plt_dur), "min"))+ theme(plot.title = element_text(face = "bold", size = 10)) + coord_cartesian(clip = "off") #plot(plt) assign(paste0("Figure7_",i), plt) } # combine part 1-9 to final Figure 7 TDIR_fill <- scale_fill_manual( name = "Trend direction and significance:", values = c( negative_sign="#754C19", negative="#D6B991", positive="#95C4C0", positive_sign="#185A50" ), labels = c( negative_sign = "negative significant", negative = "negative non-significant", positive = "positive non-significant", positive_sign = "positive significant" ), drop = FALSE ) # add dashed red "caution" frames to duration 5 & 30 minutes caution_frame <- geom_rect( aes(xmin = -Inf, xmax = Inf,ymin = -Inf, ymax = Inf ), inherit.aes = FALSE,fill = NA,colour = "red3", linewidth = 2, linetype = "dashed") Figure7_1 <- Figure7_1 + caution_frame Figure7_2 <- Figure7_2 + caution_frame Figure7 <- ( Figure7_1 + Figure7_2 + Figure7_3 + Figure7_4 + Figure7_5 + Figure7_6 + Figure7_7 + Figure7_8 + Figure7_9) + 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(Figure7) ggsave( Figure7, file = paste0(home, "Figure7.png"), width = 36, height = 36, dpi = 300, units = "cm" ) ################################################################################ # Figures 8-10 ################################################################# ################################################################################ # 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","grey", "grey", "#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(dur_cols)){ plt_dur <- dur_cols [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>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>0.05 & plt_input_syear_dur$AIC=="Jump") plt_input_sign_notjump <- subset(plt_input_syear_dur, plt_input_syear_dur$PVALUE<=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<=0.05 & plt_input_syear_dur$AIC=="Jump") plt <- ggplot()+ # 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="Sen's slope [%]") + ggtitle(paste("Duration:",gsub("DUR_","",plt_dur),"min")) #plot(plt) assign(paste0("Figure8_",j), plt) } # add dashed red frames (panel border) to DUR_5 and DUR_30 maps caution_border <- theme(panel.border = element_rect(colour = "red3", fill = NA, linewidth = 2, linetype = "dashed")) Figure8_1 <- Figure8_1 + caution_border Figure8_2 <- Figure8_2 + caution_border Sen_levels <- LabelsDefined_SensSlope # same order as your cut labels SenSlope_fill <- scale_fill_manual( name = "Sen's slope [%]:", 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 = "Akaike Information Criterion:", 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") ) ) Figure8 <- ( Figure8_1 + Figure8_2 + Figure8_3 + Figure8_4 + Figure8_5 + Figure8_6 + Figure8_7 + Figure8_8 + Figure8_9 ) + 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(Figure8, file=paste0(home, "Figure8.png"), width = 32 , height = 38, dpi = 300, units = "cm") } else if(plt_SYEAR==1971){ ggsave(Figure8, file=paste0(home, "Figure9.png"), width = 32 , height = 38, dpi = 300, units = "cm") }else if(plt_SYEAR==1991){ ggsave(Figure8, file=paste0(home, "Figure10.png"), width = 32 , height = 38, dpi = 300, units = "cm") } } ################################################################################ # Figure 11 #################################################################### ################################################################################ # compute loess l_min <- 70 # predefined minimum number of yrs to compute LOESS loess_swnd <- 42 # predefined number of years for LOESS-smoothing window loess_results <- as.data.frame(matrix(ncol=5,nrow=0)) colnames(loess_results) <- c("CODE","YYYY","DURATION","AMPTs","LOESS") for (i in 1:length(dur_cols)){ loess_dur <- dur_cols[i] loess_input_dur <- ampts[,c("CODE","YYYY",loess_dur)] colnames(loess_input_dur)[3] <- "AMPTs" for(j in 1:length(unique(loess_input_dur$CODE))){ loess_code <- unique(loess_input_dur$CODE)[j] loess_input_dur_code <- subset(loess_input_dur,loess_input_dur$CODE==loess_code) loess_input_dur_code <- loess_input_dur_code[!is.na(loess_input_dur_code$AMPTs),] if(nrow(loess_input_dur_code)>0){ loess_first_year_avl <- min(loess_input_dur_code$YYYY) loess_last_year_avl <- max(loess_input_dur_code$YYYY) loess_no_years_expected <- loess_last_year_avl - loess_first_year_avl +1 loess_no_years_real <- length(unique(loess_input_dur_code$YYYY)) loess_no_years_missing <- loess_no_years_expected - loess_no_years_real if(loess_no_years_expected>=70 & loess_no_years_missing==0){ loess_span <- loess_swnd/loess_no_years_expected loess_fit <- loess(loess_input_dur_code$AMPTs ~ loess_input_dur_code$YYYY, span=loess_span) loess_input_dur_code$LOESS <- predict(loess_fit) loess_input_dur_code$DURATION <- loess_dur loess_input_dur_code <- loess_input_dur_code[,c("CODE","YYYY","DURATION","AMPTs","LOESS")] loess_results <- rbind(loess_results,loess_input_dur_code) } } } cat("Duration:", gsub("DUR","",loess_dur), "min - ready\n") } # add metadata and plot loess_results <- left_join(loess_results,metadata,by=c("CODE")) loess_results$ALT_cut <- factor(cut(as.numeric(loess_results$ALT), breaks = BreaksDefined_Alt, labels = LabelsDefined_Alt, include.lowest = TRUE, right = FALSE), levels = LabelsDefined_Alt) for (i in 1:length(dur_cols)){ loess_dur <- dur_cols[i] loess_input_dur <- subset(loess_results,loess_results$DURATION==loess_dur) median_per_year <- loess_input_dur %>% filter(YYYY >= 1951, YYYY <= 2020) %>% group_by(YYYY) %>% summarise(MEDIAN = median(LOESS, na.rm = TRUE)) LimitDown_def <- 0 if(loess_dur == "DUR5"){ LimitUp_def <- 16 } else if (loess_dur == "DUR30"){ LimitUp_def <- 32 } else if (loess_dur == "DUR60"){ LimitUp_def <- 40 } else if (loess_dur == "DUR180"){ LimitUp_def <- 48 } else if (loess_dur == "DUR420"){ LimitUp_def <- 56 } else if (loess_dur == "DUR720"){ LimitUp_def <- 64 } else if (loess_dur == "DUR1440"){ LimitUp_def <- 160 } else if (loess_dur == "DUR4320"){ LimitUp_def <- 240 } else if (loess_dur == "DUR10080"){ LimitUp_def <- 320 } Break_def <- LimitUp_def/4 plt <- ggplot()+ geom_line(data = loess_input_dur, aes(x=YYYY,y=LOESS,group=CODE, colour=ALT_cut))+ geom_line(data = median_per_year, aes(x = YYYY, y =MEDIAN), colour = "black", size = 0.8) + theme_bw()+ labs(x ="\nYears", y="LOESS-smoothed AMPTs [mm] \n", colour = "Altitude above sea level [m]")+ scale_x_continuous(limits = c(1901, 2020), expand = c(0, 0))+ scale_y_continuous( limits = c(LimitDown_def, LimitUp_def), breaks = seq(LimitDown_def, LimitUp_def, by = Break_def), labels = scales::number_format(accuracy = 0.1)) + scale_colour_manual(values = ColorsDefined_Alt, labels = c("1750.0" = "≥ 1750.0"), drop=FALSE, na.value="red", na.translate=TRUE)+ ggtitle(paste("Duration:",gsub("DUR","",loess_dur),"min"))+ theme(plot.title = element_text(face = "bold", size = 10), legend.position = "bottom") # <- IMPORTANT: hide subplot legends #plot(plt) assign(paste0("Figure11_",i), plt) } Alt_colour <- scale_colour_manual( name = "Altitude above sea level [m]", values = ColorsDefined_Alt, limits = LabelsDefined_Alt, breaks = LabelsDefined_Alt, drop = FALSE, na.translate = TRUE, na.value = "red", guide = guide_legend( nrow = 2, byrow = TRUE, reverse = FALSE, override.aes = list(shape = 21, size = 4) ) ) Figure11 <- ( Figure11_1 + Figure11_2 + Figure11_3 + Figure11_4 + Figure11_5 + Figure11_6 + Figure11_7 + Figure11_8 + Figure11_9) + plot_layout(ncol = 3, guides = "collect") & Alt_colour & theme(legend.position = "bottom") #plot(Figure11) ggsave(Figure11, file=paste0(home, "Figure11.png"), width = 36 , height = 36, dpi = 500, units = "cm")