##### https://sites.google.com/view/climate-access-cooperative/code?authuser=0 ##### Help regarding opening netCDF files library(ncdf4) library(ncdf4.helpers) library(tidyverse) dev.off() rm(list=ls(all=T)) setwd() # set working directory # Load CMEMS SST L4 product for Storfjärden Pixel ncin <- nc_open("DMI_BAL_SST_L4_REP_OBSERVATIONS_010_016_1634224072972.nc") ncin ncin$dim$time #get longitutes and latitudes lon <- ncvar_get(ncin,"lon") summary(lon) nlon <- dim(lon) lat <- ncvar_get(ncin,"lat", verbose=F) nlat <- dim(lat) which.min(abs(lat)) sst <- ncvar_get(ncin, "analysed_sst") lon_index <- which.min(abs(lon - 23.26)) #get pixel closest to storfjärden lon_index lat_index <- which.min(abs(lat - 59.86)) #get pixel closest to storfjärden sst_time <- nc.get.time.series(ncin, v = "analysed_sst", time.dim.name = "time") list(X = lon_index, Y = lat_index) sst <- nc.get.var.subset.by.axes(ncin, "analysed_sst", axis.indices = list(X = lon_index, Y = lat_index)) sst Sat_sst <- data_frame(time = sst_time, sst = as.vector(sst)) %>% mutate(time = as.Date(format(time, "%Y-%m-%d"))) Sat_sst$temp <- Sat_sst$sst - 273.15 #convert to degrees celsius #################### Progress with Correlation and KS-test library(plyr) library(dplyr) library(lubridate) library(RColorBrewer) library(gridExtra) library(ggpubr) library(zoo) filtered_Sat_sst<-Sat_sst%>% mutate(Date = as.Date(time, format = "%Y-%m-%d"))%>% mutate(year = lubridate::year(Date), month = lubridate::month(Date), day = lubridate::day(Date), doy = lubridate::yday(Date)) # Create climatological baseline and thresholds for CMEMS SST for 1982-2011 RefP_COP<- filtered_Sat_sst%>% filter(Date >= "1982-01-01")%>% filter(Date <= "2011-12-31")%>% mutate(Date = as.Date(time, format = "%Y-%m-%d")) RefP_COP_Average<-ddply(RefP_COP,c("doy"), summarise, RefP_COP_Temp=mean(temp, na.rm=TRUE), sd = sd(temp, na.rm=TRUE), n = sum(!is.na(temp)), se = sd/sqrt(n)) names(RefP_COP_Average) Quantiles_RefP_COP_10_90<-ddply(RefP_COP, .(doy), function(x) quantile(x$temp, c(0.1,0.9), na.rm = TRUE)) names(Quantiles_RefP_COP_10_90) Sat_Clim<-filtered_Sat_sst%>% left_join(Quantiles_RefP_COP_10_90, by=c("doy"))%>% left_join(RefP_COP_Average)%>% select(Date, doy, temp,RefP_COP_Temp, '10%', '90%') names(Sat_Clim) colnames(Sat_Clim)[5]<-"RefP_COP_Q10" colnames(Sat_Clim)[6]<-"RefP_COP_Q90" # Get smoothed climatological mean and thresholds through moving average of 31 day range centered on desired doy RefP_COP_values <- Sat_Clim %>% select(doy, RefP_COP_Temp, RefP_COP_Q10, RefP_COP_Q90)%>% distinct(doy, .keep_all = TRUE) Smoothing_RefP_COP_start <- RefP_COP_values%>% filter (doy >= 352 & doy <=366) Smoothing_RefP_COP_end <- RefP_COP_values%>% filter (doy >= 1 & doy <= 15) Smoothing_RefP_COP <- rbind(Smoothing_RefP_COP_start, RefP_COP_values,Smoothing_RefP_COP_end ) %>% mutate (smooth_RefCOPtemp = rollmean(RefP_COP_Temp, k=31, fill = NA), smooth_RefCOPQ90 = rollmean(RefP_COP_Q90, k=31, fill = NA), smooth_RefCOPQ10 = rollmean(RefP_COP_Q10, k=31, fill = NA))%>% drop_na() ####################################### Get in situ surface data and keep only measured data for correlation with all data after 1982-01-01 (until 2019-08-30, end of L4) in_situ <- read.table("20220513_Storfjarden_interpol.txt",sep=" ",dec=".",header=TRUE, stringsAsFactors = TRUE) filtered_in_situ <- in_situ%>% filter(factored_depth == "surface")%>% select(Date, factored_depth, MeanTemp,value)%>% mutate(Date = as.Date(Date), doy = lubridate::yday(Date))%>% filter (Date > "1981-12-31")%>% filter (value == "real") ######################################## Check Correlation SST_vs_insitu <- Sat_Clim%>% left_join(filtered_in_situ, by = c("Date", "doy"))%>% select(Date,doy, MeanTemp, temp) SST_vs_insitu <- SST_vs_insitu[complete.cases(SST_vs_insitu),] ### check max deviation SST_vs_insitu$Max <- SST_vs_insitu$MeanTemp - SST_vs_insitu$temp ggplot(SST_vs_insitu, aes(Date, Max))+ geom_point() Max <- SST_vs_insitu%>% filter(Max > 10) ### continue correlation Fig7 <- ggscatter(SST_vs_insitu,x = "temp", y = "MeanTemp", add = "reg.line", conf.int = TRUE,font.label = c(20, "bold"), cor.coef = FALSE, cor.method = "pearson", xlab = "RS-SST [°C]", ylab = "in situ SST [°C]", title = "Fig.7: Correlation of satellite derived SST and in situ SST")+ stat_cor(method = "pearson", size = 6, label.x.npc = "left", label.y.npc = "top")+ font("title", size = 22,face = "bold")+ font("subtitle", size = 10, color = "orange")+ font("caption", size = 10, color = "orange")+ font("xlab", size = 20)+ font("ylab", size = 20)+ font("xy.text", size = 16, 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=20)) ggsave(file="Fig7.svg", plot=Fig7, width=9, height=6) ####################################### Prepare datasets for KS-Test of 1982-2011 Temp_RefP_insitu <- in_situ%>% filter(factored_depth == "surface")%>% mutate(Date = as.Date(Date, format = "%Y-%m-%d"))%>% filter (Date >= "1982-01-01")%>% filter (Date <= "2011-12-31")%>% mutate (In_situ_Temp = MeanTemp)%>% select (Date, In_situ_Temp) Temp_RefP_COP <- Sat_Clim %>% mutate(Date = as.Date(Date, format = "%Y-%m-%d"))%>% filter (Date >= "1982-01-01")%>% filter (Date <= "2011-12-31")%>% mutate (Cop_SST = temp)%>% select (Date, Cop_SST) KS_COP_insitu <- Temp_RefP_insitu %>% left_join( Temp_RefP_COP)%>% mutate(doy = lubridate::yday(Date))%>% mutate(doy = as.factor(doy)) which(is.na(KS_COP_insitu$Cop_SST)) # one data point missing in CMEMS data KS_COP_insitu[10296,3] <- -0.36 # average of 10295 and 10297 ############################################################### KS Test KS_results <- KS_COP_insitu %>% select(doy, In_situ_Temp, Cop_SST)%>% group_by(doy)%>% mutate( Dmax= unlist(ks.test(In_situ_Temp,Cop_SST)[1]), p_value=ks.test(In_situ_Temp, Cop_SST)[2])%>% select(doy, Dmax,p_value) KS_results$p_value <- as.numeric(KS_results$p_value) KS_results$doy <- as.numeric(KS_results$doy) KS_results_summary <- ddply(KS_results,c("doy"), summarise, Dmax=mean(Dmax, na.rm=TRUE), p_value = mean(p_value)) coeff <- 1.25 # for second y axis SupplFig2 <- ggplot()+ geom_line(data = KS_results, size=1, aes( doy, Dmax, color = "Dmax"))+ geom_line(data = KS_results, 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 in situ and remotely sensed SST for RP82-11 ")+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 = "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)) ggsave(file="SupplementaryFig2.svg", plot=SupplFig2, width=18, height=11) KS_results_summary$Significance <- cut(KS_results_summary$p_value, breaks = c(0,0.05, 1)) levels(KS_results_summary$Significance) = c("p below 5%", "p above 5%") ############# Load climatological mean and thresholds of each period for plot in_situ_clim <- read.table("20220516_Reference Period 3_Sat - 1982-2011.txt",sep=" ",dec=".",header=TRUE, stringsAsFactors = TRUE) filtered_in_situ_clim <- in_situ_clim %>% filter(factored_depth == "surface")%>% select( doy, smooth_Reftemp, smooth_RefQ90) Summary <- KS_results_summary%>% left_join(Smoothing_RefP_COP, by =c("doy"))%>% left_join(filtered_in_situ_clim)%>% select (doy, smooth_Reftemp,smooth_RefQ90, smooth_RefCOPtemp, smooth_RefCOPQ90, Significance) Summary$Signi_y <- 0 Signi_colors <- c('p above 5%' ="#F98400", "p below 5%" = "#5BBCD6") Fig8<- ggplot(data = Summary, aes (doy))+ geom_line(size=1.5, aes( y=smooth_RefCOPtemp, color = "satellite-derived baseline and threshold"))+ geom_line(size=1.5, aes( y=smooth_RefCOPQ90, color = "satellite-derived baseline and threshold"), linetype = "dashed")+ geom_line(size=1.5, aes( y=smooth_RefQ90, color = "in situ-surface baseline and threshold"), linetype = "dashed")+ geom_line(size=1.5, aes( y= smooth_Reftemp, color = "in situ-surface baseline and threshold"))+ #geom_point(aes( y= Signi_y, color = Significance), size = 1)+ geom_tile(aes( y = -1, fill = Significance, height = 1.5))+ scale_colour_manual(name = "Colour", values = c("satellite-derived baseline and threshold" = "grey", "in situ-surface baseline and threshold" = "black", "satellite-derived baseline and threshold" = "grey", "in situ-surface baseline and threshold" = "black")) + scale_fill_manual(values = Signi_colors)+ ylab("Temperature [°C]")+ xlab("Day of the Year")+ ggtitle("Fig.8: Kolmogorov-Smirnov-Test of in situ and satellite derived SST for RP82-11")+ theme(legend.position = "bottom")+ scale_x_continuous(expand = c(0,0), breaks = c(50, 100,150, 200, 250, 300, 350))+ font("title", size = 28,face = "bold")+ font("subtitle", size = 10, color = "orange")+ font("caption", size = 10, color = "orange")+ font("xlab", size = 24)+ font("ylab", size = 24)+ font("xy.text", size = 18, face = "bold")+ theme(legend.position = "bottom", legend.title = element_text(size = 20), legend.text = element_text(size = 18))+ 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=22)) ggsave(file="Fig8.svg", plot=Fig8, width=24, height=11) #### How many doys are different? names(Summary) Summary %>% group_by(Significance)%>% summarise(no_rows = length(Significance)) count(Summary$Significance) #1 p above 5% 174 #2 p below 5% 192