Merge pull request #207 from lauramarques/master

Corrected forcing data script and added function to run biomee multiple times

R/cost_likelihood_biomee.R

@@ -10,7 +10,7 @@
 #' Make sure that
 #' the order of the error terms in \code{par} coincides with the order provided in
 #' the \code{targets} argument.
-#' @param obs A nested data frame of observations, following the structure of \code{biomee_validation_2},
+#' @param obs A nested data frame of observations, following the structure of \code{biomee_validation},
 #' for example.
 #' @param drivers A nested data frame of driver data, for example \code{biomee_gs_leuning_drivers}.
 #' @param targets A character vector indicating the target variables for which the
@@ -36,7 +36,7 @@
 #' cost_likelihood_biomee(
 #'  par = c(3.5, 3.5, 1, 1,    # model params
 #'          0.5),              # err_GPP
-#'  obs = biomee_validation_2,
+#'  obs = biomee_validation,
 #'  drivers = biomee_gs_leuning_drivers,
 #'  targets = c("GPP")
 #' )

R/cost_rmse_biomee.R

@@ -9,7 +9,7 @@
 #' 
 #' @param par A vector containing parameter values for \code{'phiRL',
 #' 'LAI_light', 'tf_base', 'par_mort'} in that order.
-#' @param obs A nested data frame of observations, following the structure of \code{biomee_validation_2},
+#' @param obs A nested data frame of observations, following the structure of \code{biomee_validation},
 #' for example.
 #' @param drivers A nested data frame of driver data, for example \code{biomee_gs_leuning_drivers}.
 #' 
@@ -30,7 +30,7 @@
 #' # and example data
 #' cost_rmse_biomee(
 #'  par = c(3.5, 3.5, 1, 1),
-#'  obs = biomee_validation_2,
+#'  obs = biomee_validation,
 #'  drivers = biomee_gs_leuning_drivers
 #' )
 #' }

R/data.R

@@ -242,7 +242,7 @@
 #' }
 "biomee_gs_leuning_drivers"
 
-#' rsofun BiomeE GPP validation data
+#' rsofun BiomeE targets validation data
 #'
 #' Small tests dataset to validate 
 #' calibration routines
@@ -254,19 +254,6 @@
 #' }
 "biomee_validation"
 
-#' rsofun BiomeE GPP validation data
-#'
-#' Small tests dataset to validate 
-#' calibration routines
-#'
-#' @format A tibble of driver data:
-#' \describe{
-#'   \item{sitename}{site name}
-#'   \item{data}{validation dta}
-#' }
-"biomee_validation_2"
-
-
 #' rsofun BiomeE GPP validation data
 #'
 #' Small tests dataset to validate 
@@ -279,7 +266,6 @@
 #' }
 "biomee_p_model_output"
 
-
 #' rsofun BiomeE GPP validation data
 #'
 #' Small tests dataset to validate 

R/run_biomee_f_bysite.R

@@ -19,6 +19,10 @@
 #'     LAImax according to mineral N in soil.}
 #'   \item{do_closedN_run}{A logical value indicating whether doing N closed 
 #'     runs to recover N balance.}
+#'   \item{do_reset_veg}{A logical value indicating whether reseting vegetation 
+#'     for disturbance runs.}
+#'   \item{dist_frequency}{Value indicating the frequency of the disturbance event (in years) 
+#'     (e.g. 100 indicates a disturbance event occurs every 100 years, i.e., at a rate of 0.01)}
 #'   \item{code_method_photosynth}{String specifying the method of photosynthesis 
 #'     used in the model, either "pmodel" or "gs_leuning".}
 #'   \item{code_method_mortality}{String indicating the type of mortality in the 
@@ -511,6 +515,8 @@ run_biomee_f_bysite <- function(
       do_U_shaped_mortality = as.logical(params_siml$do_U_shaped_mortality),
       update_annualLAImax   = as.logical(params_siml$update_annualLAImax),
       do_closedN_run        = as.logical(params_siml$do_closedN_run),
+      do_reset_veg          = as.logical(params_siml$do_reset_veg),
+      dist_frequency        = as.integer(params_siml$dist_frequency),
       code_method_photosynth= as.integer(code_method_photosynth),
       code_method_mortality = as.integer(code_method_mortality),
 
@@ -557,7 +563,8 @@ run_biomee_f_bysite <- function(
       n                = as.integer(nrow(forcing)), # n here is for hourly (forcing is hourly), add n for daily and annual outputs
       n_daily          = as.integer(n_daily), # n here is for hourly (forcing is hourly), add n for daily and annual outputs
       n_annual         = as.integer(runyears), # n here is for hourly (forcing is hourly), add n for daily and annual outputs
-      n_annual_cohorts = as.integer(params_siml$nyeartrend), # n here is for hourly (forcing is hourly), add n for daily and annual outputs
+      #n_annual_cohorts = as.integer(params_siml$nyeartrend), # n here is for hourly (forcing is hourly), add n for daily and annual outputs
+      n_annual_cohorts = as.integer(runyears), # n here is for hourly (forcing is hourly), add n for daily and annual outputs
       forcing          = as.matrix(forcing)
       )
 
@@ -693,13 +700,13 @@ run_biomee_f_bysite <- function(
     annual_values <- c(
       "year","cID",
       "PFT","layer","density",
-      "f_layer","dDBH","dbh",
-      "height","age","Acrown",
-      "wood","nsc","NSN","NPPtr",
-      "seed","NPPL","NPPR","NPPW",
-      "GPP_yr","NPP_yr","Rauto",
-      "N_uptk","N_fix","maxLAI",
-      "Volume","n_deadtrees",
+      "flayer","DBH","dDBH","height",
+      "age","BA","dBA","Acrown","Aleaf",
+      "nsc","seedC","leafC","rootC",
+      "sapwC","woodC","nsn","treeG",
+      "fseed","fleaf","froot","fwood",
+      "GPP","NPP","Rauto",
+      "Nupt","Nfix","n_deadtrees",
       "c_deadtrees","deathrate"
     )
 

analysis/biome_segmentation_fault.R

@@ -0,0 +1,65 @@
+
+library(rsofun)
+library(dplyr)
+library(ggplot2)
+library(tibble)
+library(patchwork)
+library(devtools)
+
+load("data/biomee_gs_leuning_drivers.rda")
+biomee_gs_leuning_drivers$params_siml
+head(biomee_gs_leuning_drivers$forcing)
+
+load("data/biomee_p_model_drivers.rda")
+biomee_p_model_drivers$params_siml
+head(biomee_p_model_drivers$forcing)
+
+set.seed(2023)
+
+# run the model
+# build()
+# install()
+out <- runread_biomee_f(
+  biomee_p_model_drivers,
+  makecheck = TRUE,
+  parallel = FALSE
+)
+
+# split out the annual data
+biomee_pmodel_output <- out$data[[1]]$output_annual_tile
+
+# plot
+cowplot::plot_grid(
+  biomee_pmodel_output |>
+    ggplot() +
+    geom_line(aes(x = year, y = GPP)) +
+    theme_classic()+labs(x = "Year", y = "GPP"),
+  biomee_pmodel_output |>
+    ggplot() +
+    geom_line(aes(x = year, y = plantC)) +
+    theme_classic()+labs(x = "Year", y = "plantC")
+)
+
+# function to run biomee multiple times and detect the segmentation fault error
+biomeextimes <- function(drivers,n){
+  results <-  list()
+  for(i in 1:n) {
+    set.seed(2023)
+    out <- runread_biomee_f(
+    drivers,
+    makecheck = TRUE,
+    parallel = FALSE
+  )$data[[1]]$output_annual_tile['plantC'] %>% 
+    slice(tail(row_number(), 1)) 
+  if(is.na(out)|out==0) message("Error: Simulation failed")
+  results <- append(results, out)
+  }
+  save(results, file= "data/results.rda")
+  return(results)
+}
+
+# run the model n times 
+n = 500
+biomeextimes(biomee_p_model_drivers,n) 
+load("data/results.rda")
+results