In computationales/sofunCalVal: Calibrate and validate rsofun

This workflow runs the P-model as implemented in {rsofun} with its latest version.

if ("rsofun" %in% installed.packages()) remove.packages("rsofun")
devtools::install_github(
  "geco-bern/rsofun",
  ref = "HEAD",
  upgrade = "never",
  force = TRUE
  )

# Load packages
library(rpmodel)
library(rsofun)
library(dplyr)
library(tidyr)
library(readr)
library(ggplot2)
library(lubridate)
library(knitr)
library(ggthemes)
library(here)

# Load functions
source(here::here("R/eval_sofun.R"))
source(here::here("R/get_stats.R"))
source(here::here("R/analyse_modobs2.R"))
source(here::here("R/create_data_split.R"))
source(here::here("R/calibrate_rsofun.R"))
source(here::here("R/align_events.R"))
source(here::here("R/eval_droughtresponse.R"))

# site selection: good-quality sites from analysis of Stocker et al. (2018) New Phyt.
df_flue <- readr::read_csv(here("data/flue_stocker18nphyt.csv"))

# Set seed
set.seed(1982)

Calibrate

Before compiling this vignette, run

# for multiple folds
# source(here::here("analysis/01_submit_calibration_current_run.R"))

# for single fold
source(here::here("analysis/rscript_calibrate_current.R"))

to calibrate parameters for five folds, split by site for spatial cross-validation. This creates five files containing calibrated parameters and the {BayesianTools} model object returned from the MCMC. Read the files.

# source(here::here("analysis/01_submit_calibration_current_run.R"))
par_calib <- readRDS(here::here("data/par_calib_light_current.rds"))

knitr::kable(t(par_calib))

Predict and test

# load driver data for all sites
drivers <- readRDS(here::here("data/drivers.rds"))

# Combine fixed and calibrated model parameters
# NOTE: this must be consistent with values specified in calibrate_rsofun().
par_fixed <- list(
  beta_unitcostratio = 146.0,
  kc_jmax            = 0.41,
  rd_to_vcmax        = 0.014,
  tau_acclim         = 20.0
)

par <- c(par_calib, par_fixed)

output <- runread_pmodel_f(
  drivers = drivers,
  par = par
) 

Run evaluation

This runs the eval_sofun() routine included in the sofunCalVal package, and outputs daily, monthly and annual summary statistics comparing observed with simulated values. We'll retake this same routine when running on the latest release (i.e. the workflow is equivalent but for the release of the rsofun package)

obs_eval <- read_rds(here::here("data/obs_eval.rds"))

evalsites <- output |>
  mutate(ntsteps = purrr::map_dbl(data, ~nrow(.))) |>
  dplyr::filter(ntsteps > 0) |>
  pull(sitename)

settings_eval <- list(
  benchmark = list( gpp = c("fluxnet"), le = c("fluxnet") ),
  sitenames = evalsites,
  agg       = 8
)

##  Evaluate model ----
out_eval <- eval_sofun(
  output,
  settings_eval,
  obs_eval = obs_eval,
  overwrite = TRUE,
  light = FALSE
)

Results

GPP

out_eval$gpp$fluxnet$metrics |> 
  bind_rows(.id = "Level") |> 
  kable()

8-daily, spatial and annual

out_eval$gpp$fluxnet$plot$gg_modobs_xdaily
out_eval$gpp$fluxnet$plot$gg_modobs_spatial_annual

Mean seasonal cycle

out_eval$gpp$fluxnet$data$meandoydf_byclim |> 
  dplyr::filter(climatezone %in% c(
    "Aw south", "Am south", "BSk north", "Bsk north",
    "Cfa north", "Cfb north", "Cfb south", 
    "Csa north", "Csb north", "Dfb north",
    "Dfc north", "ET north" )
    ) |>
  dplyr::filter(koeppen_code != "-") |> 
  pivot_longer(
    c(obs_mean, mod_mean),
    names_to = "source",
    values_to = "gpp"
    ) |> 
  ggplot() +
  geom_ribbon(
    aes(x = doy, ymin = obs_min, ymax = obs_max), 
    fill = "black", 
    alpha = 0.2
    ) +
  geom_line(aes(x = doy, y = gpp, color = source), size = 0.4) +
  labs(y = expression( paste("Simulated GPP (g C m"^-2, " d"^-1, ")" ) ), 
       x = "DOY") +
  facet_wrap( ~climatezone ) +
  theme_gray() +
  theme(legend.position = "bottom") +
  scale_color_manual(
    name="Setup: ",
    values=c("red", "black")
    )

Drought response

df_dday_agg <- eval_droughtresponse( 
  df = out_eval$gpp$fluxnet$data$ddf |> 
    rename(site = sitename), 
  df_flue,
  before = 20,
  after = 105,
  leng_threshold = 10, 
  nbins = 10, 
  do_norm = TRUE
  )

df_dday_agg |> 
  ggplot() +
  geom_hline(
    yintercept = 0,
    color = "black",
    linetype = "dotted"
    ) +
  geom_vline(
    xintercept = 0,
    color = "black",
    linetype = "dotted"
    ) +
  geom_line(
    aes(
      x = dday,
      y = median
      ),
    size = 0.9) +
  geom_ribbon(
    aes(
      x = dday,
      ymin = q33,
      ymax = q66
      ), 
    alpha = 0.3) +
  scale_color_manual(
    values = c(
      "BRC" = "black",
      "FULL" = "royalblue"
      ),
    name = "Setup"
    ) +
  scale_fill_manual(
    values = c(
      "BRC" = "black",
      "FULL" = "royalblue"
      ),
    name = "Setup") +
  ylim(-1.2, 2.2) + 
  xlim(-20, 105) +
  scale_x_continuous(
    expand = c(0,0)
    ) +
  scale_y_continuous(
    expand = c(0,0)
    ) +
  labs(
    x = "Days after drought onset",
    y = expression( paste( "Bias (g C m"^{-1}, " d"^{-1}, ")")) 
    ) +
  theme_classic()

LE

out_eval$le$fluxnet$metrics |> 
  bind_rows(.id = "Level") |> 
  kable()

8-daily, spatial and annual

out_eval$le$fluxnet$plot$gg_modobs_xdaily
out_eval$le$fluxnet$plot$gg_modobs_spatial_annual

Mean seasonal cycle

out_eval$le$fluxnet$data$meandoydf_byclim |> 
  dplyr::filter(climatezone %in% c(
    "Aw south", "Am south", "BSk north", "Bsk north",
    "Cfa north", "Cfb north", "Cfb south", 
    "Csa north", "Csb north", "Dfb north",
    "Dfc north", "ET north" )
    ) |>
  dplyr::filter(koeppen_code != "-") |> 
  pivot_longer(
    c(obs_mean, mod_mean),
    names_to = "source",
    values_to = "le"
    ) |> 
  ggplot() +
  geom_ribbon(
    aes(x = doy, ymin = obs_min, ymax = obs_max), 
    fill = "black", 
    alpha = 0.2
    ) +
  geom_line(aes(x = doy, y = le, color = source), size = 0.4) +
  labs(y = expression( paste("Simulated GPP (g C m"^-2, " d"^-1, ")" ) ), 
       x = "DOY") +
  facet_wrap( ~climatezone ) +
  theme_gray() +
  theme(legend.position = "bottom") +
  scale_color_manual(
    name="Setup: ",
    values=c("royalblue", "black")
    )

Drought response

df_dday_agg <- eval_droughtresponse( 
  df = out_eval$le$fluxnet$data$ddf |> 
    rename(site = sitename), 
  df_flue,
  before = 20,
  after = 105,
  leng_threshold = 10, 
  nbins = 10, 
  do_norm = TRUE
  )

df_dday_agg |> 
  ggplot() +
  geom_hline(
    yintercept = 0,
    color = "black",
    linetype = "dotted"
    ) +
  geom_vline(
    xintercept = 0,
    color = "black",
    linetype = "dotted"
    ) +
  geom_line(
    aes(
      x = dday,
      y = median
      ),
    size = 0.9) +
  geom_ribbon(
    aes(
      x = dday,
      ymin = q33,
      ymax = q66
      ), 
    alpha = 0.3) +
  scale_color_manual(
    values = c(
      "BRC" = "black",
      "FULL" = "royalblue"
      ),
    name = "Setup"
    ) +
  scale_fill_manual(
    values = c(
      "BRC" = "black",
      "FULL" = "royalblue"
      ),
    name = "Setup") +
  ylim(-1.2, 2.2) + 
  xlim(-20, 105) +
  scale_x_continuous(
    expand = c(0,0)
    ) +
  scale_y_continuous(
    expand = c(0,0)
    ) +
  labs(
    x = "Days after drought onset",
    y = expression( paste( "Bias (g C m"^{-1}, " d"^{-1}, ")")) 
    ) +
  theme_classic()

Consistency with rpmodel

Since rsofun takes time series forcing, overwrite forcing with constant values corresponding to the arguments provided to rpmodel::rpmodel().

Setup ORG

out_pmodel <- rpmodel( 
  tc             = 20,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 30,           # mol/m2/d,
  elv            = 0,            # m.a.s.l.,
  kphio          = par$kphio,    # quantum yield efficiency as calibrated
  beta           = par$beta_unitcostratio, # unit cost ratio a/b,
  c4             = FALSE,
  method_jmaxlim = "wang17",
  do_ftemp_kphio = FALSE,        # corresponding to setup ORG
  do_soilmstress = FALSE,        # corresponding to setup ORG
  verbose        = TRUE
  ) |> 
  as_tibble()

## overwrite forcing with constant conditions (for one site)
drivers_rpmodel <- drivers |> 
  slice(1) |> 
  mutate(forcing = purrr::map(
    forcing, 
    ~mutate(
      ., 
      temp = 20,
      vpd = 1000,
      ppfd = 30 / (60*60*24),
      patm = 101325,
      fapar = 1.0,
      co2 = 400,
      tmin = 20,
      tmax = 20
      )
    )) |> 
  mutate(params_siml = purrr::map(
    params_siml, 
    ~mutate(
      .,
      soilmstress = FALSE,
      tempstress = FALSE
      ))
    )

out_rsofun <- runread_pmodel_f(
  drivers_rpmodel,
  par = par
  ) |> 
  dplyr::select(data) |>
  unnest(data) |> 
  slice(1)

print("Are values equivalent for:")
paste("- ci:ca:", all.equal(out_pmodel$chi, out_rsofun$chi, tolerance = 1e-5))
paste("- GPP:", all.equal(out_pmodel$gpp, out_rsofun$gpp, tolerance = 1e-5))
paste("- Vcmax:", all.equal(out_pmodel$vcmax / (60*60*24), out_rsofun$vcmax, tolerance = 1e-5))
paste("- Vcmax25:", all.equal(out_pmodel$vcmax25 / (60*60*24), out_rsofun$vcmax25, tolerance = 1e-5))
paste("- Jmax:", all.equal(out_pmodel$jmax / (60*60*24), out_rsofun$jmax, tolerance = 1e-5))
paste("- Jmax25:", all.equal(out_pmodel$jmax25 / (60*60*24), out_rsofun$jmax25, tolerance = 1e-5))

Appendix

sessionInfo()

computationales/sofunCalVal documentation built on June 12, 2025, 6:59 a.m.