---
title: "Merge taxon metrics — validity flags + range / rarity / endemism / pct_marine"
msens:
target_name: merge_taxon
workflow_type: merge
dependency: [merge_models, build_cell_grid, ingest_worms] # reads worms for taxonomy-based sp_cat
output: data/manifests/merge_taxon.json
editor_options:
chunk_output_type: console
---
Per-taxon attributes on the merged `model_cell`, joined to the `cell` table (ocean
cells with `area_km2` / `in_usa` / `in_pra`):
- **3 nested validity flags** (`is_valid_global ⊇ is_valid_usa ⊇ is_valid_pra`):
a taxon is *valid_global* if it has a mapped distribution reaching the ocean, isn't
extinct, and isn't flagged non-marine (`worms_is_marine`); *valid_usa* / *valid_pra*
add "has cells in the US study area / Program Areas".
- **range / rarity / endemism / pct_marine**: `range_km2` = Σ `area_km2` over the
taxon's ocean cells (v7-comparable marine range); `pct_marine` = ocean-cell share of
the whole (land+ocean) range; `us_endemism` = US share of the marine range; rarity
class from the marine `range_km2` distribution.
- **`is_marine`** (marine-relevance filter): keeps genuinely marine/coastal **birds** in the
score and culls terrestrial ones (v8 captures each bird's whole, largely-terrestrial range).
## Design
```{mermaid}
%%| label: fig-design
%%| fig-cap: "Per-taxon attributes on the merged model_cell: validity flags, range/rarity, sp_cat, is_marine"
flowchart LR
mc["dist_merged model_cell<br/>∩ cell (area/in_usa/in_pra)"] --> tc["taxon_cell tallies<br/>n_cells/n_ocean/range_km2"]
tc --> flg["validity flags +<br/>rarity + pct_marine"]
wrm["WoRMS class/phylum<br/>(+ API supplement)"] --> cat["sp_cat by taxonomy"]
bl["raw bl ranges +<br/>marine families / curation"] --> mar["is_marine (birds)"]
flg --> tx[("taxon<br/>merge.duckdb")]
cat --> tx
mar --> tx
tx --> mf["hash_query(taxon)<br/>→ content-addressed manifest"]
```
## Setup + compute
```{r}
#| label: compute
librarian::shelf(DBI, dplyr, duckdb, fs, glue, here, jsonlite, logger, readr, stringr,
tibble, ropensci/worrms, MarineSensitivity/msens, quiet = T)
source(here("libs/paths.R"))
source(here("libs/vars.R"))
dir_atlas <- glue("{dir_big_v}/marine-atlas")
merged <- glue("{dir_atlas}/dist_merged/dataset=ms_merge")
merged_global <- glue("{dir_atlas}/dist_merged_global/dataset=ms_merge") # whole-range → is_valid_global
dir_dist <- glue("{dir_atlas}/dist") # raw per-dataset (global) ranges, for bird pct_marine
merge_db <- glue("{dir_atlas}/merge.duckdb")
manifest <- here("data/manifests/merge_taxon.json")
dir_create(path_dir(manifest))
stopifnot("run merge_models first" = dir_exists(merged),
"run build_cell_grid first" = file_exists(sdm_db))
con <- dbConnect(duckdb(merge_db))
dbExecute(con, glue("ATTACH '{sdm_db}' AS grid (READ_ONLY)"))
dbExecute(con, glue("ATTACH '{spp_db}' AS spp (READ_ONLY)")) # WoRMS class/phylum for sp_cat
# per-taxon cell tallies against the ocean cell table
dbExecute(con, glue("CREATE OR REPLACE TABLE taxon_cell AS
WITH mc AS (SELECT mdl_key, cell_id FROM read_parquet('{merged}/**/*.parquet')),
j AS (SELECT mc.mdl_key, c.area_km2, c.in_usa, c.in_pra,
(c.cell_id IS NOT NULL) AS is_ocean
FROM mc LEFT JOIN grid.cell c USING (cell_id))
SELECT mdl_key,
count(*) AS n_cells,
count(*) FILTER (WHERE is_ocean) AS n_ocean,
count(*) FILTER (WHERE in_usa) AS n_usa,
count(*) FILTER (WHERE in_pra) AS n_pra,
round(sum(area_km2) FILTER (WHERE is_ocean), 1) AS range_km2,
round(sum(area_km2) FILTER (WHERE in_usa), 1) AS range_usa_km2
FROM j GROUP BY mdl_key"))
# whole-range (global) ocean-cell tally → is_valid_global, so a species valid OUTSIDE the US (e.g.
# Sotalia guianensis, IUCN range wholly non-US) is still is_valid_global — it appears in the app's
# "US waters"-UNCHECKED picker and renders via its merged whole-range COG. am-only taxa aren't in the
# global surface (they reuse am COGs), so those fall back to the US n_ocean in is_valid_global below.
if (dir_exists(merged_global)) {
dbExecute(con, glue("CREATE OR REPLACE TABLE taxon_cell_global AS
SELECT mc.mdl_key, count(*) FILTER (WHERE c.cell_id IS NOT NULL) AS n_global
FROM read_parquet('{merged_global}/**/*.parquet') mc LEFT JOIN grid.cell c USING (cell_id)
GROUP BY mc.mdl_key"))
} else dbExecute(con, "CREATE OR REPLACE TABLE taxon_cell_global (mdl_key VARCHAR, n_global BIGINT)")
```
## Flags + range/rarity/endemism → taxon
```{r}
#| label: flags
# idempotent: drop any derived columns a prior run added so `x.*` stays base taxon (+ er)
for (col in c("n_cells","n_ocean","n_global","n_usa","n_pra","range_km2","range_usa_km2",
"is_valid_global","is_valid_usa","is_valid_pra","pct_marine","us_endemism",
"rarity","is_marine","pct_marine_bl","sp_cat"))
try(dbExecute(con, glue("ALTER TABLE taxon DROP COLUMN IF EXISTS {col}")), silent = TRUE)
# rarity thresholds on the marine range (km2); tune vs v7
dbExecute(con, "CREATE OR REPLACE TABLE taxon AS
WITH t AS (
SELECT x.*, tc.n_cells, tc.n_ocean, tc.n_usa, tc.n_pra, tc.range_km2, tc.range_usa_km2,
tcg.n_global
FROM taxon x LEFT JOIN taxon_cell tc ON x.ms_merge_key = tc.mdl_key
LEFT JOIN taxon_cell_global tcg ON x.ms_merge_key = tcg.mdl_key)
SELECT *,
-- 3 nested validity flags (global superset of usa superset of pra); is_valid_global uses the
-- WHOLE-RANGE tally (n_global) so non-US species count, falling back to US n_ocean for am-only taxa
((coalesce(n_global, 0) > 0 OR coalesce(n_ocean, 0) > 0)
AND coalesce(worms_is_extinct, FALSE) = FALSE
AND coalesce(worms_is_marine, TRUE) = TRUE) AS is_valid_global,
(n_usa > 0) AS is_valid_usa,
(n_pra > 0) AS is_valid_pra,
round(100.0 * n_ocean / nullif(n_cells, 0), 1) AS pct_marine,
round(100.0 * n_usa / nullif(n_ocean, 0), 1) AS us_endemism,
CASE
WHEN range_km2 IS NULL THEN NULL
WHEN range_km2 < 10000 THEN 'very_rare'
WHEN range_km2 < 100000 THEN 'rare'
WHEN range_km2 < 1000000 THEN 'uncommon'
WHEN range_km2 < 10000000 THEN 'common'
ELSE 'widespread' END AS rarity
FROM t")
smry <- dbGetQuery(con, "SELECT
count(*) n_taxa,
count(*) FILTER (WHERE is_valid_global) valid_global,
count(*) FILTER (WHERE is_valid_usa) valid_usa,
count(*) FILTER (WHERE is_valid_pra) valid_pra,
round(median(pct_marine),1) med_pct_marine FROM taxon")
knitr::kable(smry)
```
## Species category (`sp_cat`) by taxonomy
`sp_cat` is assigned by **taxonomy** (WoRMS class/phylum/kingdom + BirdLife→bird), not a v7
crosswalk — so new taxa land in the right class instead of a catch-all `other` (eliminated).
v7's six scored classes (bird, mammal, turtle, fish, coral, invertebrate) **plus
`primary_producer`** (marine plants/algae/diatoms — distinct from the carbon `primprod`
metric). **Reptiles/amphibians are categorized but EXCLUDED from scoring** (sea turtles stay
`turtle`; `score_cell_metrics` scores only the seven classes). **AlgaeBase-sourced taxa
(macroalgae/seagrass) are license-excluded from the downloadable WoRMS DwC-A**, so
worms-authority AphiaIDs absent from the local snapshot are fetched from the **WoRMS API**
(`worrms::wm_record`) and cached in `data/worms_taxonomy_supplement.csv` — self-maintaining
(only new ids are fetched; offline runs use the cache).
```{r}
#| label: sp_cat
# Self-maintaining WoRMS supplement: AlgaeBase-sourced taxa (macroalgae/seagrass) are
# excluded from the downloadable WoRMS DwC-A by license, so any worms-authority taxon used
# here but absent from the local snapshot is fetched from the WoRMS API and cached. Only
# *new* unmatched AphiaIDs are fetched; offline runs fall back to the existing cache.
supp_csv <- here("data/worms_taxonomy_supplement.csv")
supp <- if (file_exists(supp_csv))
read_csv(supp_csv, show_col_types = FALSE, col_types = cols(taxon_id = "c")) else
tibble(taxon_id = character(), scientificName = character(),
kingdom = character(), phylum = character(), class = character(), order = character())
need <- setdiff(dbGetQuery(con, "
SELECT DISTINCT t.taxon_id FROM taxon t
LEFT JOIN spp.worms w ON t.taxon_authority='worms' AND CAST(w.scientificNameID AS VARCHAR)=t.taxon_id
WHERE t.taxon_authority='worms' AND w.scientificNameID IS NULL")$taxon_id, supp$taxon_id)
if (length(need)) {
log_info("WoRMS API: fetching {length(need)} AphiaIDs missing from the snapshot (AlgaeBase-excluded etc.)")
fetched <- tryCatch({
ids <- as.integer(need); ids <- ids[!is.na(ids)]
bind_rows(lapply(split(ids, ceiling(seq_along(ids) / 50)), \(b) tryCatch({
r <- worrms::wm_record(b)
tibble(taxon_id = as.character(r$AphiaID), scientificName = r$scientificname,
kingdom = r$kingdom, phylum = r$phylum, class = r$class, order = r$order)
}, error = \(e) NULL)))
}, error = \(e) { log_warn("WoRMS API fetch failed ({conditionMessage(e)}); using cached supplement"); NULL })
if (!is.null(fetched) && nrow(fetched)) {
supp <- bind_rows(supp, fetched) |> distinct(taxon_id, .keep_all = TRUE)
write_csv(supp, supp_csv)
log_info("supplement now {nrow(supp)} taxa ({nrow(fetched)} newly fetched + cached)")
}
}
supp <- supp |> transmute(taxon_id = as.character(taxon_id),
s_kingdom = kingdom, s_phylum = phylum, s_class = class)
turtle_keys <- dbGetQuery(con,
"SELECT DISTINCT ms_merge_key FROM taxon_model WHERE mdl_key LIKE 'rng_turtle_swot_dps|%'")$ms_merge_key
tx <- dbGetQuery(con, "
SELECT t.ms_merge_key, t.taxon_authority AS auth, t.taxon_id,
w.kingdom, w.phylum, w.\"class\" AS cls
FROM taxon t
LEFT JOIN spp.worms w ON t.taxon_authority='worms' AND CAST(w.scientificNameID AS VARCHAR)=t.taxon_id") |>
left_join(supp, by = "taxon_id") |>
mutate(kingdom = coalesce(kingdom, s_kingdom),
phylum = coalesce(phylum, s_phylum),
cls = coalesce(cls, s_class))
# the rule itself is msens::sp_cat_from_taxonomy() (unit-tested, one row per branch) -- the
# single source of truth, also used by ingest_aquax.qmd to report species by component
# before any merge exists
tx <- tx |> mutate(sp_cat = msens::sp_cat_from_taxonomy(
cls, kingdom, is_botw = auth == "botw", is_turtle = ms_merge_key %in% turtle_keys))
duckdb_register(con, "tmp_spcat", tx |> select(ms_merge_key, sp_cat))
dbExecute(con, "CREATE OR REPLACE TABLE taxon AS
SELECT t.*, s.sp_cat FROM taxon t LEFT JOIN tmp_spcat s USING (ms_merge_key)")
duckdb_unregister(con, "tmp_spcat")
knitr::kable(dbGetQuery(con, "SELECT sp_cat, count(*) n, count(*) FILTER(WHERE is_valid_usa) valid_usa
FROM taxon GROUP BY 1 ORDER BY n DESC"))
```
## Marine relevance (birds) — `is_marine`
A bird's whole global range is largely terrestrial (v8 keeps the entire range, no land
mask), so only genuinely **marine/coastal** birds should drive the marine-sensitivity
score. `is_marine` (bird) = member of a marine/coastal **family**
(`data/marine_bird_families.csv`) **and** whole-range `pct_marine ≥ 5%`, **or** on the
curated **include** list — minus the curated **exclude** list
(`data/marine_birds_curation.csv`); both are GitHub-tracked + reviewable. The share is
computed from the **raw global** `bl` ranges (the merged `model_cell` is US-scoped, so its
`pct_marine` is ~100% and unusable here). Non-bird taxa are marine by dataset → `TRUE`.
```{r}
#| label: marine_relevance
marine_fam <- read_csv(here("data/marine_bird_families.csv"), show_col_types = FALSE)$family
cur <- read_csv(here("data/marine_birds_curation.csv"), show_col_types = FALSE)
inc_sci <- clean_sci_name(cur$scientific_name[cur$action == "include"])
exc_sci <- clean_sci_name(cur$scientific_name[cur$action == "exclude"])
# whole-range ocean share per bird from the RAW global bl ranges ∩ ocean cell table
bl_pm <- dbGetQuery(con, glue("
WITH b AS (SELECT mc.mdl_key, (c.cell_id IS NOT NULL) AS is_ocean
FROM read_parquet('{dir_dist}/dataset=bl/*.parquet') mc
LEFT JOIN grid.cell c USING (cell_id))
SELECT mdl_key, round(100.0 * count(*) FILTER (WHERE is_ocean) / count(*), 1) AS pct_marine_bl
FROM b GROUP BY mdl_key"))
# family + cleaned name per bl model; map to the merged taxon via taxon_model
bl_mod <- read_csv(glue("{dir_dist}/model_bl.csv"), show_col_types = FALSE) |>
transmute(mdl_key, family, sci = clean_sci_name(sci_name))
tmod <- dbGetQuery(con, "SELECT ms_merge_key, mdl_key FROM taxon_model WHERE mdl_key LIKE 'bl|%'")
bird <- bl_pm |>
left_join(bl_mod, by = "mdl_key") |>
left_join(tmod, by = "mdl_key") |>
mutate(is_marine_bird =
((family %in% marine_fam & pct_marine_bl >= 5) | sci %in% inc_sci) & !sci %in% exc_sci) |>
group_by(ms_merge_key) |>
summarise(pct_marine_bl = max(pct_marine_bl),
is_marine_bird = any(is_marine_bird), .groups = "drop")
duckdb_register(con, "tmp_bird_marine", bird)
dbExecute(con, "CREATE OR REPLACE TABLE taxon AS
SELECT t.*,
CASE WHEN t.taxon_authority = 'botw'
THEN coalesce(b.is_marine_bird, FALSE) ELSE TRUE END AS is_marine,
b.pct_marine_bl
FROM taxon t LEFT JOIN tmp_bird_marine b USING (ms_merge_key)")
duckdb_unregister(con, "tmp_bird_marine")
mr <- dbGetQuery(con, "SELECT
count(*) FILTER (WHERE taxon_authority='botw') n_birds,
count(*) FILTER (WHERE taxon_authority='botw' AND is_marine) n_marine,
count(*) FILTER (WHERE taxon_authority='botw' AND is_valid_usa AND is_marine) n_marine_usa
FROM taxon")
log_info("marine birds: {mr$n_marine} of {mr$n_birds} ({mr$n_marine_usa} valid_usa)")
knitr::kable(mr)
# content fingerprint of the taxon output table (before disconnect)
taxon_hash <- msens::hash_query(con, "taxon")
dbDisconnect(con, shutdown = TRUE)
```
## Taxon extinction-risk
The governing `extrisk_code` / `er_score` (US-national/floored, via
`msens::compute_er_score()`) are computed upstream in `merge_models_prep` — so
`merge_models` can apply them to the range-cell values — and carried through here on the
`taxon` table (via `x.*` in the flags chunk above). No recompute needed.
## Outputs + manifest
```{r}
#| label: manifest
msens::report_table(smry, caption = "merge_taxon: validity flags")
# content-addressed: fingerprint of the taxon output table (no wall-clock)
msens::write_manifest(
manifest, target = "merge_taxon", content_hash = taxon_hash,
stats = list(ver = ver, n_taxa = smry$n_taxa,
valid_global = smry$valid_global, valid_usa = smry$valid_usa,
valid_pra = smry$valid_pra),
force = msens::force_target("merge_taxon"))
```