---
title: "Score cell metrics — extrisk_{spcat} → ecoregion rescale → primprod → composite"
msens:
target_name: score_cell_metrics
workflow_type: score
dependency: [score_zones]
output: data/manifests/score_cell_metrics.json
editor_options:
chunk_output_type: console
---
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).
## Design
```{mermaid}
%%| label: fig-design
%%| fig-cap: "Per-cell scoring: extrisk_{spcat} → ecoregion rescale, + primprod → equal-weight composite"
flowchart LR
mt["model_cell × taxon<br/>(er_score, sp_cat, is_marine)"] --> ex["extrisk_{spcat}<br/>Σ er·val/100"]
ex --> rs["ecoregion min/max<br/>→ _ecoregion_rescaled"]
vg["VGPM npp_avg → cells"] --> pp["primprod → rescaled"]
rs --> cm[("cell_metric<br/>sdm.duckdb")]
pp --> cm
cm --> co["equal-weight composite"] --> cm
cm --> mf["hash_query(cell_metric)<br/>→ content-addressed manifest"]
```
## Setup + metric registry
```{r}
#| label: setup
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)")
dbExecute(con, "CREATE TABLE cell_metric (cell_id INTEGER, metric_seq INTEGER, val DOUBLE)")
dbExecute(con, "CREATE TABLE zone_metric (zone_seq INTEGER, metric_seq INTEGER, val DOUBLE)")
# 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}")
```
## Stage 1a — extrisk_{spcat} per cell
```{r}
#| label: extrisk
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")
```
## Stage 1b/1c + primprod — ecoregion min/max + rescale
```{r}
#| label: 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"))
duckdb_unregister(con, "tmp_pp_df")
rescale_ecoregion("primprod")
log_info("stage 1b/1c + primprod done")
```
## Stage 3 — equal-weight composite per cell
```{r}
#| label: composite
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"))
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")
# content fingerprint of the full cell_metric output table (before disconnect)
cm_hash <- msens::hash_query(con, "cell_metric")
dbDisconnect(con, shutdown=TRUE)
```
## Manifest
```{r}
#| label: 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"))
```