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

Published

2026-08-10

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")
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.
redo_mc <- nzchar(Sys.getenv("REDO_SCORE_ZONES"))
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) {
  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)")
  log_info("model_cell (re)built: {n_pq} rows")
} else log_info(
  "model_cell up-to-date ({n_mc} rows == merged parquet) — skipping ~30-min rebuild (REDO_SCORE_ZONES=1 to force)")

3 Taxon prep — sp_cat + 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 is_er_spatial (turtles: SWOT+DPS ER baked into the merged cell val ->
# er_score=100 at scoring so it isn't double-counted). 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

dbExecute(con, "CREATE OR REPLACE TABLE v7_cat AS
  SELECT taxon_authority, CAST(CAST(taxon_id AS BIGINT) AS VARCHAR) AS taxon_id,
         bool_or(is_ok) AS is_ok
  FROM v7.taxon WHERE taxon_id IS NOT NULL GROUP BY 1, 2")

[1] 17561

Code

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

[1] 37051

Code

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	9410
fish	6286
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
ecoregion_key	12	630944
programarea_key	20	350938
subregion_key	5	1259994

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, "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"))