Merge pull request #3788 from infotroph/rothc-opts

Rothc: Soil file support + read options from settings file

Author: infotroph

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

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
 
 ## Fixed

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/man/netcdf2df.Rd

@@ -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

base/utils/tests/testthat/test-netcdf2df.R

@@ -1,3 +1,4 @@
 Dockerfile
 model_info.json
 inst/workflow_example
+^docs$

models/rothc/.Rbuildignore

@@ -0,0 +1 @@
+docs/

models/rothc/.gitignore

@@ -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/Dockerfile

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