From 2b6df9fec9c4ba85b098301970f539e483de5894 Mon Sep 17 00:00:00 2001
From: Beni Stocker <benjamin.stocker@gmail.com>
Date: Tue, 16 Jan 2024 15:06:03 +0100
Subject: [PATCH] Added plotting output for rsofun manuscript as comments;
 removed plotting of prior flux estimates.

---
 analysis/01-sensitivity-analysis.R   |  7 +++++--
 analysis/02-bayesian-calibration.R   |  9 ++++++---
 analysis/03-uncertainty-estimation.R | 14 ++++++--------
 3 files changed, 17 insertions(+), 13 deletions(-)

diff --git a/analysis/01-sensitivity-analysis.R b/analysis/01-sensitivity-analysis.R
index da090029..822e30f9 100644
--- a/analysis/01-sensitivity-analysis.R
+++ b/analysis/01-sensitivity-analysis.R
@@ -97,7 +97,7 @@ morrisOut.df <- data.frame(
   arrange( mu.star )
 
 # Create barplot to show sensitivity analysis output
-morrisOut.df |>
+gg <- morrisOut.df |>
   tidyr::pivot_longer( -parameter, names_to = "variable", values_to = "value") |>
   ggplot(aes(
     reorder(parameter, value),
@@ -116,4 +116,7 @@ morrisOut.df |>
     axis.title = element_blank(),
     legend.position = c(0.9, 0.1), legend.justification = c(0.95, 0.05)
   ) +
-  coord_flip()    # make horizontal
\ No newline at end of file
+  coord_flip()    # make horizontal
+
+# ggsave("./analysis/paper_results_files/morris.pdf", plot = gg, width = 5, height = 3)
+# ggsave("./analysis/paper_results_files/morris.png", plot = gg, width = 5, height = 3)
diff --git a/analysis/02-bayesian-calibration.R b/analysis/02-bayesian-calibration.R
index 0a45f8b5..a57f7243 100644
--- a/analysis/02-bayesian-calibration.R
+++ b/analysis/02-bayesian-calibration.R
@@ -56,7 +56,7 @@ plot_prior_posterior_density <- function(
   df_plot$distrib <- as.factor(df_plot$distrib)
   
   # Plot with facet wrap
-  df_plot |>
+  gg <- df_plot |>
     tidyr::gather(variable, value, kphio:err_gpp) |>
     ggplot(
       aes(x = value, fill = distrib)
@@ -66,9 +66,9 @@ plot_prior_posterior_density <- function(
     facet_wrap( ~ variable , nrow = 2, scales = "free") +
     theme(legend.position = "bottom",
           axis.title.x = element_text("")) +
-    ggtitle("Marginal parameter uncertainty") +
     scale_fill_manual(values = c("#29a274ff", t_col("#777055ff"))) # GECO colors
   
+  return(gg)
 }
 
 # Set random seed for reproducibility
@@ -114,4 +114,7 @@ par_calib <- calib_sofun(
 toc() # Stop measuring time
 
 # Plot prior and posterior distributions
-plot_prior_posterior_density(par_calib$mod)
+gg <- plot_prior_posterior_density(par_calib$mod)
+
+# ggsave("./analysis/paper_results_files/prior_posterior.pdf", plot = gg, width = 6, height = 5)
+# ggsave("./analysis/paper_results_files/prior_posterior.png", plot = gg, width = 6, height = 5)
diff --git a/analysis/03-uncertainty-estimation.R b/analysis/03-uncertainty-estimation.R
index 3c73fd86..9491a792 100644
--- a/analysis/03-uncertainty-estimation.R
+++ b/analysis/03-uncertainty-estimation.R
@@ -1,6 +1,7 @@
 # Script getting the uncertainty estimation, using the
 # output from 02-bayesian-calibration.R
 getwd()
+
 # Load libraries
 library(rsofun)
 library(dplyr)
@@ -79,12 +80,6 @@ plot_gpp_error <- ggplot(data = pmodel_runs |>
                              # Merge GPP validation data (first year)
                              p_model_validation$data[[1]][1:365, ] |>
                                dplyr::rename(gpp_obs = gpp),
-                             by = "date") |>
-                           dplyr::left_join(
-                             # Merge GPP before calibration
-                             output$data[[1]][1:365, ] |>
-                               dplyr::select(date, gpp) |>
-                               dplyr::rename(gpp_no_calib = gpp),
                              by = "date")
 ) +             # Plot only first year
   geom_ribbon(
@@ -121,7 +116,7 @@ plot_gpp_error <- ggplot(data = pmodel_runs |>
   )
 
 # Include observations in the plot
-plot_gpp_error +  
+plot_gpp_error <- plot_gpp_error +  
   scale_color_manual(name = "",
                      breaks = c("Observations",
                                 "Predictions",
@@ -133,4 +128,7 @@ plot_gpp_error +
                     breaks = c("Model uncertainty",
                                "Parameter uncertainty"),
                     values = c(t_col("#E69F00", 60),
-                               t_col("#009E73", 40)))
\ No newline at end of file
+                               t_col("#009E73", 40)))
+
+# ggsave("./analysis/paper_results_files/gpp_predictions_observations.pdf", plot = plot_gpp_error, width = 6, height = 5)
+# ggsave("./analysis/paper_results_files/gpp_predictions_observations.png", plot = plot_gpp_error, width = 6, height = 5)