Merge pull request #4079 from divine7022/n_ncc_tweaks

Make NCC compost timing crop specific and align both event workflows on the latest version monitoring

Author: mdietze

workflows/fertilization-statewide/01-build-parcel-design.R

@@ -88,13 +88,11 @@ code_lookup <- code_map |>
 
 PEcAn.logger::logger.info(sprintf("Resolved %d CADWR codes via crosswalk", nrow(code_lookup)))
 
-# the event date is anchored to green-up (leafonday) from the gap-filled
-# phenology product, observed where the satellite retrieval succeeded and
-# crop-calendar filled otherwise, so this covers the full ~600k ag universe
-# instead of the ~377k strict-matched subset. crop class per season comes
-# from the CADWR Land Use crops product. the crops product's own emergence
-# date is empty statewide, so the gap-filled green-up is the only populated
-# anchor available.
+# the event date is the anchor itself; this workflow applies no offset. the
+# anchor transition is chosen per PFT so annuals are timed to planting and
+# perennials to leaf-on, matching the split the monitoring event products use.
+# both come from the gap-filled LandIQ to MSLSP match, which keys every
+# transition by (parcel_id, year, season), so cycles join on a real season key.
 
 years <- config[["years"]]
 PEcAn.logger::logger.info("Reading crops and gap-filled phenology for years: ",
@@ -117,50 +115,78 @@ crops <- DBI::dbGetQuery(con, sprintf(
   dplyr::rename(year = "yr") |>
   dplyr::mutate(code = paste0(.data$CLASS, .data$SUBCLASS))
 
-# the phenology product has no season key, but from 2018 on it carries a second
-# green-up for most double-crop parcels, so rank green-ups within a parcel-year
-# and match the nth crop cycle to the nth green-up rather than collapsing to the
-# earliest. phenology_source is carried through for audit.
-phen_anchor_col <- config[["phen_anchor_col"]]
-phen_raw <- DBI::dbGetQuery(con, sprintf(
-  "SELECT * FROM read_parquet('%s') WHERE \"year\" IN (%s)",
+pft_anchor <- unlist(config[["pft_anchor"]])
+anchor_cols <- sort(unique(pft_anchor))
+
+# "**" is the LandIQ sentinel for subclass not specified. it becomes NA on both
+# sides of the join, so those rows land on the class-level fallback
+pft_lookup <- readr::read_csv(config[["pft_lookup_path"]],
+                              show_col_types = FALSE) |>
+  dplyr::filter(!is.na(.data$PFT))
+pft_by_code <- pft_lookup |>
+  dplyr::transmute(CLASS = .data$CLASS,
+                   SUBCLASS = as.integer(dplyr::na_if(as.character(.data$SUBCLASS), "**")),
+                   pft_group = .data$PFT) |>
+  dplyr::distinct()
+pft_by_class <- pft_lookup |>
+  dplyr::count(.data$CLASS, .data$PFT) |>
+  dplyr::slice_max(.data$n, n = 1, by = "CLASS", with_ties = FALSE) |>
+  dplyr::transmute(CLASS = .data$CLASS, pft_group_class = .data$PFT)
+
+phen <- DBI::dbGetQuery(con, sprintf(
+  "SELECT CAST(parcel_id AS INTEGER) AS parcel_id, CAST(\"year\" AS INTEGER) AS year,
+          CAST(season AS INTEGER) AS season, gapfill_date_source, %s
+   FROM read_parquet('%s') WHERE \"year\" IN (%s)",
+  paste(anchor_cols, collapse = ", "),
   file.path(config[["phen_dir"]], config[["phen_glob"]]), yr_list))
-phen_id_col <- if ("parcel_id" %in% names(phen_raw)) "parcel_id" else "site_id"
-phen_source_col <- if ("phenology_source" %in% names(phen_raw)) {
-  "phenology_source"
-} else if ("gapfill_date_source" %in% names(phen_raw)) {
-  "gapfill_date_source"
-} else {
-  NA_character_
-}
-phen <- phen_raw |>
-  dplyr::transmute(
-    parcel_id = as.integer(.data[[phen_id_col]]),
-    year = as.integer(.data$year),
-    date = as.Date(.data[[phen_anchor_col]]),
-    phenology_source = if (is.na(phen_source_col)) {
-      NA_character_
-    } else {
-      as.character(.data[[phen_source_col]])
-    }
-  ) |>
-  dplyr::filter(!is.na(.data$date)) |>
-  dplyr::arrange(.data$parcel_id, .data$year, .data$date) |>
-  dplyr::mutate(phen_rank = dplyr::row_number(), .by = c("parcel_id", "year"))
 
-# where a parcel-year has fewer green-ups than crop cycles, the later cycles
-# reuse the last available one
-phen_max <- phen |>
-  dplyr::summarize(max_rank = max(.data$phen_rank), .by = c("parcel_id", "year"))
+missing_cols <- setdiff(anchor_cols, names(phen))
+if (length(missing_cols) > 0) {
+  PEcAn.logger::logger.severe(
+    "phenology product has no column(s) named in pft_anchor: ",
+    paste(missing_cols, collapse = ", "))
+}
 
+# a cycle with no matched phenology row has no anchor, and is dropped rather
+# than given a substitute date
 plant <- crops |>
-  dplyr::mutate(season_rank = dplyr::dense_rank(.data$season),
-                .by = c("parcel_id", "year")) |>
-  dplyr::inner_join(phen_max, by = c("parcel_id", "year")) |>
-  dplyr::mutate(phen_rank = pmin(.data$season_rank, .data$max_rank)) |>
-  dplyr::inner_join(phen, by = c("parcel_id", "year", "phen_rank"))
-PEcAn.logger::logger.info(sprintf("Loaded %d cycles across %d parcels (phenology anchored)",
-                                  nrow(plant), dplyr::n_distinct(plant$parcel_id)))
+  dplyr::inner_join(phen, by = c("parcel_id", "year", "season")) |>
+  dplyr::left_join(pft_by_code, by = c("CLASS", "SUBCLASS")) |>
+  dplyr::left_join(pft_by_class, by = "CLASS") |>
+  dplyr::mutate(pft_group = dplyr::coalesce(.data$pft_group, .data$pft_group_class))
+PEcAn.logger::logger.info(sprintf(
+  "Anchored %d of %d crop cycles (%.1f%%) across %d parcels",
+  nrow(plant), nrow(crops), 100 * nrow(plant) / nrow(crops),
+  dplyr::n_distinct(plant$parcel_id)))
+
+# non-crop pfts have no anchor rule. report them so a crop type missing a rule
+# is visible rather than silently absent
+dropped <- plant |>
+  dplyr::filter(!.data$pft_group %in% names(pft_anchor)) |>
+  dplyr::count(.data$pft_group, sort = TRUE)
+if (nrow(dropped) > 0) {
+  PEcAn.logger::logger.info(sprintf(
+    "Dropping %d cycles whose pft has no anchor rule:", sum(dropped$n)))
+  for (i in seq_len(nrow(dropped))) {
+    PEcAn.logger::logger.info(sprintf("  %s: %d cycles",
+                                      dropped$pft_group[i], dropped$n[i]))
+  }
+}
+plant <- plant |> dplyr::filter(.data$pft_group %in% names(pft_anchor))
+
+# index a numeric matrix so the anchor stays config driven, not a branch per pft
+anchor_idx <- cbind(seq_len(nrow(plant)),
+                    match(pft_anchor[plant$pft_group], anchor_cols))
+anchor_num <- do.call(cbind, lapply(plant[anchor_cols], as.numeric))
+plant$date <- as.Date(anchor_num[anchor_idx], origin = "1970-01-01")
+
+# a matched row is expected to carry every transition, so a NULL anchor means the
+# product changed rather than a cycle being legitimately undated
+no_anchor <- sum(is.na(plant$date))
+if (no_anchor > 0) {
+  PEcAn.logger::logger.severe(sprintf(
+    "%d cycles have a NULL anchor in the gap-filled product", no_anchor))
+}
 
 ## subsample
 # parcel set is sampled once and applied to all years so the same parcels
@@ -226,16 +252,16 @@ if (nrow(zero_env) > 0) {
 }
 
 kept <- design |> dplyr::filter(.data$rate_source == "crosswalk")
-src <- kept |> dplyr::count(.data$phenology_source, sort = TRUE)
-PEcAn.logger::logger.info("Anchor provenance (phenology_source):")
+src <- kept |> dplyr::count(.data$gapfill_date_source, sort = TRUE)
+PEcAn.logger::logger.info("Anchor provenance (gapfill_date_source):")
 for (i in seq_len(nrow(src))) {
   PEcAn.logger::logger.info(sprintf("  %s: %d cycles (%.1f%%)",
-                                    src$phenology_source[i], src$n[i],
+                                    src$gapfill_date_source[i], src$n[i],
                                     100 * src$n[i] / nrow(kept)))
 }
 
 design <- kept |>
-  dplyr::select("parcel_id", "year", "season", "date", "code",
+  dplyr::select("parcel_id", "year", "season", "date", "code", "pft_group",
                 "min_n_lbs_acre", "max_n_lbs_acre") |>
   # fixed row order so the per row draws in 02 are reproducible under the seed
   dplyr::arrange(.data$parcel_id, .data$year, .data$season)

workflows/fertilization-statewide/README.md

@@ -9,13 +9,14 @@ Source: California Department of Water Resources. (2016-2023). Statewide Crop Ma
 Configuration parameters live in `config.yml`. Most setups only need:
 
 - `crops_path`: the harmonized CADWR Land Use crops parquet. Override with `CCMMF_CROPS_PATH`
-- `phen_dir`: the gap-filled phenology (green-up) directory. Override with `CCMMF_PHEN_DIR`
-- `phen_glob`: file glob under `phen_dir` (default `phenology_statewide_*.parquet`). Override with `CCMMF_PHEN_GLOB`
-- `phen_anchor_col`: phenology date column (default `leafonday`). Override with `CCMMF_PHEN_ANCHOR_COL`
+- `phen_dir`: the gap-filled LandIQ to MSLSP match directory, the same product the ncc workflow reads. Override with `CCMMF_PHEN_DIR`
+- `phen_glob`: file glob under `phen_dir` (default `assigned_year=*_gapfilled.parquet`). Override with `CCMMF_PHEN_GLOB`
+- `pft_lookup_path`: crop code to PFT table, the same one the monitoring products use. Override with `CCMMF_PFT_LOOKUP`
+- `pft_anchor`: phenology transition used as the anchor, per PFT. See Application timing
 - `crosswalk_path`: the CADWR to FREP to UC ANR crop name crosswalk TSV, versioned in this folder
 - `output_dir`: output directory for parquet shards
 - `n_parcels`, `n_ensemble`, `batch_size`, `workers`: settings per profile
-- `nh4_fraction`: share of total synthetic N going to ammonium; the rest goes to nitrate (default 0.5 for a 50/50 split)
+- `nh4_fraction`: share of total synthetic N going to ammonium; the rest goes to nitrate (default 1, all ammonium)
 
 # Run
 
@@ -38,15 +39,32 @@ Events carry no organic C or N: these are synthetic mineral fertilizer applicati
 
 # Known limitations
 
-- The phenology product has no season key, so crop cycles are matched to green-ups by
-  rank: the nth cycle of a parcel-year takes the nth green-up. From 2018 on the product
-  carries a second green-up for most double-crop parcels, so this resolves the majority
-  of them. Where a parcel-year has fewer green-ups than cycles, the later cycles reuse
-  the last available green-up. 2016 is the exception, carrying one green-up per parcel
-  year, so its multi-season cycles all share an anchor.
-- The event date is the green-up date itself; this workflow applies no offset. Dates
-  outside the configured crop years occur because the phenology product itself assigns
-  some green-ups to the previous calendar year.
+- A crop cycle with no matched phenology row gets no anchor and is dropped. Against the
+  gap-filled LandIQ table this is about 99% of cycles.
+- The event date is the anchor itself; this workflow applies no offset. Dates outside
+  the configured crop years occur because the phenology product itself assigns some
+  transitions to an adjacent calendar year.
+
+# Application timing
+
+The anchor transition is chosen per PFT, so annuals are timed to planting and perennials
+to leaf-on. This matches the split the monitoring event products use: they report planting
+for annuals and leaf-on for perennials, not both for both.
+
+| PFT | anchor | transition means |
+|---|---|---|
+| row | `mslsp_OGI` | onset of greenness increase, 15%, used as planting |
+| rice | `mslsp_OGI` | as above |
+| hay | `mslsp_50PCGI` | 50% greenness increase, leaf-on |
+| woody | `mslsp_50PCGI` | as above |
+
+The event date is the anchor itself. Unlike the compost workflow there is no offset window,
+because `ca_n_application_rate` is an annual total per crop rather than an application
+schedule.
+- `ca_n_application_rate` is an annual total per crop, so the whole season's N budget
+  is applied as one event rather than split across pre-plant, side-dress and
+  fertigation. The date is therefore a single anchored placeholder for a schedule, not
+  a measured application date.
 - Only crop codes present in the crosswalk resolve to an N rate envelope. Cycles whose
   code does not resolve are dropped and reported at run time.
 

workflows/fertilization-statewide/config.yml

@@ -8,15 +8,33 @@ default:
   nh4_fraction: 1
   seed: 42
   years: [2016, 2018, 2019, 2020, 2021, 2022, 2023]
+  # crops is the gap-filled LandIQ table (v4.1.2), the same inventory the
+  # phenology match is built against. v4.1 is the ungapfilled table and covers
+  # parcels the match does not.
   # large inputs are not distributed with the repo. the defaults resolve on SCC;
   # set the matching environment variable to run anywhere else. final home for
   # these is pending ccmmf/organization#257
-  crops_path: !expr path.expand(Sys.getenv("CCMMF_CROPS_PATH", "/projectnb/dietzelab/ccmmf/LandIQ-harmonized-v4.1/crops_all_years.parq"))
-  phen_dir: !expr path.expand(Sys.getenv("CCMMF_PHEN_DIR", "/projectnb/dietzelab/ccmmf/usr/akash/phen_v2"))
-  phen_glob: !expr Sys.getenv("CCMMF_PHEN_GLOB", "phenology_statewide_*.parquet")
-  phen_anchor_col: !expr Sys.getenv("CCMMF_PHEN_ANCHOR_COL", "leafonday")
+  crops_path: !expr path.expand(Sys.getenv("CCMMF_CROPS_PATH", "/projectnb/dietzelab/ccmmf/LandIQ-harmonized-v4.1.2/crops_all_years.parq"))
+  # gap-filled LandIQ to MSLSP match, the same product the ncc workflow reads
+  phen_dir: !expr path.expand(Sys.getenv("CCMMF_PHEN_DIR", "/projectnb/dietzelab/ccmmf/management/phenology/matched_landiq_mslsp_v4.1.2/gapfill_dates"))
+  phen_glob: !expr Sys.getenv("CCMMF_PHEN_GLOB", "assigned_year=*_gapfilled.parquet")
+
+  # crop code to PFT. same table the monitoring products use
+  pft_lookup_path: !expr path.expand(Sys.getenv("CCMMF_PFT_LOOKUP", "/projectnb/dietzelab/ccmmf/management/LandIQ_cropCode_lookup_table.csv"))
+
+  # anchor transition per PFT. annuals are timed to planting and perennials to
+  # leaf-on, matching the annual/perennial split the monitoring event products
+  # use: they report planting for annuals and leaf-on for perennials, not both
+  # for both. the event date is the anchor itself, this workflow applies no offset
+  pft_anchor:
+    row: mslsp_OGI
+    rice: mslsp_OGI
+    hay: mslsp_50PCGI
+    woody: mslsp_50PCGI
   crosswalk_path: workflows/fertilization-statewide/crop_name_crosswalk.tsv
-  output_dir: !expr path.expand(Sys.getenv("CCMMF_FERT_OUT", "/projectnb/dietzelab/ccmmf/usr/akash/event_files/fertilization"))
+  # versioned alongside ncc: v2.0 is the pass built on the v4.1.2 monitoring
+  # products. bump rather than writing over
+  output_dir: !expr path.expand(Sys.getenv("CCMMF_FERT_OUT", "/projectnb/dietzelab/ccmmf/usr/akash/event_files/fertilization/v2.0"))
   batch_size: 100
   workers: 1