made removal of negative obs. GPP an option (was forced before), upda…

…ted test_implementations.Rmd.

R/calib_sofun.R

@@ -552,6 +552,7 @@ get_obs_bysite <- function( sitename, settings_calib, settings_sims, settings_in
         getvars = getvars, 
         getswc = FALSE,
         threshold_GPP = settings_calib$threshold_GPP, 
+        remove_neg = FALSE,
         verbose = TRUE
         ) %>% 
         dplyr::mutate( gpp_obs = case_when("NT" %in% datasource & !("DT" %in% datasource) ~ GPP_NT_VUT_REF,
@@ -634,6 +635,7 @@ get_obs_bysite <- function( sitename, settings_calib, settings_sims, settings_in
         getvars = c("LE_F_MDS", "LE_RANDUNC"), 
         getswc = FALSE,
         threshold_LE = 0.6, 
+        remove_neg = FALSE,
         verbose = TRUE
         ) %>% 
         dplyr::right_join( ddf, by = "date" )

R/get_obs_bysite_fluxnet2015.R

@@ -95,6 +95,8 @@
 #' meaning no data is excluded.
 #' @param return_qc A logical specifying whether quality control variables
 #' should be returned.
+#' @param remove_neg A logical specifying whether negative GPP values are to
+#' be removed (replaces with NA).
 #'
 #' @return A data frame (tibble) containing cleaned observational data,
 #' named and in units corresponding to rsofun standard.
@@ -105,7 +107,8 @@
 get_obs_bysite_fluxnet2015 <- function( sitename, path_fluxnet2015, path_fluxnet2015_hh=NULL,
   timescale, getvars, getswc=TRUE,
   threshold_GPP=0.0, threshold_LE=0.0, threshold_H=0.0, threshold_SWC=0.0, 
-  threshold_WS=0.0, threshold_USTAR=0.0, threshold_T=0.0, threshold_NETRAD=0.0, return_qc=FALSE, verbose=TRUE ){
+  threshold_WS=0.0, threshold_USTAR=0.0, threshold_T=0.0, threshold_NETRAD=0.0, return_qc=FALSE, 
+  remove_neg = FALSE, verbose=TRUE ){
 
   if (verbose) print(paste("Getting FLUXNET data for", sitename, "..."))
 
@@ -390,7 +393,7 @@ get_obs_bysite_fluxnet2015 <- function( sitename, path_fluxnet2015, path_fluxnet
     # df$GPP_DT_VUT_REF <- out_clean$gpp_dt
 
     df <- df %>% 
-      clean_fluxnet_gpp(threshold = threshold_GPP)
+      clean_fluxnet_gpp(threshold = threshold_GPP, remove_neg = remove_neg)
   }
 
   ## clean energy data (sensible and latent heat flux) data - often has spuriously equal values
@@ -692,7 +695,7 @@ convert_energy_fluxnet2015 <- function( le ){
   return(le_converted)
 }
 
-clean_fluxnet_gpp <- function(df, nam_gpp_nt, nam_gpp_dt, nam_nt_qc, nam_dt_qc, threshold){
+clean_fluxnet_gpp <- function(df, nam_gpp_nt, nam_gpp_dt, nam_nt_qc, nam_dt_qc, threshold, remove_neg = FALSE){
   ##--------------------------------------------------------------------
   ## Cleans daily data using criteria 1-4 as documented in Tramontana et al., 2016
   ## gpp_nt: based on nighttime flux decomposition ("NT")
@@ -728,21 +731,26 @@ clean_fluxnet_gpp <- function(df, nam_gpp_nt, nam_gpp_dt, nam_nt_qc, nam_dt_qc,
   q975 <- quantile( df$res, probs = 0.975, na.rm=TRUE )
 
   df <- df %>% 
+
+    ## remove data outside the quartiles of the residuals between the DT and NT estimates
     mutate(GPP_NT_VUT_REF = replace_with_na_res(GPP_NT_VUT_REF, res, q025, q975),
            GPP_DT_VUT_REF = replace_with_na_res(GPP_DT_VUT_REF, res, q025, q975)
            ) %>% 
 
-    # xxx don't do this - distorts error distribution
-    # ## remove negative GPP
-    # mutate(GPP_NT_VUT_REF = replace_with_na_neg(GPP_NT_VUT_REF),
-    #        GPP_DT_VUT_REF = replace_with_na_neg(GPP_DT_VUT_REF)
-    #        ) %>% 
-
     ## remove outliers
     mutate(GPP_NT_VUT_REF = remove_outliers(GPP_NT_VUT_REF, coef = 1.5),
            GPP_DT_VUT_REF = remove_outliers(GPP_DT_VUT_REF, coef = 1.5)
            )
 
+  ## remove negative GPP
+  if (remove_neg){
+    df <- df %>% 
+      mutate(GPP_NT_VUT_REF = replace_with_na_neg(GPP_NT_VUT_REF),
+             GPP_DT_VUT_REF = replace_with_na_neg(GPP_DT_VUT_REF)
+      )
+  }
+
+
   return(df)
 }
 

R/get_obs_eval.R

@@ -61,6 +61,7 @@ get_obs_eval <- function( settings_eval, settings_sims, overwrite = TRUE, light
                        getvars = getvars,
                        getswc=!light,
                        threshold_GPP=0,
+                       remove_neg = TRUE,
                        verbose=FALSE
                        ) %>%
                        mutate( year =lubridate::year(date),
@@ -110,6 +111,7 @@ get_obs_eval <- function( settings_eval, settings_sims, overwrite = TRUE, light
                        getvars = getvars,
                        getswc=!light,
                        threshold_GPP=0.5,
+                       remove_neg = TRUE,
                        verbose=FALSE
                        ) %>%
                        mutate( year =lubridate::year(date),
@@ -156,6 +158,7 @@ get_obs_eval <- function( settings_eval, settings_sims, overwrite = TRUE, light
                        getvars = getvars, 
                        getswc=!light,                              
                        threshold_GPP=0.5, 
+                       remove_neg = FALSE,
                        verbose=TRUE
                        ) %>%
                        mutate( year =lubridate::year(date),

vignettes/test_implementations.Rmd

@@ -28,14 +28,15 @@ vignette: >
 Load the `rsofun` package. This contains all the necessary wrapper functions to set up and run SOFUN and read its output. 
 ```{r message=FALSE, echo=FALSE, warning=FALSE}
 library(rsofun)
+library(rpmodel)
 
 ## load all rsofun dependencies
 load_dependencies_rsofun()
 
 ## other crap
 knitr::opts_knit$set( root.dir = rprojroot::find_rstudio_root_file() ) # does not work properly
 # if (is.null(options()$rsofun.dir.sofun)) rlang::abort( "Option rsofun.dir.sofun not set. Do so by `options( list( rsofun.dir.sofun=string_path_where_sofun_is ))`" )
-options( list( rsofun.dir.sofun="~/sofun/trunk/" ))
+options( list( rsofun.dir.sofun="~/sofun/" ))
 ```
 
 ## Run implementations for a temperature range
@@ -58,17 +59,17 @@ fapar <- 1
 
 ## R, Vcmax based on Wang Han's formulation
 pmodel_stdrd_R <- purrr::map( as.list( seq( 0, 35, length.out = 100 ) ), 
-  ~rpmodel( tc = ., vpd = vpd, co2 = co2, elv = elv, kphio = 0.05, fapar = fapar, ppfd = ppfd, method_optci="prentice14", method_jmaxlim = "wang17", do_ftemp_kphio = FALSE ) 
+  ~rpmodel::rpmodel( tc = ., vpd = vpd, co2 = co2, elv = elv, kphio = 0.05, fapar = fapar, ppfd = ppfd, method_optci="prentice14", method_jmaxlim = "wang17", do_ftemp_kphio = FALSE ) 
   )
 
 pmodel_smith <- purrr::map( as.list( seq( 0, 35, length.out = 100 ) ), 
-  ~rpmodel( tc = ., vpd = vpd, co2 = co2, elv = elv, kphio = 0.257, fapar = fapar, ppfd = ppfd, method_optci="prentice14", method_jmaxlim = "smith19", do_ftemp_kphio = FALSE ) 
+  ~rpmodel::rpmodel( tc = ., vpd = vpd, co2 = co2, elv = elv, kphio = 0.257, fapar = fapar, ppfd = ppfd, method_optci="prentice14", method_jmaxlim = "smith19", do_ftemp_kphio = FALSE ) 
   )
 
 ## Fortran, Vcmax basedon Wang Han's formulation
 ## update quantum yield parameter in file
-params_opt <- readr::read_csv( paste0( path.package("rsofun"), "/extdata/params_opt_kphio_soilm_global.csv" ) )
-params_opt$kphio <- 0.05
+#params_opt <- readr::read_csv( paste0( path.package("rsofun"), "/extdata/params_opt_kphio_soilm_global.csv" ) )
+params_opt <- tibble(kphio = 0.05)
 nothing <- update_params( params_opt, options()$rsofun.dir.sofun )
 
 pmodel_fortran <- purrr::map( as.list( seq( 0, 35, length.out = 100 ) ), 
@@ -123,7 +124,7 @@ legend( "topright", c("rsofun standard (R)", "rsofun standard (Fortran)", "Smith
 ```{r, eval=TRUE, message=FALSE, warning=FALSE}
 vcmax_stdrd_R   <- pmodel_stdrd_R   %>% purrr::map_dbl("vcmax")
 vcmax_fortran   <- pmodel_fortran   %>% purrr::map_dbl("vcmax")
-vcmax_smith     <- pmodel_smith     %>% purrr::map_dbl("vcmax_star")
+# vcmax_smith     <- pmodel_smith     %>% purrr::map_dbl("vcmax_star")
 
 plot(  seq( 0, 35, length.out = 100 ), vcmax_stdrd_R, type = "l", xlab = "Temperature (deg C)", ylab = "Vcmax", lwd=6, ylim=c(0,max(vcmax_stdrd_R)) )
 lines( seq( 0, 35, length.out = 100 ), vcmax_fortran, col="green", lwd=3 )