Score cell metrics — extrisk_{spcat} → ecoregion rescale → primprod → composite

Per-cell scoring (replicates v7 calc_scores.qmd exactly):

1. extrisk_{spcat} = SUM(er_score × val / 100) over valid species in the category (is_er_spatial/turtles → er_score = 100; MMPA floor already baked into er_score).
2. ecoregion min/max of that metric per ecoregion; _ecoregion_rescaled per cell = (x − eco_min)/(eco_max − eco_min) × 100, coverage-normalized across ecoregions.
3. primprod = VGPM npp_avg bilinear→cells; same ecoregion min/max + rescale.
4. composite score_extriskspcat_primprod_ecoregionrescaled_equalweights = ROUND(AVG(present rescaled components)) (7 extrisk + primprod, equal weight).

1 Design

Figure 1: Per-cell scoring: extrisk_{spcat} → ecoregion rescale, + primprod → equal-weight composite

2 Setup + metric registry

librarian::shelf(DBI, dplyr, duckdb, fs, glue, here, jsonlite, logger, terra,
                 MarineSensitivity/msens, quiet = T)
source(here("libs/paths.R")); source(here("libs/vars.R"))
manifest <- here("data/manifests/score_cell_metrics.json"); dir_create(path_dir(manifest))
vgpm_tif <- glue("{dir_data}/raw/oregonstate.edu/vgpm.r2022.v.chl.v.sst.2160x4320_2014-2023.avg.sd.tif")
stopifnot(all(file_exists(c(sdm_db, cellid_tif, vgpm_tif))))

con <- dbConnect(duckdb(sdm_db))
# fresh recreate (a prior run may have left these with a `value` column, not `val`)
for (t in c("cell_metric", "zone_metric", "metric")) dbExecute(con, paste("DROP TABLE IF EXISTS", t))
dbExecute(con, "CREATE TABLE metric (metric_seq INTEGER PRIMARY KEY, metric_key VARCHAR, description VARCHAR)")

[1] 0

dbExecute(con, "CREATE TABLE cell_metric (cell_id INTEGER, metric_seq INTEGER, val DOUBLE)")

[1] 0

dbExecute(con, "CREATE TABLE zone_metric (zone_seq INTEGER, metric_seq INTEGER, val DOUBLE)")

[1] 0

# idempotent metric_seq: fetch or insert, wiping any prior values for a clean rebuild
mseq <- function(key, desc = "") {
  ex <- dbGetQuery(con, glue("SELECT metric_seq FROM metric WHERE metric_key = '{key}'"))
  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}, '{key}', '{gsub(\"'\",\"\",desc)}')")); s0 }
  dbExecute(con, glue("DELETE FROM cell_metric WHERE metric_seq = {s}"))
  dbExecute(con, glue("DELETE FROM zone_metric WHERE metric_seq = {s}"))
  s
}

# 7 scoring categories (taxonomy-based sp_cat from merge_taxon). primary_producer replaces
# the old catch-all 'other'; reptiles/amphibians are NOT in this list, so they're excluded.
sp_cats <- c("bird","coral","fish","invertebrate","mammal","primary_producer","turtle")

# species set: default = US-valid AND marine-relevant (is_marine culls terrestrial birds;
# non-birds are TRUE). SCORE_V7COMMON=1 -> restrict to v7's scored set (in_v7) for an
# apples-to-apples v7 comparison that isolates grid/algorithm from the v8 species-set change.
# SCORE_ALLBIRDS=1 -> drop the marine filter (keep every bird) for diagnostics.
taxa_where <- if (nzchar(Sys.getenv("SCORE_V7COMMON"))) "t.is_valid_usa AND t.in_v7" else "t.is_valid_usa"
if (!nzchar(Sys.getenv("SCORE_ALLBIRDS"))) taxa_where <- paste(taxa_where, "AND t.is_marine")
log_info("scoring taxa filter: {taxa_where}")

3 Stage 1a — extrisk_{spcat} per cell

for (cat in sp_cats) {
  s <- mseq(glue("extrisk_{cat}"), glue("Extinction-risk-weighted suitability, {cat}"))
  dbExecute(con, glue("INSERT INTO cell_metric (cell_id, metric_seq, val)
    SELECT mc.cell_id, {s}, ROUND(SUM(
        (CASE WHEN t.is_er_spatial THEN 100 ELSE t.er_score END) * mc.val) / 100.0, 2) AS val
    FROM model_cell mc JOIN taxon t ON mc.mdl_key = t.ms_merge_key
    WHERE {taxa_where} AND t.sp_cat = '{cat}'
    GROUP BY mc.cell_id"))
}
log_info("stage 1a extrisk: {dbGetQuery(con,'SELECT count(*) n FROM cell_metric')$n} cell-metric rows")

4 Stage 1b/1c + primprod — ecoregion min/max + rescale

# ecoregion min/max of a base metric -> zone_metric; then per-cell rescale -> cell_metric
rescale_ecoregion <- function(base_key) {
  bs   <- dbGetQuery(con, glue("SELECT metric_seq FROM metric WHERE metric_key='{base_key}'"))$metric_seq[1]
  smin <- mseq(glue("{base_key}_ecoregion_min")); smax <- mseq(glue("{base_key}_ecoregion_max"))
  srs  <- mseq(glue("{base_key}_ecoregion_rescaled"))
  for (fx in list(c(smin,"min"), c(smax,"max")))
    dbExecute(con, glue("INSERT INTO zone_metric (zone_seq, metric_seq, val)
      SELECT z.zone_seq, {fx[1]}, {fx[2]}(cm.val)
      FROM zone z JOIN zone_cell zc USING(zone_seq) JOIN cell_metric cm ON zc.cell_id=cm.cell_id
      WHERE z.fld='ecoregion_key' AND cm.metric_seq={bs} GROUP BY z.zone_seq"))
  dbExecute(con, glue("INSERT INTO cell_metric (cell_id, metric_seq, val)
    WITH cell_ecoregion AS (
      SELECT zc.cell_id, zc.zone_seq,
             zc.pct_covered * 100.0 / SUM(zc.pct_covered) OVER (PARTITION BY zc.cell_id) AS norm_pct
      FROM zone_cell zc JOIN zone z USING(zone_seq) WHERE z.fld='ecoregion_key'),
    mm AS (SELECT mn.zone_seq, mn.val min_v, mx.val max_v
           FROM (SELECT zone_seq,val FROM zone_metric WHERE metric_seq={smin}) mn
           JOIN (SELECT zone_seq,val FROM zone_metric WHERE metric_seq={smax}) mx USING(zone_seq))
    SELECT cm.cell_id, {srs},
      SUM((cm.val - mm.min_v)/(mm.max_v - mm.min_v) * (ce.norm_pct/100.0)) * 100 AS val
    FROM cell_metric cm JOIN cell_ecoregion ce ON cm.cell_id=ce.cell_id JOIN mm ON ce.zone_seq=mm.zone_seq
    WHERE cm.metric_seq={bs} AND mm.min_v < mm.max_v
    GROUP BY cm.cell_id"))
}
for (cat in sp_cats) rescale_ecoregion(glue("extrisk_{cat}"))

# primprod: VGPM npp_avg bilinear -> cells (restricted to US cells) -> cell_metric, then rescale
sp <- mseq("primprod", "Primary productivity VGPM/VIIRS npp_avg (mg C/m2/day)")
pp <- cells_from_raster(subset(rast(vgpm_tif), "npp_avg"), cellid_tif,
                        method = "bilinear", min_value = 0, zero_fill = FALSE)
duckdb_register(con, "tmp_pp_df", pp)
dbExecute(con, glue("INSERT INTO cell_metric (cell_id, metric_seq, val)
  SELECT p.cell_id, {sp}, p.val FROM tmp_pp_df p JOIN cell c ON p.cell_id = c.cell_id WHERE c.in_usa"))

[1] 623616

duckdb_unregister(con, "tmp_pp_df")
rescale_ecoregion("primprod")

[1] 617578

log_info("stage 1b/1c + primprod done")

5 Stage 3 — equal-weight composite per cell

comp_keys <- c(glue("extrisk_{sp_cats}_ecoregion_rescaled"), "primprod_ecoregion_rescaled")
sc <- mseq("score_extriskspcat_primprod_ecoregionrescaled_equalweights", "Equal-weight composite")
in_list <- paste(sprintf("'%s'", comp_keys), collapse=",")
dbExecute(con, glue("INSERT INTO cell_metric (cell_id, metric_seq, val)
  SELECT cm.cell_id, {sc}, ROUND(AVG(cm.val))
  FROM cell_metric cm JOIN metric m USING(metric_seq)
  WHERE m.metric_key IN ({in_list}) GROUP BY cm.cell_id"))

[1] 623212

smry <- dbGetQuery(con, glue("SELECT count(*) n_cells, round(min(val),1) v_min, round(max(val),1) v_max,
  round(avg(val),1) v_avg FROM cell_metric WHERE metric_seq={sc}"))
msens::report_table(smry, caption = "composite score per cell")

composite score per cell

n_cells	v_min	v_max	v_avg
623212	0	96	22.1

# content fingerprint of the full cell_metric output table (before disconnect)
cm_hash <- msens::hash_query(con, "cell_metric")
dbDisconnect(con, shutdown=TRUE)

6 Manifest

# content-addressed: fingerprint of the cell_metric output table (no wall-clock)
msens::write_manifest(
  manifest, target = "score_cell_metrics", content_hash = cm_hash,
  stats = list(ver = ver, composite_cells = smry$n_cells,
               v_min = smry$v_min, v_max = smry$v_max),
  force = msens::force_target("score_cell_metrics"))