SIPNET workflow for restarting with events

CHANGELOG.md

@@ -39,6 +39,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha
   `MCMCpack`, `mvtnorm`, `neonUtilities`, `neonstore`, `PEcAn.benchmark`,
   `PEcAn.visualization`, `rjags`, `sirt`, and `sp` from `Imports` to
   `Suggests` (@omkarrr2533, #3599).
+- Management events specified via `events.json` are now required to specify a crop code for each planting event, so that models can know when to restart with a different PFT (#3828, #3836).
 
 
 

base/remote/R/check_model_run.R

@@ -10,9 +10,9 @@ check_model_run <- function(out, stop.on.error = TRUE) {
     success <- FALSE
     msg <- paste0("Model run aborted with the following error:\n", out)
     if (stop.on.error) {
-      PEcAn.logger::logger.severe(msg)
+      PEcAn.logger::logger.severe(msg, wrap = FALSE)
     } else {
-      PEcAn.logger::logger.error(msg)
+      PEcAn.logger::logger.error(msg, wrap = FALSE)
     }
   } else {
     success <- TRUE

book_source/03_topical_pages/02_pecan_standards.Rmd

@@ -97,6 +97,86 @@ See the [Soil Data] section on more into on creating a standard soil data file.
 
 See the [Vegetation Data] section on more info on creating a standard vegetation data file
 
+### Management events
+
+Externally imposed changes in system state, such as planting and harvest or tillage of crop systems, are represented as a sequence of "events" that are each defined by a date, event type, and one or more event-specific properties that describe the kind and magnitude of change imposed. PEcAn passes the same representation of a given event to all models, but interpretation of all properties -- including dates and event types -- is necessarily model-specific and could include ignoring an event entirely.
+
+Event types are defined by the PEcAn events spec, which is currently provided as a JSON schema file bundled into the PEcAn.data.land package (TODO consider exposing this in a more linkable format?). Valid events files are stored as JSON, with each events file containing an array of sites. Each site contains an id, a schema version, and an array of events, usually recorded in date order. Each event has a date, an event type, and one or more values from the required and optional properties that are documented below for each event type. Any event is also allowed to contain arbitrary additional properties not mentioned in the spec; PEcAn will pass these on to models with no further validation.
+
+The currently supported event types are: planting, harvest, irrigation, fertilization, and tillage. Each of these is described briefly below. You will note that all of them are focused on agronomic management of croplands; future work will extend this framework to include disturbance of unmanaged systems such as fire, inundation or insect outbreaks.
+
+
+#### planting
+
+Addition of live biomass that is expected to grow according to the model's plant growth rules. The current spec expects a "planting" to be a transplantation of leafed-out seedlings; models which represent planting as addition of unsprouted seeds might choose to adjust the date to account for germination time.
+
+Required properties:
+
+* `crop_code`: a machine-readable identifier that specifies what type of plant was added. Should preferably be short and use only characters that are safe in filenames.
+* `leaf_c_kg_m2`: Amount of leaf biomass added. 
+
+Optional properties:
+
+* `wood_c_kg_m2`: Biomass added as aboveground stems. Assumed zero if not specified.
+* `fine_root_c_kg_m2`: Biomass added to belowground roots. Assumed zero if not specified.
+* `cultivar`: Identifier for finer-grainer plant type identifiers that will be meaningful to models with detailed growth parameterizations but likely ignored by PFT-based models.
+* `crop_display`: A human-readable name for the added plant type.
+
+#### harvest
+
+Removal of biomass from the living pool, either by transferring it to dead (litter) biomass or by removing it from the system. Specified as a fraction of live biomass.
+
+Required properties:
+* `frac_above_removed_0to1`: Fraction of existing aboveground biomass removed from the system.
+
+Optional properties:
+
+* `frac_below_removed_0to1`: Fraction of existing beloweground biomass removed from the system. Assumed zero if not specified.
+* `frac_above_to_litter_0to1`, `frac_below_to_litter_0to1`: fraction of existing above- or below-ground live biomass converted to dead litter. Assumed zero if not specified.
+* `crop_code`: A machine_readable identifier that specifies what type of plant was harvested.
+* `cultivar`: Identifier for finer-grainer plant type identifiers that will be meaningful to models with detailed growth parameterizations but likely ignored by PFT-based models.
+* `crop_display`: A human-readable name for the added plant type.
+
+#### irrigation
+
+Addition of water to the system beyond that specified as precipitation in met files. Note that only the total amount of water is specified and not the application rate; assumptions about the duration of the event will be model-specific.
+
+Required properties:
+
+* `method`: How water was applied. One of "soil", "canopy", "flood".
+* `amount_mm`: Amount of water applied.
+
+Optional properties:
+
+* `immed_evap_frac_0to1`: fraction of the applied water that evaporates before reaching the soil water pool.
+
+#### fertilization
+
+Addition of exogenous nutrients, including non-crop carbon (e.g. manure, compost). Nutrients other than N and C are not currently included but may be added in future schemas.
+
+Required properties:
+
+* at least one of `org_c_kg_m2`, `nh4_n_kg_m2`, `no3_n_kg_m2`, which are respectively the amount of organic C, ammonium N, and nitrate N added.
+
+Optional properties:
+
+* any of `org_c_kg_m2`, `nh4_n_kg_m2`, `no3_n_kg_m2` not already counted as required. All are assumed zero if not specified.
+* `org_n_kg_m2`: the amount of N added in organic compounds. Assumed zero if not specified. If specifying, you probably want to specify `org_c_kg_m2` as well but the schema does not enforce this.
+
+#### tillage
+
+Physical disturbance of the soil and litter pools. Models differ greatly in their representation of tillage effects and this spec intentionally limits its representation to a simple intensity score whose maximum should be interpreted as "the most complete possible mixing/disturbance." Whether this is sufficient detail will be re-evaluated in future revisions.
+
+Required properties:
+
+* `tillage_eff_0to1`: Relative tillage intensity, with 0 meaning no disturbance at all and 1 meaning the most complete disturbance possible.
+
+Optional properties:
+
+* `intensity_category`: string giving additional information on the type or amount of tillage performed.
+* `depth_m`: Depth of soil affected by the event.
+
+
 ## Output Standards
 
 * created by `model2netcdf` functions

models/sipnet/NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(combine_sipnet_out)
 export(mergeNC)
 export(met2model.SIPNET)
 export(model2netcdf.SIPNET)

models/sipnet/R/combine_sipnet_out.R

@@ -0,0 +1,36 @@
+#!/usr/bin/env Rscript
+
+#' Combine a bunch of `sipnet.out` files into a single file
+#'
+#' @param directory Parent directory to search for `sipnet.out` files. You must
+#' provide either this or an explicit vector of `files`.
+#' @param outfile File to which combined `sipnet.out` will be written
+#' @param files Optional vector of paths to files to combine. All must be
+#' readable with [read_sipnet_out()]. If `NULL` (default), looks for files
+#' named `sipnet.out` in `directory` and its subdirectories, recursively.
+#'
+#' @return `outfile` (path to output file), invisibly
+#' @export
+combine_sipnet_out <- function(directory, outfile, files = NULL) {
+  if (missing(directory) && is.null(files)) {
+    PEcAn.logger::logger.severe("Must provide either `directory` or `files`")
+  }
+  if (is.null(files)) {
+    # NOTE that this expects file paths (including parent directories) to be
+    # lexicographically sorted. For the common case of segmented SIPNET runs,
+    # the parent directories are named `segment_001`, `segment_002`, etc., so
+    # this will work automatically.
+    # If you don't want to make this assumption, or have a custom sort order,
+    # pass `files` directly to this function.
+    files <- sort(list.files(directory, "sipnet\\.out", full.names = TRUE, recursive = TRUE))
+  }
+  if (length(files) == 0) {
+    PEcAn.logger::logger.severe("No files provided; nothing to combine.")
+  }
+  flist <- lapply(files, read_sipnet_out)
+  combined <- do.call(rbind.data.frame, flist)
+  # Mimic the SIPNET fixed-width right-aligned format
+  combined_fwf <- format(combined, justify = "right")
+  utils::write.table(combined_fwf, outfile, row.names = FALSE, col.names = TRUE, quote = FALSE, sep = "  ")
+  invisible(outfile)
+}

models/sipnet/R/model2netcdf.SIPNET.R

@@ -13,6 +13,7 @@
 #' }
 #' @export mergeNC
 #' @name mergeNC
+#' @importFrom rlang .data
 #' @source https://github.com/RS-eco/processNC/blob/main/R/mergeNC.R
 mergeNC <- function(
     ##title<< Aggregate data in netCDF files
@@ -61,38 +62,7 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
                                 overwrite = FALSE, conflict = FALSE) {
   ### Read in model output in SIPNET format
   sipnet_out_file <- file.path(outdir, prefix)
-  # SIPNET v1 had a "Notes" comment line before the header; v2 removed it.
-  # if the first line starts with "year", there is no Notes line.
-  first_line <- readLines(sipnet_out_file, n = 1)
-  skip_n <- if (grepl("^year", first_line)) 0 else 1
-  # Temporary workaround until
-  # https://github.com/PecanProject/sipnet/issues/304 is resolved.
-  sipnet_output <- tryCatch({
-    utils::read.table(sipnet_out_file, header = TRUE, skip = skip_n, sep = "")
-  }, error = function(err) {
-    PEcAn.logger::logger.warn(
-      "Failed to read using `read.table`. ",
-      "Trying to parse output manually."
-    )
-    raw_lines <- readLines(sipnet_out_file)
-    raw_header <- raw_lines[[1 + skip_n]]
-    raw_body <- utils::tail(raw_lines, -(1 + skip_n))
-    # SIPNET output is right-aligned with the column names in the header.
-    # We use this to figure out where the numbers end if there are no spaces.
-    token_matches <- gregexpr("\\S+", raw_header, perl = TRUE)
-    proc_header <- regmatches(raw_header, token_matches)[[1]]
-    col_ends <- token_matches[[1]] + attr(token_matches[[1]], "match.length") - 1
-    col_starts <- c(1, utils::head(col_ends, -1) + 1)
-    col_widths <- col_ends - col_starts + 1
-    result <- utils::read.fwf(
-      textConnection(raw_body),
-      widths = col_widths,
-      col.names = proc_header,
-      na.strings = c("nan", "-nan")
-    )
-    result[] <- lapply(result, as.numeric)
-    result
-  })
+  sipnet_output <- read_sipnet_out(sipnet_out_file)
   #sipnet_output_dims <- dim(sipnet_output)
 
   ### Determine number of years and output timestep

models/sipnet/R/read_sipnet_out.R

@@ -0,0 +1,40 @@
+#' Read `sipnet.out` file to `data.frame`
+#'
+#' @param sipnet_out_file Path to `sipnet.out` file
+#'
+#' @return `data.frame` of SIPNET output
+read_sipnet_out <- function(sipnet_out_file) {
+  # SIPNET v1 had a "Notes" comment line before the header; v2 removed it.
+  # if the first line starts with "year", there is no Notes line.
+  first_line <- readLines(sipnet_out_file, n = 1)
+  skip_n <- if (grepl("^year", first_line)) 0 else 1
+  # Temporary workaround until
+  # https://github.com/PecanProject/sipnet/issues/304 is resolved.
+  sipnet_output <- tryCatch({
+    utils::read.table(sipnet_out_file, header = TRUE, skip = skip_n, sep = "")
+  }, error = function(err) {
+    PEcAn.logger::logger.warn(
+      "Failed to read using `read.table`. ",
+      "Trying to parse output manually."
+    )
+    raw_lines <- readLines(sipnet_out_file)
+    raw_header <- raw_lines[[1 + skip_n]]
+    raw_body <- utils::tail(raw_lines, -(1 + skip_n))
+    # SIPNET output is right-aligned with the column names in the header.
+    # We use this to figure out where the numbers end if there are no spaces.
+    token_matches <- gregexpr("\\S+", raw_header, perl = TRUE)
+    proc_header <- regmatches(raw_header, token_matches)[[1]]
+    col_ends <- token_matches[[1]] + attr(token_matches[[1]], "match.length") - 1
+    col_starts <- c(1, utils::head(col_ends, -1) + 1)
+    col_widths <- col_ends - col_starts + 1
+    result <- utils::read.fwf(
+      textConnection(raw_body),
+      widths = col_widths,
+      col.names = proc_header,
+      na.strings = c("nan", "-nan")
+    )
+    result[] <- lapply(result, as.numeric)
+    result
+  })
+  sipnet_output
+}