---
title: "Merge prep — taxon crosswalk (each model → WORMS / BOTW taxon)"
msens:
target_name: merge_models_prep
workflow_type: merge
dependency: [ingest_worms, ingest_aquamaps, ingest_aquax, ingest_dps_nmfs, ingest_rng_iucn, ingest_birdlife_botw, ingest_rng_fws, ingest_ch_nmfs, ingest_ch_fws, ingest_turtles_swot_dps, ingest_ca_nmfs, ingest_listings] # reads worms for the WORMS/BOTW taxon crosswalk; add ingest_sdm_nc/ingest_sdm_gm here when their density is folded into the merge
output: data/manifests/merge_models_prep.json
editor_options:
chunk_output_type: console
---
Resolve every per-dataset model to a **common taxon**, so models of the same species
across datasets merge together. Two taxon namespaces (as in v7):
- **Birds → `BOTW:{sisid}`** — native BirdLife id. WoRMS has poor bird coverage, so
birds are NOT resolved through WoRMS.
- **Non-birds → `WORMS:{worms_id}`** — resolved from `scientific_name` by **reusing
v7's taxon table** (`scientific_name → worms_id`, cheap + already curated) and
`msens::match_taxa()` (ITIS→WoRMS crosswalk, exact WoRMS match, then WoRMS REST) for
species v7 didn't have.
Emits `taxon` (one row per taxon: id, authority, name, `worms_is_marine`…) and
`taxon_model` (`mdl_key ↔ taxon_id`), keyed for merge by
`mdl_key_merged(authority, taxon_id)` → `ms_merge|WORMS:…` / `ms_merge|BOTW:…`.
## Design
```{mermaid}
%%| label: fig-design
%%| fig-cap: "Crosswalk each model → common taxon (BOTW birds / WORMS rest) + governing extinction-risk floors"
flowchart LR
csv["dist/model_{ds}.csv<br/>∩ ingested mdl_keys"] --> res["resolve taxon:<br/>bird→BOTW:sisid, else WORMS:id<br/>(reuse v7 + match_taxa)"]
res --> tm[("taxon_model + taxon<br/>merge.duckdb")]
lst["listing + WoRMS Mammalia<br/>(MMPA/MBTA floors)"] --> er["governing er_score<br/>compute_er_score()"]
er --> tm
tm --> mf["hash_query(taxon)+hash_query(taxon_model)<br/>→ content-addressed manifest"]
```
## Setup + read all model crosswalks
```{r}
#| label: setup
librarian::shelf(
arrow, DBI, dplyr, duckdb, fs, glue, here, jsonlite, logger, purrr, readr, stringr,
tibble, MarineSensitivity/msens, quiet = T)
source(here("libs/paths.R")) # dir_big_v, spp_db, ver, ver_prev, dir_big
source(here("libs/vars.R"))
options(readr.show_col_types = F)
dir_atlas <- glue("{dir_big_v}/marine-atlas")
dir_dist <- glue("{dir_atlas}/dist")
v7_sdm_db <- glue("{dir_big}/{ver_prev}/sdm.duckdb") # reuse v7 taxon resolution
manifest <- here("data/manifests/merge_models_prep.json")
dir_create(path_dir(manifest))
# model_{ds}.csv is the dataset's species crosswalk (may list MORE than were
# ingested, e.g. FWS lists all 2,196 but ingests 106) — so key off the mdl_keys that
# actually have model_cell Parquet, joined to the CSVs for scientific_name.
sci_cols <- c("scientific_name", "sci_name", "sciname", "SCIENAME", "taxa")
lookup <- dir_ls(dir_dist, glob = "*model_*.csv") |>
set_names(\(f) str_remove(path_ext_remove(path_file(f)), "^model_")) |>
imap(\(f, ds) {
m <- read_csv(f)
# only a CROSSWALK belongs in this namespace: a file without mdl_key/scientific name is
# something else that landed on the glob, and it must not crash the merge prep
if (!"mdl_key" %in% names(m) || !length(intersect(sci_cols, names(m)))) {
log_warn("skipping {path_file(f)}: not a model crosswalk (no mdl_key / scientific name)"); return(NULL) }
sci <- intersect(sci_cols, names(m))[1]
iucn <- intersect(c("category", "code", "iucn"), names(m))[1] # IUCN Red List code (rng_iucn/bl)
tibble(
mdl_key = m$mdl_key,
ds_key = ds,
scientific_name = m[[sci]],
iucn_code = if (is.na(iucn)) NA_character_ else as.character(m[[iucn]]),
sisid = if ("sisid" %in% names(m)) as.character(m$sisid) else NA_character_,
# a dataset keyed by WoRMS AphiaID natively (AquaX `ax`) resolves by id, never by name
worms_id = if ("worms_id" %in% names(m)) as.character(m$worms_id) else NA_character_) }) |>
bind_rows() |> distinct(mdl_key, .keep_all = TRUE)
# read ingested mdl_keys per dataset so a mid-transfer / corrupt dataset (e.g. am
# actively rsyncing, with a truncated "no magic bytes" parquet) is skipped with a
# warning rather than crashing the whole glob — its taxa resolve on the next run.
con0 <- dbConnect(duckdb())
ds_dirs <- dir_ls(dir_dist, type = "directory", regexp = "dataset=")
ingested <- unlist(lapply(ds_dirs, \(d) tryCatch(
dbGetQuery(con0, glue("SELECT DISTINCT mdl_key FROM read_parquet('{d}/*.parquet')"))$mdl_key,
error = \(e) { log_warn("skip {path_file(d)}: {conditionMessage(e)}"); character(0) })))
dbDisconnect(con0, shutdown = TRUE)
models <- lookup |> filter(mdl_key %in% ingested)
# STABILITY: a model the previous version already crosswalked keeps its taxon by EXACT mdl_key.
# Re-resolving every model by cleaned scientific name against ver_prev's taxon table re-keyed 12
# models on the v9 control run (v8's table keeps ONE name per taxon where v7's kept one per model,
# so a name can now land on a synonym id), moved one taxon's ER and broke the plumbing check by 3
# taxa. Name / native-id resolution below is for models ver_prev never had (ax, new ranges).
prev_tm_pq <- path.expand(glue("{dir_big}/{ver_prev}/marine-atlas/tables/taxon_model.parquet"))
prev_tm <- if (file_exists(prev_tm_pq)) {
arrow::read_parquet(prev_tm_pq) |>
transmute(mdl_key, prev_authority = taxon_authority, prev_taxon_id = as.character(taxon_id))
} else { log_warn("no ver_prev taxon_model at {prev_tm_pq} — every model resolves afresh"); tibble(mdl_key = character(), prev_authority = character(), prev_taxon_id = character()) }
models <- models |> left_join(prev_tm, by = "mdl_key")
log_info("crosswalk reuse: {sum(!is.na(models$prev_taxon_id))} models keep {ver_prev}'s taxon by mdl_key; {sum(is.na(models$prev_taxon_id))} resolve afresh")
# the CONTROL run (AX_SUPERSEDE=0) merges EXACTLY ver_prev's model set, so its hash is binary
if (!ax_supersede) {
n0 <- nrow(models); models <- models |> filter(!is.na(prev_taxon_id))
log_info("control run: {n0 - nrow(models)} model(s) {ver_prev} never crosswalked are EXCLUDED so the surface can reproduce its checkpoint")
}
# density datasets (gm, nc) are ingested + published natively, but are NOT yet folded into the merge:
# merge treats every non-`am` dataset as a binary range (max-merge), which would misrepresent graded
# density and silently change scoring (a composite change to validate with `pra_score_delta`). So
# exclude them here even when their `dist/dataset=` Parquet is present on disk, until MERGE_FOLD_DENSITY=1.
merge_exclude_ds <- if (nzchar(Sys.getenv("MERGE_FOLD_DENSITY"))) character(0) else c("gm", "nc")
# the v9 CONTROL run (AX_SUPERSEDE=0): AquaX is left OUT of the merge entirely, so the inputs are
# exactly ver_prev's and the merged surface must reproduce its checkpoint hash. (Leaving ax in but
# unsuperseded would still max() it with AquaMaps -- that is not a control, it is a third design.)
if (!ax_supersede) merge_exclude_ds <- c(merge_exclude_ds, "ax")
if (length(merge_exclude_ds)) {
n0 <- nrow(models); models <- models |> filter(!ds_key %in% merge_exclude_ds)
if (nrow(models) < n0) log_info(
"excluded {n0 - nrow(models)} density model(s) [{paste(merge_exclude_ds, collapse='/')}] from merge — set MERGE_FOLD_DENSITY=1 to fold in")
}
count(models, ds_key)
```
## Non-bird taxa → worms_id (reuse v7, then match_taxa)
```{r}
#| label: worms
nonbird <- models |> filter(ds_key != "bl", is.na(worms_id), is.na(prev_taxon_id)) |> distinct(scientific_name) |>
mutate(sci_clean = clean_sci_name(scientific_name), # drop (=synonym), ssp. markers
sci_binom = clean_sci_name(scientific_name, binomial = TRUE)) # Genus species fallback
# reuse v7's taxon resolution, matched on the cleaned name
con_v7 <- dbConnect(duckdb(v7_sdm_db, read_only = TRUE))
# v1-v7 carried `worms_id`; v8+ carry the id as `taxon_id` under taxon_authority = 'worms'.
# Resolved by introspecting the column set, never by version string.
v7_tax <- if ("worms_id" %in% dbListFields(con_v7, "taxon")) {
tbl(con_v7, "taxon") |> filter(taxon_authority == "worms") |>
select(scientific_name, worms_id, worms_is_marine, worms_is_extinct) |> collect()
} else {
tbl(con_v7, "taxon") |> filter(taxon_authority == "worms") |>
select(scientific_name, worms_id = taxon_id, worms_is_marine, worms_is_extinct) |> collect() |>
mutate(worms_id = as.numeric(worms_id))
}
dbDisconnect(con_v7, shutdown = TRUE)
v7_tax <- v7_tax |> mutate(sci_clean = clean_sci_name(scientific_name)) |>
distinct(sci_clean, .keep_all = TRUE) |> select(-scientific_name)
resolved <- nonbird |> left_join(v7_tax, by = "sci_clean")
log_info("non-bird taxa: {nrow(resolved)} | reused from v7: {sum(!is.na(resolved$worms_id))}")
# match_taxa (ITIS->WoRMS, exact, REST) for the rest: cleaned name, then binomial fallback
match_names <- function(v) {
con_spp <- dbConnect(duckdb(spp_db, read_only = TRUE)); on.exit(dbDisconnect(con_spp, shutdown = TRUE))
match_taxa(tibble(scientific_name = unique(v)), con_spp) |>
filter(!is.na(worms_id)) |> distinct(scientific_name, .keep_all = TRUE)
}
for (col in c("sci_clean", "sci_binom")) {
todo <- resolved |> filter(is.na(worms_id)) |> pull(.data[[col]]) |> unique()
if (length(todo) == 0) next
m <- match_names(todo) |> transmute(!!col := scientific_name, wid = worms_id)
resolved <- resolved |> left_join(m, by = col) |>
mutate(worms_id = coalesce(worms_id, wid)) |> select(-wid)
}
log_info("non-bird taxa resolved: {sum(!is.na(resolved$worms_id))} of {nrow(resolved)}")
```
## Build taxon + taxon_model
```{r}
#| label: taxon
# BOTW-sisid <-> WoRMS-aphia crosswalk (birds): BirdLife is the bird authority, so BOTW:{sisid}
# is a bird's canonical taxon. Map cleaned scientific_name -> sisid so a bird's NON-bl datasets
# (FWS ranges / critical habitat, AquaMaps — which resolve to WORMS by name) merge into that SAME
# BOTW model instead of splitting the species across BOTW + WORMS (the duplicate-picker bug). Match
# on the cleaned name, with a binomial fallback for authorship/subspecies differences.
bird_xw <- models |> filter(ds_key == "bl", !is.na(sisid)) |>
transmute(sci_clean = clean_sci_name(scientific_name),
sci_binom = clean_sci_name(scientific_name, binomial = TRUE),
sisid = as.character(sisid))
# distinct maps (s_cln/s_bin names avoid colliding with models' own all-NA `sisid` column)
sisid_by_clean <- bird_xw |> distinct(sci_clean, .keep_all = TRUE) |> transmute(sci_clean, s_cln = sisid)
sisid_by_binom <- bird_xw |> distinct(sci_binom, .keep_all = TRUE) |> transmute(sci_binom, s_bin = sisid)
# non-bl models: resolve to worms by name, THEN override to the bird's BOTW:{sisid} when the
# species is a BirdLife bird (so all of a bird's datasets land on one BOTW model).
tm_nonbl <- models |> filter(ds_key != "bl") |>
mutate(sci_clean = clean_sci_name(scientific_name),
sci_binom = clean_sci_name(scientific_name, binomial = TRUE)) |>
left_join(resolved |> select(scientific_name, worms_id_name = worms_id), by = "scientific_name") |>
# a native WoRMS id (ax) wins over the name resolution; `format()` not `as.character()`, which
# renders 1e+06-style ids for large AphiaIDs
mutate(worms_id = coalesce(worms_id, ifelse(is.na(worms_id_name), NA_character_,
format(worms_id_name, scientific = FALSE, trim = TRUE)))) |>
left_join(sisid_by_clean, by = "sci_clean") |>
left_join(sisid_by_binom, by = "sci_binom") |>
mutate(bird_sisid = coalesce(s_cln, s_bin), # bird? (clean match, else binomial)
taxon_authority = if_else(!is.na(bird_sisid), "botw", "worms"),
taxon_id = if_else(!is.na(bird_sisid), bird_sisid, worms_id),
# ver_prev's crosswalk wins for a model it already had (stability across versions)
taxon_authority = coalesce(prev_authority, taxon_authority),
taxon_id = coalesce(prev_taxon_id, taxon_id)) |>
filter(!is.na(taxon_id)) |>
transmute(mdl_key, taxon_authority, taxon_id)
# bl models -> botw taxon (native sisid; ver_prev's key wins where it exists)
tm_botw <- models |> filter(ds_key == "bl", !is.na(sisid) | !is.na(prev_taxon_id)) |>
transmute(mdl_key, taxon_authority = coalesce(prev_authority, "botw"),
taxon_id = coalesce(prev_taxon_id, as.character(sisid)))
taxon_model <- bind_rows(tm_nonbl, tm_botw) |>
mutate(ms_merge_key = mdl_key_merged(taxon_authority, taxon_id))
log_info("bird crosswalk: {sum(tm_nonbl$taxon_authority=='botw')} non-bl models re-keyed BOTW (would split otherwise)")
taxon <- taxon_model |>
left_join(models |> select(mdl_key, scientific_name, iucn_code), by = "mdl_key") |>
group_by(taxon_authority, taxon_id, ms_merge_key) |>
summarise(
scientific_name = first(scientific_name),
iucn_code = first(c(iucn_code[!is.na(iucn_code)], NA_character_)), # IUCN category if any model has it
n_models = n_distinct(mdl_key),
n_datasets = n_distinct(str_extract(mdl_key, "^[^|]+")), .groups = "drop") |>
left_join(resolved |> select(scientific_name, worms_is_marine, worms_is_extinct),
by = "scientific_name")
con <- dbConnect(duckdb(glue("{dir_atlas}/merge.duckdb")))
dbWriteTable(con, "taxon_model", taxon_model, overwrite = TRUE)
dbWriteTable(con, "taxon", taxon, overwrite = TRUE)
dbDisconnect(con, shutdown = TRUE)
glue("{nrow(taxon)} taxa ({sum(taxon$taxon_authority=='worms')} WORMS + ",
"{sum(taxon$taxon_authority=='botw')} BOTW) from {nrow(taxon_model)} models")
```
## Governing extinction-risk score (`er_score`)
Each **ingest** layer holds the ER *as its own dataset knew it* — IUCN category for
`rng_iucn`/BOTW, FWS status for `rng_fws`/`ch_fws`, NMFS for SWOT/critical-habitat —
each via `msens::compute_er_score()`. Those source values can be **out of date** (an
older IUCN category) or **out of context** (a US national NMFS/FWS listing, or an
MMPA/MBTA floor, that should override the IUCN international score). So here we compute
each taxon's **governing** `er_score`: the most-protective across its datasets, US
national overriding IUCN, with MMPA/MBTA floors — again via `compute_er_score()`.
`merge_models` applies this to the **range-cell values** (the fitting point where a
species' presence should carry its governing ER, not the raw source code); v7 followed
the same model. Needs the `listing` table (`ingest_listings`).
The two statutory floors differ in how they are assigned:
- **MMPA (`is_mmpa`, +20 floor)** — the Marine Mammal Protection Act protects *all*
marine mammals, so we assign it by **taxonomy**: every WoRMS class **Mammalia** taxon
(matched by AphiaID = `taxon_id`), not the incomplete NMFS directory (which lists only
~70 managed stocks vs 137 marine-mammal taxa here).
- **MBTA (`is_mbta`, +10 floor)** — the Migratory Bird Treaty Act protects only
**migratory birds native to the US / US territories**, an explicit species list, so we
assign it from the authoritative **FWS CFR 50 §10.13** list (`ingest_listings`, matched
by `clean_sci_name`). It is **not** all Aves — a bird absent from 10.13 (e.g. a
non-US / non-migratory species) gets **no** MBTA floor. Most 10.13 birds are terrestrial;
the marine-relevance filter (`score_zones`) is what keeps terrestrial birds out of the
score, while this floor only sets the *value* of those that remain.
```{r}
#| label: er_score
con <- dbConnect(duckdb(glue("{dir_atlas}/merge.duckdb")))
tx <- dbGetQuery(con, "SELECT ms_merge_key, taxon_authority, taxon_id, scientific_name, iucn_code FROM taxon") |>
mutate(sci = clean_sci_name(scientific_name))
# US ESA + MBTA + BCC stay species-specific from the listing. is_mbta is the FWS CFR
# 50 §10.13 list (migratory birds native to the US) — an explicit species list, NOT all
# Aves; a bird absent from 10.13 gets no MBTA floor.
if ("listing" %in% dbListTables(con)) {
tx <- tx |> left_join(dbGetQuery(con, "SELECT sci, nmfs_esa, fws_esa, is_mbta, is_bcc FROM listing"), by = "sci")
} else {
# a HARD stop, not a warning: this used to fall back to IUCN-only, which silently dropped every
# US ESA listing and every MBTA floor from the governing er_score. It fired on v9's first merge --
# `listing` lives in merge.duckdb, which a new version starts empty (bootstrap-release skill).
stop("no `listing` table in merge.duckdb — render ingest_nmfs-fws-listings.qmd (ingest_listings) for this version first")
}
# MMPA is assigned by TAXONOMY (all marine mammals are protected): every WoRMS class
# Mammalia taxon (AphiaID = taxon_id), more complete than the ~70-stock NMFS directory.
con_spp <- dbConnect(duckdb(spp_db, read_only = TRUE))
mammalia <- dbGetQuery(con_spp, "SELECT DISTINCT CAST(scientificNameID AS VARCHAR) aphia FROM worms WHERE class='Mammalia'")$aphia
dbDisconnect(con_spp, shutdown = TRUE)
esa_rank <- function(x) match(coalesce(x, "LC"), c("LC", "TN", "EN"))
tx <- tx |> mutate(
is_bcc = coalesce(is_bcc, FALSE),
is_mmpa = taxon_authority == "worms" & taxon_id %in% mammalia, # MMPA: all marine mammals
is_mbta = coalesce(is_mbta, FALSE) & taxon_authority == "botw", # MBTA: CFR 10.13, birds only
us_code = c("LC", "TN", "EN")[pmax(esa_rank(nmfs_esa), esa_rank(fws_esa))],
is_us = us_code != "LC" | is_mmpa | is_mbta,
extrisk_code = case_when(
is_us & is_mmpa ~ paste0("NMFS:", us_code), # US national + MMPA overrides IUCN
is_us ~ paste0("FWS:", us_code),
iucn_code %in% c("CR","EN","VU","NT","LC","DD") ~ paste0("IUCN:", iucn_code),
TRUE ~ NA_character_))
tx$er_score <- 1L
ok <- !is.na(tx$extrisk_code)
tx$er_score[ok] <- compute_er_score(tx$extrisk_code[ok], is_mmpa = tx$is_mmpa[ok], is_mbta = tx$is_mbta[ok])
dbWriteTable(con, "taxon_er",
tx |> select(ms_merge_key, extrisk_code, er_score, is_mmpa, is_mbta, is_bcc), overwrite = TRUE)
dbExecute(con, "CREATE OR REPLACE TABLE taxon AS
SELECT t.*, e.extrisk_code, e.er_score, e.is_mmpa, e.is_mbta, e.is_bcc
FROM taxon t LEFT JOIN taxon_er e USING (ms_merge_key)")
log_info("taxonomy floors: {sum(tx$is_mmpa)} Mammalia (MMPA), {sum(tx$is_mbta)} Aves (MBTA); {sum(tx$is_us)} US-listed of {nrow(tx)}")
dbDisconnect(con, shutdown = TRUE)
```
## Outputs + manifest
```{r}
#| label: manifest
# content-addressed: fingerprint of the taxon + taxon_model output tables (reopen
# read-only since the er_score chunk already disconnected)
con <- dbConnect(duckdb(glue("{dir_atlas}/merge.duckdb"), read_only = TRUE))
h <- paste0(msens::hash_query(con, "taxon"), msens::hash_query(con, "taxon_model"))
dbDisconnect(con, shutdown = TRUE)
smry <- tibble(n_taxa = nrow(taxon), n_models = nrow(taxon_model),
n_worms = sum(taxon$taxon_authority == "worms"),
n_botw = sum(taxon$taxon_authority == "botw"))
msens::report_table(smry, caption = "merge_models_prep: taxa + models")
msens::write_manifest(
manifest, target = "merge_models_prep", content_hash = h,
stats = list(ver = ver, n_models = smry$n_models, n_taxa = smry$n_taxa,
n_worms = smry$n_worms, n_botw = smry$n_botw),
force = msens::force_target("merge_models_prep"))
```