Score zone metrics — program-area composite + v7 equivalence gate

Published

2026-08-10

Aggregate the composite components to zones and compute the program-area score (v7’s distinct path, not a roll-up of the cell composite):

component → zone = SUM(cm.val × pct_covered) / SUM(pct_covered)
PA pct_area down-weight: val × (data_coverage / total_coverage) (backs up _prepctareaweighting)
PA score = SUM(val)/COUNT(val) (plain mean) of the pct_area-weighted components

Then the pra_score_delta gate: compare the v8 PA composite to v7’s per program area.

1 Design

Code

flowchart LR
  cm["cell_metric components"] --> za["component → PA<br/>Σ(val·cov)/Σcov"]
  za --> pw["pct_area down-weight"]
  pw --> sc["PA composite score"]
  sc --> zm[("zone_metric (PA rows)<br/>sdm.duckdb")]
  zm --> gate["pra_score_delta<br/>v8 vs v7 gate"]
  zm --> mf["hash_query(zone_metric)<br/>→ content-addressed manifest"]

flowchart LR
  cm["cell_metric components"] --> za["component → PA<br/>Σ(val·cov)/Σcov"]
  za --> pw["pct_area down-weight"]
  pw --> sc["PA composite score"]
  sc --> zm[("zone_metric (PA rows)<br/>sdm.duckdb")]
  zm --> gate["pra_score_delta<br/>v8 vs v7 gate"]
  zm --> mf["hash_query(zone_metric)<br/>→ content-addressed manifest"]

Figure 1: Aggregate components to program areas (coverage-weighted, pct_area down-weight) → PA composite + v7 gate

2 Setup

Code

librarian::shelf(DBI, dplyr, duckdb, fs, glue, here, jsonlite, logger, sf, quiet = T)
source(here("libs/paths.R")); source(here("libs/vars.R"))
manifest <- here("data/manifests/score_zone_metrics.json"); dir_create(path_dir(manifest))
pa_gpkg  <- glue("{dir_data}/derived/{ver_prev}/ply_programareas_2026_{ver_prev}.gpkg")
stopifnot(file_exists(sdm_db), file_exists(pa_gpkg))

con <- dbConnect(duckdb(sdm_db))
mseq_of <- function(k) dbGetQuery(con, glue("SELECT metric_seq FROM metric WHERE metric_key='{k}'"))$metric_seq[1]
mseq_new <- function(k, desc="") { ex<-dbGetQuery(con,glue("SELECT metric_seq FROM metric WHERE metric_key='{k}'"))
  s <- if(nrow(ex)) ex$metric_seq[1] else { s0<-dbGetQuery(con,"SELECT coalesce(max(metric_seq),0)+1 s FROM metric")$s
    dbExecute(con,glue("INSERT INTO metric VALUES ({s0},'{k}','{desc}')")); s0 }
  dbExecute(con, glue("DELETE FROM zone_metric WHERE metric_seq={s} AND zone_seq IN (SELECT zone_seq FROM zone WHERE fld IN ({flds_sql}))")); s }

sp_cats   <- c("bird","coral","fish","invertebrate","mammal","primary_producer","turtle")
comp_keys <- c(glue("extrisk_{sp_cats}_ecoregion_rescaled"), "primprod_ecoregion_rescaled")
comp_seqs <- sapply(comp_keys, mseq_of)

# Aggregate over EVERY zone set in this release, not only its Program Areas.
# zone now carries zone_set_key, so a release is scored over each registered
# spatial unit and one unit becomes comparable across releases. Program-Area
# results are unaffected: the aggregation is per zone_seq and independent
# across zones, so widening the scope only ADDS rows.
zone_flds <- dbGetQuery(con, "SELECT DISTINCT fld FROM zone ORDER BY 1")$fld
flds_sql  <- paste(sprintf("'%s'", zone_flds), collapse = ", ")
log_info("scoring zone fields: {paste(zone_flds, collapse=', ')}")

3 Component → program-area (coverage-weighted mean) + pct_area down-weight

Code

for (k in comp_keys) {
  bs <- mseq_of(k)
  bk <- mseq_new(glue("{k}_prepctareaweighting"), glue("pre-pct_area backup of {k}"))
  # 2a coverage-weighted mean over PA cells, into zone_metric (this metric_seq, PA zones)
  dbExecute(con, glue("DELETE FROM zone_metric WHERE metric_seq={bs} AND zone_seq IN (SELECT zone_seq FROM zone WHERE fld IN ({flds_sql}))"))
  dbExecute(con, glue("INSERT INTO zone_metric (zone_seq, metric_seq, val)
    SELECT zc.zone_seq, {bs}, SUM(cm.val*zc.pct_covered)/SUM(zc.pct_covered)
    FROM zone_cell zc JOIN zone z USING(zone_seq) JOIN cell_metric cm ON zc.cell_id=cm.cell_id
    WHERE z.fld IN ({flds_sql}) AND cm.metric_seq={bs} AND cm.val IS NOT NULL
    GROUP BY zc.zone_seq"))
  # backup pre-weight val
  dbExecute(con, glue("INSERT INTO zone_metric (zone_seq, metric_seq, val)
    SELECT zone_seq, {bk}, val FROM zone_metric WHERE metric_seq={bs}
      AND zone_seq IN (SELECT zone_seq FROM zone WHERE fld IN ({flds_sql}))"))
  # 2b pct_area = (coverage of cells with data for this metric) / (coverage of all cells in zone); val *= pct_area
  dbExecute(con, glue("UPDATE zone_metric zm SET val = val * (
      SELECT SUM(zc.pct_covered) FILTER (WHERE cm.val IS NOT NULL) * 1.0 / SUM(zc.pct_covered)
      FROM zone_cell zc LEFT JOIN cell_metric cm ON zc.cell_id=cm.cell_id AND cm.metric_seq={bs}
      WHERE zc.zone_seq = zm.zone_seq)
    WHERE zm.metric_seq={bs} AND zm.zone_seq IN (SELECT zone_seq FROM zone WHERE fld IN ({flds_sql}))"))
}
log_info("zone components + pct_area done for {length(comp_keys)} metrics")

4 Program-area composite score

Code

sc <- mseq_new("score_extriskspcat_primprod_ecoregionrescaled_equalweights", "Equal-weight composite (PA)")
seqs_sql <- paste(comp_seqs, collapse=",")
dbExecute(con, glue("INSERT INTO zone_metric (zone_seq, metric_seq, val)
  SELECT zone_seq, {sc}, SUM(val)/COUNT(val)
  FROM zone_metric
  WHERE metric_seq IN ({seqs_sql}) AND val IS NOT NULL
    AND zone_seq IN (SELECT zone_seq FROM zone WHERE fld IN ({flds_sql}))
  GROUP BY zone_seq"))

[1] 37

Code

# PA-only: the v7 equivalence gate compares Program Areas, and the composite is
# now computed for every zone set
v8_pa <- dbGetQuery(con, glue("SELECT z.val AS programarea_key, round(zm.val,2) AS score_v8
  FROM zone z JOIN zone_metric zm USING(zone_seq)
  WHERE zm.metric_seq={sc} AND z.fld='programarea_key' ORDER BY 1"))
# content fingerprint of the zone_metric output table (before disconnect)
zm_hash <- msens::hash_query(con, "zone_metric")
dbDisconnect(con, shutdown=TRUE)

5 Zone × taxon table (what the apps read)

zone_taxon is the per-zone species list the scores app shows (and downloads as CSV). It is precomputed here, not aggregated in the app: the server holds only the KB-sized serve.duckdb, whose model_cell is a view over S3 Parquet partitioned by mdl_id for per-model point reads (titiler tiles). A zone-wide aggregation there would have to list and scan the whole ~580M-row dataset over HTTPS, and in practice fails outright (IO Error: ... HTTP GET .../serve/model_cell/). v7 shipped such a table; v8 dropped it, which is what left the app’s Table of Species broken.

Rules live in msens::species_for_zone() / build_zone_taxon() and are unit-tested against synthetic v7 and v8 fixtures, so this notebook and the tests cannot drift.

Code

# the pa_score chunk closes `con`; reopen for this step
con <- dbConnect(duckdb(sdm_db))
n_zt <- msens::build_zone_taxon(con)
d_zt <- tbl(con, "zone_taxon") |>
  group_by(zone_fld) |>
  summarize(n_zones = n_distinct(zone_value), n_rows = n(), .groups = "drop") |>
  collect()
log_info("zone_taxon: {n_zt} rows across {sum(d_zt$n_zones)} zones")
dbDisconnect(con, shutdown = TRUE)
d_zt

# A tibble: 3 × 3
  zone_fld        n_zones n_rows
  <chr>             <dbl>  <dbl>
1 ecoregion_key        12  53876
2 programarea_key      20  39067
3 subregion_key         4  22757

6 v7 equivalence gate (pra_score_delta)

Code

v7_pa <- st_read(pa_gpkg, quiet=TRUE) |> st_drop_geometry() |>
  transmute(programarea_key, score_v7 = round(score_extriskspcat_primprod_ecoregionrescaled_equalweights, 2))
cmp <- v8_pa |> inner_join(v7_pa, by="programarea_key") |>
  mutate(delta = round(score_v8 - score_v7, 2)) |> arrange(desc(abs(delta)))
msens::report_table(cmp, caption = "program-area score: v8 vs v7")

program-area score: v8 vs v7
programarea_key	score_v8	score_v7	delta
SHU	19.39	29.7	-10.31
KOD	22.19	32.1	-9.91
GAB	18.14	25.7	-7.56
GOA	32.08	39.5	-7.42
SOC	22.58	29.7	-7.12
NOC	14.33	20.7	-6.37
GEO	28.94	33.9	-4.96
CEC	17.01	21.9	-4.89
NAV	22.60	27.3	-4.70
COK	47.93	52.6	-4.67
HOP	48.43	52.7	-4.27
ALB	14.94	12.3	2.64
CHU	23.11	24.9	-1.79
HAR	9.28	7.5	1.78
ALA	30.01	28.3	1.71
BFT	12.38	11.3	1.08
BOW	19.08	18.4	0.68
NOR	29.35	29.0	0.35
GAA	32.87	33.1	-0.23
MAT	24.85	24.9	-0.05

Code

gate <- with(cmp, list(n=nrow(cmp), mean_abs=round(mean(abs(delta)),2), max_abs=round(max(abs(delta)),2),
  rmse=round(sqrt(mean(delta^2)),2), cor=round(cor(score_v8, score_v7),3)))
str(gate)

List of 5
 $ n       : int 20
 $ mean_abs: num 4.12
 $ max_abs : num 10.3
 $ rmse    : num 5.18
 $ cor     : num 0.937

Code

# content-addressed: fingerprint of the zone_metric output table + deterministic gate stats
msens::write_manifest(
  manifest, target = "score_zone_metrics", content_hash = zm_hash,
  stats = c(list(ver = ver), gate),
  force = msens::force_target("score_zone_metrics"))