Rothc: Soil file support + read options from settings file

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.
 - Add function `qsub_sda()` for submitting SDA batch jobs by splitting a large number of sites into multiple small groups of sites (#3634).
 - 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.

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)

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

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
+}

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

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 \

models/rothc/R/read_soil_physics.R

@@ -0,0 +1,83 @@
+#' 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_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 |>
+    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) |>
+    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

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
+)