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

Published

2026-07-15 17:54:42

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

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

2 Setup — assemble the scoring DB

Code 
librarian::shelf(DBI, dplyr, duckdb, exactextractr, 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)
Spherical geometry (s2) switched off
Code 
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 
zone_specs <- list(
  list(tbl = glue("ply_programareas_2026_{ver}"), fld = "programarea_key", gpkg = pa_gpkg),
  list(tbl = "ply_ecoregions_2025",               fld = "ecoregion_key",   gpkg = er_gpkg))

r_cell <- rast(cellid_tif)
zone_rows <- list(); zc_rows <- list(); zseq <- 0L
for (zs in zone_specs) {
  ply <- st_read(zs$gpkg, quiet = TRUE) |> st_make_valid() |>
    st_transform(4326) |> select(all_of(zs$fld))
  for (i in seq_len(nrow(ply))) {
    zseq <- zseq + 1L
    zval <- ply[[zs$fld]][i]
    zone_rows[[length(zone_rows) + 1]] <- tibble(
      zone_seq = zseq, tbl = zs$tbl, fld = zs$fld, val = as.character(zval))
    # coverage fraction of this zone polygon over each cell -> pct_covered (0-100)
    df <- exact_extract(r_cell, ply[i, ], progress = FALSE)[[1]]   # exactextractr's own `value` col
    df <- df[!is.na(df$value) & df$coverage_fraction > 0, , drop = FALSE]
    if (nrow(df) > 0) zc_rows[[length(zc_rows) + 1]] <- tibble(
      zone_seq = zseq, cell_id = as.integer(df$value),
      pct_covered = as.integer(round(df$coverage_fraction * 100)))
  }
  log_info("{zs$fld}: {nrow(ply)} zones")
}

# subregion zones (USA / AK / GA / PA) = unions of their member Program Areas (the v7 grouping,
# `{ver_prev}/subregion_programareas.csv`). No polygon source — build from the PA zones' cells so
# the apps' "Study area" filter (subregion_key) works.
sr_map   <- read_csv(glue("{dir_data}/derived/{ver_prev}/subregion_programareas.csv"), show_col_types = FALSE)
zone_df  <- bind_rows(zone_rows); zc_df <- bind_rows(zc_rows)
pa_zones <- zone_df |> filter(fld == "programarea_key") |> select(pa_seq = zone_seq, programarea_key = val)
for (sr in sort(unique(sr_map$subregion_key))) {
  zseq <- zseq + 1L
  member_seqs <- pa_zones$pa_seq[pa_zones$programarea_key %in% sr_map$programarea_key[sr_map$subregion_key == sr]]
  sr_cells <- zc_df |> filter(zone_seq %in% member_seqs) |>
    group_by(cell_id) |> summarise(pct_covered = max(pct_covered), .groups = "drop")
  zone_rows[[length(zone_rows) + 1]] <- tibble(
    zone_seq = zseq, tbl = glue("ply_subregions_2026_{ver}"), fld = "subregion_key", val = sr)
  if (nrow(sr_cells)) zc_rows[[length(zc_rows) + 1]] <- tibble(
    zone_seq = zseq, cell_id = sr_cells$cell_id, pct_covered = sr_cells$pct_covered)
}
log_info("subregion_key: {n_distinct(sr_map$subregion_key)} zones")

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 4 697280
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"))