From 8a0f49bd56e532e4a2c920a0c90b7e5b8c622171 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Wed, 7 Jan 2026 13:34:48 -0800 Subject: [PATCH 01/19] read soil water options from settings, increment model version --- models/rothc/Dockerfile | 11 +---------- models/rothc/R/write.config.RothC.R | 7 ++++--- models/rothc/inst/RothC_input_template.dat | 2 +- models/rothc/inst/workflow_example/template.xml | 8 +++++--- 4 files changed, 11 insertions(+), 17 deletions(-) diff --git a/models/rothc/Dockerfile b/models/rothc/Dockerfile index a9b5af2dd2..2811e0a0f3 100644 --- a/models/rothc/Dockerfile +++ b/models/rothc/Dockerfile @@ -7,21 +7,12 @@ ARG IMAGE_VERSION="latest" # ---------------------------------------------------------------------- FROM pecan/models:${IMAGE_VERSION} -# ---------------------------------------------------------------------- -# INSTALL MODEL SPECIFIC PIECES -# ---------------------------------------------------------------------- - -#RUN apt-get update \ -# && apt-get install -y --no-install-recommends \ -# python3.9 \ -# && rm -rf /var/lib/apt/lists/* - # ---------------------------------------------------------------------- # SETUP FOR SPECIFIC MODEL # ---------------------------------------------------------------------- # Some variables that can be used to set control the docker build -ARG MODEL_VERSION=2.1.0 +ARG MODEL_VERSION=2.1.1 RUN mkdir -p /tmp/rothc \ && curl -sSL https://github.com/Rothamsted-Models/RothC_Code/archive/refs/tags/v${MODEL_VERSION}.tar.gz \ diff --git a/models/rothc/R/write.config.RothC.R b/models/rothc/R/write.config.RothC.R index 452396cfad..26f586bd5d 100644 --- a/models/rothc/R/write.config.RothC.R +++ b/models/rothc/R/write.config.RothC.R @@ -125,18 +125,19 @@ write.config.RothC <- function(defaults, trait.values, settings, run.id) { config.text <- gsub("@ENSNAME@", run.id, config.text) config.text <- gsub("@OUTFILE@", paste0("out", run.id), config.text) - # TODO make these editable -- hard-coding for MVP # OPT_RMMOIST: soil water parameterization. # 1: Standard RothC soil water parameters # 2: Van Genuchten soil properties and soil is allowed to be drier # (ie hygroscopic / capillary water, -1000bar) # 3: Van Genuchten soil properties, but uses the Standard RothC # soil water function - config.text <- gsub("@OPT_RMMOIST@", "1", config.text) + rmmoist <- settings$model$opt_RMmoist %||% 1 + config.text <- gsub("@OPT_RMMOIST@", rmmoist, config.text) # Bare SMD: wilting point configuration # 1: Standard RothC bareSMD # 2: bareSMD is set to wilting point -15bar (could be better for dry soils) - config.text <- gsub("@OPT_SDDBARE@", "1", config.text) + smdbare <- settings$model$opt_RMmoist %||% 1 + config.text <- gsub("@OPT_SMDBARE@", smdbare, config.text) ## Climate data # (read here to use length in soil params, remainder of processing happens below) diff --git a/models/rothc/inst/RothC_input_template.dat b/models/rothc/inst/RothC_input_template.dat index 761ee96086..d02a6c45c2 100644 --- a/models/rothc/inst/RothC_input_template.dat +++ b/models/rothc/inst/RothC_input_template.dat @@ -2,7 +2,7 @@ # @START_YEAR@-@START_MONTH@-@START_DAY@ to @END_YEAR@-@END_MONTH@-@END_DAY@ # # opt_RMmoist opt_SMDbare - @OPT_RMMOIST@ @OPT_SDDBARE@ + @OPT_RMMOIST@ @OPT_SMDBARE@ # (%) (cm) (t C/ha) (-) (%) (g/m3) (%) (-) # clay depth iom nsteps siltper BD OC minRM_Moist @SOIL_PARAMS@ diff --git a/models/rothc/inst/workflow_example/template.xml b/models/rothc/inst/workflow_example/template.xml index 3f19846611..b05cc0b4d3 100644 --- a/models/rothc/inst/workflow_example/template.xml +++ b/models/rothc/inst/workflow_example/template.xml @@ -38,11 +38,13 @@ - rothc_2_1_0 + rothc_2_1_1 RothC - 2.1.0 + 2.1.1 FALSE - /usr/local/bin/RothC_v2.1.0 + /usr/local/bin/RothC_v2.1.1 + 1 + 1 From fd8a5b5cf0fcf1d7331a6060c8853a1f99c1af97 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Wed, 7 Jan 2026 13:35:19 -0800 Subject: [PATCH 02/19] changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b05e4780c8..45e79bff9b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha ## Unreleased ### Added +- New package `PEcAn.RothC` runs the RothC soil carbon model. - Add function `PEcAn.MA::meta_analysis_standalone` to run meta-analysis without database or file IO. - Added Demo 03: Meta Analysis Quarto notebook (`documentation/tutorials/Demo_03_Meta_Analysis/meta_analysis.qmd`) to demonstrate how to perform Bayesian meta-analysis and visualize posterior distributions using pre-generated trait data. From 10b724d0ff47427996e412d9a6a233e22cf9eff6 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Wed, 7 Jan 2026 13:35:55 -0800 Subject: [PATCH 03/19] formatting --- models/rothc/R/write.config.RothC.R | 73 ++++++++++++++--------------- 1 file changed, 36 insertions(+), 37 deletions(-) diff --git a/models/rothc/R/write.config.RothC.R b/models/rothc/R/write.config.RothC.R index 26f586bd5d..21189453d4 100644 --- a/models/rothc/R/write.config.RothC.R +++ b/models/rothc/R/write.config.RothC.R @@ -1,15 +1,15 @@ -##' Writes a RothC config file. -##' -##' Requires a pft xml object, a list of trait values for a single model run, -##' and the name of the file to create -##' -##' @param defaults list of defaults to process -##' @param trait.values vector of samples for a given trait -##' @param settings list of settings from pecan settings file -##' @param run.id id of run -##' @return configuration file for MODEL for given run -##' @export -##' @author Chris Black +#' Writes a RothC config file. +#' +#' Requires a pft xml object, a list of trait values for a single model run, +#' and the name of the file to create +#' +#' @param defaults list of defaults to process +#' @param trait.values vector of samples for a given trait +#' @param settings list of settings from pecan settings file +#' @param run.id id of run +#' @return configuration file for MODEL for given run +#' @export +#' @author Chris Black write.config.RothC <- function(defaults, trait.values, settings, run.id) { # find out where to write run/ouput @@ -140,7 +140,8 @@ write.config.RothC <- function(defaults, trait.values, settings, run.id) { config.text <- gsub("@OPT_SMDBARE@", smdbare, config.text) ## Climate data - # (read here to use length in soil params, remainder of processing happens below) + # (we read it here to use its length in soil params, + # remainder of processing happens below) met_path <- settings$run$inputs$met$path met_in <- utils::read.table(met_path, header = TRUE) n_met <- nrow(met_in) @@ -177,31 +178,29 @@ write.config.RothC <- function(defaults, trait.values, settings, run.id) { OA_HUM_f = 0.02 ) -input_rows <- met_in |> - dplyr::bind_cols(inputs) |> - dplyr::select( - "year", "month", - "modern_pct", - "Tmp_C", "Rain_mm", "Evap_mm", - "C_inp_tC_ha", "FYM_tC_ha", "PC", - "PL_DPM_f", "PL_RPM_f", - "OA_DPM_f", "OA_RPM_f", "OA_BIO_f", "OA_HUM_f" - ) |> - # Duplicate first year as the equilibrium block - # TODO we probably want a more principled approach here - duplicate_first_year() |> - dplyr::mutate( - dplyr::across(dplyr::where(is.double), zapsmall) - ) |> - # Kinda ugly: Convert to one string to cram it into the template via gsub - format() |> - apply(1, paste, collapse = " ") |> - paste(collapse = "\n") - + input_rows <- met_in |> + dplyr::bind_cols(inputs) |> + dplyr::select( + "year", "month", + "modern_pct", + "Tmp_C", "Rain_mm", "Evap_mm", + "C_inp_tC_ha", "FYM_tC_ha", "PC", + "PL_DPM_f", "PL_RPM_f", + "OA_DPM_f", "OA_RPM_f", "OA_BIO_f", "OA_HUM_f" + ) |> + # Duplicate first year as the equilibrium block + # TODO we probably want a more principled approach here + duplicate_first_year() |> + dplyr::mutate( + dplyr::across(dplyr::where(is.double), zapsmall) + ) |> + # Kinda ugly: Convert to one string to cram it into the template via gsub + format() |> + apply(1, paste, collapse = " ") |> + paste(collapse = "\n") config.text <- gsub("@CLIM_DATA@", input_rows, config.text) - config.file.name <- "RothC_input.dat" - writeLines(config.text, con = file.path(rundir, config.file.name)) + writeLines(config.text, con = file.path(rundir, "RothC_input.dat")) invisible(config.text) } @@ -211,4 +210,4 @@ duplicate_first_year <- function(df) { dplyr::slice_min(.data$year) |> dplyr::mutate(year = 1) |> dplyr::bind_rows(df) -} \ No newline at end of file +} From 9f65c88aa423192cdf510b7cda5641aea960c4d7 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Fri, 9 Jan 2026 13:26:45 -0800 Subject: [PATCH 04/19] read soil water options from settings, soil params from soil_physics file --- models/rothc/R/rothc_varname_map.R | 11 ++++- models/rothc/R/write.config.RothC.R | 62 +++++++++++++++++++++++++++-- 2 files changed, 69 insertions(+), 4 deletions(-) diff --git a/models/rothc/R/rothc_varname_map.R b/models/rothc/R/rothc_varname_map.R index 3adf9f6660..e8162e0856 100644 --- a/models/rothc/R/rothc_varname_map.R +++ b/models/rothc/R/rothc_varname_map.R @@ -19,5 +19,14 @@ rothc_varname_map <- dplyr::tribble( "Hum_t_C_ha", "Humified organic matter", "t C ha-1", "slow_soil_pool_carbon_content", "IOM_t_C_ha", "Inert organic matter", "t C ha-1", "structural_soil_pool_carbon_content", "SOC_t_C_ha", "Total soil organic carbon", "t C ha-1", "TotSoilCarb", - "CO2_t_C_ha", "Accumulated CO2", "t C ha-1", NA # needs time dimension to match to "TotalResp" which is kg C m-2 sec-1 + "CO2_t_C_ha", "Accumulated CO2", "t C ha-1", NA, # needs time dimension to match to "TotalResp" which is kg C m-2 sec-1 + ## The remaining lines are parameters not variables, but mapping them here for "convenience" + "clay_pct", "Clay content of the soil", "pct", "fraction_of_clay_in_soil", + "depth_cm", "Depth of soil layer sampled", "cm", "depth", + "iom_tC_ha", "inert organic matter", "t C ha-1", NA, + # Note: these 4 params only used when opt_rmmoist = 2 + "silt_pct", "Silt content of the soil", "pct", "fraction_of_silt_in_soil", + "bulkdens_g_cm3", "Bulk density of the soil", "g cm-3", "soil_bulk_density", + "org_C_pct", "Organic carbon content of the soil", "pct", NA, # need to back-convert from soil_organic_carbon_stock? + "min_RM_moist", "min value for soil moisture rate modifier, reached at -1500 kPa", "1", NA ) diff --git a/models/rothc/R/write.config.RothC.R b/models/rothc/R/write.config.RothC.R index 21189453d4..218c1e712b 100644 --- a/models/rothc/R/write.config.RothC.R +++ b/models/rothc/R/write.config.RothC.R @@ -133,12 +133,26 @@ write.config.RothC <- function(defaults, trait.values, settings, run.id) { # soil water function rmmoist <- settings$model$opt_RMmoist %||% 1 config.text <- gsub("@OPT_RMMOIST@", rmmoist, config.text) + # Bare SMD: wilting point configuration # 1: Standard RothC bareSMD # 2: bareSMD is set to wilting point -15bar (could be better for dry soils) smdbare <- settings$model$opt_RMmoist %||% 1 config.text <- gsub("@OPT_SMDBARE@", smdbare, config.text) + # min_RM_moist: Lowest value taken by the soil moisture rate modifier, + # reached when soil water potential is -1500 kPa. + # Note: Only used if smdbare == 2. + # 0.2 is compatible with RothC 1.0; 0.15 or 0.1 may be better for dry systems. + # see Farina et al 2013 10.1016/j.geoderma.2013.01.021 + min_rmmoist <- settings$model$min_RM_moist %||% 0.2 + # not added to config.text here bc written as part of @SOIL_PARAMS@ below + + # Soil depth to simulate. + # Caution: Currently gets overridden by layer boundaries in soil file, + # so it's really "simulate no deeper than this". TODO: FIXME. + model_depth_cm <- settings$model$soil_depth_cm %||% 23 + ## Climate data # (we read it here to use its length in soil params, # remainder of processing happens below) @@ -147,12 +161,54 @@ write.config.RothC <- function(defaults, trait.values, settings, run.id) { n_met <- nrow(met_in) ## Soil parameters + soil_list <- PEcAn.data.land::pool_ic_netcdf2list( + settings$run$inputs$soil_physics$path + ) + soil_params <- soil_list$vals |> + as.data.frame() |> + dplyr::mutate(depth_cm = soil_list$dims$depth) |> + # TODO this drops layers that extend past bottom + # (eg with depth=23 and 0-10/10-30 layering, would use only 0-10) + # Least-code/copout approach: + # Make everyone generate their soil files with layers that match model depth + dplyr::filter(.data$depth_cm <= model_depth_cm) |> + dplyr::summarize( + # TODO consider weighting by layer thickness? + depth_cm = max(.data$depth_cm), + clay_pct = .data$fraction_of_clay_in_soil |> + mean() |> + PEcAn.utils::ud_convert("1", "%"), + silt_pct = .data$fraction_of_silt_in_soil |> + mean() |> + PEcAn.utils::ud_convert("1", "%"), + bulkdens_g_cm3 = .data$bulk_density |> + mean() |> + PEcAn.utils::ud_convert("kg m-3", "g cm-3"), + org_C_pct = .data$soil_organic_carbon_stock |> + sum() |> + PEcAn.utils::ud_convert("kg m-2", "g cm-2") |> + (\(x) x / (.data$depth_cm * .data$bulkdens_g_cm3))() |> + PEcAn.utils::ud_convert("1", "%"), + iom_tC_ha = .data$soil_organic_carbon_stock |> + sum() |> + PEcAn.utils::ud_convert("kg m-2", "t ha-1") |> + # Approximation from Falloon et al. 1998, 10.1016/S0038-0717(97)00256-3 + (\(x) 0.049 * x^1.139)() + ) + + ## Build string from soil params ## (plus number of timesteps, weirdly snuck into the middle) - # TODO: read from run$inputs$soil_physics soil_param_string <- paste( - "23.4 23.0 3.0041", # clay_pct, depth_cm, iom_tC_ha + soil_params$clay_pct, + soil_params$depth_cm, + soil_params$iom_tC_ha, n_met + 12, # nsteps (includes extra year for spinup) - "58.6 1.27 0.94 0.2" # silt_pct, bulkdens_g_m3, org_C_pct, min_RM_moist + # NB these last 4 params are always written, + # but only used by model if smdbar == 2 + soil_params$silt_pct, + soil_params$bulkdens_g_m3, + soil_params$org_C_pct, + min_rmmoist ) config.text <- gsub("@SOIL_PARAMS@", soil_param_string, config.text) From 3f9c57080140e66b35a357a9e75263799a4690f1 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Sat, 10 Jan 2026 22:22:15 -0800 Subject: [PATCH 05/19] typos + move input params to own map --- models/rothc/R/rothc_varname_map.R | 11 +++++++++-- models/rothc/R/write.config.RothC.R | 2 +- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/models/rothc/R/rothc_varname_map.R b/models/rothc/R/rothc_varname_map.R index e8162e0856..1b153f8420 100644 --- a/models/rothc/R/rothc_varname_map.R +++ b/models/rothc/R/rothc_varname_map.R @@ -19,8 +19,15 @@ rothc_varname_map <- dplyr::tribble( "Hum_t_C_ha", "Humified organic matter", "t C ha-1", "slow_soil_pool_carbon_content", "IOM_t_C_ha", "Inert organic matter", "t C ha-1", "structural_soil_pool_carbon_content", "SOC_t_C_ha", "Total soil organic carbon", "t C ha-1", "TotSoilCarb", - "CO2_t_C_ha", "Accumulated CO2", "t C ha-1", NA, # needs time dimension to match to "TotalResp" which is kg C m-2 sec-1 - ## The remaining lines are parameters not variables, but mapping them here for "convenience" + "CO2_t_C_ha", "Accumulated CO2", "t C ha-1", NA # needs time dimension to match to "TotalResp" which is kg C m-2 sec-1 +) + + +# Mapping from RothC _input_ parameters to PEcAn standard names. +# Not used anywhere in the code yet! +# Just writing them down for reference and this file seemed like a good place +rothc_param_map <- dplyr::tribble( + ~rothc_name, ~rothc_description, ~rothc_unit, ~pecan_name, "clay_pct", "Clay content of the soil", "pct", "fraction_of_clay_in_soil", "depth_cm", "Depth of soil layer sampled", "cm", "depth", "iom_tC_ha", "inert organic matter", "t C ha-1", NA, diff --git a/models/rothc/R/write.config.RothC.R b/models/rothc/R/write.config.RothC.R index 218c1e712b..3d45bb9bf0 100644 --- a/models/rothc/R/write.config.RothC.R +++ b/models/rothc/R/write.config.RothC.R @@ -181,7 +181,7 @@ write.config.RothC <- function(defaults, trait.values, settings, run.id) { silt_pct = .data$fraction_of_silt_in_soil |> mean() |> PEcAn.utils::ud_convert("1", "%"), - bulkdens_g_cm3 = .data$bulk_density |> + bulkdens_g_cm3 = .data$soil_bulk_density |> mean() |> PEcAn.utils::ud_convert("kg m-3", "g cm-3"), org_C_pct = .data$soil_organic_carbon_stock |> From 11f64f91f92088bef00c53497908e4821c2c41e7 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Sat, 10 Jan 2026 22:29:44 -0800 Subject: [PATCH 06/19] add soil_physics input to rothc demo --- models/rothc/inst/workflow_example/README.md | 23 +++++-- .../inst/workflow_example/fetch_soil_data.R | 67 +++++++++++++++++++ models/rothc/inst/workflow_example/run.sh | 5 ++ .../rothc/inst/workflow_example/template.xml | 8 ++- .../rothc/inst/workflow_example/xml_build.R | 19 +++++- 5 files changed, 111 insertions(+), 11 deletions(-) create mode 100755 models/rothc/inst/workflow_example/fetch_soil_data.R diff --git a/models/rothc/inst/workflow_example/README.md b/models/rothc/inst/workflow_example/README.md index f4e9028785..adf4a664d6 100644 --- a/models/rothc/inst/workflow_example/README.md +++ b/models/rothc/inst/workflow_example/README.md @@ -11,21 +11,30 @@ this set simple enough to be easy to (1) understand as demos, and ## Important caveats This package is still under development and many inputs are still hard-coded. -As I write this on 2025-12-04, only met data is read from site-specific files; -all soil and management is hard-coded. Do not interpret the outputs as -meaningful predictions yet. +As I write this on 2026-01-09, met and soil data are read from site-specific +files, and all management and initial conditions are hard-coded. +Do not interpret the outputs as meaningful predictions yet. ## Required but not yet provided here * ERA5 weather data in PEcAn standard netcdf format (instructions TK) -* a RothC binary, available from https://github.com/Rothamsted-Models/RothC_Code +* a RothC binary, compiled from https://github.com/Rothamsted-Models/RothC_Code ## To run -* put weather data in `data_raw/ERA5_CA_nc/` (or update paths to where your weather already is) -* Update line 45 of `template.xml` to the path where you installed your copy of RothC - (or put RothC at `/usr/local/bin/RothC_v2.1.0`) +* put weather data in `data_raw/ERA5_CA_nc/` (or update paths to where your + weather already is) +* Update line 45 of `template.xml` to the path where you installed your copy + of RothC (or put RothC at `/usr/local/bin/RothC_v2.1.1`) * update `site_info.csv` with your sites of interest * TK: Add site-specific management and soil information once implemented * `./run.sh` + +## Troubleshooting + +* `fetch_soil_data.R` (called in `run.sh`) uses + `PEcAn.data.land::extract_soil_gssurgo()`, which sometimes fails for sites + with only one soil type. A more robust version is under test at + https://github.com/PecanProject/pecan/pull/3643/files + If some of your sites fail, try updating PEcAn.data.land from that branch. diff --git a/models/rothc/inst/workflow_example/fetch_soil_data.R b/models/rothc/inst/workflow_example/fetch_soil_data.R new file mode 100755 index 0000000000..3dd7e9f949 --- /dev/null +++ b/models/rothc/inst/workflow_example/fetch_soil_data.R @@ -0,0 +1,67 @@ +#!/usr/bin/env Rscript + +# creates ensembles of soil condition files for each requested location. + +## run time option parsing +options <- list( + optparse::make_option("--site_info_path", + default = "site_info.csv", + help = "Path to a csv file with at least columns `id`, `lat`, `lon`", + ), + optparse::make_option("--out_dir", + default = "data/soil/", + help = paste( + "Output path: Will contain one subdirectory per site_id,", + "containing n_ensemble netCDF files each named", + "gSSURGO_soil_.dat" + ) + ), + optparse::make_option("--n_ens", + default = 20, + help = "number of files to generate per site" + ), + optparse::make_option("--overwrite", + default = FALSE, + help = paste( + "Replace files for existing sites, or skip them?", + "Note: when FALSE, refuses to overwrite even empty directories" + ) + ), + optparse::make_option("--n_cores", + default = Sys.getenv("NCPUS", 1L), + help = "number of CPUs to use in parallel", + ), + optparse::make_option("--parallel_strategy", + default = "multisession", + help = "Strategy for parallel conversion, passed to future::plan()", + ) +) |> + # Show default values in help message + purrr::modify(\(x) { + x@help <- paste(x@help, "[default: %default]") + x + }) + +args <- optparse::OptionParser(option_list = options) |> + optparse::parse_args() + +## end options + + +future::plan(args$parallel_strategy, workers = as.numeric(args$n_cores)) +# TODO option call not working -- do I need to assign it somewhere? +furrr::furrr_options(seed = TRUE) + +site_info <- + read.csv( + args$site_info_path, + colClasses = c(id = "character") + ) |> + dplyr::mutate(outdir = file.path(args$out_dir, id)) |> + dplyr::select(lat, lon, outdir) +if (!args$overwrite) { + site_info <- site_info |> + dplyr::filter(!dir.exists(outdir)) +} +site_info |> + furrr::future_pwalk(PEcAn.data.land::extract_soil_gssurgo, size = args$n_ens) diff --git a/models/rothc/inst/workflow_example/run.sh b/models/rothc/inst/workflow_example/run.sh index 53cf9fee89..8bebb5fba1 100755 --- a/models/rothc/inst/workflow_example/run.sh +++ b/models/rothc/inst/workflow_example/run.sh @@ -15,6 +15,11 @@ if [[ ! -d data/ERA5_CA_RothC ]]; then ./ERA5_nc_to_RothC.R --n_cores="$NCPUS" fi + +echo "================= Retrieving soil data from gSSURGO =================" +# fetches soil data for all sites +./fetch_soil_data.R --n_cores="$NCPUS" + echo "================= Building XML settings =================" # generates settings.xml ./xml_build.R --output_dir_name=output_"$start_time" diff --git a/models/rothc/inst/workflow_example/template.xml b/models/rothc/inst/workflow_example/template.xml index b05cc0b4d3..b124ee466b 100644 --- a/models/rothc/inst/workflow_example/template.xml +++ b/models/rothc/inst/workflow_example/template.xml @@ -29,10 +29,10 @@ + sampling - --> + @@ -45,6 +45,8 @@ /usr/local/bin/RothC_v2.1.1 1 1 + 0.1 + 30 diff --git a/models/rothc/inst/workflow_example/xml_build.R b/models/rothc/inst/workflow_example/xml_build.R index 8e0389b59e..d02918fb2b 100755 --- a/models/rothc/inst/workflow_example/xml_build.R +++ b/models/rothc/inst/workflow_example/xml_build.R @@ -32,7 +32,13 @@ options <- list( "Should contain subdirs named by site id" ) ), - + optparse::make_option("--soil_dir", + default = "data/soil", + help = paste( + "Directory containing netCDFs of soil data.", + "Should contain subdirs named by site id" + ) + ), optparse::make_option("--site_file", default = "site_info.csv", help = paste( @@ -137,6 +143,17 @@ settings <- settings |> # path = args$ic_dir, # path_template = "{path}/{id}/IC_site_{id}_{n}.nc" # ) |> + setEnsemblePaths( + n_reps = args$n_ens, + input_type = "soil_physics", + path = args$soil_dir, + # n+1 bc current implementation of extract_soil_gSSURGO saves a + # gSSURGO_soil_1.nc containing pseudo-layers from _every_ + # soil type found at site, then at least n more files after that + # ("at least" bc it also writes at least 1 file per soil type) + # TODO + path_template = "{path}/{id}/gSSURGO_soil_{n+1}.nc" + ) |> papply(add_soil_pft) # Update just the first component of the output directory, From d04530b3f76095291924291922f557aa2499758b Mon Sep 17 00:00:00 2001 From: Chris Black Date: Tue, 13 Jan 2026 09:55:49 -0800 Subject: [PATCH 07/19] move furrr opts --- models/rothc/inst/workflow_example/fetch_soil_data.R | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/models/rothc/inst/workflow_example/fetch_soil_data.R b/models/rothc/inst/workflow_example/fetch_soil_data.R index 3dd7e9f949..82d4215e9f 100755 --- a/models/rothc/inst/workflow_example/fetch_soil_data.R +++ b/models/rothc/inst/workflow_example/fetch_soil_data.R @@ -49,8 +49,6 @@ args <- optparse::OptionParser(option_list = options) |> future::plan(args$parallel_strategy, workers = as.numeric(args$n_cores)) -# TODO option call not working -- do I need to assign it somewhere? -furrr::furrr_options(seed = TRUE) site_info <- read.csv( @@ -64,4 +62,9 @@ if (!args$overwrite) { dplyr::filter(!dir.exists(outdir)) } site_info |> - furrr::future_pwalk(PEcAn.data.land::extract_soil_gssurgo, size = args$n_ens) + furrr::future_pwalk( + .f = PEcAn.data.land::extract_soil_gssurgo, + size = args$n_ens, + grid_size = 5, + .options = furrr::furrr_options(seed = TRUE) + ) From 35784efc71de7d9eca21ffbe103fefa5ea396abd Mon Sep 17 00:00:00 2001 From: Chris Black Date: Tue, 13 Jan 2026 20:23:00 -0800 Subject: [PATCH 08/19] round soil params, align with headers, fix hulk density typo g_m3->g_cm3 --- models/rothc/R/write.config.RothC.R | 15 ++++++++------- models/rothc/inst/RothC_input_template.dat | 4 ++-- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/models/rothc/R/write.config.RothC.R b/models/rothc/R/write.config.RothC.R index 3d45bb9bf0..9ed91df3d2 100644 --- a/models/rothc/R/write.config.RothC.R +++ b/models/rothc/R/write.config.RothC.R @@ -199,16 +199,17 @@ write.config.RothC <- function(defaults, trait.values, settings, run.id) { ## Build string from soil params ## (plus number of timesteps, weirdly snuck into the middle) soil_param_string <- paste( - soil_params$clay_pct, - soil_params$depth_cm, - soil_params$iom_tC_ha, + round(soil_params$clay_pct, 1), + round(soil_params$depth_cm, 1), + round(soil_params$iom_tC_ha, 2), n_met + 12, # nsteps (includes extra year for spinup) # NB these last 4 params are always written, # but only used by model if smdbar == 2 - soil_params$silt_pct, - soil_params$bulkdens_g_m3, - soil_params$org_C_pct, - min_rmmoist + round(soil_params$silt_pct, 1), + round(soil_params$bulkdens_g_cm3, 3), + round(soil_params$org_C_pct, 1), + min_rmmoist, + sep = " " # just to align with headers, RothC doesn't care ) config.text <- gsub("@SOIL_PARAMS@", soil_param_string, config.text) diff --git a/models/rothc/inst/RothC_input_template.dat b/models/rothc/inst/RothC_input_template.dat index d02a6c45c2..d6bb7ca2c0 100644 --- a/models/rothc/inst/RothC_input_template.dat +++ b/models/rothc/inst/RothC_input_template.dat @@ -3,8 +3,8 @@ # # opt_RMmoist opt_SMDbare @OPT_RMMOIST@ @OPT_SMDBARE@ -# (%) (cm) (t C/ha) (-) (%) (g/m3) (%) (-) -# clay depth iom nsteps siltper BD OC minRM_Moist +# (%) (cm) (t C/ha) (-) (%) (g/cm3) (%) (-) +# clay depth iom nsteps siltper BD OC minRM_Moist @SOIL_PARAMS@ # (-) (-) (%) (C) (mm) (mm) (t C/ha) (t C/ha) (-) (-) (-) (-) (-) (-) (-) # year month modern Tmp Rain Evap C_inp FYM PC PL_DPM_f PL_RPM_f OA_DPM_f OA_RPM_f OA_BIO_f OA_HUM_f From f76e8f71d02a4e2de57ad3ed7a904ab470c64dc5 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Tue, 13 Jan 2026 20:24:19 -0800 Subject: [PATCH 09/19] typo in soil depth param --- models/rothc/inst/workflow_example/template.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/rothc/inst/workflow_example/template.xml b/models/rothc/inst/workflow_example/template.xml index b124ee466b..f22607e8fb 100644 --- a/models/rothc/inst/workflow_example/template.xml +++ b/models/rothc/inst/workflow_example/template.xml @@ -46,7 +46,7 @@ 1 1 0.1 - 30 + 30 From acf88eda669a032ddaa737c20fc888a2e36f7a83 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Thu, 15 Jan 2026 11:19:33 -0800 Subject: [PATCH 10/19] notes for later --- models/rothc/R/write.config.RothC.R | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/models/rothc/R/write.config.RothC.R b/models/rothc/R/write.config.RothC.R index 9ed91df3d2..7c22758da8 100644 --- a/models/rothc/R/write.config.RothC.R +++ b/models/rothc/R/write.config.RothC.R @@ -166,11 +166,17 @@ write.config.RothC <- function(defaults, trait.values, settings, run.id) { ) soil_params <- soil_list$vals |> as.data.frame() |> + # netCDF metadata and PEcAn.data.land::soil.units("soil_depth") both + # say depth should be meters, but the gSSURGO files I've checked are in cm. + # TODO: fix either code or unit labels upstream, + # and check that they are consistent with other sources + # TODO 2: Assumes depth is given to bottom of layer -- is that correct? dplyr::mutate(depth_cm = soil_list$dims$depth) |> # TODO this drops layers that extend past bottom # (eg with depth=23 and 0-10/10-30 layering, would use only 0-10) - # Least-code/copout approach: - # Make everyone generate their soil files with layers that match model depth + # Consider rescaling partial layers + # (Or throwing an error on mismatch and making everyone generate their soil + # files with layers that match model depth?) dplyr::filter(.data$depth_cm <= model_depth_cm) |> dplyr::summarize( # TODO consider weighting by layer thickness? From d05b247252c81032af7f1df489d536c3e7d30e60 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Sat, 17 Jan 2026 23:58:45 -0800 Subject: [PATCH 11/19] move soil file reading to own function --- models/rothc/R/read_soil_physics.R | 70 +++++++++++++++++ models/rothc/R/write.config.RothC.R | 42 +--------- models/rothc/man/read_soil_physics.Rd | 29 +++++++ .../tests/testthat/test-read_soil_physics.R | 78 +++++++++++++++++++ 4 files changed, 178 insertions(+), 41 deletions(-) create mode 100644 models/rothc/R/read_soil_physics.R create mode 100644 models/rothc/man/read_soil_physics.Rd create mode 100644 models/rothc/tests/testthat/test-read_soil_physics.R diff --git a/models/rothc/R/read_soil_physics.R b/models/rothc/R/read_soil_physics.R new file mode 100644 index 0000000000..64521577ef --- /dev/null +++ b/models/rothc/R/read_soil_physics.R @@ -0,0 +1,70 @@ +#' Read soil parameters from a PEcAn soil physics file +#' +#' Reads values from a netCDF that follows the PEcAn soil standard, +#' aggregates layers across the depth to be simulated, +#' and converts to names and units expected by RothC. +#' +#' TODO: Currently drops all of any layer that extends below `model_depth`. +#' It would be better instead to weight all layers by their contribution +#' to the requested depth. +#' +#' @param path filepath to a netcdf, +#' probably read from `settings$run$inputs$soil_physics$path` +#' @param model_depth Soil depth to be simulated, in cm +#' +#' @return one-row dataframe with columns +#' `depth_cm`, `clay_pct`, `silt_pct`, +#' `bulkdens_g_cm3`, `org_C_pct`, `iom_tC_ha` +#' +#' @importFrom rlang .data .env +#' +read_soil_physics <- function(path, model_depth = 23) { + soil_list <- PEcAn.data.land::pool_ic_netcdf2list(path) + + # Input is documented to be in meters, + # but providing cm instead seems to be a _very_ common error; + # might as well try to handle it here + if (all(soil_list$dims$depth < 10)) { + depth_cm <- PEcAn.utils::ud_convert(soil_list$dims$depth, "m", "cm") + } else { + PEcAn.logger::logger.warn( + "Soil depths should be in meters, but found values >= 10", + "in file", path, + "Assuming these are mislabeled cm and treating them as such." + ) + depth_cm <- soil_list$dims$depth + } + soil_list$vals |> + as.data.frame() |> + # TODO: Assumes depth is given to bottom of layer -- is that correct? + dplyr::mutate(depth_cm = .env$depth_cm) |> + # TODO this drops layers that extend past bottom + # (eg with depth=23 and 0-10/10-30 layering, would use only 0-10) + # Consider rescaling partial layers + # (Or throwing an error on mismatch and making everyone generate their soil + # files with layers that match model depth?) + dplyr::filter(.data$depth_cm <= model_depth) |> + dplyr::summarize( + # TODO consider weighting by layer thickness? + depth_cm = max(.data$depth_cm), + clay_pct = .data$fraction_of_clay_in_soil |> + mean() |> + PEcAn.utils::ud_convert("1", "%"), + silt_pct = .data$fraction_of_silt_in_soil |> + mean() |> + PEcAn.utils::ud_convert("1", "%"), + bulkdens_g_cm3 = .data$soil_bulk_density |> + mean() |> + PEcAn.utils::ud_convert("kg m-3", "g cm-3"), + org_C_pct = .data$soil_organic_carbon_stock |> + sum() |> + PEcAn.utils::ud_convert("kg m-2", "g cm-2") |> + (\(x) x / (.data$depth_cm * .data$bulkdens_g_cm3))() |> + PEcAn.utils::ud_convert("1", "%"), + iom_tC_ha = .data$soil_organic_carbon_stock |> + sum() |> + PEcAn.utils::ud_convert("kg m-2", "t ha-1") |> + # Approximation from Falloon et al. 1998, 10.1016/S0038-0717(97)00256-3 + (\(x) 0.049 * x^1.139)() + ) +} diff --git a/models/rothc/R/write.config.RothC.R b/models/rothc/R/write.config.RothC.R index 7c22758da8..d645ecb012 100644 --- a/models/rothc/R/write.config.RothC.R +++ b/models/rothc/R/write.config.RothC.R @@ -160,47 +160,7 @@ write.config.RothC <- function(defaults, trait.values, settings, run.id) { met_in <- utils::read.table(met_path, header = TRUE) n_met <- nrow(met_in) - ## Soil parameters - soil_list <- PEcAn.data.land::pool_ic_netcdf2list( - settings$run$inputs$soil_physics$path - ) - soil_params <- soil_list$vals |> - as.data.frame() |> - # netCDF metadata and PEcAn.data.land::soil.units("soil_depth") both - # say depth should be meters, but the gSSURGO files I've checked are in cm. - # TODO: fix either code or unit labels upstream, - # and check that they are consistent with other sources - # TODO 2: Assumes depth is given to bottom of layer -- is that correct? - dplyr::mutate(depth_cm = soil_list$dims$depth) |> - # TODO this drops layers that extend past bottom - # (eg with depth=23 and 0-10/10-30 layering, would use only 0-10) - # Consider rescaling partial layers - # (Or throwing an error on mismatch and making everyone generate their soil - # files with layers that match model depth?) - dplyr::filter(.data$depth_cm <= model_depth_cm) |> - dplyr::summarize( - # TODO consider weighting by layer thickness? - depth_cm = max(.data$depth_cm), - clay_pct = .data$fraction_of_clay_in_soil |> - mean() |> - PEcAn.utils::ud_convert("1", "%"), - silt_pct = .data$fraction_of_silt_in_soil |> - mean() |> - PEcAn.utils::ud_convert("1", "%"), - bulkdens_g_cm3 = .data$soil_bulk_density |> - mean() |> - PEcAn.utils::ud_convert("kg m-3", "g cm-3"), - org_C_pct = .data$soil_organic_carbon_stock |> - sum() |> - PEcAn.utils::ud_convert("kg m-2", "g cm-2") |> - (\(x) x / (.data$depth_cm * .data$bulkdens_g_cm3))() |> - PEcAn.utils::ud_convert("1", "%"), - iom_tC_ha = .data$soil_organic_carbon_stock |> - sum() |> - PEcAn.utils::ud_convert("kg m-2", "t ha-1") |> - # Approximation from Falloon et al. 1998, 10.1016/S0038-0717(97)00256-3 - (\(x) 0.049 * x^1.139)() - ) + soil_params <- read_soil_physics(settings$run$inputs$soil_physics$path) ## Build string from soil params ## (plus number of timesteps, weirdly snuck into the middle) diff --git a/models/rothc/man/read_soil_physics.Rd b/models/rothc/man/read_soil_physics.Rd new file mode 100644 index 0000000000..9987663485 --- /dev/null +++ b/models/rothc/man/read_soil_physics.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/read_soil_physics.R +\name{read_soil_physics} +\alias{read_soil_physics} +\title{Read soil parameters from a PEcAn soil physics file} +\usage{ +read_soil_physics(path, model_depth = 23) +} +\arguments{ +\item{path}{filepath to a netcdf, +probably read from `settings$run$inputs$soil_physics$path`} + +\item{model_depth}{Soil depth to be simulated, in cm} +} +\value{ +one-row dataframe with columns + `depth_cm`, `clay_pct`, `silt_pct`, + `bulkdens_g_cm3`, `org_C_pct`, `iom_tC_ha` +} +\description{ +Reads values from a netCDF that follows the PEcAn soil standard, + aggregates layers across the depth to be simulated, + and converts to names and units expected by RothC. +} +\details{ +TODO: Currently drops all of any layer that extends below `model_depth`. + It would be better instead to weight all layers by their contribution + to the requested depth. +} diff --git a/models/rothc/tests/testthat/test-read_soil_physics.R b/models/rothc/tests/testthat/test-read_soil_physics.R new file mode 100644 index 0000000000..e3917b4ea6 --- /dev/null +++ b/models/rothc/tests/testthat/test-read_soil_physics.R @@ -0,0 +1,78 @@ + +# Helper: Send test values to a soil netcdf, +# read back using read_soil_physics(). +# Note that it recycles all inputs to common length +nc_roundtrip <- function(depth = c(15, 30, 60), # cm + clay = 1 / 3, # frac + silt = 1 / 3, # frac + bulk = 1350, # kg m-3 + oc = 2, # kg m-2 + model_depth = 30) { # cm + nc <- withr::local_tempfile() + + PEcAn.data.land::soil2netcdf( + soil.data = data.frame( + soil_depth = PEcAn.utils::ud_convert(depth, "cm", "m"), + fraction_of_clay_in_soil = clay, + fraction_of_silt_in_soil = silt, + soil_bulk_density = bulk, + soil_organic_carbon_stock = oc + ), + new.file = nc + ) + + read_soil_physics(nc, model_depth = model_depth) +} + +test_that("one layer", { + res <- nc_roundtrip(depth = 30) + expect_equal(res$clay_pct, 33.3333, tolerance = 1e-3) + expect_equal(res$bulkdens_g_cm3, 1.35) + # (stock in g cm-3) / (bulk * thickness) * 100% + expect_equal(res$org_C_pct, (2 / 10) / (1.35 * 30) * 100, tolerance = 1e-6) + + expect_equal( + nc_roundtrip(depth = 30, model_depth = 30) |> + dplyr::select(-"depth_cm", -"org_C_pct"), + nc_roundtrip(depth = 20, model_depth = 20) |> + dplyr::select(-"depth_cm", -"org_C_pct") + ) + expect_equal( + nc_roundtrip(depth = 30, model_depth = 30) |> + dplyr::select(-"depth_cm", -"iom_tC_ha"), + nc_roundtrip(depth = 20, model_depth = 20, oc = 2 * 20 / 30) |> + dplyr::select(-"depth_cm", -"iom_tC_ha"), + tolerance = 1e-6 + ) +}) + +test_that("two layers", { + res <- nc_roundtrip(depth = c(15, 30), clay = c(0.1, 0.2), oc = c(2, 4)) + expect_equal(res$clay_pct, 15, tolerance = 1e-3) + expect_equal(res$bulkdens_g_cm3, 1.35) + expect_equal(res$org_C_pct, ((2 + 4) / 10) / (1.35 * 30) * 100) +}) + +test_that("n layers", { + expect_equal( + nc_roundtrip(depth = 30), + nc_roundtrip(depth = 1:30, oc = 2 / 30), + tolerance = 1e-6 + ) +}) + +test_that("deeper layers ignored", { + expect_equal( + nc_roundtrip(depth = 30, clay = 0.5, oc = 2), + nc_roundtrip(depth = c(30, 60), clay = c(0.5, 0.1), oc = c(2, 20)), + tolerance = 1e-6 + ) +}) + + +test_that("detects unit errors", { + expect_output( + nc_roundtrip(depth = c(1500, 3000, 6000)), + "Assuming these are mislabeled cm" + ) +}) From bd37829326c02ed04feb32969d790cd0cd0af78c Mon Sep 17 00:00:00 2001 From: Chris Black Date: Thu, 22 Jan 2026 13:14:29 -0800 Subject: [PATCH 12/19] add netcdf2df --- CHANGELOG.md | 2 + base/utils/NAMESPACE | 1 + base/utils/NEWS.md | 6 ++ base/utils/R/netcdf2df.R | 48 ++++++++++++++++ base/utils/tests/testthat/test-netcdf2df.R | 65 ++++++++++++++++++++++ models/rothc/R/read_soil_physics.R | 2 +- 6 files changed, 123 insertions(+), 1 deletion(-) create mode 100644 base/utils/R/netcdf2df.R create mode 100644 base/utils/tests/testthat/test-netcdf2df.R diff --git a/CHANGELOG.md b/CHANGELOG.md index 45e79bff9b..15fc516641 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,8 @@ For more information about this file see also [Keep a Changelog](http://keepacha ## Unreleased ### Added +- New function `PEcAn.utils::netcdf2df()` flattens all dims and vars of a netCDF into a dataframe, + with units attached as an attribute. - New package `PEcAn.RothC` runs the RothC soil carbon model. - Add function `PEcAn.MA::meta_analysis_standalone` to run meta-analysis without database or file IO. - Added Demo 03: Meta Analysis Quarto notebook (`documentation/tutorials/Demo_03_Meta_Analysis/meta_analysis.qmd`) to demonstrate how to perform Bayesian meta-analysis and visualize posterior distributions using pre-generated trait data. diff --git a/base/utils/NAMESPACE b/base/utils/NAMESPACE index 6f50b6ee55..985a0b0b4e 100644 --- a/base/utils/NAMESPACE +++ b/base/utils/NAMESPACE @@ -36,6 +36,7 @@ export(n_leap_day) export(nc_merge_all_sites_by_year) export(nc_write_varfiles) export(need_packages) +export(netcdf2df) export(paste.stats) export(r2bugs.distributions) export(read.output) diff --git a/base/utils/NEWS.md b/base/utils/NEWS.md index e6effc30ee..1f017ca54d 100644 --- a/base/utils/NEWS.md +++ b/base/utils/NEWS.md @@ -1,3 +1,9 @@ +# PEcAn.utils 1.8.2.9000 + +## Added + +* New function `netcdf2df()` reads a netCDF as a dataframe, with units attached as an attribute. It is intended so far for files that are already "basically a rectangle" (dense grids where all variables have the same dimensions and all dimensions are scalars or fully-crossed vectors), but if you try it on files with more complex dimensions please report how it goes. + # PEcAn.utils 1.8.2 ## Added diff --git a/base/utils/R/netcdf2df.R b/base/utils/R/netcdf2df.R new file mode 100644 index 0000000000..df3e6ae544 --- /dev/null +++ b/base/utils/R/netcdf2df.R @@ -0,0 +1,48 @@ +#' Read all variables from a netCDF into a single data frame +#' +#' Reads all dimensions and variables from a netCDF and returns them as +#' a single data frame with one row per cell of the source file's +#' dimension array. +#' Units are also read and are attached to the result as attribute "units" +#' +#' Written mostly for files where all dimensions are 1d vectors, +#' e.g. single-site soil or met files. +#' Many files with more complex dimensions should work +#' (at least for cases where it's clear how to rectangle them), +#' but they are not yet well tested. +#' +#' @param path path to a netcdf file +#' +#' @return data frame with columns for each dim and var of the input file. +#' Units for all of these are attached as an attribute. +#' +#' @author Chris Black +#' @export +#' +netcdf2df <- function(path) { + nc <- ncdf4::nc_open(path) + on.exit(ncdf4::nc_close(nc), add = TRUE) + + dim_vals <- nc$dim |> + sapply(\(x) c(x[["vals"]]), simplify = FALSE) |> + # Load-bearing assumption: + # Dimension listed first in nc$dim is fastest-varying (as for expand.grid) + # TODO verify whether this is reliably true. + expand.grid() + + var_vals <- nc$var |> + names() |> + sapply(\(v) c(ncdf4::ncvar_get(nc, v)), simplify = FALSE) + + if (any(lengths(var_vals) != nrow(dim_vals))) { + PEcAn.logger::logger.error("Not all variables have same length") + } + + res <- cbind(dim_vals, as.data.frame(var_vals)) + attr(res, "units") <- c( + sapply(nc$dim, `[[`, "units"), + sapply(nc$var, `[[`, "units") + ) + + res +} diff --git a/base/utils/tests/testthat/test-netcdf2df.R b/base/utils/tests/testthat/test-netcdf2df.R new file mode 100644 index 0000000000..c646234079 --- /dev/null +++ b/base/utils/tests/testthat/test-netcdf2df.R @@ -0,0 +1,65 @@ +test_that("one dim", { + local_edition(3) + + res <- c("a", "b") |> + example_netcdf(file_path = withr::local_tempfile()) |> + netcdf2df() + + expect_equal(dim(res), c(365, 3)) + expect_equal(colnames(res), c("time", "a", "b")) + expect_equal( + attr(res, "units"), + c(time = "days since 2001-01-01", a = "kg", b = "kg") + ) + expect_equal(res$time, (0L:364L)) +}) + + + +test_that("multiple dims, all but one scalar", { + local_edition(3) + + # TODO does this work if pkg not installed (eg check time)? + # I think testthat may shim system.file + res <- netcdf2df( + system.file("test-data/CRUNCEP.2000.nc", package = "PEcAn.utils") + ) + + expect_equal(dim(res), c(366 * 4, 11)) + expect_equal( + colnames(res), + c( + "latitude", "longitude", "time", + "air_temperature", "surface_downwelling_longwave_flux_in_air", + "air_pressure", "surface_downwelling_shortwave_flux_in_air", + "eastward_wind", "northward_wind", + "specific_humidity", "precipitation_flux" + ) + ) + + # Check units. NB some of these are nonstandard; + # Key point is they should match what the file reports + expect_equal( + attr(res, "units"), + c( + latitude = "degree_north", + longitude = "degree_east", + time = "days since 2000-01-01T00:00:00Z", + air_temperature = "Kelvin", + surface_downwelling_longwave_flux_in_air = "W/m2", + air_pressure = "Pascal", + surface_downwelling_shortwave_flux_in_air = "W/m2", + eastward_wind = "m/s", + northward_wind = "m/s", + specific_humidity = "g/g", + precipitation_flux = "kg/m2/s" + ) + ) + # dims + expect_equal(res$time, seq(0, 365 + 3 / 4, 1 / 4) + 1 / 8, tolerance = 1e-12) + expect_equal(unique(res$latitude), 45.25, tolerance = 1e-6) + expect_equal(unique(res$longitude), -84.75, tolerance = 1e-6) + expect_equal(mean(res$air_temperature), 278.798636, tolerance = 1e-6) +}) + +# TODO test with more than one non-degenerate dimension diff --git a/models/rothc/R/read_soil_physics.R b/models/rothc/R/read_soil_physics.R index 64521577ef..6a6dc43465 100644 --- a/models/rothc/R/read_soil_physics.R +++ b/models/rothc/R/read_soil_physics.R @@ -19,7 +19,7 @@ #' @importFrom rlang .data .env #' read_soil_physics <- function(path, model_depth = 23) { - soil_list <- PEcAn.data.land::pool_ic_netcdf2list(path) + soil_list <- PEcAn.utils::netcdf2df(path) # Input is documented to be in meters, # but providing cm instead seems to be a _very_ common error; From 065643aaf4750212d214b595b548c630cfb2442c Mon Sep 17 00:00:00 2001 From: Chris Black Date: Thu, 22 Jan 2026 15:03:09 -0800 Subject: [PATCH 13/19] use netcdf2df instead of import from data.land --- base/utils/man/netcdf2df.Rd | 31 +++++++++++++ models/rothc/R/read_soil_physics.R | 45 ++++++++++++------- .../tests/testthat/test-read_soil_physics.R | 33 +++++++++----- 3 files changed, 81 insertions(+), 28 deletions(-) create mode 100644 base/utils/man/netcdf2df.Rd diff --git a/base/utils/man/netcdf2df.Rd b/base/utils/man/netcdf2df.Rd new file mode 100644 index 0000000000..54829f3ad4 --- /dev/null +++ b/base/utils/man/netcdf2df.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/netcdf2df.R +\name{netcdf2df} +\alias{netcdf2df} +\title{Read all variables from a netCDF into a single data frame} +\usage{ +netcdf2df(path) +} +\arguments{ +\item{path}{path to a netcdf file} +} +\value{ +data frame with columns for each dim and var of the input file. +Units for all of these are attached as an attribute. +} +\description{ +Reads all dimensions and variables from a netCDF and returns them as +a single data frame with one row per cell of the source file's +dimension array. +Units are also read and are attached to the result as attribute "units" +} +\details{ +Written mostly for files where all dimensions are 1d vectors, +e.g. single-site soil or met files. +Many files with more complex dimensions should work +(at least for cases where it's clear how to rectangle them), +but they are not yet well tested. +} +\author{ +Chris Black +} diff --git a/models/rothc/R/read_soil_physics.R b/models/rothc/R/read_soil_physics.R index 6a6dc43465..ba7b5c4592 100644 --- a/models/rothc/R/read_soil_physics.R +++ b/models/rothc/R/read_soil_physics.R @@ -19,26 +19,29 @@ #' @importFrom rlang .data .env #' read_soil_physics <- function(path, model_depth = 23) { - soil_list <- PEcAn.utils::netcdf2df(path) + soil_vals <- netcdf2df(path) + soil_units <- as.list(attr(soil_vals, "units")) # Input is documented to be in meters, # but providing cm instead seems to be a _very_ common error; # might as well try to handle it here - if (all(soil_list$dims$depth < 10)) { - depth_cm <- PEcAn.utils::ud_convert(soil_list$dims$depth, "m", "cm") - } else { + if (tolower(soil_units$depth) %in% c("m", "meters") + && any(soil_vals$depth >= 10)) { PEcAn.logger::logger.warn( - "Soil depths should be in meters, but found values >= 10", + "Soil depths reported to be in meters, but found values >= 10", "in file", path, "Assuming these are mislabeled cm and treating them as such." ) - depth_cm <- soil_list$dims$depth + soil_units$depth <- "cm" } - soil_list$vals |> - as.data.frame() |> - # TODO: Assumes depth is given to bottom of layer -- is that correct? - dplyr::mutate(depth_cm = .env$depth_cm) |> - # TODO this drops layers that extend past bottom + + soil_vals |> + dplyr::mutate( + depth_cm = .data$depth |> + PEcAn.utils::ud_convert(soil_units$depth, "cm") + ) |> + # TODO 1: Assumes depth is given to bottom of layer -- is that correct? + # TODO 2: this drops layers that extend past bottom # (eg with depth=23 and 0-10/10-30 layering, would use only 0-10) # Consider rescaling partial layers # (Or throwing an error on mismatch and making everyone generate their soil @@ -49,22 +52,32 @@ read_soil_physics <- function(path, model_depth = 23) { depth_cm = max(.data$depth_cm), clay_pct = .data$fraction_of_clay_in_soil |> mean() |> - PEcAn.utils::ud_convert("1", "%"), + PEcAn.utils::ud_convert(soil_units$fraction_of_clay_in_soil, "%"), silt_pct = .data$fraction_of_silt_in_soil |> mean() |> - PEcAn.utils::ud_convert("1", "%"), + PEcAn.utils::ud_convert(soil_units$fraction_of_silt_in_soil, "%"), bulkdens_g_cm3 = .data$soil_bulk_density |> mean() |> - PEcAn.utils::ud_convert("kg m-3", "g cm-3"), + PEcAn.utils::ud_convert(soil_units$soil_bulk_density, "g cm-3"), org_C_pct = .data$soil_organic_carbon_stock |> sum() |> - PEcAn.utils::ud_convert("kg m-2", "g cm-2") |> + PEcAn.utils::ud_convert( + soil_units$soil_organic_carbon_stock, + "g cm-2" + ) |> (\(x) x / (.data$depth_cm * .data$bulkdens_g_cm3))() |> PEcAn.utils::ud_convert("1", "%"), iom_tC_ha = .data$soil_organic_carbon_stock |> sum() |> - PEcAn.utils::ud_convert("kg m-2", "t ha-1") |> + PEcAn.utils::ud_convert( + soil_units$soil_organic_carbon_stock, + "t ha-1" + ) |> # Approximation from Falloon et al. 1998, 10.1016/S0038-0717(97)00256-3 (\(x) 0.049 * x^1.139)() ) } + +# an internal wrapper to allow stubbing the function out under test +# without affecting code outside the package. +netcdf2df <- PEcAn.utils::netcdf2df diff --git a/models/rothc/tests/testthat/test-read_soil_physics.R b/models/rothc/tests/testthat/test-read_soil_physics.R index e3917b4ea6..13d61128ce 100644 --- a/models/rothc/tests/testthat/test-read_soil_physics.R +++ b/models/rothc/tests/testthat/test-read_soil_physics.R @@ -8,20 +8,29 @@ nc_roundtrip <- function(depth = c(15, 30, 60), # cm bulk = 1350, # kg m-3 oc = 2, # kg m-2 model_depth = 30) { # cm - nc <- withr::local_tempfile() + with_mocked_bindings( + { + read_soil_physics("path/ignored", model_depth = model_depth) + }, + netcdf2df = function(...) { + dat <- data.frame( + depth = PEcAn.utils::ud_convert(depth, "cm", "m"), + fraction_of_clay_in_soil = clay, + fraction_of_silt_in_soil = silt, + soil_bulk_density = bulk, + soil_organic_carbon_stock = oc + ) + attr(dat, "units") <- c( + depth = "meters", + fraction_of_clay_in_soil = "1", + fraction_of_silt_in_soil = "1", + soil_bulk_density = "kg m-3", + soil_organic_carbon_stock = "kg m-2" + ) - PEcAn.data.land::soil2netcdf( - soil.data = data.frame( - soil_depth = PEcAn.utils::ud_convert(depth, "cm", "m"), - fraction_of_clay_in_soil = clay, - fraction_of_silt_in_soil = silt, - soil_bulk_density = bulk, - soil_organic_carbon_stock = oc - ), - new.file = nc + dat + } ) - - read_soil_physics(nc, model_depth = model_depth) } test_that("one layer", { From 8c1ce9b01ccecd2c08734c67c0c5e975718c3778 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Tue, 23 Jun 2026 06:23:25 -0700 Subject: [PATCH 14/19] Update models/rothc/R/read_soil_physics.R --- models/rothc/R/read_soil_physics.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/rothc/R/read_soil_physics.R b/models/rothc/R/read_soil_physics.R index ba7b5c4592..ff1f0421b2 100644 --- a/models/rothc/R/read_soil_physics.R +++ b/models/rothc/R/read_soil_physics.R @@ -30,7 +30,7 @@ read_soil_physics <- function(path, model_depth = 23) { PEcAn.logger::logger.warn( "Soil depths reported to be in meters, but found values >= 10", "in file", path, - "Assuming these are mislabeled cm and treating them as such." + "Assuming the units have been mislabeled and treating all depths as cm." ) soil_units$depth <- "cm" } From d7c10389b829b27008724b362d32a0abf92d09d0 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Tue, 23 Jun 2026 06:48:03 -0700 Subject: [PATCH 15/19] error on empty result --- models/rothc/R/read_soil_physics.R | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/models/rothc/R/read_soil_physics.R b/models/rothc/R/read_soil_physics.R index ff1f0421b2..b01f0be81f 100644 --- a/models/rothc/R/read_soil_physics.R +++ b/models/rothc/R/read_soil_physics.R @@ -35,7 +35,7 @@ read_soil_physics <- function(path, model_depth = 23) { soil_units$depth <- "cm" } - soil_vals |> + soil_vals <- soil_vals |> dplyr::mutate( depth_cm = .data$depth |> PEcAn.utils::ud_convert(soil_units$depth, "cm") @@ -46,7 +46,15 @@ read_soil_physics <- function(path, model_depth = 23) { # Consider rescaling partial layers # (Or throwing an error on mismatch and making everyone generate their soil # files with layers that match model depth?) - dplyr::filter(.data$depth_cm <= model_depth) |> + dplyr::filter(.data$depth_cm <= model_depth) + + if (nrow(soil_vals) == 0) { + PEcAn.logger::logger.error( + "No soil layers <=", model_depth, "found in file", sQuote(path) + ) + } + + soil_vals |> dplyr::summarize( # TODO consider weighting by layer thickness? depth_cm = max(.data$depth_cm), From ee2503ca773c1c5f80d1f13a0e6a31abe4fa1d81 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Tue, 23 Jun 2026 07:40:34 -0700 Subject: [PATCH 16/19] define %||% --- models/rothc/R/PEcAn.RothC-package.R | 15 +++++++++++++ models/rothc/R/model2netcdf.RothC.R | 1 - models/rothc/R/read_soil_physics.R | 2 -- models/rothc/man/PEcAn.RothC-package.Rd | 28 +++++++++++++++++++++++++ 4 files changed, 43 insertions(+), 3 deletions(-) create mode 100644 models/rothc/R/PEcAn.RothC-package.R create mode 100644 models/rothc/man/PEcAn.RothC-package.Rd diff --git a/models/rothc/R/PEcAn.RothC-package.R b/models/rothc/R/PEcAn.RothC-package.R new file mode 100644 index 0000000000..3c4f9006ba --- /dev/null +++ b/models/rothc/R/PEcAn.RothC-package.R @@ -0,0 +1,15 @@ +#' PEcAn.RothC-package +#' +#' PEcAn support for the RothC soil carbon model +#' +#' @importFrom rlang .data .env +#' @keywords internal +"_PACKAGE" + + +# Define null-coalescing operator unless already available in base (R >= 4.4). +# (We could also import it from rlang, but why add a whole import for one +# line of code?) +if (!exists("%||%", envir=.BaseNamespaceEnv)) { + `%||%` <- function(x, y) if (is.null(x)) y else x +} diff --git a/models/rothc/R/model2netcdf.RothC.R b/models/rothc/R/model2netcdf.RothC.R index 57bc86ad45..219c5233b1 100644 --- a/models/rothc/R/model2netcdf.RothC.R +++ b/models/rothc/R/model2netcdf.RothC.R @@ -6,7 +6,6 @@ #' @param start_date Start time of the simulation #' @param end_date End time of the simulation #' @export -#' @importFrom rlang .data .env #' #' @author Chris Black model2netcdf.RothC <- function(outdir, sitelat, sitelon, start_date, end_date) { diff --git a/models/rothc/R/read_soil_physics.R b/models/rothc/R/read_soil_physics.R index b01f0be81f..21a0f72940 100644 --- a/models/rothc/R/read_soil_physics.R +++ b/models/rothc/R/read_soil_physics.R @@ -16,8 +16,6 @@ #' `depth_cm`, `clay_pct`, `silt_pct`, #' `bulkdens_g_cm3`, `org_C_pct`, `iom_tC_ha` #' -#' @importFrom rlang .data .env -#' read_soil_physics <- function(path, model_depth = 23) { soil_vals <- netcdf2df(path) soil_units <- as.list(attr(soil_vals, "units")) diff --git a/models/rothc/man/PEcAn.RothC-package.Rd b/models/rothc/man/PEcAn.RothC-package.Rd new file mode 100644 index 0000000000..19e93d40ca --- /dev/null +++ b/models/rothc/man/PEcAn.RothC-package.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PEcAn.RothC-package.R +\docType{package} +\name{PEcAn.RothC-package} +\alias{PEcAn.RothC} +\alias{PEcAn.RothC-package} +\title{PEcAn.RothC-package} +\description{ +PEcAn support for the RothC soil carbon model +} +\seealso{ +Useful links: +\itemize{ + \item \url{https://pecanproject.github.io} + \item Report bugs at \url{https://github.com/PecanProject/pecan/issues} +} + +} +\author{ +\strong{Maintainer}: Chris Black \email{chris@ckblack.org} + +Other contributors: +\itemize{ + \item Rothamsted Research [copyright holder] +} + +} +\keyword{internal} From 5b6910b9a1db29eb7a9b7fd9f05970f951adf848 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Tue, 23 Jun 2026 09:01:13 -0700 Subject: [PATCH 17/19] fix test wording --- models/rothc/tests/testthat/test-read_soil_physics.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/rothc/tests/testthat/test-read_soil_physics.R b/models/rothc/tests/testthat/test-read_soil_physics.R index 13d61128ce..c30a85d0ae 100644 --- a/models/rothc/tests/testthat/test-read_soil_physics.R +++ b/models/rothc/tests/testthat/test-read_soil_physics.R @@ -82,6 +82,6 @@ test_that("deeper layers ignored", { test_that("detects unit errors", { expect_output( nc_roundtrip(depth = c(1500, 3000, 6000)), - "Assuming these are mislabeled cm" + "Assuming the units have been mislabeled" ) }) From 4cf83282b92f3bd4f3341996ab23070b491afd71 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Tue, 23 Jun 2026 10:50:14 -0700 Subject: [PATCH 18/19] ignore docs folder --- models/rothc/.Rbuildignore | 1 + 1 file changed, 1 insertion(+) diff --git a/models/rothc/.Rbuildignore b/models/rothc/.Rbuildignore index 3f270fb469..b7966b40fe 100644 --- a/models/rothc/.Rbuildignore +++ b/models/rothc/.Rbuildignore @@ -1,3 +1,4 @@ Dockerfile model_info.json inst/workflow_example +^docs$ From 3f12f45c279dccfb1ab48b9fc2a39dd9d0680b2b Mon Sep 17 00:00:00 2001 From: Chris Black Date: Tue, 23 Jun 2026 23:02:42 -0700 Subject: [PATCH 19/19] gitignore --- models/rothc/.gitignore | 1 + 1 file changed, 1 insertion(+) create mode 100644 models/rothc/.gitignore diff --git a/models/rothc/.gitignore b/models/rothc/.gitignore new file mode 100644 index 0000000000..77f12ae2e5 --- /dev/null +++ b/models/rothc/.gitignore @@ -0,0 +1 @@ +docs/