update following paper revision

Updated uncertainty analysis, updated scripts for figures and updated model validation end date.

R-scripts/Figure3_Contour_plots.R

@@ -9,11 +9,15 @@ library(ggplot2)
 library(ggpubr)
 
 #Set working directory
-setwd(".../FCR-GLM-metrics")
+#setwd(".../FCR-GLM-metrics")
+#setwd('/Users/Kamilla/Compile_15_11/Bubbles')
+setwd('/Users/Kamilla/Compile_old_GLM1602/Bubbles')
+sim_folder <- getwd()
+output <- nc_open("output/output.nc")
 
 #Modelled oxy
-sim_folder <- getwd()
-output <- nc_open("Calibrated_models/Deepm2_naive/output/output.nc")
+#sim_folder <- getwd()
+#output <- nc_open("Calibrated_models/Deepm2_naive/output/output.nc")
 oxy<- ncvar_get(output, "OXY_oxy")
 depth<- ncvar_get(output, "z")
 depth[depth >= 100] <- NA
@@ -62,8 +66,8 @@ ct<- na.omit(ct)
 df <- cbind(co$Var2, co$value, ct$value)
 
 #Creating dataframe for time 
-time <- data.frame(seq(as.Date("2016-12-02"), as.Date("2019-12-31"), by="day"))
-ID <- seq.int(1:1125)
+time <- data.frame(seq(as.Date("2016-12-20"), as.Date("2019-12-30"), by="day"))
+ID <- seq.int(1:1106)
 time <- cbind(ID, time)
 colnames(time) <- c("ID", "DateTime")
 colnames(df) <- c("ID", "Oxy", "Depth")
@@ -223,8 +227,10 @@ melt_d<- na.omit(melt_d)
 df <- cbind(melt_t$Var2, melt_t$value, melt_d$value)
 
 #Creating dataframe for time 
-time <- data.frame(seq(as.Date("2016-12-02"), as.Date("2019-12-31"), by="day"))
-ID <- seq.int(1:1125)
+time <- data.frame(seq(as.Date("2016-12-20"), as.Date("2019-12-30"), by="day"))
+ID <- seq.int(1:1106)
+#time <- data.frame(seq(as.Date("2016-12-02"), as.Date("2019-12-31"), by="day"))
+#ID <- seq.int(1:1125)
 #time <- data.frame(seq(as.Date("2015-07-13"), as.Date("2016-12-01"), by="day"))
 #ID <- seq.int(1:508)
 time <- cbind(ID, time)

R-scripts/Figure5_Oxy_temp_6panel_plots.R → R-scripts/Figure5_Oxy_temp_6panel_plots_error.R

@@ -19,33 +19,24 @@ library(patchwork)
 #Set working directory
 setwd(".../FCR-GLM-metrics")
 sim_folder <- getwd()
-#weight 1
-w1 <- file.path(sim_folder, 'Calibrated_models/Deepm2_exm_weight1/output/output.nc')
-#weight 2
-w2 <- file.path(sim_folder, 'Calibrated_models/Deepm2_exm_weight2/output/output.nc')
-#weight 3
-w3 <- file.path(sim_folder, 'Calibrated_models/Deepm2_exm_weight3/output/output.nc')
-error <- read.csv("observations/error_stats.csv")
-
-#oxygen
-var="OXY_oxy"
-obs_oxy<-read.csv('observations/CleanedObsOxy.csv') %>%
-  mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d", tz="EST"))) %>%
-  filter(DateTime > "2016-12-01")
-obs_oxy$DateTime <- as.Date(obs_oxy$DateTime, format="%Y-%m-%d")
 
-depths<- c(0.1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 9.2) 
+#Models w1, w2, w3
+w1 <- file.path(sim_folder, 'Calibrated_models/Deepm2_exm_weight1/output/output.nc')
 
 oxy_w1 <- get_var(w1, "OXY_oxy", reference="surface", z_out=depths) %>%
   pivot_longer(cols=starts_with("OXY_oxy_"), names_to="Depth", names_prefix="OXY_oxy_", values_to = "OXY_oxy") %>%
   mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d", tz="EST"))) 
 oxy_w1$DateTime <- as.Date(oxy_w1$DateTime, format="%Y-%m-%d")
 
+w2 <- file.path(sim_folder, 'Calibrated_models/Deepm2_exm_weight2/output/output.nc')
+
 oxy_w2 <- get_var(w2, "OXY_oxy", reference="surface", z_out=depths) %>%
   pivot_longer(cols=starts_with("OXY_oxy_"), names_to="Depth", names_prefix="OXY_oxy_", values_to = "OXY_oxy") %>%
   mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d", tz="EST"))) 
 oxy_w2$DateTime <- as.Date(oxy_w2$DateTime, format="%Y-%m-%d")
 
+w3 <- file.path(sim_folder, 'Calibrated_models/Deepm2_exm_weight3/output/output.nc')
+
 oxy_w3 <- get_var(w3, "OXY_oxy", reference="surface", z_out=depths) %>%
   pivot_longer(cols=starts_with("OXY_oxy_"), names_to="Depth", names_prefix="OXY_oxy_", values_to = "OXY_oxy") %>%
   mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d", tz="EST"))) 
@@ -59,6 +50,24 @@ oxygen <- merge(obs_oxy, oxy_w1, by=c("DateTime","Depth")) %>%
   dplyr::rename(mod_oxy_w3 = OXY_oxy)
 oxygen$DateTime <- as.Date(oxygen$DateTime, format="%Y-%m-%d")
 
+
+#oxygen
+var="OXY_oxy"
+obs_oxy<-read.csv('Observations/CleanedObsOxy.csv') %>%
+  mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d", tz="EST"))) %>%
+  filter(DateTime > "2016-12-01")
+obs_oxy$DateTime <- as.Date(obs_oxy$DateTime, format="%Y-%m-%d")
+
+depths<- c(0.1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 9.2) 
+
+#Calculating model absolute error over time (prediction - obs)
+oxygen <- oxygen %>%
+  mutate(error_w1 = abs(mod_oxy_w1 - obsoxy),
+         error_w2 = abs(mod_oxy_w2 - obsoxy),
+         error_w3 = abs(mod_oxy_w3 - obsoxy),
+         obs_dates = 0)
+
+
 #Error calculation w1
 for (i in 1:nrow(oxygen)) {
   oxygen$MEFF_1_w1[i]<- ((oxygen$mod_oxy_w1[i]- oxygen$obsoxy[i])^2)
@@ -87,30 +96,31 @@ for (i in 1:nrow(oxygen)) {
 #Adding calculated MEF to error table
 #error[error$metric=="oxy" & error$calibration=="PEST_exm_w3", "Calibration.deepm2"] <- MEFF_oxy_w3
 
-#Plots
+#Oxygen error plot
 plot_depths <- c(1, 4, 9)
-plot <- vector('list')
+plot_error_DO <- vector('list')
 
 for(i in 1:length(plot_depths)){
 
 
-  plot[[i]] <- oxygen %>%
+  plot_error_DO[[i]] <- oxygen %>%
     filter(Depth == plot_depths[i]) %>%
-    ggplot2::ggplot(ggplot2::aes(x = DateTime, y = obsoxy, colour="Obs. DO")) +
-    ggplot2::geom_point(pch=10)+
-    geom_line(data=filter(oxy_w1, Depth==plot_depths[i]), aes(x=DateTime, y=OXY_oxy, colour="Model w1 DO"), lty=1, size=0.5)+
-    geom_line(data=filter(oxy_w2, Depth==plot_depths[i]), aes(x=DateTime, y=OXY_oxy, colour="Model w2 DO"), lty=2, size=0.5)+
-    geom_line(data=filter(oxy_w3, Depth==plot_depths[i]), aes(x=DateTime, y=OXY_oxy, colour="Model w3 DO"), lty=3, size=0.8)+
-    ggplot2::ggtitle(" ", subtitle=" ")+
+    ggplot2::ggplot(ggplot2::aes(x = DateTime, y = obs_dates, colour="Obs. dates")) +
+    ggplot2::geom_point(pch=4)+
+    geom_line(data=filter(oxygen, Depth==plot_depths[i]), aes(x=DateTime, y=error_w1, colour="Model w1 DO"), lty=2, linewidth=0.5)+
+    geom_line(data=filter(oxygen, Depth==plot_depths[i]), aes(x=DateTime, y=error_w2, colour="Model w2 DO"), lty=1, linewidth=0.5)+
+    geom_line(data=filter(oxygen, Depth==plot_depths[i]), aes(x=DateTime, y=error_w3, colour="Model w3 DO"), lty=3, linewidth=0.5)+
+    ggplot2::ggtitle(" ", subtitle = "Oxygen")+
     xlab("Date")+
-    ylab(expression(bold(Oxygen~(mmol/m^{3}))))+
-    ylim(c(0, 420))+
-    ggplot2::scale_colour_manual(name="Legend", values=c("Obs. DO"="black", "Model w1 DO"="#00BA38", "Model w2 DO"="#00BA38", "Model w3 DO"="#00BA38"), guide=guide_legend(override.aes=list(linetype=c(NA, 1, 2, 3), shape=c(10, NA, NA, NA))))+
+    ylab(expression(bold(Absolute~model~error~(mmol/m^{3}))))+
+    ylim(c(0, 250))+
+    ggplot2::scale_colour_manual(name="Legend", values=c("Obs. dates" = "black", "Model w1 DO"="#00BA38", "Model w2 DO"="#00BA38", "Model w3 DO"="#00BA38"),
+    labels=c("Obs. dates", "Model w1 DO", "Model w2 DO", "Model w3 DO"), guide=guide_legend(override.aes=list(linetype=c(NA, 2, 1, 3), shape=c(4, NA, NA, NA), colour=c("black", "#00BA38", "#00BA38", "#00BA38"))))+
     ggplot2::theme_light() +
     ggplot2::theme(
-      plot.title = ggplot2::element_text(face= "bold", size = 12),
-      axis.title.y = ggplot2::element_text(face = "bold", size= 10),
-      axis.title.x = ggplot2::element_text(face = "bold", size= 10),
+      plot.title = ggplot2::element_text(face= "bold", size = 10),
+      axis.title.y = ggplot2::element_text(face = "bold", size= 8),
+      axis.title.x = ggplot2::element_text(face = "bold", size= 8),
       legend.text = ggplot2::element_text(size= 8),
       axis.text.x = ggplot2::element_text(size=9),
       axis.text.y = ggplot2::element_text(size=9),
@@ -121,16 +131,17 @@ for(i in 1:length(plot_depths)){
 
 
 }
-plot
 
-combinedPlot <- patchwork::wrap_plots(plot, ncol = 1) +
+plot_error_DO
+
+combinedPlot_DO <- patchwork::wrap_plots(plot_error_DO, ncol = 1) +
   patchwork::plot_layout(guides = "collect") &
   theme(legend.position = "bottom")
 
-combinedPlot
+combinedPlot_DO
 
 #Temperature
-obs_temp<-read.csv('observations/CleanedObsTemp.csv') %>%
+obs_temp<-read.csv('Observations/CleanedObsTemp.csv') %>%
   mutate(DateTime = as.POSIXct(strptime(DateTime, "%Y-%m-%d", tz="EST"))) %>%
   filter(DateTime > "2016-12-01" & DateTime < "2020-01-01")
 obs_temp$DateTime <- as.Date(obs_temp$DateTime, format="%Y-%m-%d")
@@ -151,13 +162,19 @@ temp_w3 <- get_var(w3, "temp", reference="surface", z_out=depths) %>%
 temp_w3$DateTime <- as.Date(temp_w3$DateTime, format="%Y-%m-%d")
 
 temp <- merge(obs_temp, temp_w1, by=c("DateTime","Depth")) %>%
-  dplyr::rename(obtemp = temp.x, modtemp = temp.y) %>%
+  dplyr::rename(obtemp = temp.x, mod_temp_w1 = temp.y) %>%
   merge(temp_w2, by=c("DateTime","Depth")) %>%
   dplyr::rename(mod_temp_w2 = temp) %>%
   merge(temp_w3, by=c("DateTime","Depth")) %>%
   dplyr::rename(mod_temp_w3 = temp)
 temp$DateTime <- as.Date(temp$DateTime, format="%Y-%m-%d")
 
+temp <- temp %>%
+  mutate(error_w1 = abs(mod_temp_w1 - obtemp),
+         error_w2 = abs(mod_temp_w2 - obtemp),
+         error_w3 = abs(mod_temp_w3 - obtemp),
+         obs_dates = 0)
+
 #Error calculation w1
 for (i in 1:nrow(temp)) {
   temp$MEFF_1_w1[i]<- ((temp$modtemp[i]- temp$obtemp[i])^2)
@@ -187,53 +204,62 @@ for (i in 1:nrow(temp)) {
 #error[error$metric=="temp" & error$calibration=="PEST_exm_w3", "Calibration.deepm2"] <- MEFF_temp_w3
 #write.csv(error, 'observations/error_stats.csv')
 
+#Temperature error plots 
 title<- c("Epilimnion", "Metalimnion", "Hypolimnion")
 depth<- c("1 m depth", "4 m depth", "9 m depth")
+plot_title <- c("Eplimnion (1 m depth)", "Metalimnion (4 m depth)", "Hypolimnion (9 m depth)")
 
-#Plots
-plot1 <- vector('list')
+plot_depths <- c(1, 4, 9)
+plot_error_temp <- vector('list')
 
 for(i in 1:length(plot_depths)){
 
 
-  plot1[[i]] <- temp %>%
+  plot_error_temp[[i]] <- temp %>%
     filter(Depth == plot_depths[i]) %>%
-    ggplot2::ggplot(ggplot2::aes(x = DateTime, y = obtemp, colour="Obs. temp")) +
-    ggplot2::geom_point(pch=10)+
-    geom_line(data=filter(temp_w1, Depth==plot_depths[i]), aes(x=DateTime, y=temp, colour="Model w1 temp"), lty=1, size=0.5)+
-    geom_line(data=filter(temp_w2, Depth==plot_depths[i]), aes(x=DateTime, y=temp, colour="Model w2 temp"), lty=2, size=0.5)+
-    geom_line(data=filter(temp_w3, Depth==plot_depths[i]), aes(x=DateTime, y=temp, colour="Model w3 temp"), lty=3, size=0.8)+
+    ggplot2::ggplot(ggplot2::aes(x = DateTime, y = obs_dates, colour="Obs. dates")) +
+    ggplot2::geom_point(pch=4)+
+    geom_line(data=filter(temp, Depth==plot_depths[i]), aes(x=DateTime, y=error_w1, colour="Model w1 temp"), lty=2, linewidth=0.5)+
+    geom_line(data=filter(temp, Depth==plot_depths[i]), aes(x=DateTime, y=error_w2, colour="Model w2 temp"), lty=1, linewidth=0.5)+
+    geom_line(data=filter(temp, Depth==plot_depths[i]), aes(x=DateTime, y=error_w3, colour="Model w3 temp"), lty=3, linewidth=0.5)+
+    ggplot2::ggtitle(" ", subtitle=" ")+
     xlab("Date")+
-    ylab("Temperature (°C)")+
-    ylim(c(0, 30))+
-    ggplot2::ggtitle(title[i], subtitle=depth[i])+
-    ggplot2::scale_colour_manual(name="Legend", values=c("Obs. temp"="black", "Model w1 temp"="#619CFF", "Model w2 temp"="#619CFF", "Model w3 temp"="#619CFF"), guide=guide_legend(override.aes=list(linetype=c(NA, 1, 2, 3), shape=c(10, NA, NA, NA))))+
+    ylab("Absolute model error (\u00B0C)")+
+    ggplot2::ggtitle(plot_title[i], subtitle= "Temperature")+
+    ylim(c(0, 5))+
+    ggplot2::scale_colour_manual(name="Legend", values=c("Obs. dates" = "black", "Model w1 temp"="#619CFF", "Model w2 temp"="#619CFF", "Model w3 temp"="#619CFF"),
+    labels=c("Obs. dates", "Model w1 temp", "Model w2 temp", "Model w3 temp"), guide=guide_legend(override.aes=list(linetype=c(NA, 2, 1, 3), shape=c(4, NA, NA, NA), colour=c("black", "#619CFF", "#619CFF", "#619CFF"))))+
     ggplot2::theme_light() +
     ggplot2::theme(
-      plot.title = ggplot2::element_text(face= "bold", size = 12),
-      axis.title.y = ggplot2::element_text(face = "bold", size= 10),
-      axis.title.x = ggplot2::element_text(face = "bold", size= 10),
+      plot.title = ggplot2::element_text(face= "bold", size = 10),
+      axis.title.y = ggplot2::element_text(face = "bold", size= 8),
+      axis.title.x = ggplot2::element_text(face = "bold", size= 8),
       legend.text = ggplot2::element_text(size= 8),
       axis.text.x = ggplot2::element_text(size=9),
       axis.text.y = ggplot2::element_text(size=9),
       plot.subtitle=element_text(size=10),
       legend.title = ggplot2::element_blank()
+
     )
 
 
 }
-plot1
 
-combinedPlot1 <- patchwork::wrap_plots(plot1, ncol = 1) +
+plot_error_temp
+
+combinedPlot_temp <- patchwork::wrap_plots(plot_error_temp, ncol = 1) +
   patchwork::plot_layout(guides = "collect") &
   theme(legend.position = "bottom")
-combinedPlot1
 
-final <- ggarrange(combinedPlot1, combinedPlot, ncol=2, nrow=1)
+combinedPlot_temp
+
+final <- ggarrange(combinedPlot_temp, combinedPlot_DO, ncol=2, nrow=1)
 final
-ggsave("Results/Figure5.png",
+ggsave("/Results/Figure5.png",
        plot = final,
        width = 246.2, 
        height = 210, 
        units = "mm")
 
+
+