library(tidyverse) library(dplyr) library(lubridate) library(plyr) library(ggpubr) library(gridExtra) library(zoo) dev.off() rm(list=ls(all=T)) setwd() # set working directory data <- read.table("20220513_Storfjarden_interpol.txt",sep=" ",dec=".",header=TRUE, stringsAsFactors = TRUE) ############################### Prepare filtered Dataframe for creation of 5 reference periods ################# ### Filter out Surface and Bottom, converting Date 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)) ############## Set Reference Periods: RefP1 01.01.1931-31.12.1960, RefP2WMO 01.01.1961-31.12.1990, RefP3 01.01.1982-30.12.2011, RefP4 01.01.1991-30.12.2020, RefP_all 01.01.1931-2019 ############## RefP1 RefP1_filtered <- filtered_data%>% filter(Date >= "1931-01-01")%>% filter(Date <= "1960-12-31") # Get daily 10th and 90th percentile for each depth layer Quantiles_RefP1_10_90<-ddply(RefP1_filtered, .(doy, factored_depth), function(x) quantile(x$MeanTemp, c(0.1,0.9), na.rm = TRUE)) # Get daily climatological mean RefP1Average<-ddply(RefP1_filtered,c("factored_depth", "doy"), summarise, RefP1_Temp=mean(MeanTemp, na.rm=TRUE), sd = sd(MeanTemp, na.rm=TRUE), n = sum(!is.na(MeanTemp)), se = sd/sqrt(n)) # Combine climatological mean and thresholds RefP1Average<-RefP1Average%>% left_join(Quantiles_RefP1_10_90, by=c("doy", "factored_depth")) names(RefP1Average) colnames(RefP1Average)[5]<-"RefP1_n" colnames(RefP1Average)[7]<-"RefP1_Q10" colnames(RefP1Average)[8]<-"RefP1_Q90" # Get smoothed climatological mean and thresholds through moving average of 31 day range centered on desired doy Smoothing_RefP1_surface_start <- RefP1Average%>% filter (factored_depth=="surface")%>% filter (doy >= 352 & doy <=366) Smoothing_RefP1_surface_end <- RefP1Average%>% filter (factored_depth=="surface")%>% filter (doy >= 1 & doy <= 15) Smoothing_RefP1_surface <- rbind(Smoothing_RefP1_surface_start, subset(RefP1Average, factored_depth == "surface"),Smoothing_RefP1_surface_end ) %>% mutate (smooth_Reftemp = rollmean(RefP1_Temp, k=31, fill = NA), smooth_RefQ90 = rollmean(RefP1_Q90, k=31, fill = NA), smooth_RefQ10 = rollmean(RefP1_Q10, k=31, fill = NA))%>% drop_na() Smoothing_RefP1_bottom_start <- RefP1Average%>% filter (factored_depth=="bottom")%>% filter (doy >= 352 & doy <=366) Smoothing_RefP1_bottom_end <- RefP1Average%>% filter (factored_depth=="bottom")%>% filter (doy >= 1 & doy <= 15) Smoothing_RefP1_bottom <- rbind(Smoothing_RefP1_bottom_start, subset(RefP1Average, factored_depth == "bottom"),Smoothing_RefP1_bottom_end ) %>% mutate (smooth_Reftemp = rollmean(RefP1_Temp, k=31, fill = NA), smooth_RefQ90 = rollmean(RefP1_Q90, k=31, fill = NA), smooth_RefQ10 = rollmean(RefP1_Q10, k=31, fill = NA))%>% drop_na() smoothed_RefP1 <- rbind (Smoothing_RefP1_surface, Smoothing_RefP1_bottom) write.table(smoothed_RefP1, "20220516_Reference Period 1 - 1931-1960.txt") ############## RefP2WMO RefP2WMO_filtered <- filtered_data%>% filter(Date >= "1961-01-01")%>% filter(Date <= "1990-12-31") # Get daily 10th and 90th percentile for each depth layer # Get daily climatological mean # Combine climatological mean and thresholds Quantiles_RefP2WMO_10_90<-ddply(RefP2WMO_filtered, .(doy, factored_depth), function(x) quantile(x$MeanTemp, c(0.1,0.9))) RefP2WMOAverage<-ddply(RefP2WMO_filtered,c("factored_depth", "doy"), summarise, RefP2WMO_Temp=mean(MeanTemp, na.rm=TRUE), sd = sd(MeanTemp, na.rm=TRUE), n = sum(!is.na(MeanTemp)), se = sd/sqrt(n)) RefP2WMOAverage<-RefP2WMOAverage%>% left_join(Quantiles_RefP2WMO_10_90, by=c("doy", "factored_depth")) names(RefP2WMOAverage) colnames(RefP2WMOAverage)[5]<-"RefP2WMO_n" colnames(RefP2WMOAverage)[7]<-"RefP2WMO_Q10" colnames(RefP2WMOAverage)[8]<-"RefP2WMO_Q90" # Get smoothed climatological mean and thresholds through moving average of 31 day range centered on desired doy Smoothing_RefP2WMO_surface_start <- RefP2WMOAverage%>% filter (factored_depth=="surface")%>% filter (doy >= 352 & doy <=366) Smoothing_RefP2WMO_surface_end <- RefP2WMOAverage%>% filter (factored_depth=="surface")%>% filter (doy >= 1 & doy <= 15) Smoothing_RefP2WMO_surface <- rbind(Smoothing_RefP2WMO_surface_start, subset(RefP2WMOAverage, factored_depth == "surface"),Smoothing_RefP2WMO_surface_end ) %>% mutate (smooth_Reftemp = rollmean(RefP2WMO_Temp, k=31, fill = NA), smooth_RefQ90 = rollmean(RefP2WMO_Q90, k=31, fill = NA), smooth_RefQ10 = rollmean(RefP2WMO_Q10, k=31, fill = NA))%>% drop_na() Smoothing_RefP2WMO_bottom_start <- RefP2WMOAverage%>% filter (factored_depth=="bottom")%>% filter (doy >= 352 & doy <=366) Smoothing_RefP2WMO_bottom_end <- RefP2WMOAverage%>% filter (factored_depth=="bottom")%>% filter (doy >= 1 & doy <= 15) Smoothing_RefP2WMO_bottom <- rbind(Smoothing_RefP2WMO_bottom_start, subset(RefP2WMOAverage, factored_depth == "bottom"),Smoothing_RefP2WMO_bottom_end ) %>% mutate (smooth_Reftemp = rollmean(RefP2WMO_Temp, k=31, fill = NA), smooth_RefQ90 = rollmean(RefP2WMO_Q90, k=31, fill = NA), smooth_RefQ10 = rollmean(RefP2WMO_Q10, k=31, fill = NA))%>% drop_na() smoothed_RefP2WMO <- rbind (Smoothing_RefP2WMO_surface, Smoothing_RefP2WMO_bottom) write.table(smoothed_RefP2WMO, "20220516_Reference Period 2WMO - 1961-1990.txt") ############## RefP3 RefP3_filtered <- filtered_data%>% filter(Date >= "1982-01-01")%>% filter(Date <= "2011-12-31") # Get daily 10th and 90th percentile for each depth layer Quantiles_RefP3_10_90<-ddply(RefP3_filtered, .(doy, factored_depth), function(x) quantile(x$MeanTemp, c(0.1,0.9))) # Get daily climatological mean RefP3Average<-ddply(RefP3_filtered,c("factored_depth", "doy"), summarise, RefP3_Temp=mean(MeanTemp, na.rm=TRUE), sd = sd(MeanTemp, na.rm=TRUE), n = sum(!is.na(MeanTemp)), se = sd/sqrt(n)) # Combine climatological mean and thresholds RefP3Average<-RefP3Average%>% left_join(Quantiles_RefP3_10_90, by=c("doy", "factored_depth")) names(RefP3Average) colnames(RefP3Average)[5]<-"RefP3_n" colnames(RefP3Average)[7]<-"RefP3_Q10" colnames(RefP3Average)[8]<-"RefP3_Q90" # Get smoothed climatological mean and thresholds through moving average of 31 day range centered on desired doy Smoothing_RefP3_surface_start <- RefP3Average%>% filter (factored_depth=="surface")%>% filter (doy >= 352 & doy <=366) Smoothing_RefP3_surface_end <- RefP3Average%>% filter (factored_depth=="surface")%>% filter (doy >= 1 & doy <= 15) Smoothing_RefP3_surface <- rbind(Smoothing_RefP3_surface_start, subset(RefP3Average, factored_depth == "surface"),Smoothing_RefP3_surface_end ) %>% mutate (smooth_Reftemp = rollmean(RefP3_Temp, k=31, fill = NA), smooth_RefQ90 = rollmean(RefP3_Q90, k=31, fill = NA), smooth_RefQ10 = rollmean(RefP3_Q10, k=31, fill = NA))%>% drop_na() Smoothing_RefP3_bottom_start <- RefP3Average%>% filter (factored_depth=="bottom")%>% filter (doy >= 352 & doy <=366) Smoothing_RefP3_bottom_end <- RefP3Average%>% filter (factored_depth=="bottom")%>% filter (doy >= 1 & doy <= 15) Smoothing_RefP3_bottom <- rbind(Smoothing_RefP3_bottom_start, subset(RefP3Average, factored_depth == "bottom"),Smoothing_RefP3_bottom_end ) %>% mutate (smooth_Reftemp = rollmean(RefP3_Temp, k=31, fill = NA), smooth_RefQ90 = rollmean(RefP3_Q90, k=31, fill = NA), smooth_RefQ10 = rollmean(RefP3_Q10, k=31, fill = NA))%>% drop_na() smoothed_RefP3 <- rbind (Smoothing_RefP3_surface, Smoothing_RefP3_bottom) write.table(smoothed_RefP3, "20220516_Reference Period 3_Sat - 1982-2011.txt") ############## RefP4 RefP4_filtered <- filtered_data%>% filter(Date >= "1991-01-01")%>% filter(Date <= "2020-12-31") # Get daily 10th and 90th percentile for each depth layer Quantiles_RefP4_10_90<-ddply(RefP4_filtered, .(doy, factored_depth), function(x) quantile(x$MeanTemp, c(0.1,0.9), na.rm = TRUE)) # Get daily climatological mean RefP4Average<-ddply(RefP4_filtered,c("factored_depth", "doy"), summarise, RefP4_Temp=mean(MeanTemp, na.rm=TRUE), sd = sd(MeanTemp, na.rm=TRUE), n = sum(!is.na(MeanTemp)), se = sd/sqrt(n)) # Combine climatological mean and thresholds RefP4Average<-RefP4Average%>% left_join(Quantiles_RefP4_10_90, by=c("doy", "factored_depth")) names(RefP4Average) colnames(RefP4Average)[5]<-"RefP4_n" colnames(RefP4Average)[7]<-"RefP4_Q10" colnames(RefP4Average)[8]<-"RefP4_Q90" # Get smoothed climatological mean and thresholds through moving average of 31 day range centered on desired doy Smoothing_RefP4_surface_start <- RefP4Average%>% filter (factored_depth=="surface")%>% filter (doy >= 352 & doy <=366) Smoothing_RefP4_surface_end <- RefP4Average%>% filter (factored_depth=="surface")%>% filter (doy >= 1 & doy <= 15) Smoothing_RefP4_surface <- rbind(Smoothing_RefP4_surface_start, subset(RefP4Average, factored_depth == "surface"),Smoothing_RefP4_surface_end ) %>% mutate (smooth_Reftemp = rollmean(RefP4_Temp, k=31, fill = NA), smooth_RefQ90 = rollmean(RefP4_Q90, k=31, fill = NA), smooth_RefQ10 = rollmean(RefP4_Q10, k=31, fill = NA))%>% drop_na() Smoothing_RefP4_bottom_start <- RefP4Average%>% filter (factored_depth=="bottom")%>% filter (doy >= 352 & doy <=366) Smoothing_RefP4_bottom_end <- RefP4Average%>% filter (factored_depth=="bottom")%>% filter (doy >= 1 & doy <= 15) Smoothing_RefP4_bottom <- rbind(Smoothing_RefP4_bottom_start, subset(RefP4Average, factored_depth == "bottom"),Smoothing_RefP4_bottom_end ) %>% mutate (smooth_Reftemp = rollmean(RefP4_Temp, k=31, fill = NA), smooth_RefQ90 = rollmean(RefP4_Q90, k=31, fill = NA), smooth_RefQ10 = rollmean(RefP4_Q10, k=31, fill = NA))%>% drop_na() smoothed_RefP4 <- rbind (Smoothing_RefP4_surface, Smoothing_RefP4_bottom) write.table(smoothed_RefP4, "20220516_Reference Period 4 - 1991-2020.txt") ############## RefP_all RefP_all_filtered <- filtered_data%>% filter(Date >= "1931-01-01")%>% filter(Date <= "2020-12-31") # Get daily 10th and 90th percentile for each depth layer Quantiles_RefP_all_10_90<-ddply(RefP_all_filtered, .(doy, factored_depth), function(x) quantile(x$MeanTemp, c(0.1,0.9), na.rm = TRUE)) # Get daily climatological mean RefP_allAverage<-ddply(RefP_all_filtered,c("factored_depth", "doy"), summarise, RefP_all_Temp=mean(MeanTemp, na.rm=TRUE), sd = sd(MeanTemp, na.rm=TRUE), n = sum(!is.na(MeanTemp)), se = sd/sqrt(n)) # Combine climatological mean and thresholds RefP_allAverage<-RefP_allAverage%>% left_join(Quantiles_RefP_all_10_90, by=c("doy", "factored_depth")) names(RefP_allAverage) colnames(RefP_allAverage)[5]<-"RefP_all_n" colnames(RefP_allAverage)[7]<-"RefP_all_Q10" colnames(RefP_allAverage)[8]<-"RefP_all_Q90" # Get smoothed climatological mean and thresholds through moving average of 31 day range centered on desired doy Smoothing_RefP_all_surface_start <- RefP_allAverage%>% filter (factored_depth=="surface")%>% filter (doy >= 352 & doy <=366) Smoothing_RefP_all_surface_end <- RefP_allAverage%>% filter (factored_depth=="surface")%>% filter (doy >= 1 & doy <= 15) Smoothing_RefP_all_surface <- rbind(Smoothing_RefP_all_surface_start, subset(RefP_allAverage, factored_depth == "surface"),Smoothing_RefP_all_surface_end ) %>% mutate (smooth_Reftemp = rollmean(RefP_all_Temp, k=31, fill = NA), smooth_RefQ90 = rollmean(RefP_all_Q90, k=31, fill = NA), smooth_RefQ10 = rollmean(RefP_all_Q10, k=31, fill = NA))%>% drop_na() Smoothing_RefP_all_bottom_start <- RefP_allAverage%>% filter (factored_depth=="bottom")%>% filter (doy >= 352 & doy <=366) Smoothing_RefP_all_bottom_end <- RefP_allAverage%>% filter (factored_depth=="bottom")%>% filter (doy >= 1 & doy <= 15) Smoothing_RefP_all_bottom <- rbind(Smoothing_RefP_all_bottom_start, subset(RefP_allAverage, factored_depth == "bottom"),Smoothing_RefP_all_bottom_end ) %>% mutate (smooth_Reftemp = rollmean(RefP_all_Temp, k=31, fill = NA), smooth_RefQ90 = rollmean(RefP_all_Q90, k=31, fill = NA), smooth_RefQ10 = rollmean(RefP_all_Q10, k=31, fill = NA))%>% drop_na() smoothed_RefP_all <- rbind (Smoothing_RefP_all_surface, Smoothing_RefP_all_bottom) #write.table(smoothed_RefP_all, "20220516_Reference Period all - 1931 - 2020.txt")