WIP: RothC coupler

Makefile

@@ -4,7 +4,7 @@ NCPUS ?= 1
 BASE := logger utils db settings visualization qaqc remote workflow
 
 MODELS := basgra biocro clm45 dalec dvmdostem ed fates gday jules linkages \
-				ldndc lpjguess maat maespa sibcasa sipnet stics template
+				ldndc lpjguess maat maespa rothc sibcasa sipnet stics template
 
 MODULES := allometry assim.batch assim.sequential benchmark \
 				 data.atmosphere data.land data.remote \

Makefile.depends

@@ -21,6 +21,7 @@ $(call depends,models/linkages): | .install/base/db .install/base/logger .instal
 $(call depends,models/lpjguess): | .install/base/logger .install/base/remote .install/base/utils
 $(call depends,models/maat): | .install/base/logger .install/base/remote .install/base/settings .install/base/utils .install/modules/data.atmosphere
 $(call depends,models/maespa): | .install/base/logger .install/base/remote .install/base/utils .install/modules/data.atmosphere
+$(call depends,models/rothc): | .install/base/logger .install/base/utils
 $(call depends,models/sibcasa): | .install/base/logger
 $(call depends,models/sipnet): | .install/base/logger .install/base/remote .install/base/utils .install/modules/data.atmosphere .install/modules/data.land
 $(call depends,models/stics): | .install/base/logger .install/base/remote .install/base/settings .install/base/utils

docker/depends/pecan_package_dependencies.csv

@@ -168,6 +168,7 @@
 "lubridate",">= 1.6.0","models/lpjguess","Imports",FALSE
 "lubridate",">= 1.6.0","models/maat","Imports",FALSE
 "lubridate",">= 1.6.0","models/maespa","Imports",FALSE
+"lubridate",">= 1.6.0","models/rothc","Imports",FALSE
 "lubridate",">= 1.6.0","models/sipnet","Imports",FALSE
 "lubridate",">= 1.6.0","modules/assim.batch","Imports",FALSE
 "lubridate",">= 1.6.0","modules/assim.sequential","Imports",FALSE
@@ -231,6 +232,7 @@
 "ncdf4","*","models/basgra","Imports",FALSE
 "ncdf4","*","models/dvmdostem","Imports",FALSE
 "ncdf4","*","models/ldndc","Imports",FALSE
+"ncdf4","*","models/rothc","Imports",FALSE
 "ncdf4","*","models/sibcasa","Imports",FALSE
 "ncdf4","*","models/stics","Imports",FALSE
 "ncdf4","*","modules/assim.sequential","Imports",FALSE
@@ -338,6 +340,7 @@
 "PEcAn.logger","*","models/lpjguess","Imports",TRUE
 "PEcAn.logger","*","models/maat","Imports",TRUE
 "PEcAn.logger","*","models/maespa","Imports",TRUE
+"PEcAn.logger","*","models/rothc","Imports",TRUE
 "PEcAn.logger","*","models/sibcasa","Imports",TRUE
 "PEcAn.logger","*","models/sipnet","Imports",TRUE
 "PEcAn.logger","*","models/stics","Imports",TRUE
@@ -429,6 +432,7 @@
 "PEcAn.utils",">= 1.4.8","models/basgra","Imports",TRUE
 "PEcAn.utils",">= 1.4.8","models/dvmdostem","Imports",TRUE
 "PEcAn.utils",">= 1.4.8","models/ldndc","Imports",TRUE
+"PEcAn.utils",">= 1.4.8","models/rothc","Imports",TRUE
 "PEcAn.utils",">= 1.4.8","models/stics","Imports",TRUE
 "PEcAn.utils",">= 1.4.8","models/template","Imports",TRUE
 "PEcAn.visualization","*","modules/assim.sequential","Suggests",TRUE
@@ -536,6 +540,7 @@
 "roxygen2","== 7.3.2","models/lpjguess","Roxygen",FALSE
 "roxygen2","== 7.3.2","models/maat","Roxygen",FALSE
 "roxygen2","== 7.3.2","models/maespa","Roxygen",FALSE
+"roxygen2","== 7.3.2","models/rothc","Roxygen",FALSE
 "roxygen2","== 7.3.2","models/sibcasa","Roxygen",FALSE
 "roxygen2","== 7.3.2","models/sipnet","Roxygen",FALSE
 "roxygen2","== 7.3.2","models/stics","Roxygen",FALSE
@@ -571,6 +576,7 @@
 "sp","*","modules/data.land","Imports",FALSE
 "sp","*","modules/data.remote","Imports",FALSE
 "stats","*","base/qaqc","Imports",FALSE
+"stats","*","models/rothc","Imports",FALSE
 "stats","*","models/sipnet","Imports",FALSE
 "stats","*","modules/allometry","Imports",FALSE
 "stats","*","modules/assim.batch","Imports",FALSE
@@ -630,6 +636,7 @@
 "testthat",">= 2.0.0","base/utils","Suggests",FALSE
 "testthat",">= 2.0.0","models/biocro","Suggests",FALSE
 "testthat",">= 2.0.0","modules/benchmark","Suggests",FALSE
+"testthat",">= 3.0.0","models/rothc","Suggests",FALSE
 "testthat",">= 3.0.0","models/sibcasa","Suggests",FALSE
 "testthat",">= 3.0.4","base/qaqc","Suggests",FALSE
 "testthat",">= 3.1.7","modules/data.atmosphere","Suggests",FALSE
@@ -664,6 +671,7 @@
 "utils","*","models/ed","Imports",FALSE
 "utils","*","models/linkages","Imports",FALSE
 "utils","*","models/lpjguess","Imports",FALSE
+"utils","*","models/rothc","Imports",FALSE
 "utils","*","modules/allometry","Imports",FALSE
 "utils","*","modules/assim.batch","Imports",FALSE
 "utils","*","modules/assim.sequential","Suggests",FALSE
@@ -682,6 +690,7 @@
 "withr","*","base/workflow","Suggests",FALSE
 "withr","*","models/basgra","Suggests",FALSE
 "withr","*","models/ed","Suggests",FALSE
+"withr","*","models/rothc","Suggests",FALSE
 "withr","*","models/sibcasa","Suggests",FALSE
 "withr","*","models/sipnet","Suggests",FALSE
 "withr","*","modules/allometry","Suggests",FALSE

models/rothc/.Rbuildignore

@@ -0,0 +1,2 @@
+Dockerfile
+model_info.json

models/rothc/DESCRIPTION

@@ -0,0 +1,30 @@
+Package: PEcAn.RothC
+Type: Package
+Title: PEcAn Package for Integration of the RothC Model
+Version: 0.0.0.9000
+Authors@R: c(
+    person("Chris", "Black", role = c("aut", "cre"), email = "[email protected]"),
+    person("Rothamsted Research", role = c("cph"))
+            )
+Description: This module provides functions to link the Rothamstead soil carbon model "RothC" to PEcAn.
+  Uses RothC 1.0.0, made available by Rothamstead Research under the Apache License 2.0. TODO: Determine whether we are including a copy of the RothC sources in this package or connecting to externally installed versions, get copyright notice and license requirements met.
+Depends: R (>= 4.1)
+Imports:
+    dplyr,
+    lubridate (>= 1.6.0),
+    ncdf4,
+    PEcAn.logger,
+    PEcAn.utils (>= 1.4.8),
+    rlang,
+    stats,
+    utils
+Suggests:
+    testthat (>= 3.0.0),
+    withr
+SystemRequirements: RothC_Py
+OS_type: unix
+License: BSD_3_clause + file LICENSE
+Copyright: Authors
+Encoding: UTF-8
+RoxygenNote: 7.3.2
+Config/testthat/edition: 3

models/rothc/Dockerfile

@@ -0,0 +1,40 @@
+# this needs to be at the top, what version are we building
+ARG IMAGE_VERSION="latest"
+
+
+# ----------------------------------------------------------------------
+# BUILD PECAN FOR MODEL
+# ----------------------------------------------------------------------
+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
+
+RUN mkdir -p /tmp/rothc \
+    && curl -sSL https://github.com/Rothamsted-Models/RothC_Code/archive/refs/tags/v${MODEL_VERSION}.tar.gz \
+        | tar xzf - -C /tmp/rothc \
+    && cd /tmp/rothc/RothC_Code-${MODEL_VERSION} \
+    && gfortran RothC.for Shell.for -o rothc_bin \
+    && mv rothc_bin /usr/local/bin/RothC_v${MODEL_VERSION} \
+    && rm -rf /tmp/rothc
+
+# TODO hard-coding this path is probably wrong
+# Setup model_info file
+# @VERSION@ is replaced with model version in the model_info.json file
+# @BINARY@ is replaced with model binary in the model_info.json file
+COPY model_info.json /work/model.json
+RUN sed -i -e "s/@VERSION@/${MODEL_VERSION}/g" \
+           -e "s#@BINARY@#/usr/local/bin/RothC_v${MODEL_VERSION}#g" /work/model.json

models/rothc/LICENSE

@@ -0,0 +1,3 @@
+YEAR: 2024
+COPYRIGHT HOLDER: PEcAn Project
+ORGANIZATION: PEcAn Project, authors affiliations

models/rothc/NAMESPACE

@@ -0,0 +1,9 @@
+# Generated by roxygen2: do not edit by hand
+
+export(met2model.RothC)
+export(model2netcdf.RothC)
+export(read_restart.ModelName)
+export(write.config.RothC)
+export(write_restart.ModelName)
+importFrom(rlang,.data)
+importFrom(rlang,.env)

models/rothc/NEWS.md

@@ -0,0 +1,3 @@
+# PEcAn.RothC 0.0.0.9000
+
+Initial development version

models/rothc/R/met2model.RothC.R

@@ -0,0 +1,180 @@
+#' Extract monthly weather from CF file for input to RothC
+#'
+#' Input files need to be named `<in.path>/<in.prefix>.YYYY.nc`
+#'
+#' Output files are named `<outfolder>/<in.prefix>.YY-mm.YY-mm.dat`
+#' with one line per month and columns for temperature, rainfall,
+#' and evaporation.
+#'
+#' Note that the created file contains only weather data and not any of the
+#' soil or management data needed for RothC's single combined input file.
+#' See `write.config.RothC()` for assembly into a model-ready RothC_input.dat`.
+#'
+#' @param in.path path on disk where CF files live
+#' @param in.prefix prefix for each file
+#' @param outfolder location where model specific output is written.
+#' @param start_date,end_date When to start and end output.
+#'  Specify as exact dates, but output will be padded to whole months.
+#' @param overwrite logical: replace output files if they already exist?
+#' @return data frame summarizing file metadata
+#' @export
+#' @author Chris Black
+met2model.RothC <- function(in.path,
+                            in.prefix,
+                            outfolder,
+                            start_date,
+                            end_date,
+                            overwrite = FALSE) {
+
+  PEcAn.logger::logger.info("START met2model.RothC")
+
+  start_date <- as.Date(start_date)
+  end_date <- as.Date(end_date)
+  start_year <- strftime(start_date, "%Y")
+  end_year <- strftime(end_date, "%Y")
+  year_regex <- paste(start_year:end_year, collapse = "|")
+
+  if (grepl("\\.nc$", in.prefix)) {
+    # Assume it's the full filename rather than a prefix
+    # NB also means we assume it contains the whole requested date range
+    name_pattern <- in.prefix
+  } else {
+    name_pattern <- paste0(in.prefix, "\\.(", year_regex, ")\\.nc$")
+  }
+
+  nc_files <- list.files(in.path, pattern = name_pattern, full.names = TRUE)
+
+  if (length(nc_files) == 0) {
+    PEcAn.logger::logger.severe(
+      "No files found matching ", in.prefix,
+      "for years", start_year, ":", end_year,
+      "; cannot process data."
+    )
+  }
+
+  # TODO complain (fail?) here if some but not all years found
+
+  out_filename <- paste(
+    in.prefix,
+    strftime(start_date, "%Y-%m"),
+    strftime(end_date, "%Y-%m"),
+    "dat",
+    sep = "."
+  )
+  out_path <- file.path(outfolder, out_filename)
+  results <- data.frame(file = out_path,
+                        host = Sys.getenv(
+                          "FQDN",
+                          unset = Sys.info()[["nodename"]]
+                        ),
+                        mimetype = "text/tab-separated-values",
+                        formatname = "RothC.dat",
+                        startdate = start_date,
+                        enddate = end_date,
+                        dbfile.name = out_filename,
+                        stringsAsFactors = FALSE)
+  PEcAn.logger::logger.info("internal results")
+  PEcAn.logger::logger.info(results)
+
+  if (file.exists(out_path) && !overwrite) {
+    PEcAn.logger::logger.debug(
+      "File '", out_path, "' already exists, skipping to next file."
+    )
+    return(invisible(results))
+  }
+
+  if (!file.exists(outfolder)) {
+    dir.create(outfolder)
+  }
+
+  met <-  nc_files |>
+    lapply(
+      read_nc,
+      varnames = c("air_temperature", "precipitation_flux", "specific_humidity")
+    ) |>
+    do.call(what = "rbind")
+
+  # TODO probably need more care with partial months here:
+  # we check if data extends to start/end, but not whether that includes enough
+  # days to treat as a whole month.
+  # e.g. if start_date = YYYY-01-31 and data starts YYYY-01-30 ->
+  # current code will aggregate those two days as if they were all of January.
+  # Consider failing if >n days missing in any output month?
+  first_month <- lubridate::floor_date(start_date, unit = "month")
+  last_month <- lubridate::ceiling_date(end_date, unit = "month")
+  met <- met[(met$timestamp >=  first_month) & (met$timestamp < last_month), ]
+  if (as.Date(min(met$timestamp)) > start_date
+      || as.Date(max(met$timestamp)) < end_date) {
+    PEcAn.logger::logger.severe(
+      "input (",
+      paste(range(met$timestamp), collapse = " to "),
+      ") does not cover requested time window (",
+      start_date, "to", end_date, ")"
+    )
+  }
+
+  met$year <- lubridate::year(met$timestamp)
+  met$month <- lubridate::month(met$timestamp)
+  met$Tmp <- met$air_temperature |>
+    PEcAn.utils::ud_convert("K", "degC")
+  met$Rain <- 0# TODO... sum up to convert from flux to accumulation, right?
+  met$Evap <- 0# TODO... how to convert Qair to pan evaporation?
+
+  met_monthly <- merge(
+    stats::aggregate(met, Tmp ~ year + month, mean),
+    stats::aggregate(met, cbind(Rain, Evap) ~ year + month, sum),
+    sort = FALSE # would treat months as strings; sort as numbers below instead
+  )
+  met_monthly <- met_monthly[order(met_monthly$year, met_monthly$month), ]
+
+  utils::write.table(
+    # as.data.frame to write integer columns as ints not floats
+    x = format(as.data.frame(met_monthly), digits = 4),
+    file = out_path,
+    quote = FALSE,
+    sep = "\t",
+    row.names = FALSE,
+    col.names = TRUE
+  )
+
+  results
+}
+
+
+
+
+
+# slurp named vars from one PEcAn nc into a dataframe with timestamp
+#
+# TODO could read other dimensions if present too, but consider if worth it --
+# maybe this function is better for files where only the time dimension varies
+#
+# if vars = NULL, read all of them
+read_nc <- function(ncfile, varnames = NULL) {
+
+  nc <- ncdf4::nc_open(ncfile)
+  on.exit(ncdf4::nc_close(nc), add = TRUE)
+
+  timestamps <- PEcAn.utils::cf2datetime(
+    nc$dim$time$vals,
+    nc$dim$time$units
+  )
+
+  if (is.null(varnames)) {
+    varnames <- names(nc$var)
+  }
+
+  var_values <- lapply(
+    varnames,
+    ncdf4::ncvar_get,
+    nc = nc
+  )
+
+  # todo handle this case (multi-loc files?)
+  stopifnot(all(sapply(var_values, length) == length(timestamps)))
+
+  var_values |>
+    stats::setNames(varnames) |>
+    as.data.frame() |>
+    transform(timestamp = timestamps)
+}