---
title: "Ingest Spatial Boundaries"
calcofi:
target_name: ingest_spatial
workflow_type: spatial
dependency: []
output: data/parquet/spatial/manifest.json
editor_options:
chunk_output_type: console
---
## Overview
Process spatial boundary datasets into **two outputs, both uploaded to GCS**:
1. **PMTiles** for map overlay β one per dataset group, served from
`gs://calcofi-files-public/_spatial/` (which is where `db-viz-hex` reads
them; the local `data/pmtiles/` copy is a build artifact, not a channel).
2. **`spatial` + `spatial_attribute` parquet** for spatial queries β the ingest
shard `release_database.qmd` picks up like every other, from
`data/parquet/spatial/`.
Plus the **reference layers** (plan 2026-09-09, D47βD54; `libs/reference_layers.R`):
registry rows with `role = reference` are not regions β no parquet, no memberships β
but their PMTiles ship beside the boundary archives and their counts land in
`data/parquet/spatial/reference_layers.json`, which `build_spatial_layers()` reads at
release. Today: the OpenStreetMap **land mask** the Explorer draws over the sea floor
and the data (`osm_land`), and the **undersea feature names** of the IHO-IOC GEBCO
Gazetteer (`gebco_gazetteer`). A `raster` row (Esri's ocean reference) is registry-only.
This notebook used to end by loading that same parquet into a
`data/calcofi_wrangling.duckdb` at the repo root. Nothing ever read it: the
release builds `spatial` from the parquet shard directly, and a repo-wide grep
found the file referenced by this notebook alone. It was a 27 MB binary,
rewritten in full on every run and tracked in git, left over from the pre-shard
pipeline. Removed 2026-08-02.
**Source**: `{dir_data}/_spatial/` (`data-public/_spatial/` in the CalCOFI Data Folder)
([US Maritime Limits](https://nauticalcharts.noaa.gov/data/us-maritime-limits-and-boundaries.html),
CA CDFW, BOEM, NOAA ONMS, MEOW, and others)
**Registry**: `metadata/spatial_layers.csv` β one row per visual layer, with
dataset grouping, styling, filter expressions, and source file paths.
### Data Flow
```{mermaid}
graph LR
A[Source Shapefiles] --> B[Read + Reproject 4326]
B --> C["clean_names + id/layer/name"]
C --> D["spatial + spatial_attribute Parquet"]
C --> E["GeoJSON β tippecanoe β PMTiles"]
D --> F["gs://calcofi-db/ingest/spatial/<br/>β release_database.qmd"]
E --> G["gs://calcofi-files-public/_spatial/<br/>β db-viz-hex, explore"]
R["OSM land polygons Β· GEBCO gazetteer<br/>(cc_stage_path('reference'))"] --> S["crop Β· rank Β· tippecanoe"]
S --> E
S --> T["reference_layers.json<br/>β build_spatial_layers()"]
```
### Pre-processing
Each dataset is standardized:
- **CRS**: reprojected to `EPSG:4326` (WGS 84 unprojected)
- **Columns**: `janitor::clean_names()` for snake_case
- **ID**: integer `id` column (per-layer sequential)
- **layer**: readable layer name from `metadata/spatial_layers.csv`
- **name**: per-feature label from `name_field` (or `layer` for dissolved layers)
## Setup
```{r}
#| label: setup
librarian::shelf(
arcgislayers,
calcofi / calcofi4db,
DBI,
dplyr,
DT,
duckdb,
fs,
glue,
here,
janitor,
jsonlite,
purrr,
readr,
sf,
stringr,
tibble,
tidyr,
quiet = T
)
# common ingest settings (overwrite, dir_data)
source(here("libs/ingest.R"))
# the reference layers (land mask, gazetteer labels): sources, ranks, tippecanoe, manifest
source(here("libs/reference_layers.R"))
dir_reference <- cc_stage_path("reference", create = TRUE)
dir_spatial <- path_expand(glue("{dir_data}/_spatial"))
dir_pmtiles <- here("data/pmtiles")
# sidecars (manifest.json) stay in the repo; the parquet bytes stage outside it
dir_parquet <- here("data/parquet/spatial")
dir_stage <- cc_stage_path("parquet", "spatial", create = TRUE)
gcs_bucket <- "calcofi-files-public"
gcs_prefix <- "_spatial"
dir_create(dir_pmtiles)
dir_create(dir_parquet)
# archive source spatial files to GCS
# gcloud config set account bebest@ucsd.edu
# TODO: fix this since it seems to put into /archive/_spatial without a timestamp
# sync_to_gcs(
# local_dir = dir_spatial,
# gcs_prefix = "archive/_spatial",
# bucket = gcs_bucket,
# exclude = c(".DS_Store", "*.tmp")
# )
# load registry
d_layers <- read_csv(
here("metadata/spatial_layers.csv"),
show_col_types = FALSE
)
datatable(
d_layers,
caption = "Spatial layers registry from `metadata/spatial_layers.csv`"
)
# the reference rows (role = reference) build in their own section below; the loop
# here is the boundary layers only (an older registry without the column is all boundaries)
if (!"role" %in% names(d_layers)) d_layers$role <- "boundary"
d_layers_ref <- d_layers |> filter(role == "reference")
d_layers <- d_layers |> filter(is.na(role) | role != "reference")
# gazetteer rows (own `source_url` + `source_layer`) are drawn from the oceanmetrics gazetteer's
# PMTiles, so no source file and no PMTiles here; since v2026.10.14 their features join `spatial`
# (section "Gazetteer layers"), so samples get memberships and the Explorer can filter by them
d_layers_gaz <- d_layers[0, ]
if ("source_url" %in% names(d_layers)) {
d_layers_gaz <- d_layers |> filter(!is.na(source_url), source_url != "")
d_layers <- d_layers |> filter(is.na(source_url) | source_url == "")
}
# compare with folders/files in dir_spatial
spatial_dirs <- dir_ls(dir_spatial, type = "directory")
spatial_dirs_not_layers_csv <- setdiff(
basename(spatial_dirs),
unique(d_layers$source_folder)
)
if (length(spatial_dirs_not_layers_csv) > 0) {
cat(
glue(
"WARNING: The following folders exist in `{dir_spatial}` but are not ",
"referenced in `metadata/spatial_layers.csv`:\n",
paste0("- ", spatial_dirs_not_layers_csv, collapse = "\n"),
"\n\n"
)
)
}
# unique dataset groups (one pmtiles per group)
d_groups <- d_layers |>
distinct(dataset_group, source_folder, source_file, source_format) |>
mutate(
source_path = path(dir_spatial, source_folder, source_file),
file_exists = file_exists(source_path)
)
cat(
glue("{nrow(d_layers)} layers across {nrow(d_groups)} dataset groups"),
"\n"
)
```
## ArcGIS Rest layer(s) to gpkg: `CA_watershed-boundaries`
```{r}
#| label: arcgis-rest
fs_to_sf <- function(url, where = NULL) {
fs <- arcgislayers::arc_open(url)
lyr <- arcgislayers::get_layer(fs, 0)
sf <- arcgislayers::arc_select(
lyr,
where = where
) |>
tibble::tibble() |>
sf::st_as_sf() |>
janitor::clean_names()
return(sf)
}
# [USGS Watershed Boundary Dataset (WBD) HUC8, California Vicinity - Overview](https://www.arcgis.com/home/item.html?id=a98ef9ead82a4bd6aa09d46869f75a49&sublayer=0)
url_rest <- "https://services2.arcgis.com/Uq9r85Potqm3MfRV/arcgis/rest/services/WBD_HU8_wm/FeatureServer"
out_gpkg <- glue("{dir_spatial}/CA_watershed-boundaries/wbd-hu8_ca.gpkg")
if (!file_exists(out_gpkg) || overwrite) {
fs_to_sf(url_rest) |>
write_sf(out_gpkg, delete_dsn = TRUE)
}
```
## Clip CA Current LME to US EEZ
```{r}
#| label: clip-layers
cacurrent_geo <- glue(
"{dir_spatial}/NOAA_IEA-CA-Current/lme-ca-current.geojson"
)
usaeez_geo <- glue(
"{dir_spatial}/NOAA_IEA-CA-Current/usa-eez.geojson"
)
cacurrent_usaeez_geo <- glue(
"{dir_spatial}/NOAA_IEA-CA-Current/lme-ca-current_usa-eez.geojson"
)
if (!file_exists(cacurrent_usaeez_geo) || overwrite) {
cacurrent <- st_read(cacurrent_geo, quiet = TRUE)
usaeez <- st_read(usaeez_geo, quiet = TRUE)
cacurrent_usaeez <- st_intersection(cacurrent, usaeez)
st_write(
cacurrent_usaeez,
cacurrent_usaeez_geo,
delete_dsn = TRUE,
quiet = TRUE
)
}
```
## US maritime zones as polygons: `MarineRegions_maritime-zones`
The NOAA Office of Coast Survey limits (`NOAA_maritime-boundaries`) are lines, so
the 12 NM territorial sea, 24 NM contiguous zone and 200 NM EEZ matched no sample
in `sample_spatial` (workflows#169). Marine Regions publishes the same zones as
polygons; `libs/download_maritime_zones.R` writes the Pacific-coast parts into one
GeoJSON, flagged `ts` / `cz` / `eez` like the NOAA lines, so the registry's filter
expressions are unchanged.
```{r}
#| label: maritime-zones
source(here("libs/download_maritime_zones.R"))
download_maritime_zones(
glue("{dir_spatial}/MarineRegions_maritime-zones/us-maritime-zones_pacific.geojson"))
```
## Helper functions
```{r}
#| label: helpers
# rename whatever the source calls its geometry column to `geometry`
#
# NOT cosmetic. `bind_rows()` in build-spatial-parquet concatenates the per-layer
# sf objects by column NAME, so a source whose geometry sits in `geom` or `Shape`
# lands in a SECOND sfc column. `st_as_sf()` then keeps `geometry` as the active
# one, every row from that source gets an EMPTY geometry, and `st_write()` drops
# the extra column with a warning that scrolls past in the render. The layer
# reaches the release as GEOMETRYCOLLECTION EMPTY β present in `spatial`, with a
# name and an id and no shape.
#
# That is exactly how National Marine Sanctuaries (`geom`, .rda), CA Watersheds
# (`geom`, .gpkg) and Ocean Disposal Sites (`Shape`, .gpkg) shipped in
# v2026.08.02: drawable nowhere, and matching ZERO samples in the
# point-in-polygon join that feeds db-viz-hex's polygon summaries β so the
# failure looked like "no CalCOFI sampling in sanctuaries", which is false.
# Shapefiles all read as `geometry`, which is why every shp-sourced layer was
# fine and the bug looked format-specific rather than name-specific.
normalize_geom_col <- function(sf_data) {
if (is.null(sf_data)) {
return(NULL)
}
geom_col <- attr(sf_data, "sf_column")
if (is.null(geom_col) || geom_col == "geometry") {
return(sf_data)
}
names(sf_data)[names(sf_data) == geom_col] <- "geometry"
st_geometry(sf_data) <- "geometry"
sf_data
}
# read source spatial data from various formats
read_spatial_source <- function(
source_folder,
source_file,
source_format,
dir_spatial
) {
if (source_format == "rda") {
e <- new.env()
load(path(dir_spatial, source_folder, source_file), envir = e)
for (nm in ls(e)) {
if (inherits(get(nm, envir = e), "sf")) {
return(
get(nm, envir = e) |>
select(where(~ !is.list(.))) |>
normalize_geom_col()
)
}
}
return(NULL)
} else if (source_format == "gdb") {
# `source_file` must be "<file>.gdb|<layer>" β a .gdb holds many layers and
# there is no sensible default. Without the suffix `parts[2]` is NA and GDAL
# fails with "Cannot open layer NA", which the caller's tryCatch turns into a
# one-line SKIPPED note; BOEM Wind Planning Areas was absent from `spatial`
# AND from the PMTiles for exactly this reason. Fail loudly instead.
parts <- str_split_1(source_file, "\\|")
if (length(parts) < 2 || is.na(parts[2]) || parts[2] == "") {
stop(glue(
"source_format 'gdb' needs source_file '<file>.gdb|<layer>'; ",
"got '{source_file}'. Layers: ",
paste(st_layers(path(dir_spatial, source_folder, parts[1]))$name,
collapse = ", ")
))
}
st_read(
path(dir_spatial, source_folder, parts[1]),
layer = parts[2],
quiet = TRUE
) |>
normalize_geom_col()
} else {
st_read(
path(dir_spatial, source_folder, source_file),
quiet = TRUE
) |>
normalize_geom_col()
}
}
# apply a MapLibre filter expression as R filter on sf
# supports: ["==", ["get", "col"], val]
apply_filter_expr <- function(sf_data, filter_expr_json) {
if (is.na(filter_expr_json) || filter_expr_json == "") {
return(sf_data)
}
expr <- jsonlite::fromJSON(filter_expr_json, simplifyVector = FALSE)
if (
length(expr) == 3 &&
expr[[1]] == "==" &&
is.list(expr[[2]]) &&
expr[[2]][[1]] == "get"
) {
col_name <- expr[[2]][[2]]
val <- expr[[3]]
sf_data |> filter(.data[[col_name]] == val)
} else {
warning(glue("unsupported filter_expr: {filter_expr_json}"))
sf_data
}
}
# build spatial_attribute rows from an sf object
build_spatial_attr <- function(sf_data, layer_name) {
# exclude core columns and geometry
geom_col <- attr(sf_data, "sf_column")
core_cols <- c("id", "layer", "name", geom_col)
attr_cols <- setdiff(names(sf_data), core_cols)
if (length(attr_cols) == 0) {
return(tibble())
}
# pivot each attribute column into EAV rows
map_dfr(attr_cols, function(col) {
vals <- sf_data[[col]]
cls <- class(vals)[1]
base <- tibble(
# spatial_key is the JOIN key, not `id`. `id` is per-layer sequential, so
# joining on it alone silently mixes features across layers β a county and
# a sanctuary both have id 1. Namespaced the same way sample_key is, and
# named `_key` per the repo convention (string natural key; `_id` is an
# integer surrogate). `id` and `layer` are kept because both are still
# meaningful on their own.
spatial_key = paste0(layer_name, ":", sf_data$id),
layer = layer_name,
id = sf_data$id,
fld = col,
val_dbl = NA_real_,
val_int = NA_integer_,
val_chr = NA_character_,
val_date = as.Date(NA),
val_lgl = NA
)
if (cls %in% c("numeric", "double")) {
base$val_dbl <- as.double(vals)
} else if (cls == "integer") {
base$val_int <- as.integer(vals)
} else if (cls %in% c("character", "factor")) {
base$val_chr <- as.character(vals)
} else if (cls == "Date") {
base$val_date <- as.Date(vals)
} else if (cls == "logical") {
base$val_lgl <- as.logical(vals)
} else {
# fallback: coerce to character
base$val_chr <- as.character(vals)
}
base
})
}
```
## Process datasets
For each dataset group: read source, add `id`/`layer`/`name` columns,
generate PMTiles, and accumulate data for consolidated `spatial` tables.
```{r}
#| label: process-datasets
#| output: asis
# accumulators for consolidated tables
all_spatial <- list()
all_spatial_attr <- list()
for (i in seq_len(nrow(d_groups))) {
grp <- d_groups[i, ]
dataset_group <- grp$dataset_group
pmtiles_path <- path(dir_pmtiles, glue("{dataset_group}.pmtiles"))
# skip if already processed and not redoing
if (!overwrite && file_exists(pmtiles_path)) {
cat(glue("- **{dataset_group}**: skipped (already exists)\n\n"))
next
}
# read source data
sf_data <- tryCatch(
read_spatial_source(
grp$source_folder,
grp$source_file,
grp$source_format,
dir_spatial
),
error = function(e) {
cat(glue(
"- **{dataset_group}**: SKIPPED β failed to read ",
"`{grp$source_folder}/{grp$source_file}`: {e$message}\n\n"
))
NULL
}
)
if (is.null(sf_data) || nrow(sf_data) == 0) {
cat(glue("- **{dataset_group}**: SKIPPED β empty dataset\n\n"))
next
}
# reproject to EPSG:4326
src_crs <- st_crs(sf_data)$epsg
if (is.na(src_crs) || src_crs != 4326) {
sf_data <- st_transform(sf_data, 4326)
}
# clean column names
sf_data <- clean_names(sf_data)
# get layer rows for this group
grp_layers <- d_layers |> filter(dataset_group == !!dataset_group)
if (nrow(grp_layers) == 1 && is.na(grp_layers$filter_expr[1])) {
# single layer, no filter
layer_name <- grp_layers$layer[1]
name_field <- grp_layers$name_field[1]
if (is.na(name_field) || name_field == "") {
# dissolve all features into one (st_make_valid for bad geometries)
sf_data <- sf_data |>
st_make_valid() |>
summarise(geometry = st_union(geometry)) |>
mutate(
id = 1L,
layer = layer_name,
name = layer_name
) |>
relocate(id, layer, name)
} else {
sf_data <- sf_data |>
mutate(
id = row_number() |> as.integer(),
layer = layer_name,
name = as.character(.data[[name_field]])
) |>
relocate(id, layer, name)
}
# accumulate for spatial tables
geom_col_name <- attr(sf_data, "sf_column")
all_spatial[[layer_name]] <- sf_data |>
select(id, layer, name, all_of(geom_col_name))
all_spatial_attr[[layer_name]] <- build_spatial_attr(sf_data, layer_name)
} else {
# multi-layer group (e.g., noaa_maritime_boundaries with filter_exprs)
parts_list <- list()
for (j in seq_len(nrow(grp_layers))) {
lr <- grp_layers[j, ]
layer_name <- lr$layer
name_field <- lr$name_field
sf_filtered <- apply_filter_expr(sf_data, lr$filter_expr)
if (nrow(sf_filtered) == 0) {
next
}
if (is.na(name_field) || name_field == "") {
sf_filtered <- sf_filtered |>
st_make_valid() |>
summarise(geometry = st_union(geometry)) |>
mutate(id = 1L, layer = layer_name, name = layer_name) |>
relocate(id, layer, name)
} else {
sf_filtered <- sf_filtered |>
mutate(
id = row_number() |> as.integer(),
layer = layer_name,
name = as.character(.data[[name_field]])
) |>
relocate(id, layer, name)
}
parts_list[[layer_name]] <- sf_filtered
geom_col_name <- attr(sf_filtered, "sf_column")
all_spatial[[layer_name]] <- sf_filtered |>
select(id, layer, name, all_of(geom_col_name))
all_spatial_attr[[layer_name]] <- build_spatial_attr(
sf_filtered,
layer_name
)
}
# for PMTiles: use the full original sf_data with layer/name added per-feature
# tag each feature with its layer based on which filter it matches
sf_data <- sf_data |> mutate(layer = NA_character_, name = NA_character_)
for (j in seq_len(nrow(grp_layers))) {
lr <- grp_layers[j, ]
if (!is.na(lr$filter_expr) && lr$filter_expr != "") {
expr <- jsonlite::fromJSON(lr$filter_expr, simplifyVector = FALSE)
if (
length(expr) == 3 &&
expr[[1]] == "==" &&
is.list(expr[[2]]) &&
expr[[2]][[1]] == "get"
) {
col_nm <- expr[[2]][[2]]
val <- expr[[3]]
mask <- sf_data[[col_nm]] == val
mask[is.na(mask)] <- FALSE
sf_data$layer[mask] <- lr$layer
nf <- lr$name_field
if (!is.na(nf) && nf != "") {
sf_data$name[mask] <- as.character(sf_data[[nf]][mask])
} else {
sf_data$name[mask] <- lr$layer
}
}
}
}
# features not matching any filter: use group name
sf_data <- sf_data |>
mutate(
layer = if_else(is.na(layer), dataset_group, layer),
name = if_else(is.na(name), layer, name)
)
if (!"id" %in% names(sf_data)) {
sf_data <- sf_data |> mutate(id = row_number() |> as.integer())
} else {
sf_data <- sf_data |> mutate(id = as.integer(id))
}
sf_data <- sf_data |> relocate(id, layer, name)
}
# write PMTiles via tippecanoe
geojson_tmp <- tempfile(fileext = ".geojson")
st_write(sf_data, geojson_tmp, delete_dsn = TRUE, quiet = TRUE)
tippecanoe_args <- c(
"-o",
pmtiles_path,
"-z10",
"-Z0",
"--simplification=10",
# `--simplification` applies at EVERY zoom by default, including the maximum
# β so the deepest tiles were coarser than the source, and since every zoom
# above z10 overzooms from z10, no view anywhere in the app ever showed the
# real boundary. Restricting simplification to the low zooms keeps z10 (and
# therefore every overzoomed view above it) at full source resolution, while
# the far-out views stay small.
"--simplify-only-low-zooms",
# keep every polygon. tippecanoe otherwise merges or drops "tiny" polygons at
# low zoom, which quietly removes small features β a pocket MPA, a disposal
# site β at exactly the zooms where the whole layer is in view and a user
# would be counting them.
"--no-tiny-polygon-reduction",
"--no-tile-size-limit",
"-l",
dataset_group,
"--force",
geojson_tmp
)
system2("tippecanoe", args = tippecanoe_args, stdout = TRUE, stderr = TRUE)
unlink(geojson_tmp)
if (!file_exists(pmtiles_path)) {
cat(glue("- **{dataset_group}**: SKIPPED β tippecanoe failed\n\n"))
next
}
pmtiles_size <- file_size(pmtiles_path)
cat(glue(
"- **{dataset_group}**: {nrow(sf_data)} features ",
"(CRS {src_crs} -> 4326), PMTiles {pmtiles_size}\n\n"
))
}
```
## Gazetteer layers
The registry rows with their own `source_url` (BOEM wind leases and planning areas, OCS planning areas,
the 11th Program draft, Pacific oil and gas leases, NOAA SoCal AOAs, National Marine Sanctuaries) are
drawn from the [Ocean Metrics gazetteer](https://oceanmetrics.io/places/). Until v2026.10.13 they
stopped there, so no sample was ever intersected with them and the Explorer's Area filter greyed them
out. Now each one's `places.parquet` is read at the version the gazetteer manifest names, its sha256
checked (`calcofi4db::read_gazetteer_layer()`), and its features join `spatial` like any other layer.
The status columns ride along in `spatial_attribute` (`status`, `status_date`, `status_source`, the
source's own `status_upstream`, the dated `status_history`), so a release records which status each
region had when it was cut. The geometries are kept whatever the status: a rescinded wind area or a
cancelled lease stays a region you can filter by, with its status beside its name.
Gazetteer layers are national (the OCS planning areas, every sanctuary, 3,307 rescinded wind blocks), so
only the features that touch the CalCOFI study area join `spatial`, whole and unclipped, so a region's
membership stays exact. The study area is the extent of the station grid (`grid`, every pattern), the
same footprint every release intersects. A layer whose pinned version and sha256 match the previous run
is read from the local cache, not downloaded again; the table at the end says which layers changed.
```{r}
#| label: gazetteer-layers
#| output: asis
stopifnot("read_gazetteer_layer() needs calcofi4db >= 4.25.0; pull and install ../calcofi4db" =
packageVersion("calcofi4db") >= "4.25.0")
# the groups above only reach `spatial` when they were (re)built this run; appending the gazetteer
# layers alone would make build-spatial-parquet rewrite `spatial` with nothing else in it
stopifnot("gazetteer layers need every boundary group rebuilt in the same run (overwrite = TRUE)" =
overwrite || nrow(d_layers_gaz) == 0)
gaz_cols <- c("place_id", "status", "status_date", "status_source", "status_upstream", "status_history")
# the CalCOFI study area: the extent of the station grid (the release's `grid` table, v2026.10.13), lon/lat
STUDY_AREA_BBOX <- c(xmin = -135.25, ymin = 18.4, xmax = -105.75, ymax = 49.25)
study_area <- st_as_sfc(st_bbox(STUDY_AREA_BBOX, crs = 4326))
# the previous run's pins: a layer whose version + sha256 are unchanged is reported as such
pins_csv <- path(dir_parquet, "gazetteer_pins.csv")
pins_prev <- if (file_exists(pins_csv)) read_csv(pins_csv, show_col_types = FALSE) else
tibble(dataset_id = character(), version = character(), sha256 = character())
# planar test against a box: s2 rejects some national multipolygons (antimeridian splits) as invalid
s2_was <- sf_use_s2(FALSE)
gaz_pins <- list()
for (i in seq_len(nrow(d_layers_gaz))) {
lr <- d_layers_gaz[i, ]
# the collection slug is the archive's folder: <base><slug>/places.pmtiles
slug <- str_match(lr$source_url, "/([^/]+)/places\\.pmtiles$")[, 2]
stopifnot("gazetteer source_url must end in /<slug>/places.pmtiles" = !is.na(slug))
sf_g <- read_gazetteer_layer(slug, cache_dir = path(dir_stage, "gazetteer"))
pin <- attr(sf_g, "gazetteer")
sf_g <- sf_g |>
normalize_geom_col() |>
apply_filter_expr(lr$filter_expr)
name_field <- if (is.na(lr$name_field) || lr$name_field == "") "name" else lr$name_field
sf_g <- sf_g |>
filter(st_geometry_type(geometry) %in% c("POLYGON", "MULTIPOLYGON")) |>
(\(x) suppressMessages(st_filter(x, study_area, .predicate = st_intersects)))() |>
mutate(
id = row_number() |> as.integer(),
layer = lr$layer,
name = as.character(.data[[name_field]])) |>
select(id, layer, name, any_of(gaz_cols), geometry)
all_spatial[[lr$layer]] <- sf_g |> select(id, layer, name, geometry)
all_spatial_attr[[lr$layer]] <- build_spatial_attr(sf_g, lr$layer)
prev <- pins_prev |> filter(dataset_id == lr$dataset_id)
gaz_pins[[lr$dataset_id]] <- tibble(
dataset_id = lr$dataset_id, layer = lr$layer, slug = pin$slug,
version = pin$version, sha256 = pin$sha256, n_features = nrow(sf_g),
change = if (nrow(prev) == 0) "new" else if (identical(prev$sha256, pin$sha256)) "unchanged"
else glue("v{prev$version} -> v{pin$version}"))
cat(glue("- **{lr$layer}** (gazetteer `{pin$slug}` v{pin$version}): {nrow(sf_g)} features in the study area\n\n"))
}
sf_use_s2(s2_was)
# which gazetteer version each layer came from, next to the manifest
d_gaz_pins <- bind_rows(gaz_pins)
write_csv(d_gaz_pins, path(dir_parquet, "gazetteer_pins.csv"))
cat(knitr::kable(d_gaz_pins, caption = "gazetteer layers in `spatial`: the pinned version and sha256 of each places.parquet"), sep = "\n")
```
## Reference layers
The land mask and the undersea feature names (plan 2026-09-09; `libs/reference_layers.R`).
Sources download once under `r dir_reference` (the OSM zip is reused from the bathymetry
build when present); the archives rebuild with `overwrite`, and the manifest merges a
skipped layer's previous record so a partial run never blanks a count.
```{r}
#| label: reference-layers
#| output: asis
ref_manifest <- path(dir_parquet, "reference_layers.json")
ref <- build_reference_layers(dir_reference, dir_pmtiles, overwrite = overwrite)
write_reference_manifest(ref, ref_manifest)
ref <- jsonlite::fromJSON(ref_manifest, simplifyVector = FALSE)
cat(glue("Reference layers (built {ref$built}, tippecanoe {ref$tippecanoe}, box {paste(ref$bbox, collapse = ', ')}):\n\n"))
for (k in names(ref$layers)) {
l <- ref$layers[[k]]
cat(glue("- **{k}**: {format(l$n_features, big.mark = ',')} features, {round(l$bytes / 1e6, 1)} MB Β· ",
"source {l$source$url} ({l$source$version})\n\n"))
}
missing_ref <- setdiff(d_layers_ref$dataset_id[d_layers_ref$source_type != "raster"], names(ref$layers))
if (length(missing_ref)) cat(glue("\nWARNING: registry reference rows without an archive: {paste(missing_ref, collapse = ', ')}\n"))
```
## Build consolidated `spatial` + `spatial_attribute` parquet
```{r}
#| label: build-spatial-parquet
pq_spatial <- path(dir_stage, "spatial.parquet")
pq_spatial_attr <- path(dir_stage, "spatial_attribute.parquet")
# skip if parquet already exists and no new layers were processed
if (length(all_spatial) == 0 && file_exists(pq_spatial)) {
message("spatial parquet already exists, skipping rebuild")
# read back for downstream references
con_tmp <- dbConnect(duckdb())
dbExecute(con_tmp, "INSTALL spatial; LOAD spatial;")
spatial_combined <- st_read(
con_tmp,
query = glue("SELECT * FROM read_parquet('{pq_spatial}')")
)
attr_combined <- arrow::read_parquet(pq_spatial_attr)
dbDisconnect(con_tmp, shutdown = TRUE)
} else {
# combine all layers into single sf. spatial_key namespaces the per-layer `id`
# so it is unique across the whole table and safe to join on by itself β see
# the note in build_spatial_attr().
spatial_combined <- bind_rows(all_spatial) |>
st_as_sf() |>
mutate(spatial_key = paste0(layer, ":", id)) |>
relocate(spatial_key)
stopifnot(
"spatial_key must be unique across all layers" =
!any(duplicated(spatial_combined$spatial_key)))
# Every feature ships with a shape, or the ingest fails here.
#
# Both checks exist because the failure mode they catch is SILENT: a second
# sfc column means some source's geometry column escaped normalize_geom_col(),
# and its rows are about to be written with no geometry at all. That shipped
# once (3 layers, 2,304 features in v2026.08.02) and was only found from the
# consumer end, when db-viz-hex could not summarize within a sanctuary.
sfc_cols <- names(spatial_combined)[
vapply(spatial_combined, inherits, logical(1), "sfc")]
stopifnot(
"exactly one geometry column β a second means a source escaped normalize_geom_col()" =
length(sfc_cols) == 1)
empty_by_layer <- spatial_combined[st_is_empty(spatial_combined), ] |>
st_drop_geometry() |>
count(layer)
if (nrow(empty_by_layer) > 0) {
stop(glue(
"{sum(empty_by_layer$n)} features have empty geometry, in: ",
paste(glue("{empty_by_layer$layer} ({empty_by_layer$n})"),
collapse = ", ")))
}
message(glue(
"spatial: {nrow(spatial_combined)} features across ",
"{length(unique(spatial_combined$layer))} layers"
))
# combine all attributes
attr_combined <- bind_rows(all_spatial_attr)
message(glue("spatial_attribute: {nrow(attr_combined)} rows"))
# remove old per-dataset parquet files
old_parquets <- dir_ls(dir_stage, glob = "*.parquet")
old_parquets <- old_parquets[!grepl("^spatial", basename(old_parquets))]
if (length(old_parquets) > 0) {
file_delete(old_parquets)
message(glue(
"Removed {length(old_parquets)} old per-dataset parquet files"
))
}
# write spatial via DuckDB (preserves geometry as WKB in parquet)
spatial_geojson <- tempfile(fileext = ".geojson")
st_write(spatial_combined, spatial_geojson, delete_dsn = TRUE, quiet = TRUE)
con_tmp <- dbConnect(duckdb())
dbExecute(con_tmp, "INSTALL spatial; LOAD spatial;")
dbExecute(
con_tmp,
glue(
"COPY (SELECT * FROM ST_Read('{spatial_geojson}'))
TO '{pq_spatial}' (FORMAT PARQUET)"
)
)
dbDisconnect(con_tmp, shutdown = TRUE)
unlink(spatial_geojson)
# write spatial_attribute as plain parquet
arrow::write_parquet(attr_combined, pq_spatial_attr)
message(glue(
"Wrote spatial.parquet ({file_size(pq_spatial)}) ",
"and spatial_attribute.parquet ({file_size(pq_spatial_attr)})"
))
} # end else (rebuild)
```
## Upload to GCS
```{r}
#| label: upload-gcs
# upload PMTiles to calcofi-files-public bucket
cat(glue(
"Uploading PMTiles to gs://{gcs_bucket}/{gcs_prefix}/\n\n"
))
sync_to_gcs(
local_dir = dir_pmtiles,
gcs_prefix = gcs_prefix,
bucket = gcs_bucket,
pattern = "\\.pmtiles$"
)
# upload parquets to calcofi-db ingest bucket (for release_database.qmd)
dir_label <- basename(dir_parquet) # "spatial"
gcs_ingest_prefix <- glue("ingest/{dir_label}")
gcs_ingest_bucket <- "calcofi-db"
cat(glue(
"\nUploading parquets to gs://{gcs_ingest_bucket}/{gcs_ingest_prefix}/\n\n"
))
sync_to_gcs(
local_dir = dir_stage,
sidecar_dir = dir_parquet,
gcs_prefix = gcs_ingest_prefix,
bucket = gcs_ingest_bucket
)
```
## Write manifest
```{r}
#| label: manifest
write_spatial_manifest(parquet_dir = dir_stage, output_dir = dir_parquet)
```
## Output summary
```{r}
#| label: summary
#| output: asis
cat("### PMTiles\n\n")
dir_ls(dir_pmtiles, glob = "*.pmtiles") |>
tibble(path = _) |>
mutate(
file = path_file(path),
size = file_size(path)
) |>
select(file, size) |>
knitr::kable()
cat("\n### Consolidated Parquet\n\n")
dir_ls(dir_stage, glob = "*.parquet") |>
tibble(path = _) |>
mutate(
file = path_file(path),
size = file_size(path)
) |>
select(file, size) |>
knitr::kable()
cat(glue(
"\n### Database Tables\n\n",
"- `spatial`: {nrow(spatial_combined)} features, ",
"{length(unique(spatial_combined$layer))} layers\n",
"- `spatial_attribute`: {nrow(attr_combined)} attribute rows\n"
))
```
