Ingest AquaMaps → global 0.05° cells (bilinear-weight interpolation)

Published

2026-07-13 18:54:07

Reingest all AquaMaps species (native HCAF 0.5°) onto the global 0.05° cell grid ([-180,180], from build_cell_grid.qmd), writing one Parquet per species to the marine-atlas release keyed by mdl_key = "am|{sp_key}".

Method — bilinear interpolation as a precomputed weight join (fast + edge-safe). The 0.5°→0.05° bilinear resample has fixed geometric weights (independent of species), so we build them once as w05(cell_id, loiczid, w): each 0.05° ocean cell → its 4 surrounding HCAF 0.5° grid cells (loiczid) + bilinear weights (which sum to 1). Then each species is interpolated by an indexed DuckDB join — no terra, no per-species raster I/O — at ~0.1 s/species (~50× faster than the per-species terra resample on shared storage).

Edges. Absent HCAF cells simply aren’t in spp_cells, so they contribute 0 to the weighted sum — the surface fades to 0 within ~one cell of the distribution edge (the v7 “surround by zeros” behaviour). Validated cell-for-cell against the terra bilinear result: values cor > 0.999 (identical where both defined), edge fade visually indistinguishable, no bleed.

Re-run control (libs/vars.R, override via env): REDO_AM_W05=1 rebuilds the weight table; REDO_AM_INGEST=1 re-interpolates every species (otherwise the loop resumes, skipping species whose Parquet already exists).

1 Design

Figure 1: AquaMaps HCAF 0.5° suitability bilinearly interpolated onto the global 0.05° grid via a precomputed weight join

2 Setup

Code

librarian::shelf(
  DBI, dplyr, duckdb, fs, glue, here, jsonlite, logger, readr, terra, tibble,
  quiet = T)
source(here("libs/paths.R"))   # dir_derived, dir_big_v, cellid_tif, am_db, w05_db, ver, s3_atlas
source(here("libs/vars.R"))    # redo_am_w05, redo_am_ingest
options(readr.show_col_types = F)

ds_key    <- "am"
threshold <- 1                                                # drop interpolated suitability < 1
dir_dist  <- glue("{dir_big_v}/marine-atlas/dist/dataset={ds_key}")
manifest  <- here("data/manifests/ingest_aquamaps.json")
dir_create(c(dir_dist, path_dir(manifest)))

stopifnot(
  "run build_cell_grid.qmd first (need cellid_tif)" = file_exists(cellid_tif),
  "AquaMaps source am.duckdb not found"             = file_exists(am_db))

3 Bilinear weight table `w05` (built once)

w05 has 4 rows per 0.05° ocean cell — its bilinear stencil onto the HCAF grid. loiczid = row*720 + col + 1 is the HCAF grid position (verified == cells.loiczid). ~68M rows, ~30 s, indexed on loiczid so the per-species join is a fast lookup.

Code

if (redo_am_w05 || !file_exists(w05_db)) {
  log_info("building AquaMaps bilinear weight table w05 (once) ...")
  if (file_exists(w05_db)) file_delete(w05_db)
  cw <- dbConnect(duckdb(w05_db))
  # AquaMaps is marine: restrict weights to OCEAN cells (cellid_tif is the GLOBAL
  # grid — land + ocean — so mask by the Bio-Oracle depth layer to keep ocean only)
  oc <- as.data.frame(
    mask(rast(cellid_tif), rast(bio_oracle_tif, lyrs = "depth_mean")),
    xy = TRUE, na.rm = TRUE)
  names(oc) <- c("lon", "lat", "cell_id")
  duckdb_register(cw, "cell_xy", oc)
  dbExecute(cw, "CREATE TABLE w05 AS
    WITH f AS (
      SELECT CAST(cell_id AS INTEGER)           AS cell_id,
             floor((lon+179.75)/0.5)            AS cl,
             floor((89.75-lat)/0.5)             AS rt,
             (lon-(-179.75+0.5*floor((lon+179.75)/0.5)))/0.5 AS fx,
             ((89.75-0.5*floor((89.75-lat)/0.5))-lat)/0.5    AS fy
      FROM cell_xy)
    SELECT cell_id, CAST(rt*720+cl+1     AS BIGINT) AS loiczid, (1-fx)*(1-fy) AS w FROM f WHERE (1-fx)*(1-fy) > 0
    UNION ALL SELECT cell_id, CAST(rt*720+cl+2     AS BIGINT), fx*(1-fy) FROM f WHERE fx*(1-fy) > 0
    UNION ALL SELECT cell_id, CAST((rt+1)*720+cl+1 AS BIGINT), (1-fx)*fy FROM f WHERE (1-fx)*fy > 0
    UNION ALL SELECT cell_id, CAST((rt+1)*720+cl+2 AS BIGINT), fx*fy     FROM f WHERE fx*fy     > 0")
  dbExecute(cw, "CREATE INDEX w05_loiczid ON w05(loiczid)")
  n_w <- dbGetQuery(cw, "SELECT count(*) n FROM w05")$n
  dbDisconnect(cw, shutdown = TRUE)
  log_info("w05 built: {format(n_w, big.mark=',')} weight rows")
}

4 Species list + model crosswalk

Code

con_am <- dbConnect(duckdb(am_db, read_only = TRUE))
sp_all <- tbl(con_am, "spp") |> select(sp_key, genus, species, iucn_id) |> collect()
dbDisconnect(con_am, shutdown = TRUE)
sp_all <- sp_all |>
  mutate(mdl_key = glue("{ds_key}|{sp_key}"), taxa = glue("{genus} {species}"))
write_csv(sp_all |> select(mdl_key, sp_key, taxa, iucn_id),
          glue("{dir_dist}/../model_{ds_key}.csv"))
nrow(sp_all)

[1] 23699

5 Interpolate every species → one Parquet each

Per species: join its present HCAF cells (spp_cells → cells → loiczid) to w05, sum w × probability×100 per 0.05° cell, keep ≥ threshold. DuckDB COPY writes the Parquet directly. Resumable; REDO_AM_INGEST=1 starts fresh.

Code

if (redo_am_ingest && dir_exists(dir_dist)) { dir_delete(dir_dist); dir_create(dir_dist) }
done <- path_ext_remove(path_file(dir_ls(dir_dist, glob = "*.parquet")))
todo <- sp_all$sp_key[!sp_all$sp_key %in% done]
log_info("AquaMaps: {length(done)} done; interpolating {length(todo)} of {nrow(sp_all)} species ...")

con <- dbConnect(duckdb(w05_db, read_only = TRUE))
dbExecute(con, glue("ATTACH '{am_db}' AS am (READ_ONLY)"))

[1] 0

Code

dbExecute(con, glue("SET threads = {max(1L, parallel::detectCores())}"))

[1] 0

Code

t0 <- Sys.time()
for (i in seq_along(todo)) {
  k <- todo[i]
  msens::copy_atlas_parquet(con, glue("
    SELECT 'am|{k}' AS mdl_key, w.cell_id, round(SUM(w.w * sc.prob), 2) AS \"val\"
    FROM w05 w
    JOIN (SELECT c.loiczid, s.probability*100 AS prob
            FROM am.spp_cells s JOIN am.cells c USING (cell_id)
           WHERE s.sp_key = '{k}') sc  ON w.loiczid = sc.loiczid
    GROUP BY w.cell_id
    HAVING SUM(w.w * sc.prob) >= {threshold}"),
    fs::path(dir_dist, k, ext = "parquet"))
  if (i %% 1000 == 0)
    log_info("  {i}/{length(todo)} ({round(as.numeric(Sys.time()-t0, units='mins'),1)} min)")
}
dbDisconnect(con, shutdown = TRUE)
log_info("interpolation done in {round(as.numeric(Sys.time()-t0, units='mins'),1)} min")

6 Verify

Code

n_pq <- length(dir_ls(dir_dist, glob = "*.parquet"))
con  <- dbConnect(duckdb())
smp  <- dbGetQuery(con, glue(
  "SELECT count(*) n_cells, min(val) v_min, max(val) v_max
     FROM read_parquet('{fs::path(dir_dist, sp_all$sp_key[1], ext=\"parquet\")}')"))
dbDisconnect(con, shutdown = TRUE)
glue("{n_pq} species Parquet written; sample {sp_all$sp_key[1]}: ",
     "{smp$n_cells} cells, val {smp$v_min}-{smp$v_max}")

23699 species Parquet written; sample Chn-013ca97e-46a8-4c5e-b398-8bb4f340f88b: 512073 cells, val 1-100

7 Outputs + manifest

Code

pq   <- dir_ls(dir_dist, glob = "*.parquet")
con  <- dbConnect(duckdb())
smry <- dbGetQuery(con, glue(
  "SELECT count(DISTINCT mdl_key) n_models, count(*) n_cells, min(val) v_min, max(val) v_max
     FROM read_parquet('{dir_dist}/*.parquet')"))
# order-independent content fingerprint of the on-disk surface (not the ingest run)
h <- msens::hash_parquet(glue("{dir_dist}/*.parquet"), con)
dbDisconnect(con, shutdown = TRUE)
msens::report_table(smry, caption = "am model_cell surface")

am model_cell surface
n_models	n_cells	v_min	v_max
23699	12982232398	1	100

Code

# content-addressed manifest: deterministic, no wall-clock, no machine paths ->
# downstream targets re-run only when this surface's content actually changes
msens::write_manifest(
  manifest, target = "ingest_aquamaps", content_hash = h,
  stats = list(ds_key = ds_key, ver = ver, n_species = nrow(sp_all),
               n_models = smry$n_models, n_cells = smry$n_cells,
               v_min = smry$v_min, v_max = smry$v_max, n_parquet = length(pq)),
  force = msens::force_target("ingest_aquamaps"))

8 Publish to S3

Code

# aws s3 sync {dir_big_v}/marine-atlas s3://oceanmetrics.io-public/marine-atlas/v8
system(glue("aws s3 sync {dir_big_v}/marine-atlas {s3_atlas}/{ver}"))

--- title: "Ingest AquaMaps → global 0.05° cells (bilinear-weight interpolation)" msens: target_name: ingest_aquamaps workflow_type: ingest dependency: [build_cell_grid] output: data/manifests/ingest_aquamaps.json dataset: {ds_key: am, response_type: suitability, source_authority: AquaMaps, temporal_interval: static, native_format: raster} editor_options: chunk_output_type: console --- Reingest **all** AquaMaps species (native HCAF 0.5°) onto the **global 0.05° cell grid** ([-180,180], from `build_cell_grid.qmd`), writing one Parquet per species to the marine-atlas release keyed by `mdl_key = "am|{sp_key}"`. **Method — bilinear interpolation as a precomputed weight join (fast + edge-safe).** The 0.5°→0.05° bilinear resample has *fixed* geometric weights (independent of species), so we build them **once** as `w05(cell_id, loiczid, w)`: each 0.05° ocean cell → its 4 surrounding HCAF 0.5° grid cells (`loiczid`) + bilinear weights (which sum to 1). Then each species is interpolated by an indexed DuckDB join — no terra, no per-species raster I/O — at ~0.1 s/species (~50× faster than the per-species terra resample on shared storage). **Edges.** Absent HCAF cells simply aren't in `spp_cells`, so they contribute 0 to the weighted sum — the surface fades to 0 within ~one cell of the distribution edge (the v7 "surround by zeros" behaviour). Validated cell-for-cell against the terra bilinear result: values `cor > 0.999` (identical where both defined), edge fade visually indistinguishable, no bleed. Re-run control (`libs/vars.R`, override via env): `REDO_AM_W05=1` rebuilds the weight table; `REDO_AM_INGEST=1` re-interpolates every species (otherwise the loop resumes, skipping species whose Parquet already exists). ## Design ```{mermaid} %%| label: fig-design %%| fig-cap: "AquaMaps HCAF 0.5° suitability bilinearly interpolated onto the global 0.05° grid via a precomputed weight join" flowchart LR am["am.duckdb<br/>HCAF 0.5° spp_cells"] --> w05["w05 weight table<br/>(cell_id → 4 loiczid, bilinear w)"] w05 --> interp["per-species DuckDB join<br/>Σ w×prob×100, ≥ threshold"] interp --> pq[("dist/dataset=am/{sp_key}.parquet<br/>(mdl_key, cell_id, val) · V2/zstd")] pq --> mf["hash_parquet → content-addressed manifest"] ``` ## Setup ```{r} #| label: setup librarian::shelf( DBI, dplyr, duckdb, fs, glue, here, jsonlite, logger, readr, terra, tibble, quiet = T) source(here("libs/paths.R")) # dir_derived, dir_big_v, cellid_tif, am_db, w05_db, ver, s3_atlas source(here("libs/vars.R")) # redo_am_w05, redo_am_ingest options(readr.show_col_types = F) ds_key <- "am" threshold <- 1 # drop interpolated suitability < 1 dir_dist <- glue("{dir_big_v}/marine-atlas/dist/dataset={ds_key}") manifest <- here("data/manifests/ingest_aquamaps.json") dir_create(c(dir_dist, path_dir(manifest))) stopifnot( "run build_cell_grid.qmd first (need cellid_tif)" = file_exists(cellid_tif), "AquaMaps source am.duckdb not found" = file_exists(am_db)) ``` ## Bilinear weight table `w05` (built once) `w05` has 4 rows per 0.05° ocean cell — its bilinear stencil onto the HCAF grid. `loiczid = row*720 + col + 1` is the HCAF grid position (verified `== cells.loiczid`). ~68M rows, ~30 s, indexed on `loiczid` so the per-species join is a fast lookup. ```{r} #| label: build_w05 if (redo_am_w05 || !file_exists(w05_db)) { log_info("building AquaMaps bilinear weight table w05 (once) ...") if (file_exists(w05_db)) file_delete(w05_db) cw <- dbConnect(duckdb(w05_db)) # AquaMaps is marine: restrict weights to OCEAN cells (cellid_tif is the GLOBAL # grid — land + ocean — so mask by the Bio-Oracle depth layer to keep ocean only) oc <- as.data.frame( mask(rast(cellid_tif), rast(bio_oracle_tif, lyrs = "depth_mean")), xy = TRUE, na.rm = TRUE) names(oc) <- c("lon", "lat", "cell_id") duckdb_register(cw, "cell_xy", oc) dbExecute(cw, "CREATE TABLE w05 AS WITH f AS ( SELECT CAST(cell_id AS INTEGER) AS cell_id, floor((lon+179.75)/0.5) AS cl, floor((89.75-lat)/0.5) AS rt, (lon-(-179.75+0.5*floor((lon+179.75)/0.5)))/0.5 AS fx, ((89.75-0.5*floor((89.75-lat)/0.5))-lat)/0.5 AS fy FROM cell_xy) SELECT cell_id, CAST(rt*720+cl+1 AS BIGINT) AS loiczid, (1-fx)*(1-fy) AS w FROM f WHERE (1-fx)*(1-fy) > 0 UNION ALL SELECT cell_id, CAST(rt*720+cl+2 AS BIGINT), fx*(1-fy) FROM f WHERE fx*(1-fy) > 0 UNION ALL SELECT cell_id, CAST((rt+1)*720+cl+1 AS BIGINT), (1-fx)*fy FROM f WHERE (1-fx)*fy > 0 UNION ALL SELECT cell_id, CAST((rt+1)*720+cl+2 AS BIGINT), fx*fy FROM f WHERE fx*fy > 0") dbExecute(cw, "CREATE INDEX w05_loiczid ON w05(loiczid)") n_w <- dbGetQuery(cw, "SELECT count(*) n FROM w05")$n dbDisconnect(cw, shutdown = TRUE) log_info("w05 built: {format(n_w, big.mark=',')} weight rows") } ``` ## Species list + model crosswalk ```{r} #| label: model_csv con_am <- dbConnect(duckdb(am_db, read_only = TRUE)) sp_all <- tbl(con_am, "spp") |> select(sp_key, genus, species, iucn_id) |> collect() dbDisconnect(con_am, shutdown = TRUE) sp_all <- sp_all |> mutate(mdl_key = glue("{ds_key}|{sp_key}"), taxa = glue("{genus} {species}")) write_csv(sp_all |> select(mdl_key, sp_key, taxa, iucn_id), glue("{dir_dist}/../model_{ds_key}.csv")) nrow(sp_all) ``` ## Interpolate every species → one Parquet each Per species: join its present HCAF cells (`spp_cells → cells → loiczid`) to `w05`, sum `w × probability×100` per 0.05° cell, keep `≥ threshold`. DuckDB `COPY` writes the Parquet directly. Resumable; `REDO_AM_INGEST=1` starts fresh. ```{r} #| label: interpolate if (redo_am_ingest && dir_exists(dir_dist)) { dir_delete(dir_dist); dir_create(dir_dist) } done <- path_ext_remove(path_file(dir_ls(dir_dist, glob = "*.parquet"))) todo <- sp_all$sp_key[!sp_all$sp_key %in% done] log_info("AquaMaps: {length(done)} done; interpolating {length(todo)} of {nrow(sp_all)} species ...") con <- dbConnect(duckdb(w05_db, read_only = TRUE)) dbExecute(con, glue("ATTACH '{am_db}' AS am (READ_ONLY)")) dbExecute(con, glue("SET threads = {max(1L, parallel::detectCores())}")) t0 <- Sys.time() for (i in seq_along(todo)) { k <- todo[i] msens::copy_atlas_parquet(con, glue(" SELECT 'am|{k}' AS mdl_key, w.cell_id, round(SUM(w.w * sc.prob), 2) AS \"val\" FROM w05 w JOIN (SELECT c.loiczid, s.probability*100 AS prob FROM am.spp_cells s JOIN am.cells c USING (cell_id) WHERE s.sp_key = '{k}') sc ON w.loiczid = sc.loiczid GROUP BY w.cell_id HAVING SUM(w.w * sc.prob) >= {threshold}"), fs::path(dir_dist, k, ext = "parquet")) if (i %% 1000 == 0) log_info(" {i}/{length(todo)} ({round(as.numeric(Sys.time()-t0, units='mins'),1)} min)") } dbDisconnect(con, shutdown = TRUE) log_info("interpolation done in {round(as.numeric(Sys.time()-t0, units='mins'),1)} min") ``` ## Verify ```{r} #| label: verify n_pq <- length(dir_ls(dir_dist, glob = "*.parquet")) con <- dbConnect(duckdb()) smp <- dbGetQuery(con, glue( "SELECT count(*) n_cells, min(val) v_min, max(val) v_max FROM read_parquet('{fs::path(dir_dist, sp_all$sp_key[1], ext=\"parquet\")}')")) dbDisconnect(con, shutdown = TRUE) glue("{n_pq} species Parquet written; sample {sp_all$sp_key[1]}: ", "{smp$n_cells} cells, val {smp$v_min}-{smp$v_max}") ``` ## Outputs + manifest ```{r} #| label: manifest pq <- dir_ls(dir_dist, glob = "*.parquet") con <- dbConnect(duckdb()) smry <- dbGetQuery(con, glue( "SELECT count(DISTINCT mdl_key) n_models, count(*) n_cells, min(val) v_min, max(val) v_max FROM read_parquet('{dir_dist}/*.parquet')")) # order-independent content fingerprint of the on-disk surface (not the ingest run) h <- msens::hash_parquet(glue("{dir_dist}/*.parquet"), con) dbDisconnect(con, shutdown = TRUE) msens::report_table(smry, caption = "am model_cell surface") # content-addressed manifest: deterministic, no wall-clock, no machine paths -> # downstream targets re-run only when this surface's content actually changes msens::write_manifest( manifest, target = "ingest_aquamaps", content_hash = h, stats = list(ds_key = ds_key, ver = ver, n_species = nrow(sp_all), n_models = smry$n_models, n_cells = smry$n_cells, v_min = smry$v_min, v_max = smry$v_max, n_parquet = length(pq)), force = msens::force_target("ingest_aquamaps")) ``` ## Publish to S3 ```{r} #| label: publish #| eval: false # aws s3 sync {dir_big_v}/marine-atlas s3://oceanmetrics.io-public/marine-atlas/v8 system(glue("aws s3 sync {dir_big_v}/marine-atlas {s3_atlas}/{ver}")) ```

1 Design

2 Setup

3 Bilinear weight table w05 (built once)

4 Species list + model crosswalk

5 Interpolate every species → one Parquet each

6 Verify

7 Outputs + manifest

8 Publish to S3

3 Bilinear weight table `w05` (built once)