---
title: "Publish input + merged surfaces to the asset store, and the release's pointer table"
msens:
target_name: publish_native
workflow_type: release
dependency: [merge_taxon]
output: data/manifests/publish_native.json
editor_options:
chunk_output_type: console
---
A release **owns pointers, never files** (asset store, 2026-10). Every distribution file this notebook
produces — the per-model COGs, the per-model range PMTiles, the merged whole-range COGs — goes into the
**content-addressed store** shared by every release, and the release gets a **pointer table**
(`native_asset`: `mdl_key × representation → asset_url`) that names those objects. Nothing is ever
written under `{ver}/`.
```
marine-atlas/cog/global05/{hash}.tif rasters ON the analysis grid: per-input "model", merged, score
marine-atlas/native/{ds_key}/{hash}.tif source-resolution originals (AquaMaps 0.5°, AquaX as delivered)
marine-atlas/native/{ds_key}/{hash}.pmtiles source polygons of a vector range, one PMTiles per model
marine-atlas/assets.parquet the catalog: one row per stored object
{ver}/tables/native_asset.parquet the release's pointers (pushed by release_marine-atlas.qmd)
```
**The key is known before the file is built.** It is the hash of the content the file will hold — the
pixels a COG stores (`msens::pixel_hashes()`: the rows quantised as the writer quantises them, repeated
cells collapsed by `max()` as the merge consumes them) or the source geometry of a range
(`msens::native_vector_hash()`) — folded with the file family's encoding tag
(`msens::asset_store_key()`). So a surface the catalog already lists costs neither a build nor an
upload; a changed surface gets a new key, and older releases keep pointing at the old, still-present
object. The object's MD5 is an integrity column of the catalog, never a key.
**Order of operations (and why).** Keys → build only what the catalog lacks → decode-verify every new COG
against its key → upload → catalogue → **gate** (`msens::store_publish_gate()`: every pointer is a
catalogued store object; nothing bulk under `{ver}/`) → write the pointer table. A pointer is written
only after its object is catalogued, so a reader can never follow one into a 404.
Two representations per input where the source has one of its own: `native` (the original: polygons,
the 0.5° AquaMaps cells, the AquaX band as delivered) and `model` (the 0.05° gridded surface the merge
consumed). AquaX is delivered ON the grid, so its pair is *Delivered / As ingested* (`dataset.on_grid`).
The served set is every taxon 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.
**Flags:** `NATIVE_NO_S3=1` (dry run: keys, builds and the upload plan; no upload, no catalog, no
pointer table), `NATIVE_TEST_N=<n>` (first n keys per class), `NATIVE_SKIP_AM` / `NATIVE_SKIP_PMTILES`,
`NATIVE_PMT_DS=a,b` (vector datasets to tile), `NATIVE_VEC_COG=1` / `PUBLISH_MERGED_COG=1` (opt-in
classes), `REDO_MERGED_COG=1` / `REDO_MERGED_COG_KEYS` / `REDO_MERGED_COG_SPATIAL=1` (repartition the
merged surface; the key then decides whether anything is new), `REDO_IUCN_SRC=1`, `REDO_KEYS=1`
(recompute the cached content hashes), `NATIVE_NO_VERIFY=1` (skip the decode check of new COGs),
`NATIVE_LEGACY_OK=1` (before the M6 prune: v8/v9 still carry pre-store copies under `{ver}/native/`;
the gate then warns instead of stopping), `NATIVE_REGISTRY_REBUILD=1` (publish ONLY what this run
built; otherwise a class this run did not build is carried forward from the published table),
`NATIVE_WORKERS=<n>`, `NATIVE_CATALOG=<assets.parquet>` (read the catalog from a file instead of the
bucket — tests).
> **`NATIVE_SKIP_AM` cannot be combined with a merged rebuild.** The am-only *aliases* — `ms_merge`
> pointer rows at an existing AquaMaps `model` object — are derived from the am rows, so skipping am
> empties their source and the `ms_merge` class shrinks; `registry_merge()` refuses that, but the run
> has spent its time by then.
## Design
```{mermaid}
%%| label: fig-design
%%| fig-cap: "Keys from content → build only what the catalog lacks → verify → upload → catalogue → gate → pointers"
flowchart LR
src["dist / merge.duckdb / am.duckdb<br/>source polygons"] --> keys["content hash ⊕ encoding<br/>= store key"]
cat[("assets.parquet<br/>(catalog)")] --> keys
keys -->|"not catalogued"| build["paint COG / tile PMTiles<br/>into the local store mirror"]
build --> verify["decode == key?"]
verify --> up["upload cog/… native/…"]
up --> cat2[("catalog + new rows")]
cat2 --> gate{"store_publish_gate"}
keys -->|"catalogued"| gate
gate --> ptr[("native_asset<br/>{ver}/tables/")]
```
## Setup
```{r}
#| label: setup
librarian::shelf(arrow, DBI, dplyr, duckdb, fs, furrr, glue, here, jsonlite, logger, purrr,
readr, sf, stringr, terra, tibble, MarineSensitivity/msens, quiet = T)
library(msens)
stopifnot("needs msens >= 0.53.0 (asset_store_key, store_publish_gate)" = packageVersion("msens") >= "0.53.0")
source(here("libs/paths.R"))
source(here("libs/url_gate.R"))
source(here("libs/app_bundle.R")) # app_bundle_aws_retry(): one `aws` call, retried only when S3 never answered
dir_atlas <- path.expand(glue("{dir_big_v}/marine-atlas"))
merge_db <- glue("{dir_atlas}/merge.duckdb")
dir_native <- glue("{dir_atlas}/native"); dir_create(dir_native) # intermediates: partitions, tile work dirs, the IUCN gpkg
dir_as <- path.expand(glue("{dir_derived}/asset_store")) # shared by every release, like the store itself
dir_store <- glue("{dir_as}/store") # local mirror of the store: {dir_store}/{key}
test_n <- as.integer(Sys.getenv("NATIVE_TEST_N", "0"))
dir_keys <- glue("{dir_as}/keys/{ver}{if (test_n > 0) '_test' else ''}") # a test subset must never feed a full run's cache
dir_create(c(dir_store, dir_keys, glue("{dir_as}/duckdb_tmp")))
manifest <- here("data/manifests/publish_native.json"); dir_create(path_dir(manifest))
stopifnot(all(file_exists(c(sdm_db, merge_db, cellid_tif))))
do_s3 <- !nzchar(Sys.getenv("NATIVE_NO_S3"))
skip_am <- nzchar(Sys.getenv("NATIVE_SKIP_AM"))
skip_pmt <- nzchar(Sys.getenv("NATIVE_SKIP_PMTILES"))
redo_keys <- nzchar(Sys.getenv("REDO_KEYS"))
do_verify <- !nzchar(Sys.getenv("NATIVE_NO_VERIFY"))
legacy_ok <- nzchar(Sys.getenv("NATIVE_LEGACY_OK"))
n_workers <- as.integer(Sys.getenv("NATIVE_WORKERS", "6"))
base <- atlas_base_url() # https://s3.us-east-1.amazonaws.com/oceanmetrics.io-public/marine-atlas
bucket_root <- "s3://oceanmetrics.io-public/marine-atlas"
grid <- grid_spec_for("global05", path.expand(cellid_tif))
if (test_n > 0 && do_s3) { log_warn("NATIVE_TEST_N={test_n}: a subset is never published -- forcing NATIVE_NO_S3"); do_s3 <- FALSE }
safe_key <- function(k) gsub("[^A-Za-z0-9._-]", "_", k) # mdl_key -> filesystem-safe stem (work dirs only)
store_path <- function(key) { p <- as.character(glue("{dir_store}/{key}")); dir_create(unique(path_dir(p))); p }
store_url <- function(key) as.character(glue("{base}/{key}"))
# the budget is SET, not passed as config: `dbConnect(duckdb(), config = list(memory_limit = ...))` is silently
# ignored by the installed driver (the connection reported 19.1 GiB / 8 threads, DuckDB's defaults, 2026-10-07)
new_con <- function(mem = "10GB", threads = "6") {
con <- dbConnect(duckdb())
dbExecute(con, glue("SET memory_limit='{mem}'; SET threads={threads}; SET temp_directory='{dir_as}/duckdb_tmp'"))
stopifnot(dbGetQuery(con, "SELECT current_setting('memory_limit') m")$m != "19.1 GiB" || mem == "19.1 GiB")
con
}
# a content-hash table is cached per (release, class) and recomputed when its source changed (file count,
# bytes, newest mtime) or REDO_KEYS=1 -- hashing dist/am is minutes, not seconds
src_sig <- function(paths) {
fi <- file_info(dir_ls(paths[dir_exists(paths)], recurse = TRUE, type = "file"))
paste(nrow(fi), sum(as.numeric(fi$size)), as.integer(max(c(fi$modification_time, as.POSIXct(0, origin = "1970-01-01")))), sep = "-")
}
cached_keys <- function(name, src_dirs, expr) {
f <- glue("{dir_keys}/{name}.parquet"); sig <- src_sig(src_dirs)
if (!redo_keys && file_exists(f)) {
x <- as_tibble(read_parquet(f))
if (identical(x$src_sig[1], sig)) return(select(x, -src_sig))
log_info("keys {name}: source changed, recomputing")
}
t0 <- Sys.time(); x <- as_tibble(expr); x$src_sig <- sig
write_parquet(x, f, compression = "zstd")
log_info("keys {name}: {nrow(x)} models in {round(as.numeric(Sys.time() - t0, units = 'mins'), 1)} min")
select(x, -src_sig)
}
# a big source is hashed in FILE BATCHES, each its own read_parquet([...]): the GROUP BY (by, cell_id)
# then never holds more than one batch. The whole-glob form over dist/am (50 GB, 23,699 one-model files)
# spilled 138 GB of DuckDB temp and filled the laptop disk (2026-10-06), and over dist_merged_global
# (3.1 G rows) ran out of memory outright (2026-10-07). Valid only where a model's rows never straddle two
# batches -- a key seen in two batches stops the run rather than hash half a surface. `files` may be globs
# (one per bucket dir or hive partition); `from(lst)` builds the source from the quoted file list
pixel_hashes_batched <- function(files, n_batch, by = "mdl_key", from = function(lst) glue("read_parquet([{lst}])"), ...) {
b <- split(files, ceiling(seq_along(files) / n_batch))
out <- vector("list", length(b))
for (i in seq_along(b)) {
con <- new_con()
lst <- paste(sprintf("'%s'", b[[i]]), collapse = ",")
out[[i]] <- pixel_hashes(con, from(lst), by, ...)
dbDisconnect(con, shutdown = TRUE)
if (i %% 5 == 0 || i == length(b)) log_info(" hashed file batch {i}/{length(b)}")
}
h <- bind_rows(out)
stopifnot("a model's rows straddle two file batches -- batching by file is invalid for this source" =
!anyDuplicated(h[[by]]))
h
}
take_n <- function(x) if (test_n > 0) head(x, test_n) else x
log_info("publish_native {ver} | s3 {do_s3} | test_n {test_n} | store mirror {dir_store}")
```
## Served raw mdl_keys (the set the app lists)
```{r}
#| label: served
con <- new_con()
dbExecute(con, glue("ATTACH '{merge_db}' AS m (READ_ONLY)"))
# the UNION of the two validity flags: 46 taxa are US-valid but absent from the global merged surface,
# 4,473 are global-only; scoping to either flag alone dropped the other's inputs
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)))
```
## The catalog
Read anonymously from the bucket (the bucket denies anonymous *listing*, so "already stored?" is an
anti-join against this file). Every class below asks `in_cat()` before building anything.
```{r}
#| label: catalog
cat_src <- Sys.getenv("NATIVE_CATALOG")
fetch <- function(url) { f <- tempfile(fileext = ".parquet"); utils::download.file(url, f, mode = "wb", quiet = TRUE); f }
catalog <- as_tibble(asset_catalog_check(as.data.frame(read_parquet(
if (nzchar(cat_src)) cat_src else fetch(glue("{base}/assets.parquet"))))))
in_cat <- function(key) key %in% catalog$key
log_info("catalog: {nrow(catalog)} objects ({sum(catalog$store == 'cog')} cog, {sum(catalog$store == 'native')} native)")
# plan rows accumulate here: one per (mdl_key, representation) this run can point at
plan <- list()
add_plan <- function(d, class) { # d: mdl_key, key, hash, content_hash, enc, file (local, may be NA when catalogued)
stopifnot(all(c("mdl_key", "key", "hash", "content_hash", "enc") %in% names(d)))
d$class <- class; d$catalogued <- in_cat(d$key)
plan[[class]] <<- d
log_info("{class}: {nrow(d)} models -> {n_distinct(d$key)} keys; {sum(d$catalogued)} catalogued, {sum(!d$catalogued)} new")
invisible(d)
}
```
## Vector ranges → per-model PMTiles (the *original* representation)
Each dataset's source polygons are read, `mdl_key` derived from the documented source id, inner-joined
to the served set, and hashed per model (`native_vector_hash()`: the same readers and normalisation M1
used, so an unchanged range resolves to the key already in the store). Only models whose key the
catalog lacks are tiled — one PMTiles per model (z0–10; packing thousands of ranges into one archive
overflowed maplibre's per-tile budget at low zoom).
```{r}
#| label: pmtiles-specs
source(here("libs/native_sources.R")) # readers shared with map_asset_store.qmd (same geometry -> same key)
iucn_gpkg <- glue("{dir_native}/src/rng_iucn.gpkg")
pmt_specs <- native_pmt_specs(dir_raw, dir_big, iucn_gpkg, bl_sisids = served_by_ds[["bl"]])
```
### IUCN indexed ranges GeoPackage
Built once from the ingest's own filtered output (presence 1:3 / marine / category already applied),
one feature set per served `id_no`, indexed on `mdl_key`, so every per-species read is an index hit.
```{r}
#| label: iucn-src
dir_create(path_dir(iucn_gpkg))
# REDO_IUCN_SRC: this gpkg is filtered to the SERVED ids, so widening the served set silently leaves the
# new species out of it -- the per-model reader then finds no features and skips them, with no error
redo_src <- nzchar(Sys.getenv("REDO_IUCN_SRC"))
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). `geom` MUST be in the SELECT list: a GPKG without a geometry
# column still reads back as a plain data.frame with the right row count, and fails much later.
ids_sql <- paste(ids_srv, collapse = ",")
if (file_exists(iucn_gpkg)) file_delete(iucn_gpkg)
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)")))
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)
})
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")
}
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")
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 }
work <- glue("{dir_native}/pmtiles/{ds}"); dir_create(work) # per-model work files, {slug}.pmtiles (resumable)
# 1. keys from the SOURCE geometry (not from the built file: a re-tile of the same range is the same object)
if (ds == "rng_iucn") {
h <- cached_keys(glue("pmtiles_{ds}"), path_dir(iucn_gpkg), {
map_dfr(keys, function(k) {
part <- st_read(iucn_gpkg, quiet = TRUE, query = sprintf("SELECT * FROM \"rng_iucn\" WHERE mdl_key = '%s'", k))
if (!nrow(part)) return(NULL)
part$mdl_key <- k; part <- st_zm(part, drop = TRUE)
tibble(mdl_key = k, n_features = nrow(part), content_hash = native_vector_hash(part))
})
})
x <- NULL
} 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 }
h <- as_tibble(native_vector_hashes(x))
}
h <- take_n(h)
k <- asset_store_key("native_pmtiles", ds, h$content_hash) |> mutate(mdl_key = h$mdl_key)
k$file <- store_path(k$key)
todo <- k |> filter(!in_cat(key), !file_exists(file))
# 2. tile only the models whose key is neither catalogued nor already staged
todo <- todo[!duplicated(todo$key), ] # two models with identical geometry = one object
if (nrow(todo)) {
log_info("{ds}: tiling {nrow(todo)} of {nrow(k)} models ({n_workers} workers)")
built <- if (ds == "rng_iucn")
publish_pmtiles_from_gpkg(iucn_gpkg, "rng_iucn", todo$mdl_key, work, workers = n_workers, minzoom = 0, maxzoom = 10)
else
publish_pmtiles_models(x[x$mdl_key %in% todo$mdl_key, ], work, layer = ds, workers = n_workers, minzoom = 0, maxzoom = 10)
if (!is.null(built) && nrow(built)) {
j <- match(built$mdl_key, todo$mdl_key)
for (i in which(!file_exists(todo$file[j]))) file_copy(built$file[i], todo$file[j][i])
}
miss <- todo$mdl_key[!file_exists(todo$file)]
if (length(miss)) log_warn("{ds}: {length(miss)} models produced no PMTiles (e.g. {miss[1]})")
}
k <- k |> filter(in_cat(key) | file_exists(file))
pmt_rows[[ds]] <- add_plan(k |> mutate(ds_key = ds), glue("pmtiles_{ds}"))
}
pmt_reg <- bind_rows(pmt_rows)
if (nrow(pmt_reg)) knitr::kable(pmt_reg |> count(ds_key, catalogued, name = "n_models"))
```
## AquaMaps → per-model COGs (the *model* representation, 0.05°)
Keys first: the pixel hash of each served model's `dist` rows (truncated into INT1U, sub-1 pixels
dropped, repeated cells collapsed by `max`). Then consolidate the served cells by **integer** key
(`mid`) via `PARTITION_BY` — one small Parquet dir per model, no global sort (a plain `ORDER BY` would
spill ~500 GB) — and paint only the models whose key is new.
```{r}
#| label: am-keys
am_keys <- if (length(served_by_ds[["am"]])) sort(served_by_ds[["am"]]) else character(0)
am_take <- take_n(am_keys)
am_reg <- tibble()
if (!skip_am && length(am_take)) {
# one file per model, so batches of files are batches of models (the straddle assertion holds it to that);
# every model is hashed, served or not, so the cache is complete for any later served set
h <- cached_keys("cog_am", glue("{dir_atlas}/dist/dataset=am"), {
pixel_hashes_batched(sort(dir_ls(glue("{dir_atlas}/dist/dataset=am"), glob = "*.parquet")), n_batch = 800,
quant = "trunc", dedup = "max")
})
h <- h |> filter(mdl_key %in% am_take)
am_reg <- asset_store_key("cog_model", "global05", h$content_hash) |> mutate(mdl_key = h$mdl_key, ds_key = "am")
am_reg$file <- store_path(am_reg$key)
}
```
```{r}
#| label: am-partition
# `mid` is POSITIONAL -- seq_along over the sorted served keys -- so it is a function of the whole served
# SET. Widening the set once moved 16,815 of 16,818 positions while the partitions were kept because the
# DIRECTORY existed, and 1,054 COGs were painted from another species' cells. 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.
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")
am_todo <- if (nrow(am_reg)) am_reg |> filter(!in_cat(key), !file_exists(file)) else tibble()
am_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 (!skip_am && length(am_take) && nrow(am_todo)) {
if (!am_map_ok) {
if (dir_exists(am_parts)) {
log_warn("_am_parts does not match the current served set — repartitioning (mid is positional)")
dir_delete(am_parts)
}
log_info("partitioning {length(am_take)} served am models by mid -> _am_parts …")
con <- new_con()
dbExecute(con, "SET partitioned_write_max_open_files=400")
dbWriteTable(con, "am_map", am_map, overwrite = TRUE)
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); am_map_ok <- TRUE
log_info("_am_parts: {length(dir_ls(am_parts))} partitions")
} else log_info("_am_parts matches the served set ({nrow(am_map)} models) — skip repartition")
} else if (!skip_am && length(am_take)) log_info("am: every served model's key is catalogued or staged — no partition needed")
```
```{r}
#| label: am-cogs
# a painter shared by every gridded class: rows (mdl_key, out, sql) -> publish_cog() into the store mirror.
# `sql` selects (cell_id, val) for the model; `round_first` applies the writer's own rounding (ax, merged-suit)
paint_rows <- function(rows, grid, round_first = FALSE, workers = n_workers) {
rows <- rows[!file_exists(rows$out), ]
rows <- rows[!duplicated(rows$out), ] # two models with identical pixels = one object
if (!nrow(rows)) return(invisible(0L))
one_batch <- function(b) {
suppressMessages({library(DBI); library(duckdb); library(terra)})
con <- dbConnect(duckdb()); dbExecute(con, "SET memory_limit='3GB'; SET threads=2") # SET: the config list is ignored
on.exit(dbDisconnect(con, shutdown = TRUE))
for (i in seq_len(nrow(b))) {
if (file.exists(b$out[i])) next
d <- dbGetQuery(con, b$sql[i])
if (round_first) d$val <- round(d$val)
if (nrow(d)) msens::publish_cog(d$cell_id, d$val, b$out[i], grid)
}
nrow(b)
}
if (workers > 1 && nrow(rows) > 1) {
plan(multisession, workers = workers); on.exit(plan(sequential), add = TRUE)
future_map(split(rows, (seq_len(nrow(rows)) - 1) %% workers), one_batch,
.options = furrr_options(packages = c("msens"), seed = TRUE))
} else one_batch(rows)
invisible(sum(file_exists(rows$out)))
}
if (nrow(am_todo) && dir_exists(am_parts)) {
todo <- am_todo |> inner_join(am_map, by = "mdl_key")
# 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, so this is the only cheap moment to catch it
conv <- new_con("2GB", "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 {n_part} cells, dist {n_dist}") }
}
dbDisconnect(conv, shutdown = TRUE)
if (bad > 0) stop(sprintf("_am_parts does not address the models am_map names (%d of %d sampled disagree); delete %s and re-run",
bad, nrow(smp), am_parts), call. = FALSE)
log_info("am partition addressing verified on {nrow(smp)} sampled models; painting {nrow(todo)} model COGs")
# max(val) per cell: what the merge consumes (and what the key hashed); a no-op where cells do not repeat
todo$sql <- sprintf("SELECT cell_id, max(val) AS val FROM read_parquet('%s/mid=%d/*.parquet') GROUP BY 1", am_parts, todo$mid)
todo$out <- todo$file
paint_rows(todo, grid)
}
if (nrow(am_reg)) {
am_reg <- am_reg |> filter(in_cat(key) | file_exists(file))
add_plan(am_reg, "cog_am")
}
```
## AquaMaps native 0.5° COGs (the *original* representation)
Each served model at its true native resolution: the HCAF `probability × 100` painted by `loiczid` onto
the 720 × 360 grid (`loiczid` is its row-major top-left index, verified against the ingest). Keyed by the
same pixel rule (`trunc` into INT1U) as every COG `publish_cog()` writes.
```{r}
#| label: am-native-cogs
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"))
h <- cached_keys("native_am", path_dir(am_db), {
con <- dbConnect(duckdb(path.expand(am_db), read_only = TRUE))
dbExecute(con, glue("SET memory_limit='10GB'; SET threads=6; SET temp_directory='{dir_as}/duckdb_tmp'"))
on.exit(dbDisconnect(con, shutdown = TRUE))
duckdb_register(con, "served_k", tibble(sp_key = sub("^am\\|", "", am_take))) # a VIEW over the df: the db is read-only
pixel_hashes(con, "(SELECT 'am|' || s.sp_key AS mdl_key, c.loiczid AS cell_id, s.probability * 100 AS val
FROM spp_cells s JOIN cells c USING (cell_id) SEMI JOIN served_k USING (sp_key)) z",
"mdl_key", quant = "trunc", dedup = "max")
})
h <- h |> filter(mdl_key %in% am_take)
am_native_reg <- asset_store_key("native_am", "am", h$content_hash) |> mutate(mdl_key = h$mdl_key, ds_key = "am")
am_native_reg$file <- store_path(am_native_reg$key)
todo <- am_native_reg |> filter(!in_cat(key), !file_exists(file))
if (nrow(todo)) {
log_info("am native 0.5° COGs: {nrow(am_native_reg)} served, {nrow(todo)} to paint")
todo$sql <- sprintf("SELECT c.loiczid AS cell_id, max(s.probability * 100) AS val FROM spp_cells s JOIN cells c USING (cell_id) WHERE s.sp_key = '%s' GROUP BY 1",
sub("^am\\|", "", todo$mdl_key))
todo <- todo[!duplicated(todo$key), ]
todo$out <- todo$file
paint_native <- function(rows, grid) {
suppressMessages({library(DBI); library(duckdb); library(terra)})
con <- dbConnect(duckdb(path.expand(am_db), read_only = TRUE)); on.exit(dbDisconnect(con, shutdown = TRUE))
for (i in seq_len(nrow(rows))) {
if (file.exists(rows$out[i])) next
d <- dbGetQuery(con, rows$sql[i])
if (nrow(d)) msens::publish_cog(d$cell_id, d$val, rows$out[i], grid)
}
nrow(rows)
}
plan(multisession, workers = n_workers)
future_map(split(todo, (seq_len(nrow(todo)) - 1) %% n_workers), paint_native, grid = native_grid,
.options = furrr_options(packages = c("msens"), seed = TRUE))
plan(sequential)
}
am_native_reg <- am_native_reg |> filter(in_cat(key) | file_exists(file))
add_plan(am_native_reg, "native_am")
}
```
## Vector inputs → 0.05° gridded COGs (the *model* representation, opt-in)
Each vector input's `dist` surface painted as a COG, so both representations exist per `mdl_key`.
Opt-in (`NATIVE_VEC_COG=1`): thousands of objects, and the key decides what is new. A model whose rows
span several files (`rng_turtle_swot_dps`) is read across all of them with `max(val)` per cell.
```{r}
#| label: vec-grid-cogs
vec_cog_reg <- tibble()
if (nzchar(Sys.getenv("NATIVE_VEC_COG")) && nrow(pmt_reg)) {
keys_by_ds <- split(pmt_reg$mdl_key, pmt_reg$ds_key)
vec_cog_reg <- map_dfr(names(keys_by_ds), function(ds) {
dd <- glue("{dir_atlas}/dist/dataset={ds}")
if (!dir_exists(dd)) return(NULL)
h <- cached_keys(glue("cog_{ds}"), dd, {
con <- new_con(); on.exit(dbDisconnect(con, shutdown = TRUE))
pixel_hashes(con, glue("read_parquet('{dd}/*.parquet')"), "mdl_key", quant = "trunc", dedup = "max")
}) |> filter(mdl_key %in% keys_by_ds[[ds]])
k <- asset_store_key("cog_model", "global05", h$content_hash) |> mutate(mdl_key = h$mdl_key, ds_key = ds)
k$file <- store_path(k$key)
todo <- k |> filter(!in_cat(key), !file_exists(file))
if (nrow(todo)) {
# one file per model where the ingest wrote one (most datasets); every file otherwise (turtles: 14 files, 6 models)
con <- new_con("2GB", "4")
fpm <- dbGetQuery(con, glue("SELECT mdl_key, count(DISTINCT filename) n_files, min(filename) f FROM read_parquet('{dd}/*.parquet', filename = true) GROUP BY 1"))
dbDisconnect(con, shutdown = TRUE)
todo <- todo |> inner_join(fpm, by = "mdl_key")
todo$sql <- ifelse(todo$n_files == 1L,
sprintf("SELECT cell_id, max(val) AS val FROM read_parquet('%s') WHERE mdl_key = '%s' GROUP BY 1", todo$f, gsub("'", "''", todo$mdl_key)),
sprintf("SELECT cell_id, max(val) AS val FROM read_parquet('%s/*.parquet') WHERE mdl_key = '%s' GROUP BY 1", dd, gsub("'", "''", todo$mdl_key)))
todo$out <- todo$file
log_info("{ds}: painting {nrow(todo)} gridded COGs")
paint_rows(todo, grid)
}
k |> filter(in_cat(key) | file_exists(file))
})
if (nrow(vec_cog_reg)) add_plan(vec_cog_reg, "cog_vec")
}
```
## NMFS DPS extinction-risk surfaces → per-species COGs
Gridded from the start (`ingest_nmfs-dps.qmd`), so their one representation is the model grid itself.
Without these rows the species app lists the input struck through.
```{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)
dd_dps <- glue("{dir_atlas}/dist/dataset=dps_nmfs")
if (length(dps_keys) && dir_exists(dd_dps)) {
h <- cached_keys("cog_dps_nmfs", dd_dps, {
con <- new_con("2GB", "4"); on.exit(dbDisconnect(con, shutdown = TRUE))
pixel_hashes(con, glue("read_parquet('{dd_dps}/*.parquet')"), "mdl_key", quant = "trunc", dedup = "max")
}) |> filter(mdl_key %in% dps_keys)
dps_reg <- asset_store_key("cog_model", "global05", h$content_hash) |> mutate(mdl_key = h$mdl_key, ds_key = "dps_nmfs")
dps_reg$file <- store_path(dps_reg$key)
todo <- dps_reg |> filter(!in_cat(key), !file_exists(file))
if (nrow(todo)) {
todo$sql <- sprintf("SELECT cell_id, max(val) AS val FROM read_parquet('%s/*.parquet') WHERE mdl_key = '%s' GROUP BY 1", dd_dps, gsub("'", "''", todo$mdl_key))
todo$out <- todo$file
paint_rows(todo, grid, workers = 1)
}
dps_reg <- dps_reg |> filter(in_cat(key) | file_exists(file))
add_plan(dps_reg, "cog_dps_nmfs")
}
```
## AquaX: both representations from the COGs the ingest built
`ingest_aquax.qmd` (`AX_COG=1`) paints the two COGs per model locally and never uploads them (its
`AX_COG_S3` was retired with the store). `native` is the band **as delivered** (Float32, every modeled
pixel; keyed by its decoded pixels, which equal the source TIF's), `model` the surface **as ingested**
(÷10, rounded into INT1U; keyed by the `dist` rows rounded the same way). Both go to the store from here.
```{r}
#| label: ax-cogs
ax_reg <- tibble()
ax_keys <- if (length(served_by_ds[["ax"]])) sort(served_by_ds[["ax"]]) else character(0)
dd_ax <- glue("{dir_atlas}/dist/dataset=ax")
dir_axn <- glue("{dir_native}/ax_native"); dir_axm <- glue("{dir_native}/ax")
if (length(ax_keys) && dir_exists(dd_ax)) {
ax_take <- take_n(ax_keys)
# model: dist rows, rounded half-to-even then truncated (ingest_aquax rounds before publish_cog())
hm <- cached_keys("cog_ax", dd_ax, {
con <- new_con(); on.exit(dbDisconnect(con, shutdown = TRUE))
pixel_hashes(con, glue("read_parquet('{dd_ax}/*.parquet')"), "mdl_key", quant = "round_trunc", dedup = "max")
}) |> filter(mdl_key %in% ax_take)
km <- asset_store_key("cog_model", "global05", hm$content_hash) |> mutate(mdl_key = hm$mdl_key, ds_key = "ax", representation = "model")
km$file <- store_path(km$key)
todo <- km |> filter(!in_cat(key), !file_exists(file))
if (nrow(todo)) {
# paint from the rows (never copy the ingest's file blind: the key promises these pixels)
con <- new_con("2GB", "4")
fpm <- dbGetQuery(con, glue("SELECT mdl_key, min(filename) f FROM read_parquet('{dd_ax}/*.parquet', filename = true) GROUP BY 1"))
dbDisconnect(con, shutdown = TRUE)
todo <- todo |> inner_join(fpm, by = "mdl_key")
todo$sql <- sprintf("SELECT cell_id, max(val) AS val FROM read_parquet('%s') WHERE mdl_key = '%s' GROUP BY 1", todo$f, gsub("'", "''", todo$mdl_key))
todo$out <- todo$file
log_info("ax model COGs: painting {nrow(todo)}")
paint_rows(todo, grid, round_first = TRUE)
}
# native: the ingest's as-delivered COG, keyed by its decoded pixels (a bit-exact copy of the source band)
f_nat <- glue("{dir_axn}/{safe_key(ax_take)}.tif"); have <- file_exists(f_nat)
if (!all(have)) log_warn("ax native: {sum(!have)} served models have no as-delivered COG under {dir_axn} (run ingest_aquax with AX_COG=1)")
hn <- cached_keys("native_ax", dir_axn, {
fl <- f_nat[have]
plan(multisession, workers = n_workers); on.exit(plan(sequential), add = TRUE)
bind_rows(future_map(split(fl, ceiling(seq_along(fl) / 100)), function(b) {
suppressMessages({library(terra); library(msens)})
rows <- purrr::map_dfr(b, function(f) msens::native_raster_rows(terra::rast(f), grid) |> dplyr::mutate(file = f))
msens::content_hashes_df(rows, by = "file")
}, .options = furrr_options(packages = c("msens", "terra", "dplyr", "purrr"), seed = TRUE)))
}) |> mutate(mdl_key = ax_take[match(file, as.character(f_nat))]) |> filter(!is.na(mdl_key))
kn <- asset_store_key("native_ax", "ax", hn$content_hash) |> mutate(mdl_key = hn$mdl_key, ds_key = "ax", representation = "native", src = hn$file)
kn$file <- store_path(kn$key)
todo <- kn |> filter(!in_cat(key), !file_exists(file))
todo <- todo[!duplicated(todo$key), ]
if (nrow(todo)) file_copy(todo$src, todo$file)
ax_reg <- bind_rows(km, kn |> select(-src)) |> filter(in_cat(key) | file_exists(file))
add_plan(ax_reg, "cog_ax")
if (nrow(ax_reg) < 2 * length(ax_take))
log_warn("ax: {length(ax_take)} served models but {nrow(ax_reg)} pointer rows")
}
```
## Merged taxon surface → per-taxon whole-range COGs (opt-in)
The merged model itself (`ds_key = "ms_merge"`), painted from `dist_merged_global` (range-derived taxa)
and, for **suitability-only** taxa the global surface never holds (2+ AquaMaps models, or any AquaX
model), from the post-supersession merge input `mc_parts` as `round(max(val))` per cell. The am-only
taxa with ONE AquaMaps model and nothing else alias that model's own `model` object (registry chunk).
```{r}
#| label: merged-cogs
do_merged_cog <- nzchar(Sys.getenv("PUBLISH_MERGED_COG"))
redo_mcog <- nzchar(Sys.getenv("REDO_MERGED_COG"))
redo_keys_m <- trimws(strsplit(Sys.getenv("REDO_MERGED_COG_KEYS"), ",")[[1]]); redo_keys_m <- redo_keys_m[nzchar(redo_keys_m)]
if (nzchar(Sys.getenv("REDO_MERGED_COG_SPATIAL"))) {
conk <- new_con("2GB", "2"); dbExecute(conk, glue("ATTACH '{merge_db}' AS m (READ_ONLY)"))
redo_keys_m <- union(redo_keys_m, 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_m)} spatial-ER taxa")
}
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)) {
con <- new_con(); dbExecute(con, glue("ATTACH '{merge_db}' AS m (READ_ONLY)"))
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 <- take_n(sort(mkeys))
# the merge writes one taxon to one bucket dir (b1..b406; 0 of 17,131 straddle, 2026-10-07), so batches of
# buckets are batches of taxa; every taxon is hashed, served or not, so the cache serves any later set
h <- cached_keys("cog_merged", dir_gmerged, {
pixel_hashes_batched(glue("{sort(dir_ls(dir_gmerged, type = 'directory'))}/*.parquet"), n_batch = 40,
quant = "trunc", dedup = "max")
}) |> filter(mdl_key %in% mkeys)
merged_reg <- asset_store_key("cog_model", "global05", h$content_hash) |> mutate(mdl_key = h$mdl_key, ds_key = "ms_merge")
merged_reg$file <- store_path(merged_reg$key)
if (length(redo_keys_m)) { # a named repaint: drop the staged file so the key is re-examined
rm_f <- merged_reg$file[merged_reg$mdl_key %in% redo_keys_m & !in_cat(merged_reg$key)]
if (length(rm_f <- rm_f[file_exists(rm_f)])) file_delete(rm_f)
}
todo <- merged_reg |> filter(!in_cat(key), !file_exists(file))
merged_map <- tibble(mdl_key = mkeys, mmid = seq_along(mkeys))
if (nrow(todo)) {
if (!dir_exists(mparts) || redo_mcog || length(dir_ls(mparts)) != nrow(merged_map)) {
# mmid is POSITIONAL in the served set (the am `mid` lesson): rebuild the partitions whenever the set moved
if (dir_exists(mparts)) dir_delete(mparts)
dbExecute(con, "SET partitioned_write_max_open_files=400")
dbWriteTable(con, "merged_map", merged_map, overwrite = TRUE)
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")
}
todo <- todo |> inner_join(merged_map, by = "mdl_key")
todo$sql <- sprintf("SELECT cell_id, max(val) AS val FROM read_parquet('%s/mmid=%d/*.parquet') GROUP BY 1", mparts, todo$mmid)
todo$out <- todo$file
log_info("merged COGs: {nrow(merged_reg)} served, painting {nrow(todo)}")
paint_rows(todo, grid)
}
dbDisconnect(con, shutdown = TRUE)
merged_reg <- merged_reg |> filter(in_cat(key) | file_exists(file))
add_plan(merged_reg, "cog_merged")
}
```
```{r}
#| label: merged-suit
# suitability-only taxa: 2+ AquaMaps models, or any AquaX model -- absent from dist_merged_global (range-derived
# surfaces only). Their whole-range merged surface is max over the surviving models per cell of the
# post-supersession merge input (`mc_parts`, partitioned by mkey_id), rounded as the writer rounds.
mm_reg <- tibble()
dir_parts <- glue("{dir_atlas}/mc_parts")
if (do_merged_cog && dir_exists(dir_parts)) {
con <- new_con(); dbExecute(con, glue("ATTACH '{merge_db}' AS mrg (READ_ONLY)"))
suitonly <- dbGetQuery(con, "
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 ('am', 'ax'))
AND NOT (count(*) = 1 AND max(split_part(tm.mdl_key,'|',1)) = 'am')") |> as_tibble()
suitonly <- suitonly |> filter(!ms_merge_key %in% merged_reg$mdl_key) |> arrange(ms_merge_key)
if (test_n > 0) suitonly <- head(suitonly, test_n)
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)) {
# one hive partition per taxon, named by the suit-only set (explicit pruning: a JOIN over the 50 GB glob is
# not guaranteed to prune); hashed by mkey_id and mapped back to the merged key
h <- cached_keys("cog_merged_suit", dir_parts, {
hh <- pixel_hashes_batched(glue("{dir_parts}/mkey_id={sort(suitonly$mkey_id)}/*.parquet"), n_batch = 500, by = "mkey_id",
from = function(lst) glue("read_parquet([{lst}], hive_partitioning = true)"),
quant = "round_trunc", dedup = "max")
hh |> mutate(mkey_id = as.integer(mkey_id)) |>
inner_join(suitonly |> transmute(mkey_id = as.integer(mkey_id), mdl_key = ms_merge_key), by = "mkey_id") |>
select(mdl_key, n, content_hash)
}) |> filter(mdl_key %in% suitonly$ms_merge_key)
mm_reg <- asset_store_key("cog_model", "global05", h$content_hash) |> mutate(mdl_key = h$mdl_key, ds_key = "ms_merge")
mm_reg$file <- store_path(mm_reg$key)
todo <- mm_reg |> filter(!in_cat(key), !file_exists(file)) |> inner_join(suitonly |> select(mdl_key = ms_merge_key, mkey_id), by = "mdl_key")
if (nrow(todo)) {
todo$sql <- sprintf("SELECT cell_id, max(val) AS val FROM read_parquet('%s/mkey_id=%d/*.parquet') GROUP BY 1", dir_parts, todo$mkey_id)
todo$out <- todo$file
log_info("suitability-only merged COGs: painting {nrow(todo)}")
paint_rows(todo, grid, round_first = TRUE)
}
mm_reg <- mm_reg |> filter(in_cat(key) | file_exists(file))
add_plan(mm_reg, "cog_merged_suit")
}
dbDisconnect(con, shutdown = TRUE)
}
```
## Verify every new COG against its key
The key promises the pixels; before anything leaves this machine, each newly staged COG is decoded and
hashed the same way (`native_raster_rows()` + `content_hashes_df()`), and the result must equal the
key's payload hash. A mismatch is a bug in a painter, never something to upload.
```{r}
#| label: verify
plan_df <- bind_rows(plan)
stopifnot("nothing to publish: every class was skipped" = nrow(plan_df) > 0)
if (!"representation" %in% names(plan_df)) plan_df$representation <- NA_character_
new_cogs <- plan_df |> filter(!catalogued, endsWith(key, ".tif")) |> distinct(key, .keep_all = TRUE)
if (nrow(new_cogs) && do_verify) {
grids <- list(cog_model = grid,
native_am = grid_spec(rast(nrows = 360, ncols = 720, xmin = -180, xmax = 180, ymin = -90, ymax = 90, crs = "EPSG:4326")),
native_ax = grid)
fam_of <- function(key) ifelse(startsWith(key, "cog/"), "cog_model", ifelse(startsWith(key, "native/am/"), "native_am", "native_ax"))
new_cogs$family <- fam_of(new_cogs$key)
dec <- function(rows, grids) {
suppressMessages({library(terra); library(msens)})
rows$decoded <- vapply(seq_len(nrow(rows)), function(i) {
r <- msens::native_raster_rows(terra::rast(rows$file[i]), grids[[rows$family[i]]])
if (!nrow(r)) return(NA_character_)
msens::content_hashes_df(cbind(mdl_key = "x", r))$content_hash
}, "")
rows[c("key", "decoded")]
}
plan(multisession, workers = n_workers)
got <- bind_rows(future_map(split(new_cogs, (seq_len(nrow(new_cogs)) - 1) %% max(1, min(n_workers, nrow(new_cogs)))), dec, grids = grids,
.options = furrr_options(packages = c("msens", "terra"), seed = TRUE)))
plan(sequential)
new_cogs <- new_cogs |> left_join(got, by = "key") # by key, never by position: the split interleaves rows
bad <- new_cogs |> filter(is.na(decoded) | decoded != content_hash)
if (nrow(bad)) {
print(bad |> select(class, mdl_key, key, content_hash, decoded) |> head(20))
stop(sprintf("%d of %d new COGs do not decode to the pixels their key promises — not uploading", nrow(bad), nrow(new_cogs)), call. = FALSE)
}
log_info("verified {nrow(new_cogs)} new COGs: decoded pixels == key")
} else if (nrow(new_cogs)) log_warn("NATIVE_NO_VERIFY: {nrow(new_cogs)} new COGs uploaded without the decode check")
```
## Pointer rows: bbox, taxon key, aliases
```{r}
#| label: registry
rep_of <- c(cog_am = "model", native_am = "native", cog_vec = "model", cog_dps_nmfs = "model", cog_merged = "model", cog_merged_suit = "model")
rows <- plan_df |> mutate(
asset_type = if_else(endsWith(key, ".pmtiles"), "pmtiles", "cog"),
representation = coalesce(representation, if_else(startsWith(class, "pmtiles_"), "native", unname(rep_of[class]))),
asset_url = store_url(key),
content_hash = hash, # the key's 16 hex (what every published pointer table carries)
rescale_min = if_else(ds_key == "ax" & representation == "native", 0L, 1L),
rescale_max = if_else(ds_key == "ax" & representation == "native", 1000L, 100L),
colormap = "spectral_r",
source_layer = if_else(asset_type == "pmtiles", ds_key, NA_character_),
source_key = if_else(asset_type == "pmtiles", mdl_key, NA_character_)) # the key the tile's FEATURES carry
stopifnot(!anyNA(rows$representation))
# per-model geographic bbox so the app fits to each input's extent, not the world. am from _am_parts (mid ->
# mdl_key), vector/dps/ax from their dist partitions, merged from _merged_parts, suit-only from the painted raster
bbox_cols <- c("xmin","xmax","ymin","ymax")
bb_sql <- function(src, by) glue("
WITH rc AS (SELECT {by}, (cell_id-1) % 7200 AS col, (cell_id-1) // 7200 AS row FROM {src})
SELECT {by}, -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 {by}")
con <- new_con()
bb <- list()
am_from_parts <- nrow(am_reg) && exists("am_map_ok") && isTRUE(am_map_ok)
if (am_from_parts)
bb$am <- dbGetQuery(con, bb_sql(glue("read_parquet('{am_parts}/**/*.parquet', hive_partitioning=true)"), "mid")) |>
inner_join(am_map, by = "mid") |> select(-mid)
vec_ds <- unique(c(pmt_reg$ds_key, dps_reg$ds_key, ax_reg$ds_key, if (nrow(am_reg) && !am_from_parts) "am"))
if (length(vec_ds)) {
globs <- paste0("'", dir_atlas, "/dist/dataset=", vec_ds, "/*.parquet'", collapse = ", ")
dbWriteTable(con, "vk", tibble(mdl_key = unique(c(pmt_reg$mdl_key, dps_reg$mdl_key, ax_reg$mdl_key, if (!am_from_parts) am_reg$mdl_key))), overwrite = TRUE)
bb$vec <- dbGetQuery(con, bb_sql(glue("(SELECT p.mdl_key, p.cell_id FROM read_parquet([{globs}]) p SEMI JOIN vk USING(mdl_key))"), "mdl_key"))
}
if (nrow(merged_reg) && dir_exists(mparts))
bb$merged <- dbGetQuery(con, bb_sql(glue("read_parquet('{mparts}/**/*.parquet', hive_partitioning = true)"), "mmid")) |>
inner_join(merged_map, by = "mmid") |> select(-mmid)
if (nrow(mm_reg)) { # the painted window = the partition's cells (publish_cog crops to them)
dbWriteTable(con, "suit_bb", suitonly |> filter(ms_merge_key %in% mm_reg$mdl_key) |> select(mdl_key = ms_merge_key, mkey_id), overwrite = TRUE)
bb$mm <- dbGetQuery(con, bb_sql(glue("(SELECT k.mdl_key, p.cell_id FROM read_parquet('{dir_parts}/mkey_id=*/*.parquet', hive_partitioning = true) p JOIN suit_bb k USING (mkey_id))"), "mdl_key"))
}
dbDisconnect(con, shutdown = TRUE)
bb_all <- bind_rows(bb) |> distinct(mdl_key, .keep_all = TRUE)
rows <- rows |> left_join(bb_all, by = "mdl_key")
# an am-native (0.5 deg) row shares its model's bbox; same for ax native vs model
log_info("bbox: {sum(!is.na(rows$xmin))} of {nrow(rows)} pointer rows framed")
# attach ms_merge_key (the merged taxon key): a raw model may feed >1 merged taxon -> one row each
con <- dbConnect(duckdb(sdm_db, read_only = TRUE))
dbExecute(con, glue("ATTACH '{merge_db}' AS m (READ_ONLY)"))
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')")
inputs <- rows |> filter(ds_key != "ms_merge") |> inner_join(tm, by = "mdl_key", relationship = "many-to-many")
merged <- rows |> filter(ds_key == "ms_merge") |> mutate(ms_merge_key = mdl_key)
# am-only taxa: the merged surface IS the AquaMaps surface (verified byte-identical on a 150-taxon sample) ->
# alias the model object: same asset_url, same bbox, zero new bytes. The aliases are ms_merge rows, so they
# exist only when that class is built (PUBLISH_MERGED_COG): without it the class must stay ABSENT for
# registry_merge() to carry the published rows forward -- aliases alone read as a class that shrank to its
# am-only members, which is where the first v9 run stopped (2026-10-06)
have_merged <- merged$ms_merge_key
amonly <- if (do_merged_cog) dbGetQuery(con, "
SELECT ms_merge_key, max(mdl_key) AS am_key FROM (
SELECT t.ms_merge_key, tm.mdl_key FROM m.taxon t JOIN m.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'") else tibble(ms_merge_key = character(), am_key = character())
alias_rows <- amonly |> filter(!ms_merge_key %in% have_merged) |>
inner_join(rows |> filter(ds_key == "am", representation == "model") |> rename(am_key = mdl_key), by = "am_key") |>
mutate(mdl_key = ms_merge_key, ds_key = "ms_merge") |> select(-am_key)
if (nrow(alias_rows)) log_info("am-only merged aliases: {nrow(alias_rows)} taxa point at their AquaMaps model object (0 new objects)")
if (!do_merged_cog) log_warn("PUBLISH_MERGED_COG unset: the ms_merge class is carried forward from the published table, not rebuilt")
dbDisconnect(con, shutdown = TRUE)
native_asset <- bind_rows(inputs, merged, alias_rows) |>
select(any_of(msens:::.NATIVE_ASSET_COLS), key, catalogued, class) |>
relocate(ms_merge_key)
for (cc in setdiff(msens:::.NATIVE_ASSET_COLS, names(native_asset))) native_asset[[cc]] <- NA
native_asset <- native_asset |> relocate(all_of(msens:::.NATIVE_ASSET_COLS))
# the registry is CUMULATIVE: a run may not delete a class it did not rebuild. The PUBLISHED pointer table
# is the prior (what readers see today); the local sdm.duckdb copy only when the release is not yet published.
pub_na <- tryCatch(as_tibble(read_parquet(fetch(glue("{base}/{ver}/tables/native_asset.parquet")))), error = function(e) NULL)
prior_asset <- if (!is.null(pub_na)) { log_info("prior = published {ver}/tables/native_asset.parquet ({nrow(pub_na)} rows)"); pub_na } else {
c0 <- dbConnect(duckdb(sdm_db, read_only = TRUE))
p0 <- if (dbExistsTable(c0, "native_asset")) { log_info("prior = local sdm.duckdb native_asset (release not published)"); dbReadTable(c0, "native_asset") } else NULL
dbDisconnect(c0, shutdown = TRUE); p0
}
if (!is.null(prior_asset) && "mdl_seq" %in% names(prior_asset) && !"mdl_key" %in% names(prior_asset)) prior_asset <- NULL # a legacy-shaped table is not a prior for a v8+ build
native_asset <- registry_merge(native_asset |> select(all_of(msens:::.NATIVE_ASSET_COLS)), prior_asset,
allow_shrink = nzchar(Sys.getenv("NATIVE_REGISTRY_REBUILD")) || test_n > 0) # a test subset is never written
carried <- attr(native_asset, "carried")
if (length(carried)) log_warn("registry: carried forward {length(carried)} class(es) not built this run: {paste(carried, collapse = '; ')}")
native_asset <- as_tibble(native_asset)
smry <- native_asset |> count(representation, asset_type, name = "n_rows") |> arrange(representation, asset_type)
report_table(smry, caption = "pointer rows by representation × asset type")
```
## Upload → store, then catalogue
Only keys the catalog lacks are uploaded, each with its local store-mirror path; a key that already
exists on S3 but not in the catalog is an error (something wrote to the store outside it). After every
upload, the listing is re-read to prove each object is there, and only then are the rows appended to
`assets.parquet` and the catalog pushed.
```{r}
#| label: upload
up <- plan_df |> filter(!catalogued) |> distinct(key, .keep_all = TRUE)
up$file <- store_path(up$key)
stopifnot("every new key needs its staged file" = all(file_exists(up$file)))
log_info("upload plan: {nrow(up)} new objects, {round(sum(file_size(up$file)) / 1e9, 2)} GB ({sum(endsWith(up$key, '.tif'))} COGs, {sum(endsWith(up$key, '.pmtiles'))} PMTiles)")
if (nrow(up)) knitr::kable(up |> count(class, name = "new_objects"))
if (do_s3 && nrow(up)) {
hash_of <- function(k) sub("^.*/([0-9a-f]{16})\\.[a-z]+$", "\\1", k)
have_s3 <- c(cog_store_index("global05"), native_store_index())
clash <- up$key[hash_of(up$key) %in% have_s3]
if (length(clash)) stop(sprintf("%d object(s) exist on S3 but not in the catalog (e.g. %s): the store was written outside the catalog; reconcile first (scripts/verify_asset_store.R)",
length(clash), clash[1]), call. = FALSE)
t0 <- Sys.time()
for (i in seq_len(nrow(up))) {
ct <- if (endsWith(up$key[i], ".tif")) "image/tiff" else "binary/octet-stream"
out <- app_bundle_aws_retry(c("s3", "cp", shQuote(up$file[i]), shQuote(glue("{bucket_root}/{up$key[i]}")),
"--content-type", ct, "--only-show-errors"))
if (!is.null(attr(out, "status")) && attr(out, "status") != 0) stop(sprintf("upload failed for %s: %s", up$key[i], paste(out, collapse = " ")), call. = FALSE)
if (i %% 500 == 0) log_info("uploaded {i}/{nrow(up)} ({round(as.numeric(Sys.time() - t0, units = 'mins'), 1)} min)")
}
have_s3 <- c(cog_store_index("global05"), native_store_index())
missing <- up$key[!hash_of(up$key) %in% have_s3]
stopifnot("objects missing from the store after upload" = length(missing) == 0)
log_info("uploaded {nrow(up)} objects in {round(as.numeric(Sys.time() - t0, units = 'mins'), 1)} min; all present in the listing")
new_rows <- asset_catalog_rows(up$key, up$content_hash, up$enc, up$file, first_ver = ver)
catalog <- asset_catalog_add(catalog, new_rows)
f_cat <- glue("{dir_as}/assets.parquet")
asset_catalog_write(catalog, f_cat)
out <- app_bundle_aws_retry(c("s3", "cp", shQuote(f_cat), shQuote(glue("{bucket_root}/assets.parquet")),
"--content-type", "application/octet-stream", "--only-show-errors"))
if (!is.null(attr(out, "status")) && attr(out, "status") != 0) stop(paste(out, collapse = "\n"))
log_info("catalog: +{length(attr(catalog, 'added'))} rows -> {nrow(catalog)} objects, pushed assets.parquet")
catalog <- as_tibble(catalog)
} else if (nrow(up)) log_info("NATIVE_NO_S3: nothing uploaded; {nrow(up)} objects staged under {dir_store}")
```
## Gate, then write the pointer table
```{r}
#| label: gate
url_gate(native_asset$asset_url, glue("{ver} native_asset (asset_url)"))
listed <- if (do_s3) {
o <- suppressWarnings(system2("aws", c("s3", "ls", "--recursive", shQuote(glue("{bucket_root}/{ver}/"))), stdout = TRUE, stderr = TRUE))
sub("^.*\\s+", "", o[nzchar(o)])
} else NULL
cat_gate <- if (do_s3 || !nrow(up)) catalog else # dry run: the catalog this run would have written
as_tibble(asset_catalog_add(catalog, asset_catalog_rows(up$key, up$content_hash, up$enc, up$file, first_ver = ver)))
g <- withCallingHandlers(
store_publish_gate(native_asset$asset_url, cat_gate, ver, listed = listed, legacy_ok = legacy_ok),
warning = function(w) { log_warn(conditionMessage(w)); invokeRestart("muffleWarning") })
log_info("store gate: {g$n_urls} pointer URLs -> {g$n_keys} catalogued store keys; bulk files under {ver}/: {length(g$legacy)}")
if (do_s3) {
con <- dbConnect(duckdb(sdm_db))
dbWriteTable(con, "native_asset", as.data.frame(native_asset), overwrite = TRUE)
dbDisconnect(con, shutdown = TRUE)
log_info("wrote native_asset ({nrow(native_asset)} rows) to {path_file(sdm_db)}; release_marine-atlas.qmd (RELEASE_S3_TABLES=1) publishes {ver}/tables/native_asset.parquet")
} else log_info("NATIVE_NO_S3: pointer table NOT written ({nrow(native_asset)} rows would be; {sum(!plan_df$catalogued)} of their objects are not yet catalogued)")
```
## Manifest
```{r}
#| label: manifest
if (do_s3) {
con <- dbConnect(duckdb(sdm_db, read_only = TRUE))
h <- hash_query(con, "native_asset")
dbDisconnect(con, shutdown = TRUE)
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),
am_cogs_model = nrow(am_reg),
am_cogs_native = nrow(am_native_reg),
vec_cogs_model = nrow(vec_cog_reg),
new_objects = nrow(up),
native_assets = nrow(native_asset)),
force = force_target("publish_native"))
} else log_info("NATIVE_NO_S3: manifest not written (dry run)")
```