Build global 0.05° cell grid (v8 sampling unit)

Global ocean 0.05° raster cells in [-180,180] with env covariates + US study-area / Program-Area membership

Published

2026-08-27

The v8 spatial sampling unit is the global 0.05° raster cell in [-180,180] (rolled back from H3 hex). The grid is anchored to the shared 0.05°, 3600×7200, EPSG:4326 topology of Bio-Oracle v3 AND the AquaX SDM TIFs — so cell_id is the raster cell index (1:ncell), and any SDM raster on this grid (Bio-Oracle, AquaX, resampled AquaMaps) maps to cell_id by position with no projection.

This notebook builds the canonical cell table (ocean cells only) — cell_id, the 12 Bio-Oracle env covariates, area_km2, and boolean in_usa / in_pra membership — plus the global cell-id COG that titiler paints (and that ingest_aquamaps.qmd uses to build its bilinear weight table).

The global build (631 MB Bio-Oracle → 25.9M cells → ~17M-row cell table + 47 MB COG) is skipped when the cell-id COG already exists (resume); force it with REDO_CELL_GRID=1 (libs/vars.R). A bbox prototype run (CELL_BBOX=...) always rebuilds into its own CELL_DB.

1 Design

Code

flowchart LR
  bo["Bio-Oracle v3<br/>0.05° 12-layer raster"] --> id["assign cell_id=1:ncell,<br/>area_km2, in_usa/in_pra"]
  id --> ocn["keep ocean cells<br/>(depth_mean non-NA)"]
  ocn --> cell[("cell table<br/>sdm.duckdb")]
  ocn --> cog["r_cellid_global.tif COG<br/>(titiler paints)"]
  cell --> mf["hash_query(cell)<br/>→ content-addressed manifest"]

flowchart LR
  bo["Bio-Oracle v3<br/>0.05° 12-layer raster"] --> id["assign cell_id=1:ncell,<br/>area_km2, in_usa/in_pra"]
  id --> ocn["keep ocean cells<br/>(depth_mean non-NA)"]
  ocn --> cell[("cell table<br/>sdm.duckdb")]
  ocn --> cog["r_cellid_global.tif COG<br/>(titiler paints)"]
  cell --> mf["hash_query(cell)<br/>→ content-addressed manifest"]

Figure 1: Bio-Oracle raster → global 0.05° cell grid: cell table + cell-id COG, content-addressed manifest

Code

librarian::shelf(DBI, dplyr, duckdb, glue, here, jsonlite, sf, terra, quiet = TRUE)
source(here("libs/paths.R"))   # ver, sdm_db, bio_oracle_tif, cellid_tif, ply_usa_gpkg, ply_pra_gpkg
source(here("libs/vars.R"))    # redo_cell_grid
stopifnot(file.exists(bio_oracle_tif), file.exists(ply_usa_gpkg), file.exists(ply_pra_gpkg))

# optional bbox restriction for prototype runs: CELL_BBOX="xmin,xmax,ymin,ymax"
cell_bbox <- Sys.getenv("CELL_BBOX", "")
cell_db   <- Sys.getenv("CELL_DB", sdm_db)
dir.create(dirname(cell_db), recursive = TRUE, showWarnings = FALSE)
dir.create(here("data/manifests"), recursive = TRUE, showWarnings = FALSE)

# skip the expensive global build if the cell-id COG already exists (resume),
# unless forced; a bbox prototype always rebuilds into its own CELL_DB.
#
# The cell-id COG is SHARED across versions, but the `cell` TABLE lives in this
# version's sdm.duckdb. A fresh version directory therefore has the COG and no
# table -- and every downstream `stopifnot("run build_cell_grid first" =
# file_exists(sdm_db))` passes on an empty database. Resume path for that case:
# copy `cell` from ver_prev's build DB (seconds) and assert it hashes to the
# checkpoint; rebuild only when there is nothing to copy.
has_cell <- function(db) {
  db <- path.expand(db)
  if (!file.exists(db)) return(FALSE)
  c0 <- dbConnect(duckdb(db, read_only = TRUE)); on.exit(dbDisconnect(c0, shutdown = TRUE))
  "cell" %in% dbListTables(c0)
}
prev_db    <- glue("{dir_big}/{ver_prev}/sdm.duckdb")
need_build <- redo_cell_grid || nzchar(cell_bbox) || !file.exists(cellid_tif) ||
              (!has_cell(cell_db) && !has_cell(prev_db))
copy_cell  <- !need_build && !has_cell(cell_db) && has_cell(prev_db)

2 Read Bio-Oracle, assign global cell_id, membership + area

Code

if (need_build) {
  r <- rast(bio_oracle_tif)   # global [-180,180], 3600 x 7200, 12 env layers
  r$cell_id  <- 1:ncell(r)    # global raster-position id (shared with AquaX/AquaMaps grids)
  r$area_km2 <- cellSize(r$depth_mean, unit = "km")

  # membership: rasterize the study-area + Program-Area polygons onto the same grid
  r$in_usa <- rasterize(vect(st_read(ply_usa_gpkg, quiet = TRUE)), r$depth_mean, field = 1, background = 0)
  r$in_pra <- rasterize(vect(st_read(ply_pra_gpkg, quiet = TRUE)), r$depth_mean, field = 1, background = 0)

  if (nzchar(cell_bbox)) {
    bb <- as.numeric(strsplit(cell_bbox, ",")[[1]])
    r  <- crop(r, ext(bb[1], bb[2], bb[3], bb[4]))
    message(glue("restricted to bbox [{cell_bbox}]"))
  }
}

3 Write cell table (ocean cells) + global cell-id COG

Code

if (need_build) {
  # ocean cells only (depth_mean non-NA); cell_id stays the global raster index
  d_cell <- as.data.frame(r, xy = TRUE, na.rm = FALSE) |>
    filter(!is.na(depth_mean)) |>
    transmute(
      cell_id = as.integer(cell_id),
      lon = x, lat = y,
      across(c(depth_mean, depth_min, depth_max, oxy_b_mean, oxy_mean, prim_prod_mean,
               ice_con_ann, salinity_b_mean, salinity_mean, sbt_an_mean, sst_an_mean, fao_area_m)),
      area_km2,
      in_usa = in_usa == 1,
      in_pra = in_pra == 1)

  con <- dbConnect(duckdb(cell_db))
  dbWriteTable(con, "cell", d_cell, overwrite = TRUE)
  dbExecute(con, "ALTER TABLE cell ADD PRIMARY KEY (cell_id)")
  dbDisconnect(con, shutdown = TRUE)

  # global cell-id COG in [-180,180] for titiler (uint32; no overviews to keep ids intact)
  writeRaster(r$cell_id, cellid_tif, filetype = "COG", datatype = "INT4U", overwrite = TRUE,
              gdal = c("COMPRESS=LZW", "OVERVIEWS=NONE"))
} else if (copy_cell) {
  # same grid, new version: the table is a straight copy of ver_prev's, verified
  # against the content-addressed checkpoint rather than trusted
  con <- dbConnect(duckdb(cell_db))
  dbExecute(con, glue("ATTACH '{path.expand(prev_db)}' AS prev (READ_ONLY)"))
  dbExecute(con, "CREATE TABLE cell AS SELECT * FROM prev.cell")
  dbExecute(con, "ALTER TABLE cell ADD PRIMARY KEY (cell_id)")
  dbExecute(con, "DETACH prev")
  h_copy <- msens::hash_query(con, "cell")
  dbDisconnect(con, shutdown = TRUE)
  mf <- here("data/manifests/build_cell_grid.json")
  if (file.exists(mf)) {
    h_ckpt <- jsonlite::fromJSON(mf)$content_hash
    stopifnot("copied cell table does not hash to the build_cell_grid checkpoint" = identical(h_copy, h_ckpt))
  }
  message(glue("cell table copied from {ver_prev} ({h_copy}) -> {cell_db}"))
} else {
  message(glue("cell grid up to date (cellid_tif + cell table exist); set REDO_CELL_GRID=1 to rebuild"))
}

4 Publish the cell-id COGs (gated)

The apps resolve a clicked point through titiler’s /cog/point/{lon},{lat}: one call for the layer’s value, one for the cell id. That needs the cell-id lookups on S3, not just on the server’s disk.

Cell ids are categorical, so these are written INT4U with no overviews — a resampled pyramid would average neighbouring ids into ids that do not exist. Both grids are published: usa05 serves v1–v7, global05 serves v8.

Resolving a click by reading the raster in the app was the source of three separate bugs: the band is named r_cellid on one grid and depth_mean on the other (so selecting $cell_id returned NULL), the app shifted longitude to 0–360 while both rasters are stored −180..180 (so every Americas click read outside the image), and a SpatRaster cached across Shiny sessions is a stale external pointer that segfaults the process. One HTTP call has none of those modes.

Code

if (nzchar(Sys.getenv("GRID_COG_S3"))) {
  grids <- msens::grid_registry()
  for (i in seq_len(nrow(grids))) {
    gid <- grids$grid_id[i]
    src <- path.expand(glue("{dir_derived}/{grids$cellid_tif[i]}"))
    if (!file.exists(src)) { message(glue("  {gid}: {basename(src)} absent - skip")); next }
    dst <- glue("{s3_atlas}/grid/{gid}/cellid.tif")
    out <- system2("aws", c("s3", "cp", shQuote(src), shQuote(dst),
                            "--only-show-errors", "--no-progress"),
                   stdout = TRUE, stderr = TRUE)
    if (!is.null(attr(out, "status")) && attr(out, "status") != 0)
      stop("aws s3 cp failed: ", paste(out, collapse = "\n"))
    message(glue("  {gid}: {round(file.size(src)/1024^2, 1)} MB -> {dst}"))
  }
} else {
  message("GRID_COG_S3 unset - not publishing cell-id COGs")
}

5 Summary + manifest

Code

# The cell grid is built on the laptop; the SERVER carries only serve.duckdb and
# the published Parquet, so `cell_db` does not exist there. Summarise whichever
# is present -- otherwise this notebook cannot be rendered on the server at all,
# which is exactly where the cell-id COGs above have to be published from.
cell_src <- if (file.exists(cell_db)) {
  list(kind = "db", from = glue("'{cell_db}'"))
} else {
  # tables/ live under big/{ver}; on the server dir_big_v collapses to
  # derived/{ver}, which is a different tree
  pq <- path.expand(glue("{dir_big}/{ver}/tables/cell.parquet"))
  if (file.exists(pq)) list(kind = "parquet", from = glue("read_parquet('{pq}')"))
  else NULL
}
stopifnot("no cell table: neither the grid database nor tables/cell.parquet" =
            !is.null(cell_src))
message(glue("summarising cell from {cell_src$kind}"))

con  <- if (cell_src$kind == "db") dbConnect(duckdb(cell_db, read_only = TRUE)) else dbConnect(duckdb())
tbl_expr <- if (cell_src$kind == "db") "cell" else cell_src$from
smry <- dbGetQuery(con, glue("
  SELECT count(*) AS n_cell,
         count(*) FILTER (WHERE in_usa) AS n_usa,
         count(*) FILTER (WHERE in_pra) AS n_pra,
         round(min(area_km2), 3) AS area_min_km2,
         round(max(area_km2), 3) AS area_max_km2,
         round(min(lon),1) AS lon_min, round(max(lon),1) AS lon_max,
         count(*) FILTER (WHERE depth_mean IS NULL) AS n_null_depth
  FROM {tbl_expr}"))
# content fingerprint of the cell table (before disconnect)
h <- msens::hash_query(con, tbl_expr)
dbDisconnect(con, shutdown = TRUE)
msens::report_table(smry, caption = "global 0.05° cell grid summary")

global 0.05° cell grid summary
n_cell	n_usa	n_pra	area_min_km2	area_max_km2	lon_min	lon_max	n_null_depth
17072105	634208	343576	0.014	30.773	-180	180	0

Code

# content-addressed manifest: deterministic (no wall-clock, no machine paths) ->
# downstream targets re-run only when the cell table's content actually changes
msens::write_manifest(
  here("data/manifests/build_cell_grid.json"), target = "build_cell_grid", content_hash = h,
  stats = c(as.list(smry), list(version = ver, bbox = cell_bbox,
            source = basename(bio_oracle_tif))),
  force = msens::force_target("build_cell_grid"))
message(glue("cell grid: {smry$n_cell} cells ({smry$n_usa} in US, {smry$n_pra} in PRAs) -> {cell_db}"))

--- title: "Build global 0.05° cell grid (v8 sampling unit)" subtitle: "Global ocean 0.05° raster cells in [-180,180] with env covariates + US study-area / Program-Area membership" msens: target_name: build_cell_grid workflow_type: grid dependency: [bootstrap_version] # a same-grid version bump clones the unchanged ingests first output: data/manifests/build_cell_grid.json format: html editor_options: chunk_output_type: console --- The v8 spatial sampling unit is the **global 0.05° raster cell** in **[-180,180]** (rolled back from H3 hex). The grid is anchored to the shared 0.05°, 3600×7200, EPSG:4326 topology of **Bio-Oracle v3 AND the AquaX SDM TIFs** — so `cell_id` is the raster cell index (`1:ncell`), and any SDM raster on this grid (Bio-Oracle, AquaX, resampled AquaMaps) maps to `cell_id` by position with **no projection**. This notebook builds the canonical `cell` table (ocean cells only) — `cell_id`, the 12 Bio-Oracle env covariates, `area_km2`, and boolean `in_usa` / `in_pra` membership — plus the global cell-id COG that titiler paints (and that `ingest_aquamaps.qmd` uses to build its bilinear weight table). The global build (631 MB Bio-Oracle → 25.9M cells → ~17M-row `cell` table + 47 MB COG) is **skipped when the cell-id COG already exists** (resume); force it with `REDO_CELL_GRID=1` (`libs/vars.R`). A bbox prototype run (`CELL_BBOX=...`) always rebuilds into its own `CELL_DB`. ## Design ```{mermaid} %%| label: fig-design %%| fig-cap: "Bio-Oracle raster → global 0.05° cell grid: cell table + cell-id COG, content-addressed manifest" flowchart LR bo["Bio-Oracle v3 0.05° 12-layer raster"] --> id["assign cell_id=1:ncell, area_km2, in_usa/in_pra"] id --> ocn["keep ocean cells (depth_mean non-NA)"] ocn --> cell[("cell table sdm.duckdb")] ocn --> cog["r_cellid_global.tif COG (titiler paints)"] cell --> mf["hash_query(cell) → content-addressed manifest"] ``` ```{r setup} librarian::shelf(DBI, dplyr, duckdb, glue, here, jsonlite, sf, terra, quiet = TRUE) source(here("libs/paths.R")) # ver, sdm_db, bio_oracle_tif, cellid_tif, ply_usa_gpkg, ply_pra_gpkg source(here("libs/vars.R")) # redo_cell_grid stopifnot(file.exists(bio_oracle_tif), file.exists(ply_usa_gpkg), file.exists(ply_pra_gpkg)) # optional bbox restriction for prototype runs: CELL_BBOX="xmin,xmax,ymin,ymax" cell_bbox <- Sys.getenv("CELL_BBOX", "") cell_db <- Sys.getenv("CELL_DB", sdm_db) dir.create(dirname(cell_db), recursive = TRUE, showWarnings = FALSE) dir.create(here("data/manifests"), recursive = TRUE, showWarnings = FALSE) # skip the expensive global build if the cell-id COG already exists (resume), # unless forced; a bbox prototype always rebuilds into its own CELL_DB. # # The cell-id COG is SHARED across versions, but the `cell` TABLE lives in this # version's sdm.duckdb. A fresh version directory therefore has the COG and no # table -- and every downstream `stopifnot("run build_cell_grid first" = # file_exists(sdm_db))` passes on an empty database. Resume path for that case: # copy `cell` from ver_prev's build DB (seconds) and assert it hashes to the # checkpoint; rebuild only when there is nothing to copy. has_cell <- function(db) { db <- path.expand(db) if (!file.exists(db)) return(FALSE) c0 <- dbConnect(duckdb(db, read_only = TRUE)); on.exit(dbDisconnect(c0, shutdown = TRUE)) "cell" %in% dbListTables(c0) } prev_db <- glue("{dir_big}/{ver_prev}/sdm.duckdb") need_build <- redo_cell_grid || nzchar(cell_bbox) || !file.exists(cellid_tif) || (!has_cell(cell_db) && !has_cell(prev_db)) copy_cell <- !need_build && !has_cell(cell_db) && has_cell(prev_db) ``` ## Read Bio-Oracle, assign global cell_id, membership + area ```{r build_grid} if (need_build) { r <- rast(bio_oracle_tif) # global [-180,180], 3600 x 7200, 12 env layers r$cell_id <- 1:ncell(r) # global raster-position id (shared with AquaX/AquaMaps grids) r$area_km2 <- cellSize(r$depth_mean, unit = "km") # membership: rasterize the study-area + Program-Area polygons onto the same grid r$in_usa <- rasterize(vect(st_read(ply_usa_gpkg, quiet = TRUE)), r$depth_mean, field = 1, background = 0) r$in_pra <- rasterize(vect(st_read(ply_pra_gpkg, quiet = TRUE)), r$depth_mean, field = 1, background = 0) if (nzchar(cell_bbox)) { bb <- as.numeric(strsplit(cell_bbox, ",")[[1]]) r <- crop(r, ext(bb[1], bb[2], bb[3], bb[4])) message(glue("restricted to bbox [{cell_bbox}]")) } } ``` ## Write cell table (ocean cells) + global cell-id COG ```{r write_cell} if (need_build) { # ocean cells only (depth_mean non-NA); cell_id stays the global raster index d_cell <- as.data.frame(r, xy = TRUE, na.rm = FALSE) |> filter(!is.na(depth_mean)) |> transmute( cell_id = as.integer(cell_id), lon = x, lat = y, across(c(depth_mean, depth_min, depth_max, oxy_b_mean, oxy_mean, prim_prod_mean, ice_con_ann, salinity_b_mean, salinity_mean, sbt_an_mean, sst_an_mean, fao_area_m)), area_km2, in_usa = in_usa == 1, in_pra = in_pra == 1) con <- dbConnect(duckdb(cell_db)) dbWriteTable(con, "cell", d_cell, overwrite = TRUE) dbExecute(con, "ALTER TABLE cell ADD PRIMARY KEY (cell_id)") dbDisconnect(con, shutdown = TRUE) # global cell-id COG in [-180,180] for titiler (uint32; no overviews to keep ids intact) writeRaster(r$cell_id, cellid_tif, filetype = "COG", datatype = "INT4U", overwrite = TRUE, gdal = c("COMPRESS=LZW", "OVERVIEWS=NONE")) } else if (copy_cell) { # same grid, new version: the table is a straight copy of ver_prev's, verified # against the content-addressed checkpoint rather than trusted con <- dbConnect(duckdb(cell_db)) dbExecute(con, glue("ATTACH '{path.expand(prev_db)}' AS prev (READ_ONLY)")) dbExecute(con, "CREATE TABLE cell AS SELECT * FROM prev.cell") dbExecute(con, "ALTER TABLE cell ADD PRIMARY KEY (cell_id)") dbExecute(con, "DETACH prev") h_copy <- msens::hash_query(con, "cell") dbDisconnect(con, shutdown = TRUE) mf <- here("data/manifests/build_cell_grid.json") if (file.exists(mf)) { h_ckpt <- jsonlite::fromJSON(mf)$content_hash stopifnot("copied cell table does not hash to the build_cell_grid checkpoint" = identical(h_copy, h_ckpt)) } message(glue("cell table copied from {ver_prev} ({h_copy}) -> {cell_db}")) } else { message(glue("cell grid up to date (cellid_tif + cell table exist); set REDO_CELL_GRID=1 to rebuild")) } ``` ## Publish the cell-id COGs (gated) The apps resolve a clicked point through titiler's `/cog/point/{lon},{lat}`: one call for the layer's value, one for the cell id. That needs the cell-id lookups **on S3**, not just on the server's disk. Cell ids are **categorical**, so these are written `INT4U` with **no overviews** — a resampled pyramid would average neighbouring ids into ids that do not exist. Both grids are published: `usa05` serves v1–v7, `global05` serves v8. > Resolving a click by reading the raster in the app was the source of three separate bugs: the > band is named `r_cellid` on one grid and `depth_mean` on the other (so selecting `$cell_id` > returned NULL), the app shifted longitude to 0–360 while both rasters are stored −180..180 (so > every Americas click read outside the image), and a `SpatRaster` cached across Shiny sessions is > a stale external pointer that **segfaults** the process. One HTTP call has none of those modes. ```{r publish-cellid} if (nzchar(Sys.getenv("GRID_COG_S3"))) { grids <- msens::grid_registry() for (i in seq_len(nrow(grids))) { gid <- grids$grid_id[i] src <- path.expand(glue("{dir_derived}/{grids$cellid_tif[i]}")) if (!file.exists(src)) { message(glue(" {gid}: {basename(src)} absent - skip")); next } dst <- glue("{s3_atlas}/grid/{gid}/cellid.tif") out <- system2("aws", c("s3", "cp", shQuote(src), shQuote(dst), "--only-show-errors", "--no-progress"), stdout = TRUE, stderr = TRUE) if (!is.null(attr(out, "status")) && attr(out, "status") != 0) stop("aws s3 cp failed: ", paste(out, collapse = "\n")) message(glue(" {gid}: {round(file.size(src)/1024^2, 1)} MB -> {dst}")) } } else { message("GRID_COG_S3 unset - not publishing cell-id COGs") } ``` ## Summary + manifest ```{r summary} # The cell grid is built on the laptop; the SERVER carries only serve.duckdb and # the published Parquet, so `cell_db` does not exist there. Summarise whichever # is present -- otherwise this notebook cannot be rendered on the server at all, # which is exactly where the cell-id COGs above have to be published from. cell_src <- if (file.exists(cell_db)) { list(kind = "db", from = glue("'{cell_db}'")) } else { # tables/ live under big/{ver}; on the server dir_big_v collapses to # derived/{ver}, which is a different tree pq <- path.expand(glue("{dir_big}/{ver}/tables/cell.parquet")) if (file.exists(pq)) list(kind = "parquet", from = glue("read_parquet('{pq}')")) else NULL } stopifnot("no cell table: neither the grid database nor tables/cell.parquet" = !is.null(cell_src)) message(glue("summarising cell from {cell_src$kind}")) con <- if (cell_src$kind == "db") dbConnect(duckdb(cell_db, read_only = TRUE)) else dbConnect(duckdb()) tbl_expr <- if (cell_src$kind == "db") "cell" else cell_src$from smry <- dbGetQuery(con, glue(" SELECT count(*) AS n_cell, count(*) FILTER (WHERE in_usa) AS n_usa, count(*) FILTER (WHERE in_pra) AS n_pra, round(min(area_km2), 3) AS area_min_km2, round(max(area_km2), 3) AS area_max_km2, round(min(lon),1) AS lon_min, round(max(lon),1) AS lon_max, count(*) FILTER (WHERE depth_mean IS NULL) AS n_null_depth FROM {tbl_expr}")) # content fingerprint of the cell table (before disconnect) h <- msens::hash_query(con, tbl_expr) dbDisconnect(con, shutdown = TRUE) msens::report_table(smry, caption = "global 0.05° cell grid summary") # content-addressed manifest: deterministic (no wall-clock, no machine paths) -> # downstream targets re-run only when the cell table's content actually changes msens::write_manifest( here("data/manifests/build_cell_grid.json"), target = "build_cell_grid", content_hash = h, stats = c(as.list(smry), list(version = ver, bbox = cell_bbox, source = basename(bio_oracle_tif))), force = msens::force_target("build_cell_grid")) message(glue("cell grid: {smry$n_cell} cells ({smry$n_usa} in US, {smry$n_pra} in PRAs) -> {cell_db}")) ```