diff --git a/CHANGELOG.md b/CHANGELOG.md index e8ff1d55362..323e14467be 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,9 @@ 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. - Added `inst/ilamb/` pipeline in PEcAn.benchmark to convert downscaled SDA reanalysis GeoTIFFs into ILAMB-compatible CF netCDF for carbon-cycle benchmarking (#4019). - Added PEcAn.PEPRMT model, including a demo run with example data - Add `format_try_for_ma()` and `try_trait_mapping()` to `PEcAn.data.remote` to convert trait data from the external TRY database into the tabular format required by the PEcAn meta-analysis module (#3717). diff --git a/base/utils/NAMESPACE b/base/utils/NAMESPACE index 6f50b6ee555..985a0b0b4e9 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 95b268ca041..55cf9fa82ca 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 ## Fixed diff --git a/base/utils/R/netcdf2df.R b/base/utils/R/netcdf2df.R new file mode 100644 index 00000000000..df3e6ae5443 --- /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/man/netcdf2df.Rd b/base/utils/man/netcdf2df.Rd new file mode 100644 index 00000000000..54829f3ad4c --- /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/base/utils/tests/testthat/test-netcdf2df.R b/base/utils/tests/testthat/test-netcdf2df.R new file mode 100644 index 00000000000..c646234079f --- /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/.Rbuildignore b/models/rothc/.Rbuildignore index 3f270fb469c..b7966b40fe6 100644 --- a/models/rothc/.Rbuildignore +++ b/models/rothc/.Rbuildignore @@ -1,3 +1,4 @@ Dockerfile model_info.json inst/workflow_example +^docs$ diff --git a/models/rothc/.gitignore b/models/rothc/.gitignore new file mode 100644 index 00000000000..77f12ae2e58 --- /dev/null +++ b/models/rothc/.gitignore @@ -0,0 +1 @@ +docs/ diff --git a/models/rothc/Dockerfile b/models/rothc/Dockerfile index a9b5af2dd24..2811e0a0f3d 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/PEcAn.RothC-package.R b/models/rothc/R/PEcAn.RothC-package.R new file mode 100644 index 00000000000..3c4f9006ba3 --- /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 57bc86ad45c..219c5233b12 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 new file mode 100644 index 00000000000..21a0f729404 --- /dev/null +++ b/models/rothc/R/read_soil_physics.R @@ -0,0 +1,89 @@ +#' 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` +#' +read_soil_physics <- function(path, model_depth = 23) { + 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 (tolower(soil_units$depth) %in% c("m", "meters") + && any(soil_vals$depth >= 10)) { + PEcAn.logger::logger.warn( + "Soil depths reported to be in meters, but found values >= 10", + "in file", path, + "Assuming the units have been mislabeled and treating all depths as cm." + ) + soil_units$depth <- "cm" + } + + soil_vals <- 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 + # files with layers that match 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), + clay_pct = .data$fraction_of_clay_in_soil |> + mean() |> + PEcAn.utils::ud_convert(soil_units$fraction_of_clay_in_soil, "%"), + silt_pct = .data$fraction_of_silt_in_soil |> + mean() |> + PEcAn.utils::ud_convert(soil_units$fraction_of_silt_in_soil, "%"), + bulkdens_g_cm3 = .data$soil_bulk_density |> + mean() |> + 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( + 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( + 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/R/rothc_varname_map.R b/models/rothc/R/rothc_varname_map.R index 3adf9f66602..1b153f84204 100644 --- a/models/rothc/R/rothc_varname_map.R +++ b/models/rothc/R/rothc_varname_map.R @@ -21,3 +21,19 @@ rothc_varname_map <- dplyr::tribble( "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 ) + + +# 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, + # 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 452396cfad4..d645ecb012b 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 @@ -125,32 +125,57 @@ 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) + + # 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 - # (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) - ## Soil parameters + 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) - # TODO: read from run$inputs$soil_physics soil_param_string <- paste( - "23.4 23.0 3.0041", # clay_pct, depth_cm, 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) - "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 + 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) @@ -176,31 +201,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) } @@ -210,4 +233,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 +} diff --git a/models/rothc/inst/RothC_input_template.dat b/models/rothc/inst/RothC_input_template.dat index 761ee96086b..d6bb7ca2c00 100644 --- a/models/rothc/inst/RothC_input_template.dat +++ b/models/rothc/inst/RothC_input_template.dat @@ -2,9 +2,9 @@ # @START_YEAR@-@START_MONTH@-@START_DAY@ to @END_YEAR@-@END_MONTH@-@END_DAY@ # # opt_RMmoist opt_SMDbare - @OPT_RMMOIST@ @OPT_SDDBARE@ -# (%) (cm) (t C/ha) (-) (%) (g/m3) (%) (-) -# clay depth iom nsteps siltper BD OC minRM_Moist + @OPT_RMMOIST@ @OPT_SMDBARE@ +# (%) (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 diff --git a/models/rothc/inst/workflow_example/README.md b/models/rothc/inst/workflow_example/README.md index f4e90287859..adf4a664d64 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 00000000000..82d4215e9f3 --- /dev/null +++ b/models/rothc/inst/workflow_example/fetch_soil_data.R @@ -0,0 +1,70 @@ +#!/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)) + +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( + .f = PEcAn.data.land::extract_soil_gssurgo, + size = args$n_ens, + grid_size = 5, + .options = furrr::furrr_options(seed = TRUE) + ) diff --git a/models/rothc/inst/workflow_example/run.sh b/models/rothc/inst/workflow_example/run.sh index 53cf9fee890..8bebb5fba15 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 3f198466117..f22607e8fbb 100644 --- a/models/rothc/inst/workflow_example/template.xml +++ b/models/rothc/inst/workflow_example/template.xml @@ -29,20 +29,24 @@ + sampling - --> + - 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 + 0.1 + 30 diff --git a/models/rothc/inst/workflow_example/xml_build.R b/models/rothc/inst/workflow_example/xml_build.R index 8e0389b59ed..d02918fb2bb 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, diff --git a/models/rothc/man/PEcAn.RothC-package.Rd b/models/rothc/man/PEcAn.RothC-package.Rd new file mode 100644 index 00000000000..19e93d40ca5 --- /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} diff --git a/models/rothc/man/read_soil_physics.Rd b/models/rothc/man/read_soil_physics.Rd new file mode 100644 index 00000000000..99876634856 --- /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 00000000000..c30a85d0ae6 --- /dev/null +++ b/models/rothc/tests/testthat/test-read_soil_physics.R @@ -0,0 +1,87 @@ + +# 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 + 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" + ) + + dat + } + ) +} + +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 the units have been mislabeled" + ) +})