calval workbook to events.json pipeline

.gitignore

@@ -47,3 +47,10 @@ po/*~
 
 # RStudio Connect folder
 rsconnect/
+
+# XLSX snapshots — kept in Google Docs as the data-entry UI and attached
+# to GitHub releases for record-keeping; not tracked in the repo.
+*.xlsx
+
+# Pipeline outputs
+data/

R/detect_missing_events.R

@@ -0,0 +1,86 @@
+#' Inject irrigation events for irrigated cropping systems that lack them
+#'
+#' Some cal/val datasets (notably White/Salinas 2020) record plantings,
+#' harvests, tillage and fertilization but never write a row for
+#' irrigation, even though the crops involved are clearly irrigated.
+#' Per David's plan: detect those gaps and inject placeholder
+#' `irrigation` events into the per-treatment events list before the
+#' priors stage fills in `amount_mm` and `method`.
+#'
+#' MVP heuristic — for each treatment:
+#'   1. If it already has any `irrigation` events, leave it alone.
+#'   2. Otherwise, find pairs of cash-crop plantings and harvests
+#'      (`cover_crop` not set / FALSE) and inject one placeholder
+#'      irrigation event at the midpoint of each window.
+#'   3. The placeholder carries a `prior_filled` tag so consumers can
+#'      see the field came from gap-filling, not source data.
+#'
+#' Distributing the right *count* of irrigation events across the
+#' planting/harvest window (using the `n_events` prior) is left for a
+#' follow-up that will also consume `seasonal_total_mm / n_events` for
+#' the per-event amount. For MVP, one placeholder per window is enough
+#' to make the events.json schema-valid; the priors stage fills
+#' `amount_mm` from `seasonal_total_mm` directly.
+#'
+#' @param events List of events for one treatment (from
+#'   [map_events_for_treatment()]).
+#' @return The events list, possibly with new `irrigation` entries.
+#' @export
+inject_missing_irrigations <- function(events) {
+  if (length(events) == 0) return(events)
+
+  has_irrigation <- any(vapply(events, function(e) identical(e$event_type, "irrigation"),
+                               logical(1)))
+  if (has_irrigation) return(events)
+
+  is_planting <- vapply(events, function(e) identical(e$event_type, "planting"),
+                        logical(1))
+  is_harvest  <- vapply(events, function(e) identical(e$event_type, "harvest"),
+                        logical(1))
+  is_cover    <- vapply(events, function(e) isTRUE(e$cover_crop), logical(1))
+
+  cash_plantings <- which(is_planting & !is_cover)
+  cash_harvests  <- which(is_harvest  & !is_cover)
+  if (length(cash_plantings) == 0 || length(cash_harvests) == 0) return(events)
+
+  # Pair each planting with the next harvest of the same crop, where
+  # "next" means earliest harvest after the planting date with a
+  # matching crop_display (best-effort; falls back to date order).
+  injected <- list()
+  for (pi in cash_plantings) {
+    p <- events[[pi]]
+    p_date <- as.Date(p$date)
+    p_crop <- p$crop_display %||% p$crop_code
+
+    after <- cash_harvests[vapply(cash_harvests, function(hi) {
+      h <- events[[hi]]
+      h_crop <- h$crop_display %||% h$crop_code
+      isTRUE(as.Date(h$date) > p_date) &&
+        (is.null(p_crop) || is.null(h_crop) || p_crop == h_crop)
+    }, logical(1))]
+    if (length(after) == 0) next
+
+    h <- events[[after[which.min(as.Date(vapply(after, function(i) events[[i]]$date,
+                                                character(1))))]]]
+    mid <- as.Date(floor((as.numeric(as.Date(p$date)) +
+                          as.numeric(as.Date(h$date))) / 2),
+                   origin = "1970-01-01")
+
+    injected[[length(injected) + 1L]] <- list(
+      event_type = "irrigation",
+      date       = format(mid, "%Y-%m-%d"),
+      source     = p$source %||% h$source,
+      injected_by = "detect_missing_events:cash_crop_window_midpoint",
+      cover_crop_window = FALSE
+    )
+  }
+
+  if (length(injected) == 0) return(events)
+
+  # Append in chronological order with the rest. Stable sort by date.
+  combined <- c(events, injected)
+  ord <- order(vapply(combined, function(e) as.Date(e$date), as.Date(NA)))
+  combined[ord]
+}
+
+`%||%` <- function(a, b) if (is.null(a) || (is.character(a) && !nzchar(a))) b else a

R/fill_with_priors.R

@@ -0,0 +1,90 @@
+#' Fill schema-required fields that the cal/val data does not specify.
+#'
+#' The cal/val workbook does not record values for a handful of fields
+#' that the events.json schema requires for certain event types
+#' (`leaf_c_kg_m2` for plantings, `frac_above_removed_0to1` for
+#' harvests, `tillage_eff_0to1` for tillage, and `amount_mm` / `method`
+#' for irrigations). Per David's MVP plan, this function samples those
+#' fields from prior distributions defined in Akash's
+#' `management_priors.yaml` (vendored at `inst/akash_priors/`).
+#'
+#' Sampling uses a per-treatment seed derived from the treatment name so
+#' the same input always yields the same output (matters for diffing
+#' and CI). When the priors are upgraded to ensembles, callers will
+#' override `seed` to vary per ensemble member.
+#'
+#' @param events A list of event objects (already produced by
+#'   [map_events_for_treatment()]).
+#' @param priors The list returned by Akash's `load_priors()`.
+#' @param seed Integer; seeds R's RNG before sampling so output is
+#'   reproducible.
+#' @return The events list with previously-missing required fields
+#'   populated.
+#' @export
+fill_with_priors <- function(events, priors, seed = 0L) {
+  set.seed(seed)
+  for (i in seq_along(events)) {
+    events[[i]] <- .fill_event(events[[i]], priors)
+  }
+  events
+}
+
+# Mapping from event_type to (practice_path, parameter_name).
+# practice_path is a vector of keys to descend into the priors list.
+.PRIOR_FIELD_MAP <- list(
+  planting = list(
+    field        = "leaf_c_kg_m2",
+    practice_path = c("crop_baselines", "processing_tomato",
+                      "events", "planting"),
+    param_key    = "leaf_c_kg_m2"
+  ),
+  harvest = list(
+    field        = "frac_above_removed_0to1",
+    practice_path = c("practices", "harvest_grain", "parameters"),
+    param_key    = "frac_above_removed_0to1"
+  ),
+  tillage = list(
+    field        = "tillage_eff_0to1",
+    practice_path = c("practices", "conventional_tillage", "parameters"),
+    param_key    = "tillage_eff_0to1"
+  ),
+  irrigation = list(
+    field        = "amount_mm",
+    practice_path = c("practices", "irrigation_sprinkler", "parameters"),
+    param_key    = "seasonal_total_mm"   # divided by n_events later
+  )
+)
+
+.fill_event <- function(event, priors) {
+  spec <- .PRIOR_FIELD_MAP[[event$event_type]]
+  if (is.null(spec)) return(event)            # nothing to fill
+  if (!is.null(event[[spec$field]])) return(event)   # already set
+
+  dist_spec <- .descend(priors, c(spec$practice_path, spec$param_key))
+  if (is.null(dist_spec)) return(event)
+
+  sample <- sample_distribution(dist_spec, n = 1)
+  event[[spec$field]] <- as.numeric(sample)
+
+  # Tag the event so consumers can tell which fields came from priors.
+  source_log <- if (!is.null(event$source)) event$source else ""
+  prior_note <- paste0(spec$field, "<-prior:",
+                       paste(spec$practice_path, collapse = "/"))
+  event$prior_filled <- if (is.null(event$prior_filled)) prior_note
+                        else paste(event$prior_filled, prior_note, sep = "; ")
+
+  # Irrigation also needs `method`; pick canopy for sprinkler systems.
+  if (event$event_type == "irrigation" && is.null(event$method)) {
+    event$method <- "canopy"
+  }
+  event
+}
+
+# Walk a nested list along a key path, returning NULL if any step fails.
+.descend <- function(x, keys) {
+  for (k in keys) {
+    if (is.null(x) || !k %in% names(x)) return(NULL)
+    x <- x[[k]]
+  }
+  x
+}

R/map_event.R

@@ -0,0 +1,187 @@
+#' Parse "key1=val1; key2=val2; ..." strings from attributes_keyvalue
+#'
+#' Returns a named character vector of values, or an empty named character
+#' vector if the input is missing/empty.
+#' @keywords internal
+parse_kv <- function(s) {
+  if (is.null(s) || length(s) == 0) return(stats::setNames(character(0), character(0)))
+  if (is.na(s) || !nzchar(s)) return(stats::setNames(character(0), character(0)))
+  pairs <- strsplit(s, "\\s*;\\s*")[[1]]
+  pairs <- pairs[nzchar(pairs)]
+  if (length(pairs) == 0) return(stats::setNames(character(0), character(0)))
+  parts <- strsplit(pairs, "=", fixed = TRUE)
+  has_eq <- vapply(parts, length, integer(1)) == 2
+  parts <- parts[has_eq]
+  if (length(parts) == 0) return(stats::setNames(character(0), character(0)))
+  vals <- vapply(parts, function(p) trimws(p[[2]]), character(1))
+  keys <- vapply(parts, function(p) trimws(p[[1]]), character(1))
+  stats::setNames(vals, keys)
+}
+
+# UAN-32 / CAN-17 N speciation per Nichols 2024 (paper text).
+# Fractions of *total N* (not mass fractions of solution).
+.UAN32_FRAC <- c(nh4 = 7.75 / 32, no3 = 7.75 / 32, urea = 16.5 / 32)
+.CAN17_FRAC <- c(nh4 = 5.4 / 17,  no3 = 11.6 / 17, urea = 0)
+
+# Convert reported `level + units` to kg m^-2 for fertilization events.
+# Returns NA if units are unrecognised.
+.to_kg_per_m2 <- function(level, units) {
+  if (is.na(level) || is.na(units)) return(NA_real_)
+  level <- as.numeric(level)
+  if (is.na(level)) return(NA_real_)
+  u <- tolower(trimws(units))
+  if (grepl("mg ha-1", u, fixed = TRUE) || grepl("dry t ha-1", u, fixed = TRUE)) {
+    return(level * 0.1)
+  }
+  if (grepl("kg n ha-1", u, fixed = TRUE) || grepl("kg ha-1", u, fixed = TRUE)) {
+    return(level * 1e-4)
+  }
+  NA_real_
+}
+
+# Map our internal mgmttype to the events_schema_v0.1.1 event_type enum.
+.event_type_for <- function(mgmttype) {
+  switch(
+    as.character(mgmttype),
+    planting             = "planting",
+    harvest              = "harvest",
+    tillage              = "tillage",
+    fertilization        = "fertilization",
+    compost_application  = "fertilization",
+    cover_crop_planting  = "planting",
+    NA_character_
+  )
+}
+
+#' Map one managements row to one events.json event element
+#'
+#' Translates a single row from the cal/val managements tab into a list
+#' that conforms to the per-event object in PEcAn's
+#' `events_schema_v0.1.1.json`. Required fields the data does not yet
+#' provide (`leaf_c_kg_m2`, `frac_above_removed_0to1`,
+#' `tillage_eff_0to1`) are returned as `NA` so the prior-sampling stage
+#' can fill them later.
+#'
+#' @param row Named list / single-row tibble from the managements tab.
+#' @param fertilizer Default fertilizer formulation when N events do not
+#'   specify; one of `"uan_32"` (default) or `"can_17"`. Used by the
+#'   Nichols 2024 fertigation rows where the source data lumps
+#'   "UAN_32 or CAN_17" without per-event resolution.
+#' @return A list representing one event, or `NULL` if the row's
+#'   `mgmttype` cannot be mapped.
+#' @export
+map_event <- function(row, fertilizer = c("uan_32", "can_17")) {
+  fertilizer <- match.arg(fertilizer)
+  row <- as.list(row)
+
+  ev_type <- .event_type_for(row$mgmttype)
+  if (is.na(ev_type)) return(NULL)
+
+  date <- midpoint_date(row$min_date, row$max_date)
+  attrs <- parse_kv(row$attributes_keyvalue)
+
+  out <- list(
+    event_type = ev_type,
+    date       = date
+  )
+  if (!is.null(row$citation) && !is.na(row$citation) && nzchar(row$citation)) {
+    out$source <- row$citation
+  }
+
+  # NOTE on missing required schema fields: the v0.1.1 schema requires
+  # `leaf_c_kg_m2` for planting, `frac_above_removed_0to1` for harvest,
+  # and `tillage_eff_0to1` for tillage. These are NOT in the cal/val
+  # data and will be injected by the prior-sampling stage. We omit them
+  # here entirely (rather than writing `null`) so the JSON is well-formed
+  # and the priors stage can simply add the fields. Validation against
+  # the schema will therefore surface a known set of gaps until the
+  # priors layer is wired in.
+
+  if (ev_type == "planting") {
+    out$crop_code    <- row$crop_name
+    out$crop_display <- row$crop_name
+    if (!is.null(row$cover_crop) && !is.na(row$cover_crop) &&
+        toupper(row$cover_crop) == "TRUE") {
+      out$cover_crop <- TRUE        # extra metadata; schema allows additionalProperties
+    }
+    return(out)
+  }
+
+  if (ev_type == "harvest") {
+    if (!is.null(row$crop_name) && !is.na(row$crop_name)) out$crop_display <- row$crop_name
+    return(out)
+  }
+
+  if (ev_type == "tillage") {
+    if ("implement" %in% names(attrs)) out$intensity_category <- attrs[["implement"]]
+    if (!is.na(suppressWarnings(as.numeric(row$level)))) {
+      val <- as.numeric(row$level)
+      u <- tolower(trimws(if (is.null(row$units)) "" else as.character(row$units)))
+      if (grepl("cm", u, fixed = TRUE))      out$depth_m <- val * 0.01
+      else if (grepl("\\bm\\b", u))           out$depth_m <- val
+    }
+    return(out)
+  }
+
+  if (ev_type == "fertilization") {
+    is_compost <- isTRUE(row$mgmttype == "compost_application") ||
+      grepl("compost|manure|organic", as.character(row$units), ignore.case = TRUE) ||
+      ("material" %in% names(attrs) && grepl("compost|manure",
+                                             attrs[["material"]], ignore.case = TRUE))
+    amount <- .to_kg_per_m2(row$level, row$units)
+
+    if (is_compost && !is.na(amount)) {
+      # Compost: derive C and N from attributes when present, else use
+      # Akash's defaults (carbon_fraction = 0.35, cn_ratio inferred from
+      # PEcAn fertilizer_composition_data later via priors stage).
+      if ("compost_C_pct" %in% names(attrs)) {
+        out$org_c_kg_m2 <- amount * as.numeric(attrs[["compost_C_pct"]]) / 100
+      } else {
+        out$org_c_kg_m2 <- amount * 0.35
+      }
+      if ("compost_N_pct" %in% names(attrs)) {
+        out$org_n_kg_m2 <- amount * as.numeric(attrs[["compost_N_pct"]]) / 100
+      } else if ("CN_ratio" %in% names(attrs) && !is.null(out$org_c_kg_m2)) {
+        out$org_n_kg_m2 <- out$org_c_kg_m2 / as.numeric(attrs[["CN_ratio"]])
+      } else if ("N_content_g_kg" %in% names(attrs)) {
+        # 15 g N / kg compost = 1.5% N
+        out$org_n_kg_m2 <- amount * as.numeric(attrs[["N_content_g_kg"]]) / 1000
+      }
+      return(out)
+    }
+
+    # Mineral N event. Look for explicit fertilizer_type, else fall back
+    # to function default.
+    n_total <- amount   # kg N m-2
+    if (is.na(n_total)) return(out)   # no amount; let priors handle later
+
+    fert_kind <- if ("fertilizer_type" %in% names(attrs)) attrs[["fertilizer_type"]]
+                 else fertilizer
+    fert_kind <- tolower(trimws(fert_kind))
+    if (grepl("uan", fert_kind))      frac <- .UAN32_FRAC
+    else if (grepl("can", fert_kind)) frac <- .CAN17_FRAC
+    else                              frac <- .UAN32_FRAC   # safe default
+
+    out$nh4_n_kg_m2 <- n_total * frac[["nh4"]]
+    out$no3_n_kg_m2 <- n_total * frac[["no3"]]
+    if (frac[["urea"]] > 0) out$org_n_kg_m2 <- n_total * frac[["urea"]]
+    return(out)
+  }
+
+  out
+}
+
+#' Map all rows for a single treatment
+#'
+#' @param mgmt Tibble of managements rows (already filtered to one
+#'   `treatments.name`).
+#' @inheritParams map_event
+#' @return A list of event lists, in the order the rows appeared.
+#' @export
+map_events_for_treatment <- function(mgmt, fertilizer = "uan_32") {
+  events <- vector("list", nrow(mgmt))
+  for (i in seq_len(nrow(mgmt))) {
+    events[[i]] <- map_event(mgmt[i, , drop = FALSE], fertilizer = fertilizer)
+  }
+  Filter(Negate(is.null), events)
+}