Example
Example
Simulation settings
Manually select some sites from which we’re going to use the data for evaluation and calibration.
mysites <- c("BE-Vie", "DE-Tha", "DK-Sor", "FI-Hyy", "IT-Col", "NL-Loo", "US-MMS", "US-WCr", "US-UMB", "US-Syv", "DE-Hai", "IT-MBo", "US-GLE", "FR-Fon", "NL-Hor", "US-UMd", "AU-Dry", "DE-Obe", "IT-Tor", "US-Wi4")Create a site meta info table that contains all the site-specific information that is used to force site-simulations (e.g. starting year, number of simulations years, elevation, etc.). For FLUXNET2015 data, required meta info is provided by the rsofun package (data frame rsofun::metainfo_Tier1_sites_kgclimate_fluxnet2015).
path_siteinfo <- "~/.siteinfo_example_fortran.csv"
siteinfo <- rsofun::metainfo_Tier1_sites_kgclimate_fluxnet2015 %>%
dplyr::filter(sitename %in% mysites) %>%
write_csv(path = path_siteinfo)Now specify the simulation parameters that are identical for all site-scale simulations.
params_siml <- list(
spinup = TRUE,
spinupyears = 10,
recycle = 1,
soilmstress = FALSE,
tempstress = FALSE,
calc_aet_fapar_vpd = FALSE,
in_ppfd = TRUE,
in_netrad = FALSE,
const_clim_year = -9999,
const_lu_year = -9999,
const_co2_year = -9999,
const_ndep_year = -9999,
const_nfert_year = -9999,
outdt = 1,
ltre = FALSE,
ltne = FALSE,
ltrd = FALSE,
ltnd = FALSE,
lgr3 = TRUE,
lgn3 = FALSE,
lgr4 = FALSE
)Run prepare_setup_sofun() to define the simulation settings that contain all the information specified by the two steps above (meta info, and simulation parameters).
settings_sims <- prepare_setup_sofun(siteinfo = siteinfo, params_siml = params_siml)Define model parameters
First, let’s do it by hand (calibration of parameters is shown later).
params_modl <- list(
kphio = 0.04997714009213085,
soilm_par_a = 1.0,
soilm_par_b = 0.0,
vpdstress_par_a = 0.2,
vpdstress_par_b = 0.2,
vpdstress_par_m = 5
)Define soil parameters
For now, this is implemented as an illustration. Should be made site-specific.
df_soiltexture <- bind_rows(
top = tibble(layer = "top", fsand = 0.4, fclay = 0.3, forg = 0.1, fgravel = 0.1),
bottom = tibble(layer = "bottom", fsand = 0.4, fclay = 0.3, forg = 0.1, fgravel = 0.1)
)Get input
First, define input settings.
settings_input <- list(
data = NA,
temperature = "fluxnet2015",
precipitation = "fluxnet2015",
vpd = "fluxnet2015",
ppfd = "fluxnet2015",
netrad = "fluxnet2015", # c("fluxnet2015", "watch_wfdei"),
patm = "fluxnet2015",
netrad = NA,
cloudcover = "cru",
path_input = "~/sofun_inputs/example/",
path_watch_wfdei = "~/data/watch_wfdei/",
path_cru = "~/data/cru/ts_4.01/",
path_MODIS_FPAR_MCD15A3H = "~/data/fluxnet_subsets/fapar_MODIS_FPAR_MCD15A3H_gee_MCD15A3H_fluxnet2015_gee_subset/",
path_MODIS_EVI_MOD13Q1 = "~/data/fluxnet_subsets/fapar_MODIS_EVI_MOD13Q1_gee_MOD13Q1_fluxnet2015_gee_subset/",
path_co2 = "~/data/co2/cCO2_rcp85_const850-1765.csv",
path_fluxnet2015 = "~/data/FLUXNET-2015_Tier1/20191024/DD/",
path_fluxnet2015_hh = "~/data/FLUXNET-2015_Tier1/20191024/HH/",
get_from_remote = FALSE,
settings_gee = get_settings_gee(
bundle = "fpar",
python_path = "/Users/benjaminstocker/Library/Enthought/Canopy_64bit/User/bin/python",
gee_path = "~/gee_subset/gee_subset/"
),
fapar = "MODIS_FPAR_MCD15A3H",
splined_fapar = TRUE
)Then, get the input data.
## [1] "Creating raster brick from file ~/data/cru/ts_4.01/cru_ts4.01.1901.2016.cld.dat.nc"Run the model
Run the model for all the sites specified in the first step.
df_drivers <- collect_drivers_sofun(
settings = settings_sims,
forcing = ddf_input,
df_soiltexture = df_soiltexture
)
## run for a single site
mod <- run_sofun_f_bysite(
df_drivers$sitename[1],
df_drivers$params_siml[[1]],
df_drivers$siteinfo[[1]],
df_drivers$forcing[[1]],
df_drivers$df_soiltexture[[1]],
params_modl = params_modl,
makecheck = TRUE
)
## Run for the full set of sites
ptm <- proc.time()
df_output <- runread_sofun_f(
df_drivers,
params_modl = params_modl,
makecheck = TRUE,
parallel = FALSE
)
print(ptm)## user system elapsed
## 6.000 0.803 9.246# microbenchmark::microbenchmark(
# runread_sofun_f(
# df_drivers,
# params_modl = params_modl,
# makecheck = TRUE,
# parallel = TRUE,
# ncores = 4
# ),
# runread_sofun_f(
# df_drivers,
# params_modl = params_modl,
# makecheck = TRUE,
# parallel = FALSE
# ),
# times = 5,
# units = 's'
# )
df_output$out_sofun[[1]] %>%
ggplot(aes(x=date, y=gpp)) +
geom_line() +
labs(title = df_output$sitename[[1]], subtitle = "SOFUN output")
Calibrate
Define calibration settings.
## Define calibration settings common for all setups
settings_calib <- list(
method = "gensa",
targetvars = c("gpp"),
timescale = list( gpp = "d" ),
path_fluxnet2015 = "~/data/FLUXNET-2015_Tier1/20191024/DD/",
path_fluxnet2015_hh = "~/data/FLUXNET-2015_Tier1/20191024/HH/",
threshold_GPP = 0.5,
path_gepisat = "~/data/gepisat/v3_fluxnet2015/daily_gpp/",
maxit = 5, # (5 for gensa) (30 for optimr) #
sitenames = mysites,
filter_temp_max = 35.0,
filter_drought = FALSE,
metric = "rmse",
dir_results = "./",
name = "ORG",
par = list( kphio = list( lower=0.03, upper=0.07, init=0.0496 ) ),
datasource = list( gpp = "fluxnet2015_NT" ),
filter_temp_min = NA,
filter_soilm_min = NA
)Get calibration target data.
ddf_obs_calib <- get_obs_calib(
settings_calib = settings_calib,
dplyr::select(df_drivers, sitename, siteinfo) %>% tidyr::unnest(siteinfo),
settings_input
)## [1] "getting obs for BE-Vie"
## [1] "Getting FLUXNET data for BE-Vie ..."
## [1] "getting obs for DE-Tha"
## [1] "Getting FLUXNET data for DE-Tha ..."
## [1] "getting obs for DK-Sor"
## [1] "Getting FLUXNET data for DK-Sor ..."
## [1] "getting obs for FI-Hyy"
## [1] "Getting FLUXNET data for FI-Hyy ..."
## [1] "getting obs for IT-Col"
## [1] "Getting FLUXNET data for IT-Col ..."
## [1] "getting obs for NL-Loo"
## [1] "Getting FLUXNET data for NL-Loo ..."
## [1] "getting obs for US-MMS"
## [1] "Getting FLUXNET data for US-MMS ..."
## [1] "getting obs for US-WCr"
## [1] "Getting FLUXNET data for US-WCr ..."
## [1] "getting obs for US-UMB"
## [1] "Getting FLUXNET data for US-UMB ..."
## [1] "getting obs for US-Syv"
## [1] "Getting FLUXNET data for US-Syv ..."
## [1] "getting obs for DE-Hai"
## [1] "Getting FLUXNET data for DE-Hai ..."
## [1] "getting obs for IT-MBo"
## [1] "Getting FLUXNET data for IT-MBo ..."
## [1] "getting obs for US-GLE"
## [1] "Getting FLUXNET data for US-GLE ..."
## [1] "getting obs for FR-Fon"
## [1] "Getting FLUXNET data for FR-Fon ..."
## [1] "getting obs for NL-Hor"
## [1] "Getting FLUXNET data for NL-Hor ..."
## [1] "getting obs for US-UMd"
## [1] "Getting FLUXNET data for US-UMd ..."
## [1] "getting obs for AU-Dry"
## [1] "Getting FLUXNET data for AU-Dry ..."
## [1] "getting obs for DE-Obe"
## [1] "Getting FLUXNET data for DE-Obe ..."
## [1] "getting obs for IT-Tor"
## [1] "Getting FLUXNET data for IT-Tor ..."
## [1] "getting obs for US-Wi4"
## [1] "Getting FLUXNET data for US-Wi4 ..."Calibrate the model.
set.seed(1982)
settings_calib <- calib_sofun(
settings_calib,
df_drivers,
ddf_obs = ddf_obs_calib
)## kphio
## 0.05231893
## [1] "writing output from GenSA function to .//out_gensa_ORG.Rdata"
## [1] "writing calibrated parameters to .//params_opt_ORG.csv"The calibrated parameters are returned by calib_sofun() as part of the list:
print(settings_calib$par_opt)## kphio
## 0.05231893Evaluate
Run the model once again with these parameters and evaluate results.
mylist <- readr::read_csv("~/eval_pmodel/myselect_fluxnet2015.csv") %>%
dplyr::filter( use==1 ) %>%
dplyr::pull( Site )
settings_eval <- list(
sitenames = settings_sims$sitename,
sitenames_siteplots = mylist,
agg = 8,
path_fluxnet2015_d = "~/data/FLUXNET-2015_Tier1/20160128/point-scale_none_1d/original/unpacked/",
path_fluxnet2015_w = "~/data/FLUXNET-2015_Tier1/20160128/point-scale_none_7d/original/unpacked/",
path_fluxnet2015_m = "~/data/FLUXNET-2015_Tier1/20160128/point-scale_none_1m/original/unpacked/",
path_fluxnet2015_y = "~/data/FLUXNET-2015_Tier1/20160128/point-scale_none_1y/original/unpacked/",
path_gepisat_d = "~/data/gepisat/v3_fluxnet2015/daily_gpp/",
benchmark = list( gpp = c("fluxnet2015_NT") ),
remove_premodis = TRUE
)Get evaluation data (benchmarking data).
filn <- "./obs_eval.Rdata"
if (file.exists(filn)){
load(filn)
} else {
obs_eval <- get_obs_eval(
settings_eval = settings_eval,
settings_sims = settings_sims,
overwrite = TRUE,
light = TRUE,
add_forcing = FALSE
)
save(obs_eval, file = filn)
} Now run the model with calibrated parameters.
params_modl <- list(
kphio = settings_calib$par_opt[["kphio"]],
soilm_par_a = 1.0,
soilm_par_b = 0.0,
vpdstress_par_a = 0.2,
vpdstress_par_b = 0.2,
vpdstress_par_m = 5
)
mod <- runread_sofun_f(
df_drivers,
params_modl = params_modl,
makecheck = TRUE
) %>%
rename(id = sitename) %>%
unnest(out_sofun)And finally do the evaluation.
out_eval <- eval_sofun(
mod,
settings_eval,
settings_sims,
obs_eval = obs_eval,
overwrite = TRUE,
light = FALSE
)## [1] "-----------------------------------"
## [1] "gpp"
## [1] "-----------------------------------"
## [1] "Aggregating model outputs..."
## [1] "Evaluate annual values..."
## [1] "Evaluate monthly values..."
## [1] "Evaluate spatial values..."
## [1] "Evaluate interannual variability..."
## [1] "Evaluate mean seasonal cycle..."
## [1] "Evaluate mean seasonal cycle by climate zones..."
## [1] "Evaluate inter-day variability..."
## [1] "Evaluate mean seasonal cycle by X-day periods..."
## [1] "Evaluate inter-X-day variability..."
## [1] "Done with eval_sofun()."Print some results.
out_eval$gpp$fluxnet2015$metrics$xdaily_pooled## # A tibble: 1 x 5
## rsq rmse slope bias nvals
## <dbl> <dbl> <dbl> <dbl> <int>
## 1 0.736 2.23 1.09 -0.172 9569out_eval$gpp$fluxnet2015$data$xdf %>% rbeni::analyse_modobs2("mod", "obs", type = "heat")## Loading required package: ggthemes## Loading required package: RColorBrewer## $df_metrics
## # A tibble: 10 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 rmse standard 2.23
## 2 rsq standard 0.736
## 3 mae standard 1.53
## 4 n standard 9569
## 5 slope standard 1.09
## 6 mean_obs standard 4.14
## 7 prmse standard 0.538
## 8 pmae standard 0.369
## 9 bias standard -0.172
## 10 pbias standard -2.89
##
## $gg
##
## $linmod
##
## Call:
## lm(formula = obs ~ mod, data = df)
##
## Coefficients:
## (Intercept) mod
## -0.1742 1.0872
##
##
## $results
## # A tibble: 1 x 6
## rsq rmse mae bias slope n
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.736 2.23 1.53 -0.172 1.09 9569out_eval$gpp$fluxnet2015$data$ddf %>%
dplyr::filter(sitename=="BE-Vie" & year(date) < 2005) %>%
ggplot(aes(x = date)) +
geom_line(aes(y = obs), col="black") +
geom_line(aes(y = mod), col="red") +
labs(title = "BE-Vie")
out_eval$gpp$fluxnet2015$data$ddf %>%
dplyr::filter(sitename=="AU-Dry" & year(date) > 2010) %>%
ggplot(aes(x = date)) +
geom_line(aes(y = obs), col="black") +
geom_line(aes(y = mod), col="red") +
labs(title = "AU-Dry")## Warning: Removed 33 rows containing missing values (geom_path).