Score setup — zones (program areas + ecoregions) + taxon prep

Published

2026-08-28

Foundation for calc_scores (v8): assemble the scoring DuckDB from the merged model_cell + taxon + cell, prep the taxon for scoring (sp_cat, is_er_spatial), and build zone + zone_cell — cells ↔︎ program areas / ecoregions with pct_covered (0–100) via exactextractr coverage of each zone polygon over the cell-id grid. Scoring formulas need zone_cell.pct_covered at every step.

1 Design

Code

flowchart LR
  mrg["dist_merged model_cell"] --> mc["model_cell (sdm.duckdb)"]
  mtx["merge.duckdb taxon<br/>+ v7 in_v7/is_er_spatial"] --> tx["taxon (scoring)"]
  ply["zones/{zone_set_key}/{grid_id}<br/>zone_cell.parquet (shared)"] --> zc["zone + zone_cell<br/>(+ zone_set_key)"]
  mc --> mf["hash_query(model_cell,taxon,zone,zone_cell)<br/>→ content-addressed manifest"]
  tx --> mf
  zc --> mf

flowchart LR
  mrg["dist_merged model_cell"] --> mc["model_cell (sdm.duckdb)"]
  mtx["merge.duckdb taxon<br/>+ v7 in_v7/is_er_spatial"] --> tx["taxon (scoring)"]
  ply["zones/{zone_set_key}/{grid_id}<br/>zone_cell.parquet (shared)"] --> zc["zone + zone_cell<br/>(+ zone_set_key)"]
  mc --> mf["hash_query(model_cell,taxon,zone,zone_cell)<br/>→ content-addressed manifest"]
  tx --> mf
  zc --> mf

Figure 1: Assemble the scoring DB: merged model_cell + scoring taxon + zone/zone_cell coverage

2 Setup — assemble the scoring DB

Code

librarian::shelf(arrow, DBI, dplyr, duckdb, fs, glue, here, jsonlite, logger, readr,
                 sf, terra, tibble, MarineSensitivity/msens, quiet = T)
source(here("libs/paths.R"))
source(here("libs/vars.R"))
sf_use_s2(FALSE)

dir_atlas <- glue("{dir_big_v}/marine-atlas")
merged    <- glue("{dir_atlas}/dist_merged/dataset=ms_merge")
merged_er <- glue("{dir_atlas}/dist_merged_er/dataset=ms_merge")   # per-cell ER beside it (dps taxa; may be absent)
merge_db  <- glue("{dir_atlas}/merge.duckdb")
pa_gpkg   <- glue("{dir_data}/derived/{ver_prev}/ply_programareas_2026_{ver_prev}.gpkg")
er_gpkg   <- glue("{dir_data}/derived/{ver_prev}/ply_ecoregions_2025.gpkg")
manifest  <- here("data/manifests/score_zones.json")
dir_create(path_dir(manifest))
stopifnot(all(file_exists(c(sdm_db, merge_db, cellid_tif, pa_gpkg, er_gpkg))),
          dir_exists(merged))

con <- dbConnect(duckdb(sdm_db))                    # sdm_db already has `cell` (build_cell_grid)
dbExecute(con, glue("ATTACH '{merge_db}' AS m (READ_ONLY)"))

[1] 0

Code

# model_cell = merged US model (mdl_key, cell_id, val); materialize for scoring joins.
# rebuilding this ~580M-row table + its index is the slow part of score_zones (~30 min), so skip it
# when already materialized AND row-count-identical to the merged parquet — i.e. an idempotent re-run
# to refresh downstream taxon flags only, with the merged surface unchanged. A changed merge (row
# count differs) forces the rebuild; REDO_SCORE_ZONES=1 forces it unconditionally.
# ...but a row count cannot see a merge that changed VALUES on the same footprint (v9.1: the DPS
# species kept their range cells and only their ER-weighted values moved), so the resume key is
# the merge's CONTENT hash (data/manifests/merge_models.json), remembered in sdm.duckdb.
redo_mc  <- nzchar(Sys.getenv("REDO_SCORE_ZONES"))
# resume key = the val surface's content hash + the ER side table's (er_hash; "none" before v9 / when absent):
# the dps taxa's ER can change with their val surface untouched, and scoring reads both
mrg_hash <- tryCatch({
  mf <- jsonlite::fromJSON(here("data/manifests/merge_models.json"))
  paste0(mf$content_hash, "+", if (is.null(mf$er_hash)) "none" else mf$er_hash)
}, error = \(e) NA_character_)
dbExecute(con, "CREATE TABLE IF NOT EXISTS pipeline_state (target VARCHAR PRIMARY KEY, content_hash VARCHAR)")

[1] 0

Code

have_hash <- dbGetQuery(con, "SELECT content_hash FROM pipeline_state WHERE target = 'model_cell'")$content_hash
n_mc <- if ("model_cell" %in% dbListTables(con))
  dbGetQuery(con, "SELECT count(*) n FROM model_cell")$n else 0
n_pq <- dbGetQuery(con, glue("SELECT count(*) n FROM read_parquet('{merged}/**/*.parquet')"))$n
if (redo_mc || n_mc != n_pq || is.na(mrg_hash) || !length(have_hash) || !identical(have_hash, mrg_hash)) {
  dbExecute(con, glue("CREATE OR REPLACE TABLE model_cell AS
    SELECT mdl_key, cell_id, val FROM read_parquet('{merged}/**/*.parquet')"))
  dbExecute(con, "CREATE INDEX IF NOT EXISTS mc_cell ON model_cell(cell_id)")
  # model_cell_er = per-cell extinction risk for the taxa whose ER is spatial but whose merged val is
  # the distribution (dps_nmfs, msens::dps_sql); empty for releases without it. score_cell_metrics
  # multiplies it in where taxon.er_mode = 'cell'.
  have_er <- dir_exists(merged_er) && length(dir_ls(merged_er, glob = "*.parquet", recurse = TRUE)) > 0
  if (have_er) {
    dbExecute(con, glue("CREATE OR REPLACE TABLE model_cell_er AS
      SELECT mdl_key, cell_id, er FROM read_parquet('{merged_er}/**/*.parquet')"))
  } else dbExecute(con, "CREATE OR REPLACE TABLE model_cell_er (mdl_key VARCHAR, cell_id INTEGER, er DOUBLE)")
  n_er <- dbGetQuery(con, "SELECT count(*) n FROM model_cell_er")$n
  dbExecute(con, glue("INSERT OR REPLACE INTO pipeline_state VALUES ('model_cell', '{mrg_hash}')"))
  log_info("model_cell (re)built: {n_pq} rows + model_cell_er {n_er} rows (merge {mrg_hash})")
} else log_info(
  "model_cell up-to-date ({n_mc} rows, merge {mrg_hash} unchanged) — skipping ~30-min rebuild (REDO_SCORE_ZONES=1 to force)")

3 Taxon prep — sp_cat + er_mode (is_er_spatial)

Code

# sp_cat is assigned by TAXONOMY in merge_taxon (WoRMS class/phylum + BirdLife->bird) and
# rides along on m.taxon; here we only add in_v7 (v7's scored set, for the same-species
# compare) and er_mode -- HOW the extinction-risk weight reaches score_cell_metrics:
#   'taxon'         the governing er_score (scalar) x merged val -- every ordinary species
#   'premultiplied' sea turtles: ER x suitability is already the merged val (msens::turtle_sql),
#                   so scoring passes 100 through and does not double-count
#   'cell'          NMFS DPS species (v9, msens::dps_sql): merged val = suitability, per-cell ER
#                   from model_cell_er x val
# is_er_spatial (= er_mode <> 'taxon') stays for the docs/app. v7 taxon_id is DOUBLE, v8 VARCHAR
# -> cast via BIGINT (R as.character gives scientific notation for large ids, won't match).
dbExecute(con, glue("ATTACH '{dir_big}/{ver_prev}/sdm.duckdb' AS v7 (READ_ONLY)"))

[1] 0

Code

# ver_prev's "scored set" column: v1-v7 `is_ok` (bakes in the marine + category cull);
# v8+ `is_valid_usa AND is_marine` (the definition changed at v8). Resolved by introspecting
# the column set, never by version string -- the column is still called `in_v7` downstream
# (SCORE_V7COMMON, schema.qmd) and means "in ver_prev's scored set".
prev_cols <- dbGetQuery(con, "SELECT column_name FROM duckdb_columns()
  WHERE database_name = 'v7' AND table_name = 'taxon'")$column_name
prev_ok <- if ("is_ok" %in% prev_cols) "is_ok" else
  if ("is_valid_usa" %in% prev_cols) "(coalesce(is_valid_usa, FALSE) AND coalesce(is_marine, TRUE))" else
  stop("ver_prev taxon has neither is_ok nor is_valid_usa")
dbExecute(con, glue("CREATE OR REPLACE TABLE v7_cat AS
  SELECT taxon_authority, CAST(TRY_CAST(taxon_id AS BIGINT) AS VARCHAR) AS taxon_id,
         bool_or({prev_ok}) AS is_ok
  FROM v7.taxon WHERE TRY_CAST(taxon_id AS BIGINT) IS NOT NULL GROUP BY 1, 2"))

[1] 37051

Code

dbExecute(con, "CREATE OR REPLACE TABLE taxon AS
  SELECT *, (er_mode <> 'taxon') AS is_er_spatial FROM (
    SELECT t.*,
      coalesce(c.is_ok, FALSE)                                 AS in_v7,   -- in v7's scored set (same-species compare)
      CASE
        WHEN t.ms_merge_key IN (SELECT DISTINCT ms_merge_key FROM m.taxon_model WHERE mdl_key LIKE 'dps_nmfs|%')            THEN 'cell'
        WHEN t.ms_merge_key IN (SELECT DISTINCT ms_merge_key FROM m.taxon_model WHERE mdl_key LIKE 'rng_turtle_swot_dps|%') THEN 'premultiplied'
        ELSE 'taxon' END                                       AS er_mode
    FROM m.taxon t
    LEFT JOIN v7_cat c ON t.taxon_authority = c.taxon_authority AND t.taxon_id = c.taxon_id)")

[1] 37067

Code

# every 'cell' taxon must have its ER beside every scored cell, or scoring would silently weight it NULL
n_gap <- dbGetQuery(con, "SELECT count(*) n FROM model_cell mc JOIN taxon t ON mc.mdl_key = t.ms_merge_key
  LEFT JOIN model_cell_er e ON e.mdl_key = mc.mdl_key AND e.cell_id = mc.cell_id
  WHERE t.er_mode = 'cell' AND e.er IS NULL")$n
stopifnot("er_mode = 'cell' taxa have scored cells without a per-cell ER (model_cell_er)" = n_gap == 0)
log_info("er_mode: {paste(capture.output(print(dbGetQuery(con, 'SELECT er_mode, count(*) n FROM taxon GROUP BY 1 ORDER BY 1'), row.names = FALSE)), collapse = ' | ')}")
dbExecute(con, "DETACH v7")

[1] 0

Code

knitr::kable(dbGetQuery(con, "SELECT sp_cat, count(*) n FROM taxon WHERE is_valid_usa GROUP BY 1 ORDER BY 2 DESC"))

sp_cat	n
invertebrate	9424
fish	6290
bird	880
coral	783
primary_producer	319
mammal	75
turtle	6
reptile	3
amphibian	1

4 Zones + zone_cell (program areas + ecoregions)

Code

# Zones now come from the SHARED zone-set store, not a per-version exactextractr
# pass. zone_cell is a function of (geometry x grid) alone, so re-extracting it
# per release was pure repetition: measured, ONE program-area geometry covers v2
# through v8 and ONE ecoregion geometry covers v1 through v8. build_zone_cells
# computes each pair once, and its gate asserts the result reproduces what this
# chunk used to build, byte for byte.
zone_sets <- readr::read_csv(here("data/zone_sets.csv"), show_col_types = FALSE)
msens::validate_zone_sets(zone_sets)
grid_id <- msens::grid_for_ver(ver)

# every zone set extracted on THIS version's grid — so a release is scored over
# each registered spatial unit, not only its contemporaneous one
zsets <- zone_sets |>
  mutate(pq = glue("{dir_derived}/zones/{zone_set_key}/{grid_id}/zone_cell.parquet")) |>
  filter(file.exists(path.expand(pq)))
stopifnot("no zone-set extracts for this grid — run build_zone_cells" = nrow(zsets) > 0)

# back-compat `tbl`/`fld`: the apps query zone by these, so they keep the names
# they have always had. `zone_set_key` is additive, and is what makes a spatial
# unit comparable across releases.
tbl_of <- c(programarea = glue("ply_programareas_2026_{ver}"),
            ecoregion   = "ply_ecoregions_2025",
            subregion   = glue("ply_subregions_2026_{ver}"),
            planarea    = "ply_planareas_2025")

zone_rows <- list(); zc_rows <- list(); zseq <- 0L
for (i in seq_len(nrow(zsets))) {
  z  <- zsets[i, ]
  zc <- arrow::read_parquet(path.expand(z$pq))
  for (k in sort(unique(zc$zone_key))) {
    zseq <- zseq + 1L
    zone_rows[[length(zone_rows) + 1]] <- tibble(
      zone_seq = zseq, zone_set_key = z$zone_set_key,
      tbl = as.character(tbl_of[[z$zone_type]]), fld = paste0(z$zone_type, "_key"),
      val = as.character(k))
    d <- zc[zc$zone_key == k, ]
    zc_rows[[length(zc_rows) + 1]] <- tibble(
      zone_seq = zseq, cell_id = d$cell_id, pct_covered = d$pct_covered)
  }
  log_info("{z$zone_set_key}: {n_distinct(zc$zone_key)} zones, {nrow(zc)} cells")
}

# `USA` stays a synthesized pseudo-zone (the union of the canonical subregions),
# because the apps use subregion_key = 'USA' for "whole study area". It is the
# absence of a filter rather than a polygon, so it is built from the cells —
# extracting a dissolved USA polygon would double-count nothing, but taking the
# max over member cells is what the previous implementation did and keeps the
# value identical.
zone_df <- bind_rows(zone_rows); zc_df <- bind_rows(zc_rows)
sr_seqs <- zone_df$zone_seq[zone_df$fld == "subregion_key"]
if (length(sr_seqs)) {
  zseq <- zseq + 1L
  usa <- zc_df |> filter(zone_seq %in% sr_seqs) |>
    group_by(cell_id) |> summarise(pct_covered = max(pct_covered), .groups = "drop")
  zone_rows[[length(zone_rows) + 1]] <- tibble(
    zone_seq = zseq, zone_set_key = zsets$zone_set_key[zsets$zone_type == "subregion"][1],
    tbl = as.character(tbl_of[["subregion"]]), fld = "subregion_key", val = "USA")
  zc_rows[[length(zc_rows) + 1]] <- tibble(
    zone_seq = zseq, cell_id = usa$cell_id, pct_covered = usa$pct_covered)
  log_info("subregion USA: {nrow(usa)} cells (union of the canonical subregions)")
}

dbWriteTable(con, "zone",      bind_rows(zone_rows), overwrite = TRUE)
dbWriteTable(con, "zone_cell", bind_rows(zc_rows) |> filter(pct_covered > 0), overwrite = TRUE)

smry <- dbGetQuery(con, "SELECT z.fld, count(DISTINCT z.zone_seq) zones, count(*) zone_cells
  FROM zone z JOIN zone_cell zc USING(zone_seq) GROUP BY 1")
msens::report_table(smry, caption = "zones + zone_cells by field")

zones + zone_cells by field
fld	zones	zone_cells
programarea_key	20	350938
subregion_key	5	1259994
ecoregion_key	12	630944

Code

# content fingerprint of everything this notebook produces for downstream scoring:
# model_cell + scoring taxon + zone + zone_cell (fold; before disconnect)
zones_hash   <- paste0(msens::hash_query(con, "model_cell"), msens::hash_query(con, "model_cell_er"), msens::hash_query(con, "taxon"),
                       msens::hash_query(con, "zone"),       msens::hash_query(con, "zone_cell"))
n_zones      <- dbGetQuery(con, "SELECT count(*) n FROM zone")$n
n_zone_cells <- dbGetQuery(con, "SELECT count(*) n FROM zone_cell")$n
dbDisconnect(con, shutdown = TRUE)

5 Manifest

Code

# content-addressed: fingerprint of model_cell + taxon + zone + zone_cell (no wall-clock)
msens::write_manifest(
  manifest, target = "score_zones", content_hash = zones_hash,
  stats = list(ver = ver, n_zones = n_zones, n_zone_cells = n_zone_cells),
  force = msens::force_target("score_zones"))