fix(spatialize): stop gridded output carrying two disagreeing territorial keys

NEWS.md

@@ -69,6 +69,29 @@
   `whep_sf_required` if `sf` is not installed instead of returning the broken
   object. No published values change: the argument-less call is untouched, and
   both in-package callers use it.
+* **`build_gridded_landuse()` and `build_gridded_livestock()` take an
+  `area_key`, and say when their output cannot join a national table.** The
+  spatialize chain allocates on the raw reporting codes its `country_areas`
+  and `country_grid` are keyed on, while whep's polity-keyed national tables
+  are aggregated on `polity_area_code`. A reporting code that is not itself a
+  bucket therefore left every output row carrying two territorial keys that
+  disagree — `area_code = 276` beside `polity_area_code = 206` — so whether a
+  consumer joined on one or the other decided whether Sudan existed in its
+  result (#582). Measured against the deployed pins, `country_grid` holds 831
+  such cells under 2 codes (276 Sudan, 277 South Sudan) and `country_areas`
+  0.64% of its harvested area; the other six codes the issue listed are no
+  longer off-bucket, because #628 gave Syria, North Macedonia, Eswatini,
+  Equatorial Guinea, New Caledonia and Palestine their own published codes.
+  The default `area_key = "grid"` is unchanged bit-for-bit and now warns
+  naming the codes that cannot join; `"polity_area"` re-keys the output on the
+  bucket before the polity columns are attached, so the two keys agree in
+  every row. **No published value changes** unless `"polity_area"` is asked
+  for: on a 2020 Sudan/South Sudan run it conserved 21,894,526 ha and 230.7 M
+  head exactly, kept the row count, and moved 13,447 crop rows and 3,671
+  livestock rows from a key no national table carries onto `206`. Under
+  `"polity_area"` the raw code is carried, not replaced, as `grid_area_code`,
+  the shape `build_cell_polity()` adopted in #579. `run_spatialize()` accepts
+  `area_key` in `overrides`.
 
 * **`get_primary_production()`, `get_wide_cbs()` and `get_processing_coefs()`
   take a `years` argument.** A scoped request now builds only that window

R/run_spatialize.R

@@ -37,6 +37,9 @@
 #'     `cft_mapping.csv` drives CFT aggregation: `cft_name`
 #'     (granular 33-class WHEP taxonomy) or `cft_lpjml`
 #'     (12 LPJmL crop CFTs + single `others` bucket).
+#'   - `area_key`: one of `"grid"` (default) or `"polity_area"`,
+#'     forwarded to both engines. See [build_gridded_landuse()]'s
+#'     *Which area code the output is keyed on*.
 #' @param paths Named list of filesystem paths. Recognised entries:
 #'   - `l_files_dir`: path to the `L_files` root, for local prepared inputs.
 #'   - `input_dir`: directory holding the prepared input parquets. If `NULL`
@@ -173,6 +176,7 @@ run_spatialize <- function(
       preset,
       years,
       resolved_years,
+      config,
       input_dir,
       out_dir
     )
@@ -237,7 +241,8 @@ run_spatialize <- function(
       multicropping = lu_inputs$multicropping,
       years = resolved_years,
       max_iterations = config$max_iterations,
-      expansion_threshold = config$expansion_threshold
+      expansion_threshold = config$expansion_threshold,
+      area_key = config$area_key
     )
   )
   list(
@@ -256,6 +261,7 @@ run_spatialize <- function(
   preset,
   years,
   resolved_years,
+  config,
   input_dir,
   out_dir
 ) {
@@ -279,7 +285,8 @@ run_spatialize <- function(
     country_grid = ls_inputs$country_grid,
     species_proxy = ls_inputs$species_proxy,
     manure_pattern = ls_inputs$manure_pattern,
-    years = resolved_years
+    years = resolved_years,
+    area_key = config$area_key
   )
   list(
     years = resolved_years,
@@ -293,13 +300,15 @@ run_spatialize <- function(
       use_type_constraint = FALSE,
       aggregate_to_cft = TRUE,
       max_iterations = 1000L,
-      expansion_threshold = 100L
+      expansion_threshold = 100L,
+      area_key = "grid"
     ),
     whep = list(
       use_type_constraint = TRUE,
       aggregate_to_cft = TRUE,
       max_iterations = 1000L,
-      expansion_threshold = 100L
+      expansion_threshold = 100L,
+      area_key = "grid"
     )
   )
 }
@@ -310,7 +319,8 @@ run_spatialize <- function(
     "aggregate_to_cft",
     "max_iterations",
     "expansion_threshold",
-    "cft_target"
+    "cft_target",
+    "area_key"
   )
 }
 

R/spatialize.R

@@ -83,6 +83,9 @@
 #'   - `expansion_threshold`: Iteration number after which crops are
 #'     allowed to expand into cells without an existing pattern.
 #'     Default: `100L`.
+#'   - `area_key`: Which area code the output is keyed on, `"grid"`
+#'     (default) or `"polity_area"`. See *Which area code the output is
+#'     keyed on*.
 #'
 #' @return A tibble with gridded crop (or CFT) harvested areas.
 #'   Columns:
@@ -92,12 +95,41 @@
 #'   - `polity_area_code`, `reporting_polity_code`,
 #'     `reporting_polity_name`, `reporting_polity_has_geometry`: Polity
 #'     metadata for `area_code`.
+#'   - `grid_area_code`: Only under `area_key = "polity_area"`; the
+#'     reporting code the engine allocated on.
 #'   - `polycell_id`, `cell_id`: Preserved when supplied in
 #'     `country_grid`.
 #'   - `crop_name` or `cft_name`: Crop or CFT identifier.
 #'   - `rainfed_ha`: Rainfed harvested area in the cell.
 #'   - `irrigated_ha`: Irrigated harvested area in the cell.
 #'
+#' @section Which area code the output is keyed on:
+#' The chain allocates *from* a national table keyed on `area_code` and
+#' *into* a `country_grid` keyed the same way, so both sides speak the raw
+#' reporting vocabulary the grid was rasterized in. WHEP's polity-keyed
+#' national tables are aggregated on `polity_area_code` instead, a bucket
+#' that a reporting code need not equal: `276` Sudan and `277` South Sudan
+#' both fall in bucket `206`. Every such output row therefore carries two
+#' territorial keys that disagree, and whether a consumer joins on
+#' `area_code` or on `polity_area_code` decides whether Sudan exists in its
+#' result (whep#582).
+#'
+#' `area_key` selects which of the two the output carries. It is not a
+#' fallback: `"grid"` is the default, reproduces today's codes
+#' bit-for-bit, and warns naming the codes that cannot join;
+#' `"polity_area"` resolves each code to its bucket through
+#' [polity_area_crosswalk] before the polity columns are attached, so
+#' `area_code` and `polity_area_code` agree in every row. It respects
+#' `options(whep.unfold_rest_of_world)` (see [folded_reporting_areas()]),
+#' so the output and the national tables agree about where a Rest-of-World
+#' member's rows belong.
+#'
+#' Under `"polity_area"` the raw reporting code is **carried, not
+#' replaced**: the output gains `grid_area_code`, joined with `+` where two
+#' reporting areas of one bucket meet in a cell and their rows collapse. So
+#' the fold stays recoverable at the join rather than baked into the output,
+#' the shape [build_cell_polity()] adopted for the same reason (whep#579).
+#'
 #' @section Methodology:
 #' This function reimplements the spatial crop allocation from the LandInG
 #' toolbox (Ostberg et al. 2023, doi:10.5194/gmd-16-3375-2023) with the
@@ -255,6 +287,10 @@ build_gridded_landuse <- function(
   }
 
   tibble::as_tibble(result) |>
+    .spatialize_apply_area_key(
+      config$area_key,
+      c("rainfed_ha", "irrigated_ha")
+    ) |>
     .add_reporting_polity_columns()
 }
 
@@ -462,7 +498,8 @@ build_gridded_landuse <- function(
     multicropping = NULL,
     max_iterations = 1000L,
     expansion_threshold = 100L,
-    n_workers = 1L
+    n_workers = 1L,
+    area_key = "grid"
   )
 }
 
@@ -483,7 +520,9 @@ build_gridded_landuse <- function(
       "i" = "Known: {.val {names(defaults)}}."
     ))
   }
-  utils::modifyList(defaults, config)
+  config <- utils::modifyList(defaults, config)
+  config$area_key <- .resolve_spatialize_area_key(config$area_key)
+  config
 }
 
 #' Validate that required columns exist.

R/spatialize_compartments.R

@@ -45,6 +45,78 @@
   unique(c(.compartment_id_cols(data), "lon", "lat"))
 }
 
+#' Validate the area key a spatialize engine was asked for.
+#' @noRd
+.resolve_spatialize_area_key <- function(area_key) {
+  if (is.null(area_key)) {
+    area_key <- "grid"
+  }
+  rlang::arg_match0(area_key, c("grid", "polity_area"), arg_nm = "area_key")
+}
+
+#' Apply the requested area key to a spatialized output.
+#'
+#' The spatialize chain allocates *from* a national table keyed on
+#' `area_code` and *into* a `country_grid` keyed the same way, so both sides
+#' speak the raw reporting vocabulary the grid was rasterized in. WHEP's
+#' polity-keyed national tables are aggregated on `polity_area_code` instead,
+#' so a reporting code that is not itself a bucket leaves the output carrying
+#' two territorial keys that disagree within one row -- `area_code = 276`
+#' beside `polity_area_code = 206` -- and a consumer's choice of join column
+#' silently decides whether Sudan exists in its result (whep#582).
+#'
+#' `"grid"` is today's behaviour, kept as the default because switching moves
+#' published values for every consumer joining gridded output on `area_code`;
+#' it only gains the diagnostic that today's silence hides. The alternative is
+#' selected, never a fallback.
+#' @noRd
+.spatialize_apply_area_key <- function(result, area_key, value_cols) {
+  if (area_key == "grid") {
+    .warn_cell_polity_off_bucket(result)
+    return(result)
+  }
+  .spatialize_to_bucket(result, value_cols)
+}
+
+#' Re-key a spatialized output on `polity_area_code`.
+#'
+#' THE RAW CODE IS CARRIED, NOT REPLACED, exactly as `build_cell_polity()`
+#' does under `area_key = "polity_area"` (whep#579): the reporting code the
+#' engine allocated on arrives as an added `grid_area_code`, so the fold this
+#' performs stays recoverable at the join instead of becoming irrecoverable in
+#' the output. Where two reporting areas of one bucket meet in a cell their
+#' rows collapse, their values are summed, and the raw codes are joined with a
+#' separator rather than one of them being picked -- picking would be the
+#' silent half of the same problem. Codes absent from the crosswalk keep their
+#' own code, so a gap stays visible rather than turning into an `NA` key.
+#' @noRd
+.spatialize_to_bucket <- function(result, value_cols) {
+  value_cols <- intersect(value_cols, names(result))
+  dt <- data.table::as.data.table(result)
+  dt[, area_code := as.integer(area_code)]
+  lookup <- data.table::as.data.table(.cell_polity_bucket_lookup())
+  dt[lookup, polity_bucket := i.polity_area_code, on = "area_code"]
+  dt[, grid_area_code := area_code]
+  dt[!is.na(polity_bucket), area_code := polity_bucket]
+  dt[, polity_bucket := NULL]
+  group_cols <- setdiff(names(dt), c(value_cols, "grid_area_code"))
+  out <- dt[,
+    c(
+      lapply(.SD, sum, na.rm = TRUE),
+      list(
+        grid_area_code = paste(sort(unique(grid_area_code)), collapse = "+")
+      )
+    ),
+    by = group_cols,
+    .SDcols = value_cols
+  ]
+  data.table::setcolorder(
+    out,
+    c(intersect(names(result), names(out)), "grid_area_code")
+  )
+  tibble::as_tibble(out)
+}
+
 #' Detect whether a country grid changes through time.
 #' @noRd
 .country_grid_is_dynamic <- function(country_grid) {

R/spatialize_livestock.R

@@ -118,13 +118,20 @@
 #'   (default), all years present in `livestock_data` are processed.
 #'   When supplied, `livestock_data`, `gridded_pasture`, and
 #'   `gridded_cropland` are filtered to this set before processing.
+#' @param area_key Which area code the output is keyed on: `"grid"`
+#'   (default, the reporting codes `livestock_data` and `country_grid` are
+#'   keyed on) or `"polity_area"` (the [polity_area_crosswalk] bucket
+#'   national tables are aggregated on). See
+#'   [build_gridded_landuse()]'s *Which area code the output is keyed on*.
 #'
 #' @return A tibble with gridded livestock data. Columns:
 #'   - `lon`, `lat`: Cell centre coordinates.
 #'   - `area_code`: WHEP polity code for this cell compartment.
 #'   - `polity_area_code`, `reporting_polity_code`,
 #'     `reporting_polity_name`, `reporting_polity_has_geometry`: Polity
 #'     metadata for `area_code`.
+#'   - `grid_area_code`: Only under `area_key = "polity_area"`; the
+#'     reporting code the engine allocated on.
 #'   - `polycell_id`, `cell_id`: Preserved when supplied in
 #'     `country_grid`.
 #'   - `year`: Integer year.
@@ -133,6 +140,8 @@
 #'   - Any additional numeric columns from `livestock_data`
 #'     (e.g. `enteric_ch4_kt`, `manure_ch4_kt`).
 #'
+#' @inheritSection build_gridded_landuse Which area code the output is keyed on
+#'
 #' @export
 #'
 #' @examples
@@ -168,8 +177,10 @@ build_gridded_livestock <- function(
   manure_pattern = NULL,
   glw_density = NULL,
   grass_productivity = NULL,
-  years = NULL
+  years = NULL,
+  area_key = c("grid", "polity_area")
 ) {
+  area_key <- rlang::arg_match(area_key)
   .validate_livestock_inputs(
     livestock_data,
     gridded_pasture,
@@ -242,6 +253,7 @@ build_gridded_livestock <- function(
     dplyr::bind_rows()
 
   result |>
+    .spatialize_apply_area_key(area_key, numeric_cols) |>
     .add_reporting_polity_columns()
 }
 

R/utils.R

@@ -1807,6 +1807,11 @@ utils::globalVariables(
     "lag_weight",
     # build_cbs.R (#187) — .canonicalise_gdp_pop_area() data.table NSE columns
     "canonical_area",
-    "key_row"
+    "key_row",
+    # spatialize_compartments.R (#582) - .spatialize_to_bucket() data.table NSE
+    # columns: the carried raw reporting code and the bucket joined onto it
+    "grid_area_code",
+    "polity_bucket",
+    "i.polity_area_code"
   )
 )