Add write.events.SIPNET function

CHANGELOG.md

@@ -1,9 +1,10 @@
 # Change Log
+
 All notable changes are kept in this file. All changes made should be added to the section called
 `Unreleased`. Once a new release is made this file will be updated to create a new `Unreleased`
 section for the next release.
 
-	For more information about this file see also [Keep a Changelog](http://keepachangelog.com/) .
+For more information about this file see also [Keep a Changelog](http://keepachangelog.com/) .
 
 ## Unreleased
 
@@ -15,6 +16,7 @@ section for the next release.
 - Added `AmeriFlux_met_ensemble()` function with ERA5 fallback for AmeriFlux meteorological data processing and ensemble generation
 - Added `all_site_nc_merge_by_year()` and `single_site_nc_merge()` functions to merge netCDF files across ensembles and sites from pecan model netCDF outputs.
 - Added parallel mode for the entire SDA workflow.
+- `write.events.SIPNET()` to generate SIPNET `events.in` files from a `events.json` file.
 - Included all relevant carbon pools (`ROOT_BIOMASS`, `AG_BIOMASS`, `SOIL_STOCK`, `LIT_BIOMASS`) in BADM-based IC extraction; excluded non-pool variables like `SOIL_CHEM`.
 - Added explicit support for `LIT_BIOMASS` to fully utilize **BADM** biomass capabilities.
 - Added `test-IC_BADM_Utilities.R` to validate BADM initial condition extraction and processing

docker/depends/pecan_package_dependencies.csv

@@ -131,6 +131,7 @@
 "IDPmisc","*","modules/assim.batch","Imports",FALSE
 "itertools","*","modules/assim.sequential","Suggests",FALSE
 "jsonlite","*","base/remote","Imports",FALSE
+"jsonlite","*","models/sipnet","Imports",FALSE
 "jsonlite","*","models/stics","Imports",FALSE
 "jsonlite","*","modules/data.atmosphere","Imports",FALSE
 "jsonlite","*","modules/data.remote","Suggests",FALSE

models/sipnet/DESCRIPTION

@@ -14,6 +14,7 @@ URL: https://pecanproject.github.io
 BugReports: https://github.com/PecanProject/pecan/issues
 Imports:
     dplyr,
+    jsonlite,
     lubridate (>= 1.6.0),
     ncdf4 (>= 1.15),
     PEcAn.data.atmosphere,

models/sipnet/NAMESPACE

@@ -9,6 +9,7 @@ export(sample.IC.SIPNET)
 export(split_inputs.SIPNET)
 export(veg2model.SIPNET)
 export(write.config.SIPNET)
+export(write.events.SIPNET)
 export(write_restart.SIPNET)
 importFrom(dplyr,"%>%")
 importFrom(rlang,.data)

models/sipnet/NEWS.md

@@ -1,3 +1,9 @@
+# Unreleased
+
+## Added
+
+* `write.events.SIPNET()` to generate SIPNET `events.in` files from a `events.json` file.
+
 # PEcAn.SIPNET 1.9.1
 
 ## Changed

models/sipnet/R/write.events.SIPNET.R

@@ -0,0 +1,118 @@
+## TODO:
+## - integrate call into write.configs.SIPNET
+## - parameterize planting allocation fractions
+## - make sure files are written in correct output directory
+## - map crops associated w/ planting and harvest --> PFTs; this will need to be handled separate from events.in
+#' Write SIPNET events.in files from a PEcAn events.json
+#'
+#' Reads a single PEcAn events.json containing one or more site objects and
+#' writes one SIPNET `events.in` file per site. Events are translated according to [SIPNET's `events.in`
+#' specification](https://pecanproject.github.io/sipnet/parameters/#agronomic-events).
+#' The writer expects inputs to already match the PEcAn MVP schema v0.1.0 naming and units where applicable.
+#'
+#' @details
+#' - Supported `event_type` values: `tillage`, `planting`, `fertilization`,
+#'   `irrigation`, `harvest`.
+#' - Units translated from PEcAn standard_vars to SIPNET events.in specification:
+#'   `kg/m^2` to `g/m^2`; irrigation `amount_mm` to `cm`.
+#' - Planting allocation uses fixed internal parameters. Future work should use the same values
+#'   that are written to `sipnet.parms` (e.g. after integrating this into `write.configs.SIPNET`)
+#'
+#' @param events_json character. Path to an `events.json` file containing an
+#'   array of site objects with `site_id`, optional `pft`, and `events`.
+#' @param outdir character. Output directory where per-site `events-<site>.in`
+#'   files are written.
+#'
+#' @return Invisibly, a vector of files written.
+#'
+#' @examples
+#' # Example with two events for a single site
+#' tmp <- withr::local_tempfile(fileext = ".json")
+#' site <- list(
+#'     site_id = "EX1",
+#'     events = list(
+#'         list(event_type = "tillage", date = "2022-02-04", tillage_eff_0to1 = 0.2),
+#'         list(event_type = "planting", date = "2022-02-19", leaf_c_kg_m2 = 0.01)
+#'     )
+#' )
+#' jsonlite::write_json(list(site), tmp, auto_unbox = TRUE)
+#' outdir <- withr::local_tempdir()
+#' files <- write.events.SIPNET(tmp, outdir)
+#' files
+#'
+#' @export
+write.events.SIPNET <- function(events_json, outdir) {
+    # TODO add overwrite argument
+    x <- jsonlite::fromJSON(events_json, simplifyVector = FALSE)
+    # allow a single site events.json that does not have a site_id 
+    site_objs <- if (!is.null(x$site_id)) list(x) else x
+    files_written <- vector()
+
+    leafAllocation <- 0.50
+    woodAllocation <- 0.15
+    fineRootAllocation <- 0.10
+    coarseRootAllocation <- 0.25
+
+    # Unit conversion helpers
+    kg2g <- as.numeric(PEcAn.utils::ud_convert(1, "kg", "g")) # 1000
+    mm2cm <- as.numeric(PEcAn.utils::ud_convert(1, "mm", "cm")) # 0.1
+
+    # For each site, build event time series and write file
+    for (site in site_objs) {
+        sid <- site$site_id
+        evs <- site$events
+        # Order by date and build lines
+        dates <- as.Date(vapply(evs, function(e) e$date, character(1)))
+        ord <- order(dates)
+        lines <- character()
+        for (e in evs[ord]) {
+            d <- as.Date(e$date)
+            year <- as.integer(format(d, "%Y"))
+            day <- as.integer(format(d, "%j"))
+            type <- e$event_type
+            if (type == "tillage") {
+                f <- if (is.null(e$tillage_eff_0to1)) 0 else e$tillage_eff_0to1
+                # TODO: consider validating up front against schema rather than here
+                lines <- c(lines, sprintf("%d  %d  till  %s", year, day, f))
+            } else if (type == "planting") {
+                # infer total planted biomass from leaf pool and allocation fraction
+                leaf_g <- as.numeric(if (is.null(e$leaf_c_kg_m2)) 0 else e$leaf_c_kg_m2) * kg2g
+                total_g <- if (leafAllocation > 0) leaf_g / leafAllocation else leaf_g
+                wood_g <- woodAllocation * total_g
+                fr_g <- fineRootAllocation * total_g
+                cr_g <- coarseRootAllocation * total_g
+                lines <- c(lines, sprintf("%d  %d  plant  %s %s %s %s", year, day, leaf_g, wood_g, fr_g, cr_g))
+            } else if (type == "fertilization") {
+                orgN_g <- as.numeric(if (is.null(e$org_n_kg_m2)) 0 else e$org_n_kg_m2) * kg2g
+                orgC_g <- as.numeric(if (is.null(e$org_c_kg_m2)) 0 else e$org_c_kg_m2) * kg2g
+                nh4_g <- as.numeric(if (is.null(e$nh4_n_kg_m2)) 0 else e$nh4_n_kg_m2) * kg2g
+                no3_g <- as.numeric(if (is.null(e$no3_n_kg_m2)) 0 else e$no3_n_kg_m2) * kg2g
+                minN_g <- nh4_g + no3_g
+                lines <- c(lines, sprintf("%d  %d  fert   %s %s %s", year, day, orgN_g, orgC_g, minN_g))
+            } else if (type == "irrigation") {
+                amt_cm <- as.numeric(if (is.null(e$amount_mm)) 0 else e$amount_mm) * mm2cm
+                method_code <- if (is.null(e$method) || e$method == "soil") 1 else 0
+                lines <- c(lines, sprintf("%d  %d  irrig  %s %s", year, day, amt_cm, method_code))
+            } else if (type == "harvest") {
+                abv_rem <- e$frac_above_removed_0to1 %||% 0
+                blw_rem <- e$frac_below_removed_0to1 %||% 0
+                abv_lit <- e$frac_above_to_litter_0to1 %||% (1.0 - abv_rem)
+                blw_lit <- e$frac_below_to_litter_0to1 %||% (1.0 - blw_rem)
+                lines <- c(
+                    lines,
+                    sprintf(
+                        "%d  %d harv   %s %s %s %s",
+                        year, day,
+                        abv_rem, blw_rem,
+                        abv_lit, blw_lit
+                    )
+                )
+            }
+        }
+        dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
+        fp <- file.path(outdir, sprintf("events-%s.in", sid))
+        writeLines(lines, fp)
+        files_written <- c(files_written, fp)
+    }
+    invisible(files_written)
+}

models/sipnet/tests/testthat/event_fixtures/events_site1.json

@@ -0,0 +1,43 @@
+[
+  {
+    "site_id": ["EX1"],
+    "events": [
+      {
+        "event_type": ["tillage"],
+        "date": ["2022-02-04"],
+        "tillage_eff_0to1": [0.2]
+      },
+      {
+        "event_type": ["tillage"],
+        "date": ["2022-02-09"],
+        "tillage_eff_0to1": [0.1]
+      },
+      {
+        "event_type": ["irrigation"],
+        "date": ["2022-02-09"],
+        "amount_mm": [50],
+        "method": ["soil"]
+      },
+      {
+        "event_type": ["fertilization"],
+        "date": ["2022-02-09"],
+        "org_n_kg_m2": [0],
+        "org_c_kg_m2": [0],
+        "nh4_n_kg_m2": [0.01]
+      },
+      {
+        "event_type": ["planting"],
+        "date": ["2022-02-19"],
+        "leaf_c_kg_m2": [0.01]
+      },
+      {
+        "event_type": ["harvest"],
+        "date": ["2022-09-07"],
+        "frac_above_removed_0to1": [0.1],
+        "frac_below_removed_0to1": [0.0],
+        "frac_above_to_litter_0to1": [0.0],
+        "frac_below_to_litter_0to1": [0.0]
+      }
+    ]
+  }
+]

models/sipnet/tests/testthat/event_fixtures/events_site1_site2.json

@@ -0,0 +1,38 @@
+[
+  {
+    "site_id": "S1",
+    "pft": "PFT",
+    "events": [
+      {
+        "event_type": "tillage",
+        "date": "2022-01-15",
+        "tillage_eff_0to1": 0.1
+      },
+      {
+        "event_type": "harvest",
+        "date": "2022-09-01",
+        "frac_above_removed_0to1": 0.2,
+        "frac_below_removed_0to1": 0.0,
+        "frac_above_to_litter_0to1": 0.0,
+        "frac_below_to_litter_0to1": 0.0
+      }
+    ]
+  },
+  {
+    "site_id": "S2",
+    "pft": "PFT",
+    "events": [
+      {
+        "event_type": "planting",
+        "date": "2022-03-01",
+        "leaf_c_kg_m2": 0.01
+      },
+      {
+        "event_type": "irrigation",
+        "date": "2022-03-10",
+        "amount_mm": 10,
+        "method": "soil"
+      }
+    ]
+  }
+]

models/sipnet/tests/testthat/test-write.events.SIPNET.R

@@ -0,0 +1,37 @@
+context("write.events.SIPNET")
+
+# Helper to remove excess whitespace
+norm <- function(x) gsub("\\s+", " ", trimws(x))
+
+testthat::test_that("write.events.SIPNET produces expected lines", {
+    ev_json1 <- testthat::test_path("event_fixtures/events_site1.json")
+    outdir <- withr::local_tempdir()
+    files <- write.events.SIPNET(ev_json1, outdir)
+    expect_length(files, 1)
+    got <- readLines(files[1])
+    expected <- c(
+        "2022  35  till  0.2",
+        "2022  40  till  0.1",
+        "2022  40  irrig  5 1",
+        "2022  40  fert   0 0 10",
+        "2022  50  plant  10 3 2 5",
+        "2022  250 harv   0.1 0 0 0"
+    )
+    expect_equal(norm(got), norm(expected))
+    # TODO determine What's generating the whitespace differences and eliminate use of norm()
+})
+
+testthat::test_that("write.events.SIPNET handles multi-site events.json (one file per site)", {
+    ev_json2 <- testthat::test_path("event_fixtures/events_site1_site2.json")
+    outdir <- withr::local_tempdir()
+    files <- write.events.SIPNET(ev_json2, outdir)
+    testthat::expect_length(files, 2)
+    testthat::expect_true(all(file.exists(files)))
+    # quick sanity checks for each site's first/last event ordering
+    got1 <- readLines(files[grepl("events-S1\\.in$", files)])
+    got2 <- readLines(files[grepl("events-S2\\.in$", files)])
+    testthat::expect_true(startsWith(norm(got1[1]), "2022 15 till"))
+    testthat::expect_true(startsWith(norm(tail(got1, 1)), "2022 244 harv"))
+    testthat::expect_true(startsWith(norm(got2[1]), "2022 60 plant"))
+    testthat::expect_true(startsWith(norm(tail(got2, 1)), "2022 69 irrig"))
+})

models/sipnet/tests/testthat/test.met2model.R

@@ -16,8 +16,11 @@ add_gaps_to_nc <- function(src_nc, gapped_nc,
 }
 
 test_that("Met conversion runs without error", {
-  nc_path <- system.file("test-data", "CRUNCEP.2000.nc",
-                         package = "PEcAn.utils")
+  nc_path <- system.file(
+    "test-data",
+    "CRUNCEP.2000.nc",
+    package = "PEcAn.utils"
+  )
   in.path <- dirname(nc_path)
   in.prefix <- "CRUNCEP"
   start_date <- "2000-01-01"