library(tidyverse) library(plyr) library(dplyr) library(lubridate) library(gridExtra) library(ggpubr) library(gridExtra) dev.off() rm(list=ls(all=T)) setwd() # set working directory data <- read.table("20220513_Storfjarden_interpol.txt",sep=" ",dec=".",header=TRUE, stringsAsFactors = TRUE) # Data preparation RefP4 <- read.table("20220516_Reference Period 4 - 1991-2020.txt",sep=" ",dec=".",header=TRUE, stringsAsFactors = TRUE) RefP4_clim_sur <- RefP4%>% filter(factored_depth == "surface") RefP4_clim_bot <- RefP4 %>% filter (factored_depth == "bottom") RefP1 <- read.table("20220516_Reference Period 1 - 1931-1960.txt",sep=" ",dec=".",header=TRUE, stringsAsFactors = TRUE) RefP1_clim_sur <- RefP1 %>% filter (factored_depth == "surface") RefP1_clim_bot <- RefP1 %>% filter (factored_depth == "bottom") filtered_data<-data%>% mutate(Date = as.Date(Date, format = "%Y-%m-%d"))%>% mutate(year = lubridate::year(Date), month = lubridate::month(Date), day = lubridate::day(Date), doy = lubridate::yday(Date)) RefP1_surface <- filtered_data%>% filter(Date >= "1931-01-01")%>% filter(Date <= "1960-12-31")%>% filter(factored_depth == "surface")%>% select(doy, MeanTemp) colnames(RefP1_surface)[2] <- "Temp_P1" RefP1_bottom <- filtered_data%>% filter(Date >= "1931-01-01")%>% filter(Date <= "1960-12-31")%>% filter(factored_depth == "bottom")%>% select(doy, MeanTemp) colnames(RefP1_bottom)[2] <- "Temp_P1" RefP4_surface <- filtered_data%>% filter(Date >= "1991-01-01")%>% filter(Date <= "2020-12-31")%>% filter(factored_depth == "surface")%>% select(doy, MeanTemp) colnames(RefP4_surface)[2] <- "Temp_P4" RefP4_bottom <- filtered_data%>% filter(Date >= "1991-01-01")%>% filter(Date <= "2020-12-31")%>% filter(factored_depth == "bottom")%>% select(doy, MeanTemp) colnames(RefP4_bottom)[2] <- "Temp_P4" merged_surface <- RefP4_surface merged_surface$Temp_P1 <- NA merged_surface$Temp_P1[seq_along(RefP1_surface$Temp_P1)] <- RefP1_surface$Temp_P1 merged_bottom <- RefP4_bottom merged_bottom$Temp_P1 <- NA merged_bottom$Temp_P1[seq_along(RefP1_bottom$Temp_P1)] <- RefP1_bottom$Temp_P1 # Applying KS test to temperature distributions of every day of the year - surface layer KS_results_surface <- merged_surface %>% group_by(doy)%>% mutate( Dmax= unlist(ks.test(Temp_P1,Temp_P4)[1]), p_value=ks.test(Temp_P1, Temp_P4)[2])%>% select(doy, Dmax,p_value) KS_results_surface$p_value <- as.numeric(KS_results_surface$p_value) KS_results_surface$doy <- as.numeric(KS_results_surface$doy) KS_results_surface_summary <- ddply(KS_results_surface,c("doy"), summarise, Dmax=mean(Dmax, na.rm=TRUE), sd_Dmax = sd(Dmax, na.rm=TRUE), p_value = mean(p_value), sd_P = sd(p_value, na.rm=TRUE)) coeff <- 1.25 # for second y-axis c <- ggplot()+ geom_line(data = KS_results_surface, size=1, aes( doy, Dmax, color = "Dmax"))+ geom_line(data = KS_results_surface, size=1, aes( doy, p_value/coeff, color = "p-value"))+ scale_colour_manual(name = "Colour", values = c("Dmax" = "black", "p-value" = "red"))+ scale_y_continuous(name = "Dmax", sec.axis = sec_axis(~.*coeff, name="p-value"))+ ggtitle("Kolmogorov-Smirnov-Test of daily temperature distribution for RP31-60 and RP91-20")+theme(legend.position = "right")+ xlab("Day of the Year")+ font("title", size = 22,face = "bold")+ font("subtitle", size = 10, color = "orange")+ font("caption", size = 10, color = "orange")+ font("xlab", size = 16)+ font("ylab", size = 16)+ font("xy.text", size = 12, face = "bold")+ theme(legend.position = "none")+ theme(panel.grid.minor=element_blank(),panel.grid.major=element_blank())+ theme(panel.background=element_rect(fill="transparent"), panel.border=element_rect(colour="black", fill=NA, size=0.2), strip.text=element_text(size=14)) KS_results_surface_summary$Significance <- cut(KS_results_surface_summary$p_value, breaks = c(0,0.05, 1)) levels(KS_results_surface_summary$Significance) = c("p below 5%", "p above 5%") Signi_sur <- RefP1_clim_sur %>% select(doy) Signi_sur$y <- 0 Signi_sur$Significance <- KS_results_surface_summary$Significance # check how many days are different Signi_sur %>% group_by(Significance)%>% summarise(no_rows = length(Significance)) #1 p above 5% 30 #2 p below 5% 336 Signi_colors <- c('p above 5%' ="#F98400", "p below 5%" = "#5BBCD6") a <- ggplot()+ geom_line(data = RefP1_clim_sur,size=1.5, aes(x=doy, y=smooth_Reftemp, color = "1931-1960 baseline and threshold"))+ geom_line(data = RefP1_clim_sur,size=1.5, aes(x=doy, y=smooth_RefQ90, color = "1931-1960 baseline and threshold"), linetype = "dashed")+ geom_line(data = RefP4_clim_sur, size=1.5, aes(x=doy, y=smooth_RefQ90, color = "1991-2020 baseline and threshold"), linetype = "dashed")+ geom_line(data = RefP4_clim_sur,size=1.5, aes(x=doy, y= smooth_Reftemp, color = "1991-2020 baseline and threshold"))+ #geom_point(data = Signi_bot, aes(x=doy, y= y, color = Significance), size = 1)+ geom_tile(data= Signi_sur, aes(x = doy, y = -1, fill = Significance, height = 1.5))+ scale_colour_manual(name = "Colour", values = c("1931-1960 baseline and threshold" = "black", "1991-2020 baseline and threshold" = "grey", "1931-1960 baseline and threshold" = "black", "1991-2020 baseline and threshold" = "grey")) + scale_fill_manual(values = Signi_colors)+ ylab("Temperature [°C]")+ xlab("Day of the Year")+ ggtitle("Fig.3: Kolmogorov-Smirnov-Test of daily temperature distribution for RP31-60 and RP91-20")+ #theme(legend.position = "right")+ scale_x_continuous(expand = c(0,0), breaks = c(50, 100,150, 200, 250, 300, 350))+ font("title", size = 40,face = "bold")+ font("subtitle", size = 10, color = "orange")+ font("caption", size = 10, color = "orange")+ font("xlab", size = 22)+ font("ylab", size = 22)+ font("xy.text", size = 16, face = "bold")+ theme(legend.position = "none", legend.title = element_text(size = 18), legend.text = element_text(size = 26))+ #removed legend for final Figure theme(panel.grid.minor=element_blank(),panel.grid.major=element_blank())+ theme(panel.background=element_rect(fill="transparent"), panel.border=element_rect(colour="black", fill=NA, size=0.2), strip.text=element_text(size=20)) # Applying KS test to temperature distributions of every day of the year - bottom layer KS_results_bottom <- merged_bottom %>% group_by(doy)%>% mutate( Dmax= unlist(ks.test(Temp_P1,Temp_P4)[1]), p_value=ks.test(Temp_P1, Temp_P4)[2])%>% select(doy, Dmax,p_value) KS_results_bottom$p_value <- as.numeric(KS_results_bottom$p_value) KS_results_bottom$doy <- as.numeric(KS_results_bottom$doy) KS_results_bottom_summary <- ddply(KS_results_bottom,c("doy"), summarise, Dmax=mean(Dmax, na.rm=TRUE), sd_Dmax = sd(Dmax, na.rm=TRUE), p_value = mean(p_value), sd_P = sd(p_value, na.rm=TRUE)) d<-ggplot()+ geom_line(data = KS_results_bottom, size=1, aes( doy, Dmax, color = "Dmax"))+ geom_line(data = KS_results_bottom, size=1, aes( doy, p_value/coeff, color = "p-value"))+ scale_colour_manual(name = "Colour", values = c("Dmax" = "black", "p-value" = "red"))+ scale_y_continuous(name = "Dmax", sec.axis = sec_axis(~.*coeff, name="p-value"))+ xlab("Day of the Year")+ font("title", size = 22,face = "bold")+ font("subtitle", size = 10, color = "orange")+ font("caption", size = 10, color = "orange")+ font("xlab", size = 16)+ font("ylab", size = 16)+ font("xy.text", size = 12, face = "bold")+ theme(legend.position = "bottom")+ theme(panel.grid.minor=element_blank(),panel.grid.major=element_blank())+ theme(panel.background=element_rect(fill="transparent"), panel.border=element_rect(colour="black", fill=NA, size=0.2), strip.text=element_text(size=14)) KS_results_bottom_summary$Significance <- cut(KS_results_bottom_summary$p_value, breaks = c(0,0.05, 1)) levels(KS_results_bottom_summary$Significance) = c("p below 5%", "p above 5%") Signi_bot <- RefP1_clim_bot %>% select(doy) Signi_bot$y <- 0 Signi_bot$Significance <- KS_results_bottom_summary$Significance Signi_bot %>% group_by(Significance)%>% summarise(no_rows = length(Significance)) #1 p above 5% 62 #2 p below 5% 304 b<-ggplot()+ geom_line(data = RefP1_clim_bot,size=1.5, aes(x=doy, y=smooth_Reftemp, color = "1931-1960 baseline and threshold"))+ geom_line(data = RefP1_clim_bot,size=1.5, aes(x=doy, y=smooth_RefQ90, color = "1931-1960 baseline and threshold"), linetype = "dashed")+ geom_line(data = RefP4_clim_bot, size=1.5, aes(x=doy, y=smooth_RefQ90, color = "1991-2020 baseline and threshold"), linetype = "dashed")+ geom_line(data = RefP4_clim_bot,size=1.5, aes(x=doy, y= smooth_Reftemp, color = "1991-2020 baseline and threshold"))+ #geom_point(data = Signi_bot, aes(x=doy, y= y, color = Significance), size = 1)+ geom_tile(data= Signi_bot, aes(x = doy, y = -0.75, fill = Significance, height = 1))+ scale_colour_manual(name = "Colour", values = c("1931-1960 baseline and threshold" = "black", "1991-2020 baseline and threshold" = "grey", "1931-1960 baseline and threshold" = "black", "1991-2020 baseline and threshold" = "grey")) + scale_fill_manual(values = Signi_colors)+ ylab("Temperature [°C]")+ xlab("Day of the Year")+ #ggtitle("Kolmogorov-Smirnov-Test Reference Periods 1931-1960 and 1991-2020 bottom layer")+ #theme(legend.position = "right")+ scale_x_continuous(expand = c(0,0), breaks = c(50, 100,150, 200, 250, 300, 350))+ font("title", size = 22,face = "bold")+ font("subtitle", size = 10, color = "orange")+ font("caption", size = 10, color = "orange")+ font("xlab", size = 22)+ font("ylab", size = 22)+ font("xy.text", size = 16, face = "bold")+ theme(legend.position = "bottom", legend.title = element_text(size = 18), legend.text = element_text(size = 16))+ theme(panel.grid.minor=element_blank(),panel.grid.major=element_blank())+ theme(panel.background=element_rect(fill="transparent"), panel.border=element_rect(colour="black", fill=NA, size=0.2), strip.text=element_text(size=20)) ### combine plots lay <- rbind(c(1), c(2)) Fig3 <- grid.arrange(a,b, layout_matrix = lay) SupplementaryFig1 <- grid.arrange(c,d, layout_matrix = lay) ggsave(file="Fig3.svg", plot=Fig3, width=26, height=11) ggsave(file="SupplementaryFig1.svg", plot=SupplementaryFig1, width=26, height=11) dev.off() #######################################################