---
title: "Merge models — combine each taxon's models per cell (max-merge + range mask; turtles ×)"
msens:
target_name: merge_models
workflow_type: merge
dependency: [merge_models_prep]
output: data/manifests/merge_models.json
editor_options:
chunk_output_type: console
---
Combine every taxon's per-dataset models into one **merged `model_cell`**
(`mdl_key = ms_merge|WORMS:…` / `BOTW:…`), reproducing v7's logic at the v8 global
scale (~13B am cell-rows), so it must stream + spill rather than materialize:
- **Default = max-merge, masked to the range footprint — on BOTH surfaces.** A taxon's
merged cells are its **range** cells if it has any (`rng_*`/`ch_*`/`ca_*`/`bl`), else its
full `am` footprint; at each kept cell `val = max` across the taxon's models. Streamed
straight to Parquet via a single `COPY` (one grouped scan), so no 13B-row table. The
global surface applies the identical mask — it is the whole *range*, not the whole am
footprint (see @sec-masking).
- **Turtles = multiplicative.** `val = pmax(1, round(er × suit / 100))` with
critical-habitat max-override — computed from **only the 6 turtle taxa's files**
(never a full am scan, which is what OOM'd the naive version).
- **v9: AquaX (`ax`) supersedes AquaMaps (`am`) per taxon, inside the AquaX mask only.** Before
the merge input is partitioned, `msens::supersede_sql()` drops a superseded taxon's `am` cells
that fall inside `ax_mask` (AquaX's own extent: the union over models of their modeled pixels,
measured against `in_usa` in `ingest_aquax.qmd`). Outside the mask — the rest of the world and any
`in_usa` cell no AquaX model reaches — `am` carries on. Which taxa supersede is `data/ax_supersedes_am.csv`
(`supersedes` column, itself driven by `AX_ABSENT_SUPERSEDES`); `AX_SUPERSEDE=0` is the control
run — AquaX left OUT of the merge (`merge_models_prep`), so the inputs are ver_prev's and this
notebook's content hash must equal ver_prev's checkpoint (v8: `bd8b2931f17bb2b3`).
Both suitability datasets then feed `merge_sql(suit_ds = c("am", "ax"))` unchanged.
## Design
```{mermaid}
%%| label: fig-design
%%| fig-cap: "ONE batched merge → TWO surfaces from the same data, BOTH masked to the range footprint. GLOBAL viz = the whole range, max(er, am-at-range) → COGs. US scoring = v7-faithful US-boundary-aware rule (range∩US max(er,am-at-range) UNION am∩US for taxa with no US range) — NOT a trim of global. Turtles multiplicative."
flowchart LR
raw["dist/*/*.parquet<br/>(raw model_cell)"] --> sup["v9: supersede_sql()<br/>drop am ∩ ax_mask<br/>for AquaX taxa"] --> mcp[("mc_parts<br/>partitioned by mkey_id (int)")]
mcp --> gm["GLOBAL: whole range,<br/>max(er, am-at-range)<br/>(am beyond range MASKED)"]
mcp --> um["US: range∩US max(er,am@range)<br/>∪ am∩US for taxa w/ no US range"]
mcp --> tt["turtles ×<br/>greatest(1, er·suit/100)"]
mcp --> dd["NMFS DPS species<br/>val = suitability, er per cell<br/>(dist_merged_er)"]
gm --> g[("dist_merged_global/**<br/>whole-range → viz COGs")]
tt --> g
dd --> g
um --> us[("dist_merged/**<br/>US model_cell → scoring")]
tt --> us
dd --> us
us --> mf["hash_parquet → 392cb9…<br/>content-addressed manifest"]
```
## Setup
```{r}
#| label: setup
librarian::shelf(DBI, dplyr, duckdb, fs, glue, here, jsonlite, logger, quiet = T)
source(here("libs/paths.R"))
source(here("libs/vars.R"))
dir_atlas <- glue("{dir_big_v}/marine-atlas")
dir_dist <- glue("{dir_atlas}/dist")
merge_db <- glue("{dir_atlas}/merge.duckdb")
# ONE global merge → two surfaces from the same data:
dir_global <- glue("{dir_atlas}/dist_merged_global/dataset=ms_merge") # whole-range (land+ocean) → viz COGs
dir_merged <- glue("{dir_atlas}/dist_merged/dataset=ms_merge") # US-study-area trim → model_cell scoring
manifest <- here("data/manifests/merge_models.json")
dir_create(path_dir(manifest))
stopifnot("run merge_models_prep first" = file_exists(merge_db),
"run build_cell_grid first" = file_exists(sdm_db))
suit_ds <- c("am", "ax") # v9: AquaX joins AquaMaps as a suitability dataset (msens::merge_sql)
turtle_ds <- "rng_turtle_swot_dps"
# SPATIAL extinction risk, two rules (both msens, both unit-tested in test-merge.R):
# * sea turtles (v4b, msens::turtle_sql): val = greatest(1, round(er x suit / 100)), critical
# habitat overriding with a max -- ER is BAKED INTO the merged value.
# * NMFS DPS species (v9, msens::dps_sql; dps_nmfs = humpback, Southern Resident killer whale,
# salmonid ESUs, ...): the merged value STAYS the distribution (suitability masked to the ER
# footprint) and the per-cell ER is written beside it (dist_merged_er) for scoring to multiply
# in. The plain rule's species-wide `max(er, suit)` painted the humpback a flat 100; the turtle
# rule then painted it 1 everywhere its ER is 1 -- the weight, not the whale.
dps_ds <- "dps_nmfs"
ch_ds <- c("ch_nmfs", "ch_fws", "ca_nmfs")
ch_sql <- paste(sprintf("'%s'", ch_ds), collapse = ",")
con <- dbConnect(duckdb(merge_db))
# spill config so the multi-billion-row aggregation goes out-of-core instead of OOM. Keep the
# memory limit WELL under physical RAM (other processes — browser, renders — compete), cap threads
# (each parallel aggregate/join holds hash tables), and allow a large temp spill.
mem_gb <- tryCatch(
max(4, min(8, floor(as.numeric(system("sysctl -n hw.memsize", intern = TRUE)) / 1e9 * 0.3))),
error = function(e) 6)
mem_gb <- as.numeric(Sys.getenv("MERGE_MEM_GB", mem_gb)) # override when nothing else is running
dbExecute(con, glue("PRAGMA memory_limit='{mem_gb}GB'"))
# streaming COPY/aggregates need no insertion order; keeping it pins buffers the OOM guard then hits
dbExecute(con, "SET preserve_insertion_order = false")
dbExecute(con, glue("PRAGMA temp_directory='{dir_atlas}/duckdb_tmp'"))
dbExecute(con, glue("PRAGMA max_temp_directory_size='400GB'"))
dbExecute(con, glue("PRAGMA threads={min(4L, max(1L, parallel::detectCores() - 1L))}"))
# ONE GLOBAL merge, two outputs from the same data: the whole-range merged surface
# (dist_merged_global → viz COGs) and, trimmed to the US study area, the model_cell scoring surface
# (dist_merged). am is DENSE (~13B rows, 50G) so it is NEVER grouped: range cells (non-am) are
# materialized once + indexed (mc_range, ~3B); am is a VIEW (mc_am) probed only at the sparse range
# cells (both/turtle) or streamed for am-only taxa. us_cells is applied AFTER the global merge.
dbExecute(con, glue("ATTACH '{sdm_db}' AS grid (READ_ONLY)"))
dbExecute(con, "CREATE OR REPLACE TABLE us_cells AS SELECT cell_id FROM grid.cell WHERE in_usa")
n_us <- dbGetQuery(con, "SELECT count(*) n FROM us_cells")$n
# PARTITION the full merge input (ranges + am, ~16B rows) by ms_merge_key — a STREAMING write that
# never holds it in memory. The merge below then reads one BATCH of taxa at a time via
# `ms_merge_key IN (...)`, so DuckDB partition-prunes to just those taxa's slices: bounded memory,
# no OOM (a bulk join over the ~3B `both` range cells otherwise blows the memory limit).
try(dbExecute(con, "DROP TABLE IF EXISTS model_cell_all"), silent = TRUE)
try(dbExecute(con, "DROP TABLE IF EXISTS mc_range"), silent = TRUE)
# ms_merge_key holds '|' and ':' which break hive partition PATHS (…/mc_parts/ms_merge_key=ms_merge|…
# can't be parsed back). Partition on an INTEGER surrogate mkey_id — exactly what the serve surface
# does with mdl_id. ms_merge_key stays a stored column.
dbExecute(con, "CREATE OR REPLACE TABLE mkey_map AS
SELECT ms_merge_key, CAST(dense_rank() OVER (ORDER BY ms_merge_key) AS INTEGER) AS mkey_id
FROM (SELECT DISTINCT ms_merge_key FROM taxon_model)")
# v9 SUPERSESSION inputs (see msens::supersede_sql): the AquaX mask + the taxa whose AquaMaps
# cells are dropped inside it. Both come from ingest_aquax.qmd; `supersedes` already encodes the
# AX_ABSENT_SUPERSEDES policy. AX_SUPERSEDE=0 empties the table -> the control run.
ax_mask_pq <- glue("{dir_dist}/ax_mask.parquet")
sup_csv <- here("data/ax_supersedes_am.csv")
has_ax <- dir_exists(glue("{dir_dist}/dataset=ax")) && file_exists(ax_mask_pq) && file_exists(sup_csv)
if (has_ax) {
dbExecute(con, glue("CREATE OR REPLACE TABLE ax_mask AS SELECT cell_id FROM read_parquet('{ax_mask_pq}')"))
sup <- readr::read_csv(sup_csv, show_col_types = FALSE)
sup_keys <- if (ax_supersede) unique(sup$ms_merge_key[sup$supersedes]) else character(0)
dbWriteTable(con, "supersede", data.frame(ms_merge_key = sup_keys), overwrite = TRUE)
# every superseded taxon must be a taxon of THIS merge, or the registry and the crosswalk disagree
n_unknown <- dbGetQuery(con, "SELECT count(*) n FROM supersede WHERE ms_merge_key NOT IN (SELECT ms_merge_key FROM taxon_model)")$n
stopifnot("ax_supersedes_am.csv names taxa absent from taxon_model" = n_unknown == 0)
log_info("supersession: {length(sup_keys)} taxa lose am inside ax_mask ({dbGetQuery(con, 'SELECT count(*) n FROM ax_mask')$n} cells)",
if (!ax_supersede) " -- AX_SUPERSEDE=0: CONTROL RUN, nothing superseded" else "")
} else {
dbExecute(con, "CREATE OR REPLACE TABLE ax_mask (cell_id INTEGER)")
dbExecute(con, "CREATE OR REPLACE TABLE supersede (ms_merge_key VARCHAR)")
log_info("no ax ingest found — no supersession (v8-style merge)")
}
dir_parts <- glue("{dir_atlas}/mc_parts")
# the partition write is the ~40-min step; skip it if the input already exists (REDO_MC_PARTS=1 forces)
if (!dir_exists(dir_parts) || nzchar(Sys.getenv("REDO_MC_PARTS"))) {
if (dir_exists(dir_parts)) dir_delete(dir_parts)
dbExecute(con, "SET partitioned_write_max_open_files=512")
# the supersession filter is applied HERE, on the partitioned merge input, so every consumer
# downstream (both surfaces, turtles, the masking check) sees the same input
msens::copy_atlas_parquet(con, glue(
"SELECT mkey_id, ms_merge_key, ds_key, cell_id, val FROM (
SELECT mm.mkey_id, tm.ms_merge_key, split_part(mc.mdl_key,'|',1) AS ds_key, mc.cell_id, mc.val
FROM read_parquet('{dir_dist}/*/*.parquet') mc JOIN taxon_model tm USING (mdl_key)
JOIN mkey_map mm USING (ms_merge_key)) mc
WHERE {msens::supersede_sql(src = 'mc')}"),
dir_parts, partition_by = "mkey_id")
} else if (nzchar(Sys.getenv("REDO_MC_PARTS_SPATIAL"))) {
# refresh ONLY the spatial-ER taxa's partitions (a re-ingested dps_nmfs / turtle surface): mkey_id
# is dense_rank over the taxon set, so refreshing in place is valid only while the set -- and
# every position -- is unchanged; assert it against the partition count before touching anything
n_map <- dbGetQuery(con, "SELECT count(*) n FROM mkey_map")$n
stopifnot("mc_parts partition count != taxa in mkey_map: the map moved, REDO_MC_PARTS=1" =
length(dir_ls(dir_parts)) == n_map)
sp_ids <- dbGetQuery(con, "SELECT DISTINCT mm.mkey_id FROM taxon_model tm JOIN mkey_map mm USING (ms_merge_key)
WHERE tm.mdl_key LIKE 'rng_turtle_swot_dps|%' OR tm.mdl_key LIKE 'dps_nmfs|%'")$mkey_id
for (id in sp_ids) { pd <- glue("{dir_parts}/mkey_id={id}"); if (dir_exists(pd)) dir_delete(pd) }
dbExecute(con, "SET partitioned_write_max_open_files=512")
msens::copy_atlas_parquet(con, glue(
"SELECT mkey_id, ms_merge_key, ds_key, cell_id, val FROM (
SELECT mm.mkey_id, tm.ms_merge_key, split_part(mc.mdl_key,'|',1) AS ds_key, mc.cell_id, mc.val
FROM read_parquet('{dir_dist}/*/*.parquet') mc JOIN taxon_model tm USING (mdl_key)
JOIN mkey_map mm USING (ms_merge_key)
WHERE mm.mkey_id IN ({paste(sp_ids, collapse = ',')})) mc
WHERE {msens::supersede_sql(src = 'mc')}"),
dir_parts, partition_by = "mkey_id")
log_info("REDO_MC_PARTS_SPATIAL: re-partitioned {length(sp_ids)} spatial-ER taxa; the other partitions are kept")
} else log_info("mc_parts exists — skipping repartition (REDO_MC_PARTS=1 to force, REDO_MC_PARTS_SPATIAL=1 for the spatial-ER taxa only)")
# VIEW over the partitioned input — used only for the per-partition COUNT below (metadata scan).
# Every DATA read goes through parts_sql(): the partitions named by PATH, never the 150k-file glob
# filtered by `mkey_id IN (...)` -- the reader buffered far beyond the 6 turtle partitions it needed
# and hit the memory limit (v9 control run, 2026-08-27).
dbExecute(con, glue("CREATE OR REPLACE VIEW mc_parts AS
SELECT mkey_id, ms_merge_key, ds_key, cell_id, val
FROM read_parquet('{dir_parts}/**/*.parquet', hive_partitioning = true)"))
parts_sql <- function(ids) glue("read_parquet([{paste0(\"'\", dir_parts, \"/mkey_id=\", ids, \"/*.parquet'\", collapse = \", \")}], hive_partitioning = true)")
all_ids <- dbGetQuery(con, "SELECT mkey_id FROM mkey_map ORDER BY mkey_id")$mkey_id
log_info("partitioned merge input: {length(all_ids)} taxa by mkey_id; US trim {n_us} in_usa cells; memory_limit={mem_gb}GB")
```
## Spatial extinction risk — sea turtles (×) and NMFS DPS species (suitability + ER per cell)
```{r}
#| label: turtles
turtle_keys <- dbGetQuery(con, "SELECT DISTINCT ms_merge_key FROM taxon_model WHERE mdl_key LIKE 'rng_turtle_swot_dps|%'")$ms_merge_key
dps_keys <- dbGetQuery(con, glue("SELECT DISTINCT ms_merge_key FROM taxon_model WHERE mdl_key LIKE '{dps_ds}|%'"))$ms_merge_key
stopifnot("a taxon cannot carry both a turtle and a dps_nmfs ER surface" = !any(dps_keys %in% turtle_keys))
ids_of <- function(keys) if (length(keys))
dbGetQuery(con, glue("SELECT mkey_id FROM mkey_map WHERE ms_merge_key IN ({paste(sprintf(\"'%s'\", keys), collapse = ',')})"))$mkey_id else integer(0)
turtle_ids <- ids_of(turtle_keys); dps_ids <- ids_of(dps_keys)
log_info("spatial-ER taxa: {length(turtle_keys)} turtles (multiplicative) + {length(dps_keys)} NMFS-DPS species (suitability + per-cell ER)")
# turtles: whole-range multiplicative merge, reading ONLY the turtle taxa's partitions
# (mkey_id IN (...) → partition-pruned, tiny). Rule = msens::turtle_sql() (unit-tested).
dbExecute(con, glue("CREATE OR REPLACE TABLE turtle_src AS
SELECT ms_merge_key, ds_key, cell_id, val FROM {parts_sql(turtle_ids)}"))
dbExecute(con, paste("CREATE OR REPLACE TABLE mc_turtle AS",
msens::turtle_sql(turtle_ds, suit_ds, ch_ds, src = "turtle_src")))
log_info("mc_turtle (whole-range): {dbGetQuery(con, 'SELECT count(*) n FROM mc_turtle')$n} cells")
# NMFS DPS species: msens::dps_sql() -- (mdl_key, cell_id, val = suitability, er = per-cell ER).
# The critical-habitat datasets are NOT an input: dps_nmfs was built from the same designations
# WITH each entity's status, while ch_* carries the species-level ER (100 for the humpback) and
# would paint the Threatened Mexico-DPS habitat back to Endangered.
# one taxon at a time over a VIEW of its partition: materializing all 19 as one table (the humpback
# alone is 15 M range cells against global AquaMaps) blew the DuckDB memory budget; per taxon each
# join is small and mc_dps just accumulates
dbExecute(con, "CREATE OR REPLACE TABLE mc_dps (mdl_key VARCHAR, cell_id INTEGER, val DOUBLE, er DOUBLE)")
for (id in dps_ids) {
dbExecute(con, glue("CREATE OR REPLACE VIEW dps_src AS
SELECT ms_merge_key, ds_key, cell_id, val FROM {parts_sql(id)} WHERE ds_key NOT IN ({ch_sql})"))
dbExecute(con, paste("INSERT INTO mc_dps", msens::dps_sql(dps_ds, suit_ds, src = "dps_src")))
}
if (length(dps_ids)) dbExecute(con, "DROP VIEW IF EXISTS dps_src")
log_info("mc_dps (whole-range): {dbGetQuery(con, 'SELECT count(*) n FROM mc_dps')$n} cells over {length(dps_ids)} taxa")
```
## Max-merge (masked) + spatial-ER taxa → streamed to Parquet
```{r}
#| label: merge_write
# taxon flags from the crosswalk (no cell scan)
suit_sql <- paste(sprintf("'%s'", suit_ds), collapse = ",")
dbExecute(con, glue("CREATE OR REPLACE TABLE taxon_flags AS
SELECT ms_merge_key,
bool_or(split_part(mdl_key,'|',1) IN ({suit_sql})) AS has_suit,
bool_or(split_part(mdl_key,'|',1) NOT IN ({suit_sql})) AS has_range
FROM taxon_model GROUP BY ms_merge_key"))
# skip the ~2h merge write when the surfaces already exist — re-renders (e.g. to refresh the HTML)
# reuse them; the summaries + hash below read the on-disk parquet either way. REDO_MERGE=1 forces
# everything; REDO_MERGE_SPATIAL=1 rewrites ONLY the turtle + dps outputs (a rule or ER-surface
# change for the ~25 spatial-ER taxa) and keeps the batched surfaces of the other ~17,750 taxa.
redo_merge <- nzchar(Sys.getenv("REDO_MERGE"))
redo_spatial <- nzchar(Sys.getenv("REDO_MERGE_SPATIAL"))
# per-cell ER of the dps taxa, US-trimmed, beside the scoring surface: score_zones materializes it as
# model_cell_er and score_cell_metrics multiplies it in (er_mode = 'cell')
dir_merged_er <- glue("{dir_atlas}/dist_merged_er/dataset=ms_merge")
have_merged <- dir_exists(dir_merged) && dir_exists(dir_global) &&
length(dir_ls(dir_merged, glob = "*.parquet", recurse = TRUE)) > 0
write_spatial <- function() {
for (d in c(glue("{dir_global}/turtle"), glue("{dir_merged}/turtle"),
glue("{dir_global}/dps"), glue("{dir_merged}/dps"), glue("{dir_merged_er}/dps")))
if (dir_exists(d)) dir_delete(d)
dir_create(dir_merged_er)
# turtles → whole-range global + US-trim (ER baked in)
msens::copy_atlas_parquet(con, "SELECT mdl_key, cell_id, val FROM mc_turtle",
glue("{dir_global}/turtle"), per_thread = TRUE)
msens::copy_atlas_parquet(con,
"SELECT mt.mdl_key, mt.cell_id, mt.val FROM mc_turtle mt JOIN us_cells u ON mt.cell_id = u.cell_id",
glue("{dir_merged}/turtle"), per_thread = TRUE)
# dps → whole-range global + US-trim (the distribution) + US-trim ER beside it
msens::copy_atlas_parquet(con, "SELECT mdl_key, cell_id, val FROM mc_dps",
glue("{dir_global}/dps"), per_thread = TRUE)
msens::copy_atlas_parquet(con,
"SELECT d.mdl_key, d.cell_id, d.val FROM mc_dps d JOIN us_cells u ON d.cell_id = u.cell_id",
glue("{dir_merged}/dps"), per_thread = TRUE)
msens::copy_atlas_parquet(con,
"SELECT d.mdl_key, d.cell_id, d.er FROM mc_dps d JOIN us_cells u ON d.cell_id = u.cell_id",
glue("{dir_merged_er}/dps"), per_thread = TRUE)
}
if (redo_merge || !have_merged) {
if (dir_exists(dir_global)) dir_delete(dir_global); dir_create(dir_global)
if (dir_exists(dir_merged)) dir_delete(dir_merged); dir_create(dir_merged)
if (dir_exists(dir_merged_er)) dir_delete(dir_merged_er); dir_create(dir_merged_er)
# survives the run so the masking check can re-run on a skipped re-render (see @sec-masking)
dbExecute(con, "CREATE OR REPLACE TABLE range_counts (ms_merge_key VARCHAR, n_range BIGINT)")
write_spatial()
# BATCH the rest so each merge reads ONLY its taxa's pruned partitions (bounded memory). Two surfaces
# come off the same batch (see glob_global_sql / us_sql below): the GLOBAL viz surface = am ∪ range
# (whole-range merged model, FULL OUTER); the US SCORING surface = v7-faithful US-boundary-aware rule
# (range footprint ∩ US valued max(er, am-at-range), UNION am ∩ US for taxa with no range in the US).
# The US rule keeps the 750 am-in-US taxa whose IUCN range lies wholly outside the US (silently
# dropped before) yet still masks am to the range for taxa that DO have US range — so it reproduces
# the prior US merge exactly and is NOT a plain trim of the global surface.
# am-only taxa have no range footprint → omitted from the GLOBAL surface (reuse am COGs), but their
# US am cells go into the scoring surface via the am-only UNION branch.
# Batch by CELL COUNT, not taxa count: a few taxa have global ranges (10-25M cells each), so a
# fixed 150-taxa batch OOMs. Pack taxa until ~40M cells/batch (a single taxon maxes at the 25.9M
# grid, so it always fits one batch). Counts come from parquet metadata (fast, no data scan).
nonturtle <- setdiff(all_ids, c(turtle_ids, dps_ids))
cnt <- dbGetQuery(con, "SELECT mkey_id, count(*) AS n FROM mc_parts GROUP BY mkey_id")
cnt <- cnt[cnt$mkey_id %in% nonturtle, ]; cnt <- cnt[order(cnt$mkey_id), ]
batches <- unname(split(cnt$mkey_id, cumsum(as.numeric(cnt$n)) %/% 40e6))
log_info("merging {length(nonturtle)} non-turtle taxa in {length(batches)} cell-bounded batches (~40M cells each)")
# The merge RULES are msens::merge_sql() — the SINGLE SOURCE OF TRUTH, unit-tested in
# msens/tests/testthat/test-merge.R (every taxon category asserted). Two surfaces from the same batch:
# GLOBAL viz ($global) = the WHOLE RANGE footprint, max(er, am-at-range) [am BEYOND the range is
# MASKED here too] → COGs. Whole *range*, not whole am footprint: this surface was briefly a FULL
# OUTER union with the entire am footprint, which painted raw AquaMaps over-prediction into every
# merged COG (MarineSensitivity/apps#8) while scoring stayed correct.
# US scoring ($us) = v7-faithful, IUCN-CONSTRAINED: range∩US max(er, am-at-range) [am BEYOND the
# range is MASKED] UNION raw am∩US for TRUE am-only taxa (global has_range=FALSE). A species whose
# IUCN range lies WHOLLY outside the US gets NO US presence — the iucn_range_outside_us_eez
# exclusion (e.g. Sotalia guianensis). See ?msens::merge_sql for the full rationale.
msq <- msens::merge_sql(suit_ds)
glob_global_sql <- msq$global
us_sql <- msq$us
for (i in seq_along(batches)) {
dbExecute(con, glue("CREATE OR REPLACE TABLE b AS
SELECT ms_merge_key, ds_key, cell_id, val FROM {parts_sql(batches[[i]])}"))
dbExecute(con, msq$b_range) # range footprint (non-am) valued by governing er_score
dbExecute(con, msq$b_am_rng) # am AT range cells only → range-footprint max(er, am-at-range)
# per-taxon range-cell tally, accumulated while b_range is in hand (free) — the masking check
# below needs it and would otherwise have to re-scan 12 GB of range parquet on every render
dbExecute(con, "INSERT INTO range_counts SELECT ms_merge_key, count(*) FROM b_range GROUP BY 1")
msens::copy_atlas_parquet(con, glob_global_sql, glue("{dir_global}/b{i}"), per_thread = TRUE)
msens::copy_atlas_parquet(con, us_sql, glue("{dir_merged}/b{i}"), per_thread = TRUE)
if (i %% 10 == 0) log_info(" batch {i}/{length(batches)}")
}
dbExecute(con, "DROP TABLE IF EXISTS b; DROP TABLE IF EXISTS b_range;
DROP TABLE IF EXISTS b_am_rng")
} else if (redo_spatial) {
write_spatial()
log_info("REDO_MERGE_SPATIAL: rewrote the turtle + dps surfaces; the batched surfaces of the other taxa are kept")
} else log_info("dist_merged + dist_merged_global exist — skipping merge write (REDO_MERGE=1 to force, REDO_MERGE_SPATIAL=1 for the spatial-ER taxa only)")
# summaries: global (viz) + US (scoring); fingerprint the US scoring surface (the tracked output)
smry_g <- dbGetQuery(con, glue("SELECT count(DISTINCT mdl_key) n_taxa, count(*) n_cells,
min(val) v_min, max(val) v_max FROM read_parquet('{dir_global}/**/*.parquet')"))
smry <- dbGetQuery(con, glue("SELECT count(DISTINCT mdl_key) n_taxa, count(*) n_cells,
min(val) v_min, max(val) v_max FROM read_parquet('{dir_merged}/**/*.parquet')"))
h <- msens::hash_parquet(glue("{dir_merged}/**/*.parquet"), con)
# the per-cell ER beside the scoring surface (dps taxa) — scored, so fingerprinted too (er_hash);
# score_zones resumes on content_hash + er_hash, since a changed ER leaves the val surface intact
have_er <- dir_exists(dir_merged_er) && length(dir_ls(dir_merged_er, glob = "*.parquet", recurse = TRUE)) > 0
smry_er <- if (have_er) dbGetQuery(con, glue("SELECT count(DISTINCT mdl_key) n_taxa, count(*) n_cells,
min(er) er_min, max(er) er_max FROM read_parquet('{dir_merged_er}/**/*.parquet')")) else
data.frame(n_taxa = 0L, n_cells = 0L, er_min = NA_real_, er_max = NA_real_)
h_er <- if (have_er) msens::hash_parquet(glue("{dir_merged_er}/**/*.parquet"), con) else "none"
msens::report_table(smry_g, caption = "GLOBAL whole-range merged surface (ms_merge) → viz COGs")
msens::report_table(smry, caption = "US-trimmed merged model_cell (ms_merge) → scoring")
msens::report_table(smry_er, caption = "US-trimmed per-cell extinction risk beside it (dps taxa; dist_merged_er) → scoring weight")
```
## Masking check {#sec-masking}
The global surface is what the species app **draws**, and it is the one surface the manifest hash
cannot see — the fingerprint is taken over `dist_merged` alone. So assert the mask here, on the
written parquet, or nothing does: for every non-turtle taxon the global surface must have **exactly
one cell per range cell**, because `$global` is `b_range LEFT JOIN b_am_rng` and `b_range` is
`DISTINCT`. Any surplus is AquaMaps painted beyond the expert range — the apps#8 regression, which
had 3,461 taxa carrying it (the walrus alone 743,749 surplus cells, reaching 9.75°N).
Turtles and the NMFS DPS species are excluded: their footprint is the extinction-risk surface, not
the union of every non-am dataset, so the identity doesn't apply — `msens::turtle_sql()` and
`msens::dps_sql()` mask them by construction and `test-merge.R` asserts it.
```{r}
#| label: masking_check
if (dbExistsTable(con, "range_counts")) {
bad <- dbGetQuery(con, glue("
WITH g AS (SELECT mdl_key, count(*) AS n_global
FROM read_parquet('{dir_global}/**/*.parquet') GROUP BY 1)
SELECT g.mdl_key, r.n_range, g.n_global, g.n_global - r.n_range AS n_surplus
FROM g JOIN range_counts r ON r.ms_merge_key = g.mdl_key
WHERE g.n_global <> r.n_range ORDER BY n_surplus DESC"))
n_ck <- dbGetQuery(con, "SELECT count(*) n FROM range_counts")$n
if (nrow(bad)) msens::report_table(head(bad, 20), caption = "UNMASKED taxa (am beyond the range)")
log_info("masking check: {n_ck - nrow(bad)}/{n_ck} taxa masked to their range footprint")
stopifnot("global surface has cells beyond the range footprint (apps#8)" = nrow(bad) == 0)
} else log_info("range_counts absent (merge skipped by a prior run) — REDO_MERGE=1 to re-check")
# v9 supersession check, on the partitioned INPUT (what both surfaces read): for every superseded
# taxon, no am cell may survive inside ax_mask. The v8 behaviour fails this by construction.
n_sup <- dbGetQuery(con, "SELECT count(*) n FROM supersede")$n
if (n_sup > 0) {
sup_ids <- dbGetQuery(con, "SELECT mkey_id FROM mkey_map WHERE ms_merge_key IN (SELECT ms_merge_key FROM supersede)")$mkey_id
leak <- dbGetQuery(con, glue("
SELECT count(*) AS n_rows, count(DISTINCT ms_merge_key) AS n_taxa
FROM {parts_sql(sup_ids)} p
WHERE p.ds_key = 'am' AND p.cell_id IN (SELECT cell_id FROM ax_mask)"))
log_info("supersession check: {leak$n_rows} am rows inside ax_mask for {leak$n_taxa} superseded taxa (expect 0)")
stopifnot("AquaMaps cells survive inside the AquaX mask for superseded taxa (supersede_sql not applied to mc_parts?)" =
leak$n_rows == 0)
# ...and the converse: those taxa DO carry ax inside the mask (the ingest reached the merge)
n_ax <- dbGetQuery(con, "SELECT count(DISTINCT ms_merge_key) n FROM mc_parts WHERE ds_key = 'ax'")$n
log_info("ax feeds {n_ax} taxa in the merge input")
}
dbDisconnect(con, shutdown = TRUE)
```
## Manifest
```{r}
#| label: manifest
# content-addressed: fingerprint of the US scoring surface (the tracked target output); the global
# viz surface (dist_merged_global) is recorded in stats but doesn't drive the hash.
msens::write_manifest(
manifest, target = "merge_models", content_hash = h,
stats = list(ver = ver, n_taxa = smry$n_taxa, n_cells = smry$n_cells,
v_min = smry$v_min, v_max = smry$v_max,
n_taxa_global = smry_g$n_taxa, n_cells_global = smry_g$n_cells,
er_hash = h_er, n_taxa_er = smry_er$n_taxa, n_cells_er = smry_er$n_cells),
force = msens::force_target("merge_models"))
```