---
title: "Publish native-format input surfaces — range PMTiles + AquaMaps COGs"
msens:
target_name: publish_native
workflow_type: release
dependency: [merge_taxon]
output: data/manifests/publish_native.json
editor_options:
chunk_output_type: console
---
Phase 4b: publish each **raw input** model in a cloud-native, pyramided format so the
species app can overlay the whole (often global) range a source contributes — surfaces
far too large (up to ~19M cells) for the dense per-cell SQL tiler:
- **vector ranges** (`rng_iucn`, `rng_fws`, `ca_nmfs`, `ch_fws`, `ch_nmfs`, `bl`,
`rng_turtle_swot_dps`) → **one PMTiles per dataset**, features tagged with `mdl_key`
(client filters by it), served from the file host `/pmtiles`.
- **AquaMaps** (`am`, raster) → **one COG per model** (bbox-cropped, overviews), served
from S3 by titiler-v8 `/cog`.
- **AquaX** (`ax`, raster, v9) → COGs are built + uploaded by `ingest_aquax.qmd`; this notebook
only **registers** them from `model_ax.csv`. AquaX is delivered ON the analysis grid, so its two
representations are not original-vs-interpolated: `native` = the band **as delivered** (Float32
0–1000, every modeled pixel), `model` = the surface **as ingested** (÷10, INT1U 1–100, the ≥ 1
threshold applied — what the merge consumes). The `dataset.on_grid` flag tells the app to label
the toggle *Delivered / As ingested*.
Both drive a `native_asset` registry (`mdl_key → asset_type + url`) that STAC + the app read.
We publish inputs for every taxon the toolkit considers valid in EITHER surface
(`(is_valid_usa OR is_valid_global) & is_marine`, non reptile/amphibian). **Neither flag alone
is a superset of the other**: v8 merges two surfaces (a global one and a US-scoped scoring one),
so 4,473 taxa are global-only and **46 are US-only** — valid and scored in the US, absent from
the global surface. Building for `is_valid_usa` alone left every non-US species' inputs
unpublished; building for `is_valid_global` alone silently dropped those 46. The union is
21,581 taxa / 25,450 input edges.
**Flags:** `REDO_NATIVE=1` (rebuild PMTiles + IUCN gpkg + re-sort am), `REDO_PMTILES=1` (rebuild
just the vector tiles), `REDO_MERGED_COG=1` (repaint just the merged whole-range COGs — what a
re-merge invalidates), `NATIVE_SKIP_AM` / `NATIVE_SKIP_PMTILES`,
`NATIVE_TEST_N=<n>` (only first n am COGs, smoke test), `NATIVE_NO_S3` (stage locally, skip upload),
`NATIVE_WORKERS=<n>`, `NATIVE_VEC_COG=1` / `PUBLISH_MERGED_COG=1` (opt-in COG classes),
`NATIVE_REGISTRY_REBUILD=1` (publish ONLY what this run built — see the registry chunk; without it
a skipped class is carried forward rather than deleted).
> **`NATIVE_SKIP_AM` cannot be combined with a merged rebuild.** The am-only *aliases* — ~14,799
> `ms_merge` registry rows pointing at an existing `native/am` COG — are derived from `am_reg`,
> so skipping am empties their source and the `ms_merge` class comes out 21,584 → 6,785 rows.
> `registry_merge()` refuses that as a shrink, which is the guard working, but the run has already
> spent its time by then. Repaint merged COGs with `REDO_MERGED_COG=1` and leave am unskipped: am
> is resumable by file existence, so an existing set costs one `file_exists` sweep.
## Design
```{mermaid}
%%| label: fig-design
%%| fig-cap: "Per-model native inputs: vector ranges → PMTiles, AquaMaps cells → COGs, unified in a native_asset registry"
flowchart LR
srv["merge.duckdb<br/>served taxa"] --> vec["vector ranges"]
srv --> am["AquaMaps"]
vec --> vn["PMTiles polygons<br/>representation=native"]
vec --> vm["0.05° gridded COG (opt-in)<br/>representation=model"]
am --> an["0.5° HCAF COG<br/>representation=native"]
am --> amm["0.05° resampled COG<br/>representation=model"]
vn --> reg[("native_asset<br/>mdl_key × representation → url, bbox")]
vm --> reg
an --> reg
amm --> reg
reg --> mf["hash_query → content-addressed manifest"]
```
## Setup
```{r}
#| label: setup
librarian::shelf(DBI, dplyr, duckdb, fs, furrr, glue, here, jsonlite, logger, purrr,
readr, sf, stringr, terra, MarineSensitivity/msens, quiet = T)
source(here("libs/paths.R"))
source(here("../msens/R/publish.R")) # override any stale install
dir_atlas <- glue("{dir_big_v}/marine-atlas")
merge_db <- glue("{dir_atlas}/merge.duckdb")
dir_native <- glue("{dir_atlas}/native"); dir_create(dir_native)
dir_cog <- glue("{dir_native}/am"); dir_create(dir_cog)
dir_pmt <- glue("{dir_native}/pmtiles"); dir_create(dir_pmt)
manifest <- here("data/manifests/publish_native.json"); dir_create(path_dir(manifest))
stopifnot(all(file_exists(c(sdm_db, merge_db, cellid_tif))))
redo <- nzchar(Sys.getenv("REDO_NATIVE"))
do_s3 <- !nzchar(Sys.getenv("NATIVE_NO_S3"))
skip_am <- nzchar(Sys.getenv("NATIVE_SKIP_AM"))
skip_pmt <- nzchar(Sys.getenv("NATIVE_SKIP_PMTILES"))
# narrower than REDO_NATIVE, which also rebuilds the ~7 GB IUCN gpkg and re-sorts am
redo_pmt <- redo || nzchar(Sys.getenv("REDO_PMTILES"))
test_n <- as.integer(Sys.getenv("NATIVE_TEST_N", "0"))
n_workers <- as.integer(Sys.getenv("NATIVE_WORKERS", "6"))
s3_ver <- glue("{s3_atlas}/{ver}")
s3_http <- glue("https://s3.us-east-1.amazonaws.com/{sub('^s3://', '', s3_ver)}")
pmt_base <- glue("https://file.marinesensitivity.org/pmtiles/{ver}")
safe_key <- function(k) gsub("[^A-Za-z0-9._-]", "_", k) # mdl_key -> filesystem/URL-safe
```
## Served raw mdl_keys (the set the app lists)
```{r}
#| label: served
con <- dbConnect(duckdb(), config = list(threads = "6", memory_limit = "10GB"))
dbExecute(con, glue("SET temp_directory='{dir_atlas}/duckdb_tmp'"))
dbExecute(con, glue("ATTACH '{merge_db}' AS m (READ_ONLY)"))
# the UNION of the two validity flags. is_valid_global is NOT a superset of is_valid_usa --
# 46 taxa are US-valid but absent from the global merged surface (n_cells == n_usa), and
# scoping to is_valid_global alone dropped their merged-COG aliases, which is what
# registry_merge() refused to publish as a shrink.
served <- dbGetQuery(con, "
SELECT DISTINCT tm.mdl_key, split_part(tm.mdl_key,'|',1) AS ds_key
FROM m.taxon_model tm
JOIN m.taxon t USING(ms_merge_key)
WHERE (t.is_valid_usa OR t.is_valid_global) AND t.is_marine
AND t.sp_cat NOT IN ('reptile','amphibian')")
dbDisconnect(con, shutdown = TRUE)
served_by_ds <- split(served$mdl_key, served$ds_key)
knitr::kable(served |> count(ds_key, name = "n_mdl") |> arrange(desc(n_mdl)))
```
## Vector ranges → per-model PMTiles
Each dataset's source polygons are read, `mdl_key` derived from the documented source id,
then **inner-joined to the served set** (reproducing each ingest's presence/marine/category
filter — only ingested-and-served models survive the join) and split into **one PMTiles per
model** at `{ds}/{sp_id}.pmtiles`. A per-model file — like the per-model AquaMaps COGs — needs
no client-side filter and its low-zoom tiles stay a few KB. Packing thousands of overlapping
global ranges into one per-dataset archive overflowed the z0–z2 world tiles past maplibre's
per-tile budget, so a selected range silently failed to draw at global zoom.
IUCN's thousands of ranges are read from a **v8 indexed GeoPackage** (`native/src/rng_iucn.gpkg`,
built once from the source shapefiles with the ingest's presence/marine/category filter,
dissolved per `id_no`, indexed) — the derived per-group gpkgs are stale v7 artifacts (410 served
species short).
```{r}
#| label: pmtiles-specs
# reader per dataset: returns an sf with a `mdl_key` column (unfiltered — the served
# inner-join happens after). paths mirror the ingest notebooks (see Explore recon).
Sys.setenv(OGR_STROKE_CURVE = "TRUE") # linearize BOTW MULTISURFACE curves
slug <- function(x) gsub("^_|_$", "", gsub("[^A-Za-z0-9]+", "_", x))
read_iucn <- function() {
# read the v8 indexed ranges gpkg (one dissolved feature per served id_no), built by the
# `iucn-src` chunk below — fast + off-Drive vs. crawling thousands of source shapefiles.
gpkg <- glue("{dir_native}/src/rng_iucn.gpkg")
stopifnot("run the iucn-src chunk first" = file_exists(gpkg))
st_read(gpkg, quiet = TRUE) |> transmute(mdl_key = as.character(mdl_key)) |> st_zm(drop = TRUE)
}
read_fws_rng <- function() {
shps <- glue("{dir_raw}/fws.gov/usfws_complete_species_current_range/usfws_complete_species_current_range_{1:2}.shp")
map_dfr(shps[file_exists(shps)], function(s) {
x <- st_read(s, quiet = TRUE); names(x) <- tolower(names(x))
x |> transmute(mdl_key = glue("rng_fws|{spcode}")) # source is NAD83 -> publish_pmtiles reprojects
})
}
read_ca_nmfs <- function() {
s <- glue("{dir_raw}/fisheries.noaa.gov/core-areas/shapefile_Rices_whale_core_distribution_area_Jun19_SERO/shapefile_Rices_whale_core_distribution_area_Jun19_SERO.shp")
st_read(s, quiet = TRUE) |> transmute(mdl_key = "ca_nmfs|Balaenoptera_ricei")
}
read_ch_fws <- function() {
x <- st_read(glue("{dir_raw}/fws.gov/crithab_all_layers/crithab_poly.shp"), quiet = TRUE)
names(x) <- tolower(names(x)); x |> transmute(mdl_key = glue("ch_fws|{spcode}"))
}
read_ch_nmfs <- function() {
x <- st_read(glue("{dir_raw}/fisheries.noaa.gov/ply.gpkg"), layer = "ply", quiet = TRUE)
x |> transmute(mdl_key = glue("ch_nmfs|{slug(SCIENAME)}"))
}
read_turtle <- function() {
codes <- c("CC","CM","DC","EI","LK","LO")
map_dfr(codes, function(cd) {
s <- glue("{dir_raw}/swot_seamap.env.duke.edu/swot_distribution/Global_Distribution_{cd}.shp")
if (!file_exists(s)) return(NULL)
st_read(s, quiet = TRUE) |> transmute(mdl_key = glue("rng_turtle_swot_dps|{cd}")) |> st_zm(drop = TRUE)
})
}
read_bl <- function(sisids) {
gpkg <- c(glue("{dir_big}/../raw/BOTW_2024_2.gpkg"),
glue("{dir_raw}/birdlife.org/BOTW_GPKG_2024_2/BOTW_2024_2.gpkg"))
gpkg <- gpkg[file_exists(gpkg)][1]; stopifnot(!is.na(gpkg))
ids <- paste(unique(sub("^bl\\|", "", sisids)), collapse = ",")
q <- glue("SELECT sisid, geom FROM all_species WHERE sisid IN ({ids}) AND presence IN (1,2,3)")
st_read(gpkg, query = q, quiet = TRUE) |> transmute(mdl_key = glue("bl|{sisid}"))
}
pmt_specs <- list(
# fast/small first (seconds each), then bl (one local gpkg query), then rng_iucn
# LAST (thousands of shapefiles off Google Drive — minutes even warm).
ca_nmfs = read_ca_nmfs,
ch_fws = read_ch_fws,
ch_nmfs = read_ch_nmfs,
rng_turtle_swot_dps = read_turtle,
rng_fws = read_fws_rng,
bl = function() read_bl(served_by_ds[["bl"]]),
rng_iucn = read_iucn)
```
### IUCN v8 indexed ranges GeoPackage
Built once from the source shapefiles with the ingest's presence/marine/category filter,
dissolved to one feature per served `id_no`, indexed on `id_no`. `read_iucn()` then reads it
in one fast, off-Drive query (the derived per-group gpkgs are stale v7 and miss 410 species).
```{r}
#| label: iucn-src
iucn_gpkg <- glue("{dir_native}/src/rng_iucn.gpkg"); dir_create(path_dir(iucn_gpkg))
# REDO_IUCN_SRC exists separately from REDO_NATIVE for the same reason REDO_PMTILES does: this
# gpkg is filtered to the SERVED ids, so widening the served set (e.g. is_valid_usa ->
# is_valid_global) silently leaves the new species out of it — publish_pmtiles_from_gpkg then
# finds no features and skips them, with no error. REDO_NATIVE would also re-sort am and rebuild
# every existing tile, which is hours of work to add species that need only this one step.
redo_src <- redo || nzchar(Sys.getenv("REDO_IUCN_SRC"))
# The INGEST's own filtered output is the source of truth: presence 1:3 / marine / category are
# already applied in it, and it is indexed on id_no. Reading the raw shapefiles reproduces that
# filter from scratch and needs ~20 GB of source this machine may not hold — keep it only as the
# fallback for a box that has the shapefiles but not the derived gpkg.
iucn_all <- glue("{dir_derived}/iucnredlist.org/rng_iucn_all.gpkg")
if (!file_exists(iucn_gpkg) || redo_src) {
iucn_ok <- c("CR","EN","VU","NT","LC","DD")
ids_srv <- as.integer(sub("^rng_iucn\\|", "", served_by_ds[["rng_iucn"]]))
if (file_exists(iucn_all)) {
log_info("building {path_file(iucn_gpkg)} for {length(ids_srv)} served ids from {path_file(iucn_all)} …")
# ogr2ogr streams: the source is ~9 GB, far too large to hold as an sf. mdl_key is derived
# here so publish_pmtiles_from_gpkg can index-seek on it, exactly as before.
ids_sql <- paste(ids_srv, collapse = ",")
if (file_exists(iucn_gpkg)) file_delete(iucn_gpkg)
# system2() goes through a shell, so every argument that contains SQL punctuation
# ( ( ) | ' ) must be shQuote'd, and `~` must be expanded here rather than left for a
# shell that will not see it inside quotes.
# `geom` MUST be in the SELECT list. ogr2ogr carries the geometry only if you ask for it,
# and a GPKG without a geometry column still reads back cleanly as a plain data.frame with
# the right row count and column names — so a check on nrow()/names() passes on a file that
# has no geometry at all. It surfaces much later, as st_zm() on a data.frame.
sql_i <- glue("SELECT id_no, ('rng_iucn|' || CAST(id_no AS TEXT)) AS mdl_key, ",
"'rng_iucn' AS ds_key, geom FROM rng_iucn WHERE id_no IN ({ids_sql})")
st <- system2("ogr2ogr", c(
"-f", "GPKG", shQuote(path.expand(iucn_gpkg)), shQuote(path.expand(iucn_all)),
"-nln", "rng_iucn", "-nlt", "PROMOTE_TO_MULTI", "-sql", shQuote(sql_i)),
stdout = TRUE, stderr = TRUE)
if (!is.null(attr(st, "status")) && attr(st, "status") != 0)
stop(paste(tail(st, 20), collapse = "\n"))
system2("ogrinfo", c(shQuote(path.expand(iucn_gpkg)), "-sql",
shQuote("CREATE INDEX idx_mdlkey ON rng_iucn(mdl_key)")))
# assert it is REALLY spatial, the same way the consumer reads it (st_read + query), not
# merely that rows exist -- see the geom note above
chk <- st_read(iucn_gpkg, quiet = TRUE, as_tibble = FALSE,
query = glue("SELECT * FROM \"rng_iucn\" WHERE id_no = {ids_srv[1]}"))
stopifnot(
"iucn src gpkg has no geometry column (geom missing from the ogr2ogr SELECT?)" =
inherits(chk, "sf"),
"iucn src gpkg geometry is empty" = nrow(chk) > 0 && !all(st_is_empty(chk)))
log_info("wrote {path_file(iucn_gpkg)} — geometry verified via st_read(query=)")
} else {
log_info("building {path_file(iucn_gpkg)} for {length(ids_srv)} served ids (source shapefiles) …")
shps <- dir_ls(glue("{dir_raw}/iucnredlist.org"), glob = "*.shp", recurse = TRUE)
gi <- map_dfr(shps, function(s) {
y <- tryCatch(st_read(s, quiet = TRUE), error = function(e) NULL)
if (is.null(y)) return(NULL); names(y) <- tolower(names(y))
if (!all(c("id_no","presence","marine","category") %in% names(y))) return(NULL)
y |> filter(presence %in% 1:3, tolower(marine) == "true", category %in% iucn_ok,
id_no %in% ids_srv) |>
transmute(id_no = as.integer(id_no)) |> st_zm(drop = TRUE)
})
# keep each presence polygon as its own feature (same mdl_key) — publish_pmtiles_models
# groups a model's features into one file; dissolving would trip s2 on invalid rings.
gi <- gi |> mutate(mdl_key = as.character(glue("rng_iucn|{id_no}")), ds_key = "rng_iucn")
gi <- gi[!st_is_empty(gi), ]
st_write(gi, iucn_gpkg, layer = "rng_iucn", quiet = TRUE, delete_dsn = TRUE)
system2("ogrinfo", c(shQuote(iucn_gpkg), "-sql",
shQuote("CREATE INDEX idx_mdlkey ON rng_iucn(mdl_key)")))
log_info("wrote {path_file(iucn_gpkg)}: {nrow(gi)} species")
}
# the src gpkg must cover the served set, or publish_pmtiles_from_gpkg silently skips the
# species it cannot find — the failure mode that left every widened run short
n_src <- as.integer(DBI::dbGetQuery(
DBI::dbConnect(RSQLite::SQLite(), path.expand(iucn_gpkg)),
"SELECT count(DISTINCT id_no) n FROM rng_iucn")$n)
if (n_src < length(ids_srv))
log_warn("iucn src gpkg holds {n_src} of {length(ids_srv)} served ids — {length(ids_srv) - n_src} will have no PMTiles")
} else log_info("{path_file(iucn_gpkg)} exists — skip (REDO_IUCN_SRC=1 to rebuild)")
```
```{r}
#| label: pmtiles-build
pmt_only <- Sys.getenv("NATIVE_PMT_DS") # optional comma-list to restrict (partial reruns/tests)
pmt_ds <- if (nzchar(pmt_only)) intersect(names(pmt_specs), strsplit(pmt_only, ",")[[1]]) else names(pmt_specs)
pmt_rows <- list()
if (!skip_pmt) for (ds in pmt_ds) {
keys <- served_by_ds[[ds]]
if (is.null(keys)) { log_warn("no served keys for {ds} — skip PMTiles"); next }
dir_ds <- glue("{dir_pmt}/{ds}"); dir_create(dir_ds)
# one PMTiles per model (resumable — skips existing unless REDO_PMTILES/REDO_NATIVE).
# Built to z10 so the deepest tiles carry full source resolution: every higher zoom
# overzooms from maxzoom, so tiles built at z6 were coarser than the source everywhere
# above it. A single range still tiles to a few KB, so there is no density blowup.
#
# REDO_PMTILES exists separately from REDO_NATIVE because the latter also rebuilds the
# ~7 GB IUCN gpkg and re-sorts AquaMaps. Forcing a tile rebuild used to mean moving the
# output directory aside by hand, which leaves no record of what was rebuilt or why.
if (ds == "rng_iucn") {
# IUCN's full-res source is ~7 GB — read each model on demand from the indexed gpkg
# (never load it all into memory), one index hit per species.
log_info("{ds}: {length(keys)} models -> per-model PMTiles from indexed gpkg ({n_workers} workers)")
built <- publish_pmtiles_from_gpkg(iucn_gpkg, "rng_iucn", keys, dir_ds,
workers = n_workers, redo = redo_pmt, minzoom = 0, maxzoom = 10)
} else {
log_info("reading source polygons for {ds} ({length(keys)} served models) …")
x <- pmt_specs[[ds]]()
x <- x |> mutate(mdl_key = as.character(mdl_key)) |>
filter(mdl_key %in% keys) |> mutate(ds_key = ds)
x <- x[!st_is_empty(x), ]
if (!nrow(x)) { log_warn("{ds}: 0 features after served-join — skip"); next }
log_info("{ds}: {nrow(x)} features, {n_distinct(x$mdl_key)} models -> per-model PMTiles ({n_workers} workers)")
built <- publish_pmtiles_models(x, dir_ds, layer = ds, workers = n_workers,
redo = redo_pmt, minzoom = 0, maxzoom = 10)
}
if (is.null(built) || !nrow(built)) { log_warn("{ds}: no PMTiles built"); next }
pmt_rows[[ds]] <- built |> transmute(
mdl_key, ds_key = ds, asset_type = "pmtiles", representation = "native", # source polygons = original
# per-model URL; ?v={mtime} busts the browser/PMTiles-protocol cache on rebuild
# (Caddy ignores the query and serves the file).
asset_url = glue("{pmt_base}/{ds}/{path_file(file)}?v={as.integer(file_info(file)$modification_time)}"),
source_layer = ds)
}
pmt_reg <- bind_rows(pmt_rows) # now ONE ROW PER MODEL (mdl_key, asset_url)
if (nrow(pmt_reg)) knitr::kable(pmt_reg |> count(ds_key, name = "n_models"))
```
## AquaMaps → per-model COGs
Consolidate the served `am` cells by **integer** key (`mid`) via `PARTITION_BY` — one small
Parquet dir per model (`_am_parts/mid={mid}/`, columns `cell_id,val` only). This drops the
40-byte `mdl_key` string from the ~12B-row intermediate and needs no global sort (a plain
ORDER BY would spill ~500GB), so each per-model read is a direct partition scan. Then paint
each model onto the global grid, cropped to its bbox, as a COG with overviews. Parallel +
resumable (skip existing COGs).
```{r}
#| label: am-partition
# test mode uses a separate parts dir + only the test subset, so it never poisons a full run
am_keys <- if (length(served_by_ds[["am"]])) sort(served_by_ds[["am"]]) else character(0)
am_take <- if (test_n > 0) head(am_keys, test_n) else am_keys
# `mid` is POSITIONAL — seq_along over the sorted served keys — so it is a function of the
# whole served SET, not of the model. Widening the served set (is_valid_usa -> is_valid_global)
# moved 16,815 of 16,818 positions: mid=4 stopped meaning `am|Chn-0e47…` and started meaning
# `am|Chn-08d9…`. The partitions were skipped because the DIRECTORY existed, so every COG built
# afterwards read another species' cells and wrote them under the right filename. 1,054 COGs
# came out that way, each internally consistent and silently wrong (am|Fis-112180: 6,469 cells
# in dist, 100,981 in its COG).
#
# So the map is PERSISTED beside the partitions and the partitions are rebuilt whenever it
# changes. Existence is not evidence of correctness; matching the map is. Same lesson as
# msens::assign_mdl_id(), which exists because dense_rank() renumbers published partitions.
am_map <- tibble(mdl_key = am_take, mid = seq_along(am_take))
am_parts<- glue("{dir_native}/_am_parts{if (test_n > 0) '_test' else ''}")
am_map_csv <- glue("{am_parts}.map.csv")
if (!skip_am && length(am_take)) {
map_ok <- file_exists(am_map_csv) && dir_exists(am_parts) && {
prev <- suppressMessages(readr::read_csv(am_map_csv, show_col_types = FALSE))
nrow(prev) == nrow(am_map) && all(prev$mdl_key == am_map$mdl_key) && all(prev$mid == am_map$mid)
}
if (!map_ok || redo) {
if (dir_exists(am_parts)) {
log_warn("_am_parts does not match the current served set — repartitioning ",
"(mid is positional; reusing it would paint the wrong species)")
dir_delete(am_parts)
}
log_info("partitioning {length(am_take)} served am models by mid -> _am_parts …")
con <- dbConnect(duckdb(), config = list(threads = "6", memory_limit = "10GB"))
dbExecute(con, glue("SET temp_directory='{dir_atlas}/duckdb_tmp'"))
dbExecute(con, "SET partitioned_write_max_open_files=400")
dbWriteTable(con, "am_map", am_map, overwrite = TRUE)
msens::copy_atlas_parquet(con, glue("
SELECT map.mid, p.cell_id, p.val
FROM read_parquet('{dir_atlas}/dist/dataset=am/*.parquet') p
JOIN am_map map ON p.mdl_key = map.mdl_key"),
am_parts, partition_by = "mid")
dbDisconnect(con, shutdown = TRUE)
readr::write_csv(am_map, am_map_csv)
log_info("_am_parts: {length(dir_ls(am_parts))} partitions, {round(sum(file_size(dir_ls(am_parts, recurse=TRUE, type='file')))/1e9,1)} GB")
} else log_info("_am_parts matches the served set ({nrow(am_map)} models) — skip repartition")
}
```
```{r}
#| label: am-cogs
am_reg <- tibble()
if (!skip_am && length(am_take) && dir_exists(am_parts)) {
todo <- am_map[!file_exists(glue("{dir_cog}/{safe_key(am_map$mdl_key)}.tif")), ]
log_info("am COGs: {nrow(am_map)} target, {nrow(todo)} to build ({nrow(am_map)-nrow(todo)} exist)")
# PROVE mid still addresses the model it names, on the DATA rather than on the map file.
# A COG painted from the wrong partition is internally valid — right filename, plausible
# raster — so nothing downstream can detect it; the only moment it is cheap to catch is
# here, by asking whether partition mid=N holds the cells `dist` records for that mdl_key.
if (nrow(todo)) {
conv <- dbConnect(duckdb(), config = list(threads = "4"))
smp <- am_map[unique(round(seq(1, nrow(am_map), length.out = min(8, nrow(am_map))))), ]
bad <- 0L
for (i in seq_len(nrow(smp))) {
pf <- glue("{am_parts}/mid={smp$mid[i]}")
if (!dir_exists(pf)) next
n_part <- dbGetQuery(conv, glue("SELECT count(*) n FROM read_parquet('{pf}/*.parquet')"))$n
n_dist <- dbGetQuery(conv, glue("
SELECT count(*) n FROM read_parquet('{dir_atlas}/dist/dataset=am/*.parquet')
WHERE mdl_key = '{smp$mdl_key[i]}'"))$n
if (n_part != n_dist) {
bad <- bad + 1L
log_error("mid={smp$mid[i]} ({smp$mdl_key[i]}): partition has {n_part} cells, dist has {n_dist}")
}
}
dbDisconnect(conv, shutdown = TRUE)
if (bad > 0)
stop(sprintf(paste0("_am_parts does not address the models am_map names (%d of %d sampled ",
"disagree). Building COGs now would publish other species' surfaces ",
"under these names — delete %s and re-run to repartition."),
bad, nrow(smp), am_parts), call. = FALSE)
log_info("am partition addressing verified on {nrow(smp)} sampled models")
}
grid <- grid_spec(rast(cellid_tif))
build_batch <- function(rows) {
suppressMessages({library(DBI); library(duckdb); library(terra)})
con <- dbConnect(duckdb())
on.exit(dbDisconnect(con, shutdown = TRUE))
for (i in seq_len(nrow(rows))) {
k <- rows$mdl_key[i]; f <- sprintf("%s/%s.tif", dir_cog, safe_key(k))
if (file.exists(f)) next
pf <- sprintf("%s/mid=%d", am_parts, rows$mid[i])
if (!dir.exists(pf)) next
d <- dbGetQuery(con, sprintf("SELECT cell_id, val FROM read_parquet('%s/*.parquet')", pf))
if (nrow(d)) msens::publish_cog(d$cell_id, d$val, f, grid)
}
nrow(rows)
}
if (nrow(todo)) {
plan(multisession, workers = n_workers)
batches <- split(todo, (seq_len(nrow(todo)) - 1) %% n_workers)
future_map(batches, build_batch,
.options = furrr_options(globals = c("am_parts","dir_cog","grid","safe_key"),
packages = c("msens"), seed = TRUE))
plan(sequential)
}
built <- glue("{dir_cog}/{safe_key(am_map$mdl_key)}.tif")
# representation = "model": these COGs are painted from dist/dataset=am (the 0.05° bilinear-
# RESAMPLED suitability actually used in scoring) — NOT AquaMaps' native 0.5° cells. The native
# surface is published separately below (representation = "native").
am_reg <- tibble(
mdl_key = am_map$mdl_key, ds_key = "am", asset_type = "cog", representation = "model",
asset_url = glue("{s3_http}/native/am/{safe_key(am_map$mdl_key)}.tif"),
rescale_min = 1L, rescale_max = 100L, colormap = "spectral_r")[file_exists(built), ]
log_info("am model COGs (0.05° resampled) on disk: {nrow(am_reg)} / {nrow(am_map)}")
}
```
## AquaMaps native 0.5° COGs (the *original* surface)
Publish each served AquaMaps model at its **true native 0.5°** resolution — the HCAF `probability`
painted by `loiczid` onto the 720×360 half-degree grid (`publish_cog` is grid-agnostic; `loiczid`
is its row-major top-left index, verified against the ingest). This is the **original**
representation, distinct from the 0.05°-resampled model COGs above; `probability×100` shares the
[0,100] rescale. Parallel + resumable.
```{r}
#| label: am-native-cogs
dir_cog_native <- glue("{dir_native}/am_native"); dir_create(dir_cog_native)
am_native_reg <- tibble()
if (!skip_am && length(am_take) && file_exists(am_db)) {
native_grid <- grid_spec(rast(nrows = 360, ncols = 720, xmin = -180, xmax = 180,
ymin = -90, ymax = 90, crs = "EPSG:4326"))
todo_n <- am_map[!file_exists(glue("{dir_cog_native}/{safe_key(am_map$mdl_key)}.tif")), ]
log_info("am native 0.5° COGs: {nrow(am_map)} target, {nrow(todo_n)} to build")
build_native <- function(rows) {
suppressMessages({library(DBI); library(duckdb); library(terra)})
con <- dbConnect(duckdb(am_db, read_only = TRUE)); on.exit(dbDisconnect(con, shutdown = TRUE))
for (i in seq_len(nrow(rows))) {
k <- rows$mdl_key[i]; sp <- sub("^am\\|", "", k)
f <- sprintf("%s/%s.tif", dir_cog_native, safe_key(k))
if (file.exists(f)) next
d <- dbGetQuery(con, sprintf(
"SELECT c.loiczid AS cell_id, s.probability*100 AS val
FROM spp_cells s JOIN cells c USING (cell_id) WHERE s.sp_key = '%s'", sp))
if (nrow(d)) msens::publish_cog(d$cell_id, d$val, f, native_grid)
}
nrow(rows)
}
if (nrow(todo_n)) {
plan(multisession, workers = n_workers)
batches <- split(todo_n, (seq_len(nrow(todo_n)) - 1) %% n_workers)
future_map(batches, build_native,
.options = furrr_options(globals = c("am_db","dir_cog_native","native_grid","safe_key"),
packages = c("msens"), seed = TRUE))
plan(sequential)
}
built_n <- glue("{dir_cog_native}/{safe_key(am_map$mdl_key)}.tif")
am_native_reg <- tibble(
mdl_key = am_map$mdl_key, ds_key = "am", asset_type = "cog", representation = "native",
asset_url = glue("{s3_http}/native/am_native/{safe_key(am_map$mdl_key)}.tif"),
rescale_min = 1L, rescale_max = 100L, colormap = "spectral_r")[file_exists(built_n), ]
log_info("am native 0.5° COGs on disk: {nrow(am_native_reg)} / {nrow(am_map)}")
}
```
## Vector inputs → 0.05° gridded COGs (the *interpolated* surface, opt-in)
Each vector input's dist surface (`(cell_id, val)` on the 0.05° grid) painted as a COG — the
**model/interpolated** representation of a vector range, complementing its native PMTiles polygons
(so both representations exist per `mdl_key`). Opt-in (`NATIVE_VEC_COG=1`) since there are
thousands; iterates each dist file once (each holds one model) and `publish_cog`s the served ones.
```{r}
#| label: vec-grid-cogs
do_vec_cog <- nzchar(Sys.getenv("NATIVE_VEC_COG"))
vec_cog_reg <- tibble()
if (do_vec_cog && nrow(pmt_reg)) {
dir_vcog <- glue("{dir_native}/vec_grid"); dir_create(dir_vcog)
keys_by_ds <- split(pmt_reg$mdl_key, pmt_reg$ds_key) # served vector models per dataset
grid <- grid_spec(rast(cellid_tif))
build_ds <- function(ds, keys) {
suppressMessages({library(DBI); library(duckdb); library(terra); library(arrow)})
files <- fs::dir_ls(file.path(dir_atlas, paste0("dist/dataset=", ds)), glob = "*.parquet")
con <- dbConnect(duckdb()); on.exit(dbDisconnect(con, shutdown = TRUE))
for (fp in files) {
k <- dbGetQuery(con, sprintf("SELECT DISTINCT mdl_key FROM read_parquet('%s') LIMIT 1", fp))$mdl_key
if (length(k) != 1 || !k %in% keys) next
f <- sprintf("%s/%s.tif", dir_vcog, safe_key(k))
if (file.exists(f)) next
d <- dbGetQuery(con, sprintf("SELECT cell_id, val FROM read_parquet('%s')", fp))
if (nrow(d)) msens::publish_cog(d$cell_id, d$val, f, grid)
}
length(files)
}
plan(multisession, workers = n_workers)
future_map2(names(keys_by_ds), keys_by_ds, build_ds,
.options = furrr_options(globals = c("dir_vcog","dir_atlas","grid","safe_key"),
packages = c("msens"), seed = TRUE))
plan(sequential)
vec_cog_reg <- pmt_reg |> distinct(mdl_key, ds_key) |> mutate(
asset_type = "cog", representation = "model",
asset_url = glue("{s3_http}/native/vec_grid/{safe_key(mdl_key)}.tif"),
rescale_min = 1L, rescale_max = 100L, colormap = "spectral_r")
vec_cog_reg <- vec_cog_reg[file_exists(glue("{dir_vcog}/{safe_key(vec_cog_reg$mdl_key)}.tif")), ]
log_info("vector grid COGs on disk: {nrow(vec_cog_reg)}")
}
```
## NMFS DPS extinction-risk surfaces → per-species COGs
The `dps_nmfs` surfaces (per-cell extinction risk: each listed entity's habitat at its status, the
species' baseline elsewhere) are *gridded from the start* — built by `ingest_nmfs-dps.qmd` from the
NMFS critical-habitat service + the IUCN range — so their one representation is the model grid
itself, painted as a COG per species. Without these rows the species app listed the input as
*NMFS DPS extinction risk* but struck it through: an input with no published surface is disabled.
Resumable by file existence; `REDO_DPS_COG=1` (or `REDO_NATIVE`) repaints.
```{r}
#| label: dps-cogs
dps_reg <- tibble()
dps_keys <- if (length(served_by_ds[["dps_nmfs"]])) sort(served_by_ds[["dps_nmfs"]]) else character(0)
if (length(dps_keys)) {
dir_dps <- glue("{dir_native}/dps_nmfs")
if ((redo || nzchar(Sys.getenv("REDO_DPS_COG"))) && dir_exists(dir_dps)) dir_delete(dir_dps)
dir_create(dir_dps)
grid <- grid_spec(rast(cellid_tif))
con <- dbConnect(duckdb())
for (k in dps_keys) {
f <- glue("{dir_dps}/{safe_key(k)}.tif"); if (file_exists(f)) next
d <- dbGetQuery(con, glue("SELECT cell_id, val FROM read_parquet('{dir_atlas}/dist/dataset=dps_nmfs/*.parquet') WHERE mdl_key = '{k}'"))
if (nrow(d)) msens::publish_cog(d$cell_id, d$val, f, grid)
}
dbDisconnect(con, shutdown = TRUE)
dps_reg <- tibble(
mdl_key = dps_keys, ds_key = "dps_nmfs", asset_type = "cog", representation = "model",
asset_url = glue("{s3_http}/native/dps_nmfs/{safe_key(dps_keys)}.tif"),
rescale_min = 1L, rescale_max = 100L, colormap = "spectral_r")
dps_reg <- dps_reg[file_exists(glue("{dir_dps}/{safe_key(dps_reg$mdl_key)}.tif")), ]
log_info("dps_nmfs COGs on disk: {nrow(dps_reg)} / {length(dps_keys)}")
}
```
## Merged taxon surface → per-taxon whole-range COGs
Publish the **merged model itself** (`ds_key = "ms_merge"`, `representation = "model"`) as a
whole-range COG per taxon, painted from `dist_merged_global` — so the species app renders the merged
layer (and non-US species' global range) via titiler `/cog` (anonymous, clickable) instead of the
creds-gated `model_cell` glob. Partition the global surface by an integer key (`mmid`) — no global
sort — then paint each per-taxon slice. Opt-in (`PUBLISH_MERGED_COG=1`); parallel + resumable.
```{r}
#| label: merged-cogs
do_merged_cog <- nzchar(Sys.getenv("PUBLISH_MERGED_COG"))
# REDO_MERGED_COG exists separately from REDO_NATIVE for the same reason REDO_PMTILES does: the
# merged COGs are a pure function of dist_merged_global, so a re-merge invalidates every one of them
# while the ~7 GB IUCN gpkg and the am re-sort REDO_NATIVE would also redo are untouched. The build
# below is resumable by file existence, so repainting means clearing the outputs -- as a FLAG, never
# a hand-`mv`, so the HTML records which run rebuilt them (apps#8 masking fix).
redo_mcog <- redo || nzchar(Sys.getenv("REDO_MERGED_COG"))
# REDO_MERGED_COG_KEYS="ms_merge|WORMS:137092,..." repaints ONLY those taxa (a rule change for the
# ~25 spatial-ER taxa need not repaint 21,000 COGs); their partitions are refreshed from the global
# surface and their tifs cleared, everything else is kept
redo_keys <- trimws(strsplit(Sys.getenv("REDO_MERGED_COG_KEYS"), ",")[[1]]); redo_keys <- redo_keys[nzchar(redo_keys)]
# REDO_MERGED_COG_SPATIAL=1 names them for you: every taxon merged through the spatial-ER branches
# (sea turtles + NMFS DPS species), resolved from the crosswalk so no key list is pasted by hand
if (nzchar(Sys.getenv("REDO_MERGED_COG_SPATIAL"))) {
conk <- dbConnect(duckdb(), config = list(threads = "2")); dbExecute(conk, glue("ATTACH '{merge_db}' AS m (READ_ONLY)"))
redo_keys <- union(redo_keys, dbGetQuery(conk, "SELECT DISTINCT ms_merge_key FROM m.taxon_model
WHERE mdl_key LIKE 'rng_turtle_swot_dps|%' OR mdl_key LIKE 'dps_nmfs|%'")$ms_merge_key)
dbDisconnect(conk, shutdown = TRUE)
log_info("REDO_MERGED_COG_SPATIAL: {length(redo_keys)} spatial-ER taxa to repaint")
}
merged_reg <- tibble()
merged_map <- tibble()
mparts <- glue("{dir_native}/_merged_parts{if (test_n > 0) '_test' else ''}")
dir_gmerged <- glue("{dir_atlas}/dist_merged_global/dataset=ms_merge")
if (do_merged_cog && dir_exists(dir_gmerged)) {
dir_mcog <- glue("{dir_native}/merged")
if (redo_mcog && dir_exists(dir_mcog)) {
log_info("REDO_MERGED_COG: clearing {length(dir_ls(dir_mcog, glob = '*.tif'))} merged COGs")
dir_delete(dir_mcog)
}
dir_create(dir_mcog)
if (length(redo_keys)) {
rm_tif <- glue("{dir_mcog}/{safe_key(redo_keys)}.tif"); rm_tif <- rm_tif[file_exists(rm_tif)]
if (length(rm_tif)) file_delete(rm_tif)
log_info("REDO_MERGED_COG_KEYS: cleared {length(rm_tif)} of {length(redo_keys)} named merged COGs")
}
con <- dbConnect(duckdb(), config = list(threads = "6", memory_limit = "10GB"))
dbExecute(con, glue("SET temp_directory='{dir_atlas}/duckdb_tmp'"))
dbExecute(con, glue("ATTACH '{merge_db}' AS m (READ_ONLY)"))
# served merged taxa (the set the app lists) present in the global surface
mkeys <- dbGetQuery(con, glue("
SELECT DISTINCT g.mdl_key FROM read_parquet('{dir_gmerged}/**/*.parquet') g
SEMI JOIN (SELECT DISTINCT ms_merge_key AS mdl_key FROM m.taxon
WHERE is_marine AND ms_merge_key IS NOT NULL AND sp_cat NOT IN ('reptile','amphibian')) s
USING (mdl_key)"))$mdl_key
mkeys <- sort(mkeys); if (test_n > 0) mkeys <- head(mkeys, test_n)
merged_map <- tibble(mdl_key = mkeys, mmid = seq_along(mkeys))
if (length(mkeys) && (!dir_exists(mparts) || redo_mcog)) {
if (dir_exists(mparts)) dir_delete(mparts)
dbExecute(con, "SET partitioned_write_max_open_files=400")
dbWriteTable(con, "merged_map", merged_map, overwrite = TRUE)
msens::copy_atlas_parquet(con, glue("
SELECT mm.mmid, g.cell_id, g.val
FROM read_parquet('{dir_gmerged}/**/*.parquet') g JOIN merged_map mm ON g.mdl_key = mm.mdl_key"),
mparts, partition_by = "mmid")
log_info("_merged_parts: {length(dir_ls(mparts))} partitions")
} else if (length(redo_keys) && dir_exists(mparts)) {
# mmid is POSITIONAL in the served set (see the am `mid` lesson above): refreshing a few
# partitions in place is only valid while the set -- and so every position -- is unchanged
stopifnot("_merged_parts partition count != served merged taxa: the map moved, REDO_MERGED_COG=1" =
length(dir_ls(mparts)) == nrow(merged_map))
rk <- merged_map[merged_map$mdl_key %in% redo_keys, ]
for (id in rk$mmid) { pd <- glue("{mparts}/mmid={id}"); if (dir_exists(pd)) dir_delete(pd) }
dbWriteTable(con, "merged_map_redo", rk, overwrite = TRUE)
msens::copy_atlas_parquet(con, glue("
SELECT mm.mmid, g.cell_id, g.val
FROM read_parquet('{dir_gmerged}/**/*.parquet') g JOIN merged_map_redo mm ON g.mdl_key = mm.mdl_key"),
mparts, partition_by = "mmid")
log_info("REDO_MERGED_COG_KEYS: re-partitioned {nrow(rk)} taxa from the global surface")
}
dbDisconnect(con, shutdown = TRUE)
grid <- grid_spec(rast(cellid_tif))
todo <- merged_map[!file_exists(glue("{dir_mcog}/{safe_key(merged_map$mdl_key)}.tif")), ]
log_info("merged COGs: {nrow(merged_map)} target, {nrow(todo)} to build")
build_merged <- function(rows) {
suppressMessages({library(DBI); library(duckdb); library(terra)})
con <- dbConnect(duckdb()); on.exit(dbDisconnect(con, shutdown = TRUE))
for (i in seq_len(nrow(rows))) {
k <- rows$mdl_key[i]; f <- sprintf("%s/%s.tif", dir_mcog, safe_key(k))
if (file.exists(f)) next
pf <- sprintf("%s/mmid=%d", mparts, rows$mmid[i]); if (!dir.exists(pf)) next
d <- dbGetQuery(con, sprintf("SELECT cell_id, val FROM read_parquet('%s/*.parquet')", pf))
if (nrow(d)) msens::publish_cog(d$cell_id, d$val, f, grid)
}
nrow(rows)
}
if (nrow(todo)) {
plan(multisession, workers = n_workers)
future_map(split(todo, (seq_len(nrow(todo)) - 1) %% n_workers), build_merged,
.options = furrr_options(globals = c("mparts","dir_mcog","grid","safe_key"),
packages = c("msens"), seed = TRUE))
plan(sequential)
}
built <- glue("{dir_mcog}/{safe_key(merged_map$mdl_key)}.tif")
merged_reg <- tibble(
mdl_key = merged_map$mdl_key, ds_key = "ms_merge", asset_type = "cog", representation = "model",
asset_url = glue("{s3_http}/native/merged/{safe_key(merged_map$mdl_key)}.tif"),
rescale_min = 1L, rescale_max = 100L, colormap = "spectral_r")[file_exists(built), ]
log_info("merged model COGs on disk: {nrow(merged_reg)} / {nrow(merged_map)}")
}
```
## native_asset registry → sdm.duckdb
```{r}
#| label: registry
# per-model asset rows across representations:
# original/native = vector PMTiles (source polygons) + am native 0.5° COGs
# model/interpolated = am 0.05° resampled COGs + vector 0.05° gridded COGs + merged whole-range COGs
# each keyed on the STABLE mdl_key; the `representation` column distinguishes them. Merged rows carry
# mdl_key == ms_merge_key, so they skip the input→taxon join and are bound in directly below.
pmt_model_rows <- if (nrow(pmt_reg)) pmt_reg else tibble()
vec_cog_rows <- if (exists("vec_cog_reg") && nrow(vec_cog_reg)) vec_cog_reg else tibble()
# AquaX (v9, `ax`): the COGs are built and uploaded by ingest_aquax.qmd (the native IS the model
# grid, so there is nothing to repaint here) and recorded in model_ax.csv with their bbox. This
# chunk only REGISTERS them, in both representations, mirroring am/am_native so the species
# app's Original/Interpolated toggle needs no change. Scoped to the served set like every other
# class. Rows are skipped (not the whole class) when a model's COG url is missing.
ax_csv <- glue("{dir_atlas}/dist/model_ax.csv")
ax_reg <- if (file_exists(ax_csv) && length(served_by_ds[["ax"]])) {
ax_m <- readr::read_csv(ax_csv, show_col_types = FALSE, col_types = readr::cols(.default = "c")) |>
filter(mdl_key %in% served_by_ds[["ax"]]) |>
mutate(across(any_of(c("xmin","xmax","ymin","ymax")), as.numeric))
bind_rows(
ax_m |> filter(!is.na(cog_url)) |>
transmute(mdl_key, ds_key = "ax", asset_type = "cog", representation = "model",
asset_url = cog_url, rescale_min = 1L, rescale_max = 100L, colormap = "spectral_r",
xmin, xmax, ymin, ymax),
ax_m |> filter(!is.na(cog_native_url)) |>
transmute(mdl_key, ds_key = "ax", asset_type = "cog", representation = "native",
asset_url = cog_native_url, rescale_min = 0L, rescale_max = 1000L, colormap = "spectral_r",
xmin, xmax, ymin, ymax))
} else tibble()
if (nrow(ax_reg)) log_info("ax COG registry rows: {nrow(ax_reg)} ({n_distinct(ax_reg$mdl_key)} models, ",
"{sum(ax_reg$representation=='model')} model / {sum(ax_reg$representation=='native')} native)")
if (length(served_by_ds[["ax"]]) && nrow(ax_reg) < 2 * length(served_by_ds[["ax"]]))
log_warn("ax: {length(served_by_ds[['ax']])} served models but only {nrow(ax_reg)} registry rows — run ingest_aquax with AX_COG=1 AX_COG_S3=1")
dps_rows <- if (exists("dps_reg") && nrow(dps_reg)) dps_reg else tibble() # NMFS DPS per-cell ER COGs
asset_info <- bind_rows(am_reg, am_native_reg, pmt_model_rows, vec_cog_rows, ax_reg, dps_rows) |>
mutate(across(everything(), ~ if (is.factor(.)) as.character(.) else .))
# per-model geographic bbox (cell_id row/col arithmetic on the global grid) so the app fits
# to each input's actual extent — NOT the whole world (lands on the globe antipode). am from
# the compact _am_parts (mid->mdl_key); vector ranges from their small dist partitions.
bbox_cols <- c("xmin","xmax","ymin","ymax")
bb_am <- tibble()
if (nrow(am_reg) && exists("am_parts") && dir_exists(am_parts)) {
con <- dbConnect(duckdb(), config = list(threads = "6"))
bb_am <- dbGetQuery(con, glue("
WITH rc AS (SELECT mid, (cell_id-1) % 7200 AS col, (cell_id-1) // 7200 AS row
FROM read_parquet('{am_parts}/**/*.parquet', hive_partitioning=true))
SELECT mid,
-180 + min(col)*0.05 AS xmin, -180 + (max(col)+1)*0.05 AS xmax,
90 - (max(row)+1)*0.05 AS ymin, 90 - min(row)*0.05 AS ymax
FROM rc GROUP BY mid")) |>
dplyr::left_join(am_map, by = "mid") |> dplyr::select(-mid)
dbDisconnect(con, shutdown = TRUE)
}
bb_vec <- tibble()
vec_ds <- unique(c(pmt_model_rows$ds_key, dps_rows$ds_key)) # the pmtiles datasets + the gridded dps surfaces
if (length(vec_ds)) {
con <- dbConnect(duckdb(), config = list(threads = "6"))
globs <- paste0("'", dir_atlas, "/dist/dataset=", vec_ds, "/*.parquet'", collapse = ", ")
dbWriteTable(con, "vk", tibble(mdl_key = unique(c(pmt_model_rows$mdl_key, dps_rows$mdl_key))), overwrite = TRUE)
bb_vec <- dbGetQuery(con, glue("
WITH rc AS (SELECT p.mdl_key, (p.cell_id-1) % 7200 AS col, (p.cell_id-1) // 7200 AS row
FROM read_parquet([{globs}]) p SEMI JOIN vk USING(mdl_key))
SELECT mdl_key,
-180 + min(col)*0.05 AS xmin, -180 + (max(col)+1)*0.05 AS xmax,
90 - (max(row)+1)*0.05 AS ymin, 90 - min(row)*0.05 AS ymax
FROM rc GROUP BY mdl_key"))
dbDisconnect(con, shutdown = TRUE)
}
bb_all <- bind_rows(bb_am, bb_vec)
if (nrow(bb_all) && "mdl_key" %in% names(asset_info)) {
# ax rows already carry their bbox (from model_ax.csv); keep it where the join has none
asset_info <- asset_info |> left_join(bb_all, by = "mdl_key", suffix = c("", ".j"))
for (b in bbox_cols) if (paste0(b, ".j") %in% names(asset_info)) {
asset_info[[b]] <- if (b %in% names(asset_info)) coalesce(asset_info[[b]], asset_info[[paste0(b, ".j")]]) else asset_info[[paste0(b, ".j")]]
asset_info[[paste0(b, ".j")]] <- NULL
}
}
# attach ms_merge_key (the merged taxon key) so the species app can join a taxon
# directly to its input surfaces. A raw model may feed >1 merged taxon -> one row each.
con <- dbConnect(duckdb(sdm_db))
dbExecute(con, glue("ATTACH '{merge_db}' AS m (READ_ONLY)"))
# same scope as `served` above — an asset built for a taxon must be able to find that taxon
tm <- dbGetQuery(con, "SELECT DISTINCT tm.mdl_key, tm.ms_merge_key
FROM m.taxon_model tm JOIN m.taxon t USING(ms_merge_key)
WHERE (t.is_valid_usa OR t.is_valid_global) AND t.is_marine
AND t.sp_cat NOT IN ('reptile','amphibian')")
dbExecute(con, "DETACH m")
native_asset <- if (nrow(asset_info) && "mdl_key" %in% names(asset_info)) {
asset_info |>
inner_join(tm, by = "mdl_key", relationship = "many-to-many") |>
relocate(ms_merge_key)
} else tibble() # no input-model rows this run (e.g. am+pmtiles skipped) — merged rows bound below
# merged whole-range COGs: mdl_key IS the taxon key -> skip the input->taxon join, bind directly
# (built for ALL valid marine taxa incl non-US, so the app can show a non-US species' global range)
if (exists("merged_reg") && nrow(merged_reg)) {
bb_merged <- tibble()
if (exists("mparts") && dir_exists(mparts)) {
conb <- dbConnect(duckdb(), config = list(threads = "6"))
bb_merged <- dbGetQuery(conb, glue("
WITH rc AS (SELECT mmid, (cell_id-1) % 7200 AS col, (cell_id-1) // 7200 AS row
FROM read_parquet('{mparts}/**/*.parquet', hive_partitioning = true))
SELECT mmid, -180 + min(col)*0.05 AS xmin, -180 + (max(col)+1)*0.05 AS xmax,
90 - (max(row)+1)*0.05 AS ymin, 90 - min(row)*0.05 AS ymax
FROM rc GROUP BY mmid")) |> left_join(merged_map, by = "mmid") |> select(-mmid)
dbDisconnect(conb, shutdown = TRUE)
}
merged_rows <- merged_reg |> mutate(ms_merge_key = mdl_key)
if (nrow(bb_merged)) merged_rows <- left_join(merged_rows, bb_merged, by = "mdl_key")
native_asset <- bind_rows(native_asset, relocate(merged_rows, ms_merge_key))
}
# am-only taxa: the merged surface IS the AquaMaps surface ----
# For a taxon whose only source is AquaMaps there is nothing to merge, so the
# pipeline never writes it to dist_merged_global and it gets no merged COG --
# 19,699 served taxa, for which the species app fell back to the custom SQL tile
# factory. Verified on a 150-taxon sample: the merged surface is byte-identical
# to the am surface (same count + bit_xor over (cell_id, val)), so the existing
# native/am COG already IS their merged raster.
#
# Alias rather than duplicate: same asset_url, same bbox, zero new bytes.
have_merged <- native_asset$ms_merge_key[native_asset$ds_key == "ms_merge"]
dbExecute(con, glue("ATTACH IF NOT EXISTS '{merge_db}' AS mrg (READ_ONLY)"))
amonly_keys <- dbGetQuery(con, "
SELECT ms_merge_key, max(mdl_key) AS am_key FROM (
SELECT t.ms_merge_key, tm.mdl_key
FROM mrg.taxon t JOIN mrg.taxon_model tm USING (ms_merge_key)
WHERE t.is_marine AND t.ms_merge_key IS NOT NULL
AND t.sp_cat NOT IN ('reptile','amphibian'))
GROUP BY 1 HAVING count(*) = 1 AND max(split_part(mdl_key,'|',1)) = 'am'")
am_src <- native_asset |>
filter(ds_key == "am", representation == "model", asset_type == "cog") |>
select(am_key = mdl_key, asset_url, any_of(bbox_cols), rescale_min, rescale_max, colormap) |>
distinct(am_key, .keep_all = TRUE)
alias_rows <- amonly_keys |>
filter(!ms_merge_key %in% have_merged) |>
inner_join(am_src, by = "am_key") |>
mutate(mdl_key = ms_merge_key, ds_key = "ms_merge",
asset_type = "cog", representation = "model") |>
select(-am_key)
if (nrow(alias_rows)) {
native_asset <- bind_rows(native_asset, relocate(alias_rows, ms_merge_key))
log_info("am-only merged aliases: {nrow(alias_rows)} taxa now have a merged COG (0 new objects)")
}
# suitability-only taxa that CANNOT be aliased: more than one model, or any AquaX model.
#
# The alias above is a pass-through, and that is only sound when a taxon has exactly one
# contributing model and it is AquaMaps. Two classes fall outside it, and both are absent from
# dist_merged_global (which holds RANGE-derived surfaces only):
# - am-only taxa with 2+ AquaMaps models (192 in v8): the merged surface is a max across them;
# - from v9, taxa with an AquaX model and no range (~8,800): inside the AquaX mask AquaMaps
# was superseded, so their whole-range merged surface is AquaMaps OUTSIDE the mask plus
# AquaX inside it -- a raster that exists nowhere else. Without this chunk the alias count
# fell 14,799 -> 5,957 and those taxa had NO merged COG (the app draws nothing).
#
# The post-supersession merge INPUT (`mc_parts`, partitioned by mkey_id) is exactly that
# surface, so paint one COG per taxon from its partition: max over the surviving models per cell.
dir_parts <- glue("{dir_atlas}/mc_parts")
suit_ds <- c("am", "ax")
suit_sql <- paste(sprintf("'%s'", suit_ds), collapse = ", ")
suitonly <- dbGetQuery(con, glue("
SELECT t.ms_merge_key, mm.mkey_id, count(*) AS n_models,
count(*) FILTER (WHERE split_part(tm.mdl_key,'|',1) = 'ax') AS n_ax
FROM mrg.taxon t JOIN mrg.taxon_model tm USING (ms_merge_key) JOIN mrg.mkey_map mm USING (ms_merge_key)
WHERE t.is_marine AND t.ms_merge_key IS NOT NULL AND t.sp_cat NOT IN ('reptile','amphibian')
GROUP BY 1, 2
HAVING bool_and(split_part(tm.mdl_key,'|',1) IN ({suit_sql}))
AND NOT (count(*) = 1 AND max(split_part(tm.mdl_key,'|',1)) = 'am')"))
suitonly <- suitonly |> filter(!ms_merge_key %in% native_asset$ms_merge_key[
native_asset$ds_key == "ms_merge" & native_asset$asset_type == "cog"])
log_info("suitability-only taxa needing a painted merged COG: {nrow(suitonly)} ({sum(suitonly$n_ax > 0)} with AquaX, {sum(suitonly$n_ax == 0)} multi-AquaMaps)")
if (nrow(suitonly)) {
stopifnot("mc_parts (the partitioned merge input) is missing — run merge_models first" = dir_exists(dir_parts))
dir_mm <- glue("{dir_native}/merged"); dir_create(dir_mm)
grid_mm <- grid_spec(rast(cellid_tif))
suitonly <- suitonly |> mutate(out = as.character(glue("{dir_mm}/{safe_key(ms_merge_key)}.tif")))
todo_mm <- suitonly |> filter(!file_exists(out))
log_info("suitability-only merged COGs: {nrow(suitonly)} target, {nrow(todo_mm)} to build")
paint_batch <- function(rows, dir_parts, grid) {
suppressMessages({library(DBI); library(duckdb); library(terra)})
cx <- dbConnect(duckdb()); on.exit(dbDisconnect(cx, shutdown = TRUE))
for (i in seq_len(nrow(rows))) {
d <- dbGetQuery(cx, sprintf("SELECT cell_id, max(val) AS val FROM read_parquet('%s/mkey_id=%d/*.parquet') GROUP BY 1",
dir_parts, rows$mkey_id[i]))
if (nrow(d)) msens::publish_cog(d$cell_id, round(d$val), rows$out[i], grid)
}
nrow(rows)
}
if (nrow(todo_mm)) {
plan(multisession, workers = n_workers)
future_map(split(todo_mm, (seq_len(nrow(todo_mm)) - 1) %% n_workers), paint_batch,
dir_parts = dir_parts, grid = grid_mm,
.options = furrr_options(seed = TRUE, packages = c("msens", "terra", "DBI", "duckdb")))
plan(sequential)
}
built <- suitonly |> filter(file_exists(out))
# bbox off the PUBLISHED raster, so the registry describes the artifact a client will fetch
bb <- do.call(rbind, lapply(built$out, function(f) { e <- terra::ext(terra::rast(f)); c(e$xmin, e$xmax, e$ymin, e$ymax) }))
mm_rows <- tibble(
ms_merge_key = built$ms_merge_key, mdl_key = built$ms_merge_key, ds_key = "ms_merge",
asset_type = "cog", representation = "model",
asset_url = glue("{s3_http}/native/merged/{safe_key(built$ms_merge_key)}.tif"),
xmin = bb[, 1], xmax = bb[, 2], ymin = bb[, 3], ymax = bb[, 4],
rescale_min = 1L, rescale_max = 100L, colormap = "spectral_r")
native_asset <- bind_rows(native_asset, mm_rows)
log_info("suitability-only merged COGs registered: {nrow(mm_rows)} taxa ({nrow(todo_mm)} newly built)")
# merged_reg drives the S3 sync of native/merged below; make sure these ride along
merged_reg <- if (exists("merged_reg") && nrow(merged_reg)) bind_rows(merged_reg, mm_rows |> select(any_of(names(merged_reg)))) else mm_rows
}
# the registry is CUMULATIVE: a run may not delete an asset class it did not rebuild ----
#
# Everything above assembles `native_asset` from what THIS run built, and the write below is
# an overwrite — so a run that skips a build chunk (NATIVE_SKIP_PMTILES, or an opt-in class
# like PUBLISH_MERGED_COG left off) used to publish a registry missing every class it never
# touched, while those assets sat published and unreferenced. That is how all 2,234 vector-range
# PMTiles rows vanished: the species app then showed AquaMaps as each taxon's only input, and the
# content-hash checkpoint below faithfully fingerprinted the truncated table (n_pmtiles_native: 0).
#
# msens::registry_merge() carries forward any class this run produced no rows for and ERRORS on a
# class that came back smaller (a partial build must not publish as a complete one). The rule lives
# in the package so testthat asserts it — see msens/tests/testthat/test-publish.R.
prior_asset <- if (dbExistsTable(con, "native_asset")) dbReadTable(con, "native_asset") else NULL
native_asset <- msens::registry_merge(
native_asset, prior_asset,
allow_shrink = nzchar(Sys.getenv("NATIVE_REGISTRY_REBUILD")))
carried <- attr(native_asset, "carried")
if (length(carried))
log_warn("registry: carried forward {length(carried)} class(es) not built this run: {paste(carried, collapse='; ')}")
dbWriteTable(con, "native_asset", native_asset, overwrite = TRUE)
smry <- dbGetQuery(con, "SELECT representation, asset_type, count(*) n_rows,
count(DISTINCT mdl_key) n_models, count(DISTINCT ms_merge_key) n_taxa, count(DISTINCT ds_key) n_ds
FROM native_asset GROUP BY 1, 2 ORDER BY 1, 2")
dbDisconnect(con, shutdown = TRUE)
msens::report_table(smry, caption = "native_asset registry by representation × asset type")
```
## Sync → S3 (COGs) + file host (PMTiles)
```{r}
#| label: sync
aws_sync <- function(local, prefix) {
st <- system2("aws", c("s3","sync", local, glue("{s3_ver}/{prefix}"),
"--only-show-errors","--no-progress"), stdout=TRUE, stderr=TRUE)
if (!is.null(attr(st,"status")) && attr(st,"status")!=0) stop(paste(st, collapse="\n"))
log_info("synced {path_file(local)} -> {s3_ver}/{prefix}")
}
if (do_s3) {
if (nrow(am_reg)) aws_sync(dir_cog, "native/am") # model (0.05° resampled)
if (nrow(am_native_reg)) aws_sync(dir_cog_native, "native/am_native") # original (0.5° native)
if (exists("vec_cog_reg") && nrow(vec_cog_reg))
aws_sync(glue("{dir_native}/vec_grid"), "native/vec_grid") # vector model grid
if (exists("merged_reg") && nrow(merged_reg))
aws_sync(glue("{dir_native}/merged"), "native/merged") # merged whole-range COGs
if (exists("dps_reg") && nrow(dps_reg))
aws_sync(glue("{dir_native}/dps_nmfs"), "native/dps_nmfs") # NMFS DPS per-cell ER surfaces
# PMTiles -> server file host (/share/data/derived/pmtiles/{ver}); also keep on S3
if (nrow(pmt_reg)) {
system2("ssh", c("msens", glue("mkdir -p /share/data/derived/pmtiles/{ver}")))
system2("rsync", c("-aq", glue("{dir_pmt}/"), glue("msens:/share/data/derived/pmtiles/{ver}/")))
st <- system2("aws", c("s3","sync", dir_pmt, glue("{s3_ver}/native/pmtiles"),
"--only-show-errors","--no-progress"), stdout=TRUE, stderr=TRUE)
log_info("synced PMTiles -> file host /pmtiles/{ver} + {s3_ver}/native/pmtiles")
}
} else log_info("NATIVE_NO_S3 set — staged locally only ({dir_native})")
```
## Manifest
```{r}
#| label: manifest
# content-addressed: fingerprint the native_asset registry table (reopen sdm.duckdb
# read-only — the registry chunk already disconnected). Deterministic: no wall-clock,
# no machine paths, no run-mode flags.
con <- dbConnect(duckdb(sdm_db, read_only = TRUE))
h <- msens::hash_query(con, "native_asset")
dbDisconnect(con, shutdown = TRUE)
msens::write_manifest(
manifest, target = "publish_native", content_hash = h,
stats = list(
ver = ver,
pmtiles_datasets = if (nrow(pmt_reg)) paste(sort(unique(pmt_reg$ds_key)), collapse = ",") else "",
n_pmtiles_native = nrow(pmt_reg), # vector original polygons
am_cogs_model = nrow(am_reg), # am 0.05° resampled
am_cogs_native = if (exists("am_native_reg")) nrow(am_native_reg) else 0L, # am 0.5° original
vec_cogs_model = if (exists("vec_cog_reg")) nrow(vec_cog_reg) else 0L, # vector 0.05° gridded
native_assets = nrow(native_asset)),
force = msens::force_target("publish_native"))
```