successful knit of vignettes/example.Rmd

DESCRIPTION

@@ -7,6 +7,6 @@ Description: Implements the SOFUN modelling framework for simulating terrestrial
   calibration and benchmarking.
 License: GPL-3
 LazyData: true
-Depends: dplyr, purrr, lubridate, tidyr, raster, lubridate, stringi, sp, gensa, stringr, rland, readr
+Depends: dplyr, purrr, lubridate, tidyr, raster, lubridate, stringi, sp, GenSA, stringr, rlang, readr
 RoxygenNote: 7.1.0
 VignetteBuilder: knitr

NAMESPACE

@@ -2,17 +2,11 @@
 
 export(calib_sofun)
 export(collect_drivers_sofun)
-export(download_from_remote)
+export(collect_obs_eval)
 export(eval_sofun)
 export(filter_days)
-export(get_obs_bysite_fluxnet2015)
-export(get_obs_bysite_gpp_gepisat)
-export(get_obs_calib)
-export(get_obs_eval)
-export(get_obs_eval2)
 export(get_stats)
 export(init_dates_dataframe)
-export(load_dependencies_rsofun)
 export(oob_calib_eval_sofun)
 export(prepare_setup_sofun)
 export(run_sofun_f_bysite)
@@ -30,6 +24,7 @@ import(readr)
 import(rlang)
 import(stringr)
 import(tidyr)
+importFrom(GenSA,GenSA)
 importFrom(grDevices,boxplot.stats)
 importFrom(grDevices,col2rgb)
 importFrom(grDevices,colorRampPalette)

R/collect_drivers_sofun.R

@@ -2,21 +2,22 @@
 #'
 #' Collect all drivers for site-level simulations into a nested data frame with one row for each site
 #'
-#' @param settings A named list containing the simulation settings (see vignette_rsofun.pdf for more information and examples).
+#' @param siteinfo A data frame containing site meta info (rows for sites). Required columns are: \code{"sitename", "date_start", "date_end", "lon", "lat", "elv"}.
 #' @param meteo A data nested data frame with rows for each site and meteo forcing data time series nested inside a column named \code{"data"}
 #' @param fapar A data nested data frame with rows for each site and fAPAR forcing data time series nested inside a column named \code{"data"}
+#' @param co2 A data nested data frame with rows for each site and CO2 forcing data time series nested inside a column named \code{"data"}
 #' @param df_soiltexture
 #'
 #' @return
 #' @export
 #'
 #' @examples mod <- collect_drivers_sofun( settings = settings, forcing, df_soiltexture )
 #' 
-collect_drivers_sofun <- function( settings, meteo, fapar, df_soiltexture ){
+collect_drivers_sofun <- function( siteinfo, meteo, fapar, co2, df_soiltexture ){
 
   ## create mega-df containing all forcing data and parameters that vary by site (not model parameters!)
-  names_metainfo <- names(settings)[-which(names(settings) %in% c("sitename", "params_siml"))]
-  df_mega <- settings %>% 
+  names_metainfo <- names(siteinfo)[-which(names(siteinfo) %in% c("sitename", "params_siml"))]
+  df_mega <- siteinfo %>% 
     tidyr::nest(siteinfo = names_metainfo) %>% 
     left_join(
       meteo %>% 
@@ -28,9 +29,36 @@ collect_drivers_sofun <- function( settings, meteo, fapar, df_soiltexture ){
         rename(fapar = data),
       by = "sitename"
     ) %>% 
+    left_join(
+      co2 %>% 
+        rename(co2 = data),
+      by = "sitename"
+    ) %>% 
     mutate(df_soiltexture = purrr::map(as.list(seq(nrow(.))), ~return(df_soiltexture)))
 
+  ## use only interpolated fapar and combine meteo data and fapar into a single nested column 'forcing'
+  df_mega <- df_mega %>% 
+    mutate(fapar = purrr::map(fapar, ~dplyr::select(., date, fapar = modisvar_interpol))) %>% 
+    mutate(co2   = purrr::map(co2  , ~dplyr::select(., date, co2))) %>% 
+    mutate(forcing = purrr::map2(meteo, fapar, ~left_join( .x, .y, by = "date"))) %>% 
+    mutate(forcing = purrr::map2(forcing, co2, ~left_join( .x, .y, by = "date"))) %>% 
+    dplyr::select(-meteo, -fapar, -co2)
+
+  ## interpolate to fill gaps in forcing time series
+  myapprox <- function(vec){
+    approx(vec, xout = 1:length(vec))$y
+  }
+  fill_na_forcing <- function(df){
+    vars <- names(df)[-which(names(df)=="date")]
+    df %>% 
+      mutate_at(vars, myapprox)
+  }
+  df_mega <- df_mega %>% 
+    mutate(forcing = purrr::map(forcing, ~fill_na_forcing(.)))
+
   return(df_mega)
 }
+
+
 
 

R/collect_obs_eval.R

@@ -0,0 +1,80 @@
+#' Get observational data for model evaluation
+#'
+#' Gets observational data for model evaluation, given the benchmarking variable, and data source used 
+#' for the evaluation. This information is specified in the evaluation settings (argument \code{settings}).
+#'
+#' @param siteinfo A data frame containing site meta info (rows for sites). Required columns are: \code{"sitename", "date_start", "date_end", "lon", "lat", "elv"}.
+#' @param settings A list specifying evaluation settings (see vignette eval_sofun.pdf for more information and examples)
+#' @param adf A nested data frame with rows for sites and time series (annual) for each site nested inside the column \code{data}. 
+#' @param mdf A nested data frame with rows for sites and time series (monthly) for each site nested inside the column \code{data}. 
+#' @param ddf A nested data frame with rows for sites and time series (daily) for each site nested inside the column \code{data}. 
+#' @param hrdf A nested data frame with rows for sites and time series (hourly) for each site nested inside the column \code{data}. 
+#' @param hhdf A nested data frame with rows for sites and time series (half-hourly) for each site nested inside the column \code{data}. 
+#'
+#' @return A list containing nested data frames of observed values aggregated to several temporal scales (ddf for daily, xdf for X-daily 
+#' (determined by argument \code{agg}), mdf for monthly, adf for annual).
+#' @export
+#'
+#' @examples \dontrun{obs <- collect_obs_eval( settings = NULL, adf = NULL, mdf = NULL, ddf = NULL, hrdf = NULL, hhdf = NULL, agg = 8 )}
+#' 
+collect_obs_eval <- function( siteinfo, settings, adf = NULL, mdf = NULL, ddf = NULL, hrdf = NULL, hhdf = NULL ){
+
+  if (is.null(adf)){
+    rlang::abort("collect_obs_eval(): object provided by argument 'adf' is missing.")
+  }
+
+  if (is.null(ddf)){
+    rlang::abort("collect_obs_eval(): object provided by argument 'ddf' is missing.")
+  }
+
+  if (is.null(mdf)){	
+    rlang::abort("collect_obs_eval(): object provided by argument 'mdf' is missing.")
+  }
+
+  ##------------------------------------------------------------
+  ## Aggregate daily data to multi-day periods
+  ## periods should start with the 1st of January each year, otherwise can't compute mean seasonal cycle
+  ##------------------------------------------------------------
+  ## Generate vector of starting dates of X-day periods, making sure the 1st of Jan is always the start of a new period
+  calc_xdf <- function(ddf, agg){
+
+    firstyear <- ddf %>% mutate(year = lubridate::year(date)) %>% pull(year) %>% min()
+    lastyear  <- ddf %>% mutate(year = lubridate::year(date)) %>% pull(year) %>% max()
+    listyears <- seq( ymd(paste0(firstyear, "-01-01")), ymd(paste0(lastyear, "-01-01")), by = "year" )    
+    breaks <- purrr::map( as.list(listyears), ~seq( from=., by=paste0( agg, " days"), length.out = ceiling(365 / agg)) ) %>% Reduce(c,.)
+
+    ## take mean across periods
+    evalvars <- names(settings$benchmark)
+    ddf %>% 
+      mutate( inbin = cut( date, breaks = breaks, right = FALSE ) ) %>%
+      tidyr::drop_na(inbin) %>% 
+      group_by( inbin ) %>%
+      summarise_at( vars(one_of(evalvars)), mean, na.rm=TRUE)
+
+  }
+  xdf <- ddf %>%
+    mutate(data = purrr::map(data, ~calc_xdf(., agg = settings$agg)))
+
+  ## complement info required for evaluation
+  ddf <- siteinfo %>% 
+    dplyr::select(sitename, lon, lat, elv, classid, c4, whc, koeppen_code, igbp_land_use, plant_functional_type) %>% 
+    right_join(ddf, by = "sitename")
+
+  xdf <- siteinfo %>% 
+    dplyr::select(sitename, lon, lat, elv, classid, c4, whc, koeppen_code, igbp_land_use, plant_functional_type) %>% 
+    right_join(xdf, by = "sitename")
+
+  mdf <- siteinfo %>% 
+    dplyr::select(sitename, lon, lat, elv, classid, c4, whc, koeppen_code, igbp_land_use, plant_functional_type) %>% 
+    right_join(mdf, by = "sitename")
+
+  adf <- siteinfo %>% 
+    dplyr::select(sitename, lon, lat, elv, classid, c4, whc, koeppen_code, igbp_land_use, plant_functional_type) %>% 
+    right_join(adf, by = "sitename")
+
+  ## get breaks
+  breaks <- xdf %>% tidyr::unnest(data) %>% pull(inbin) %>% as.character() %>% lubridate::ymd()
+
+	return( list( ddf = ddf, xdf = xdf, mdf = mdf, adf = adf, breaks = breaks ) )
+
+}

R/eval_sofun.R

@@ -25,11 +25,22 @@
 #' 
 eval_sofun <- function(mod, settings_eval, settings_sims, obs_eval = NA, overwrite = TRUE, doplot = FALSE, light = FALSE){
 
+  ## make model output a long flat table
+  mod <- mod %>% 
+    dplyr::rename(id = sitename) %>% 
+    tidyr::unnest(out_sofun)
+
   ## Evaluate daily variables
-  out <- lapply( as.list(names(settings_eval$benchmark)), 
-                 function(x) eval_sofun_byvar(x, dplyr::select(mod, sitename, date, mod = eval(x)), settings_eval, settings_sims, obs_eval = obs_eval, overwrite = TRUE, doplot = FALSE, light = light)
-                 ) %>%
-          setNames(names(settings_eval$benchmark))
+  out <- purrr::map(
+    as.list(names(settings_eval$benchmark)),
+    ~eval_sofun_byvar(., dplyr::select(mod, sitename, date, mod = {{.}}), settings_eval, settings_sims, obs_eval = obs_eval, overwrite = TRUE, doplot = FALSE, light = light)
+    ) %>% 
+    setNames(names(settings_eval$benchmark))
+
+  # out <- lapply( as.list(names(settings_eval$benchmark)), 
+  #                function(x) eval_sofun_byvar(x, dplyr::select(mod, sitename, date, mod = eval(x)), settings_eval, settings_sims, obs_eval = obs_eval, overwrite = TRUE, doplot = FALSE, light = light)
+  #                ) %>%
+  #         setNames(names(settings_eval$benchmark))
 
   return(out)
 }
@@ -82,18 +93,18 @@ eval_sofun_byvar <- function(varnam, ddf_mod, settings_eval, settings_sims, obs_
     ##------------------------------------------------------------
     ## Get observations for evaluation
     ##------------------------------------------------------------
-    if (identical(obs_eval, NA)) obs_eval <- get_obs_eval( settings_eval = settings_eval, settings_sims = settings_sims, overwrite = overwrite )
+    if (identical(obs_eval, NA)) rlang::abort("eval_sofun_byvar(): Object provided by argument 'eval_sofun_byvar' could not be identified.")
 
     ## detach
     if (varnam=="aet"){
       varnam_obs <- "latenth"
     } else {
       varnam_obs <- varnam
     }
-    adf <- obs_eval$adf %>% dplyr::select(sitename, date,  obs = {{varnam_obs}} ) # eval(varnam_obs))
-    mdf <- obs_eval$mdf %>% dplyr::select(sitename, date,  obs = {{varnam_obs}})
-    ddf <- obs_eval$ddf %>% dplyr::select(sitename, date,  obs = {{varnam_obs}})
-    xdf <- obs_eval$xdf %>% dplyr::select(sitename, inbin, obs = {{varnam_obs}})
+    adf <- obs_eval$adf %>% tidyr::unnest(data) %>% dplyr::select(sitename, date,  obs = {{varnam_obs}})
+    mdf <- obs_eval$mdf %>% tidyr::unnest(data) %>% dplyr::select(sitename, date,  obs = {{varnam_obs}})
+    ddf <- obs_eval$ddf %>% tidyr::unnest(data) %>% dplyr::select(sitename, date,  obs = {{varnam_obs}}, lat, koeppen_code)
+    xdf <- obs_eval$xdf %>% tidyr::unnest(data) %>% dplyr::select(sitename, inbin, obs = {{varnam_obs}})
 
     # obs_eval$adf <- NULL
     # obs_eval$mdf <- NULL
@@ -130,7 +141,7 @@ eval_sofun_byvar <- function(varnam, ddf_mod, settings_eval, settings_sims, obs_
     ## mean across multi-day period
     xdf <- ddf_mod %>% 
       # mutate( year = year(date), week = week(date) ) %>%
-      mutate( year = year(date), inbin = cut( date, breaks = obs_eval$breaks_xdf, right = FALSE ) ) %>%
+      mutate( year = year(date), inbin = cut( date, breaks = obs_eval$breaks, right = FALSE ) ) %>%
       tidyr::drop_na() %>% 
       group_by( sitename, inbin ) %>%
       summarise( mod_mean = mean( mod, na.rm = TRUE ), mod_min = min( mod, na.rm = TRUE ), mod_max = max( mod, na.rm = TRUE ), n_mod = sum(!is.na(mod)) ) %>%
@@ -305,7 +316,7 @@ eval_sofun_byvar <- function(varnam, ddf_mod, settings_eval, settings_sims, obs_
       if (!light){
         print("Evaluate mean seasonal cycle by climate zones...")
         meandoydf_byclim <- ddf %>% mutate( doy = yday(date) ) %>%
-          left_join( dplyr::select( metainfo_Tier1_sites_kgclimate_fluxnet2015, sitename, lat, koeppen_code ), by = "sitename" ) %>%   # 'metainfo_Tier1_sites_kgclimate_fluxnet2015' is lazy-loaded with library(rsofun)
+          # left_join( dplyr::select( metainfo_Tier1_sites_kgclimate_fluxnet2015, sitename, lat, koeppen_code ), by = "sitename" ) %>%   # 'metainfo_Tier1_sites_kgclimate_fluxnet2015' is lazy-loaded with library(rsofun)
           mutate( hemisphere = ifelse( lat>0, "north", "south" ) ) %>%
           # mutate( bias = mod - obs ) %>% 
           dplyr::select( -lat ) %>%