Ingest Spatial Boundaries

Published

2026-08-14

1 Overview

Process spatial boundary datasets into two outputs, both uploaded to GCS:

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).
spatial + spatial_attribute parquet for spatial queries — the ingest shard release_database.qmd picks up like every other, from data/parquet/spatial/.

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: ~/My Drive/projects/calcofi/data-public/_spatial/ (US Maritime Limits, 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.

1.1 Data Flow

Code

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"]

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"]

1.2 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)

2 Setup

Code

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"))

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`"
)

Code

# 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"
    )
  )
}

WARNING: The following folders exist in `/Users/bbest/My Drive/projects/calcofi/data-public/_spatial` but are not referenced in `metadata/spatial_layers.csv`:
- CA_congressional-districts
- CA_water-districts

Code

# 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"
)

19 layers across 17 dataset groups

3 ArcGIS Rest layer(s) to gpkg: `CA_watershed-boundaries`

Code

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)
}

4 Clip CA Current LME to US EEZ

Code

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
  )
}

5 Helper functions

Code

# 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
  })
}

6 Process datasets

For each dataset group: read source, add id/layer/name columns, generate PMTiles, and accumulate data for consolidated spatial tables.

Code

# 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"
  ))
}

noaa_maritime_boundaries: 260 features (CRS 4326 -> 4326), PMTiles 4.68M
ca_maritime_boundaries: 1 features (CRS 3310 -> 4326), PMTiles 120K
ca_marine_protected_areas: 155 features (CRS NA -> 4326), PMTiles 279K
noaa_onms_sanctuaries: 16 features (CRS 4326 -> 4326), PMTiles 367K
ca_cowcod_conservation: 1 features (CRS 3310 -> 4326), PMTiles 8.53K
ca_swqpa: 36 features (CRS 3310 -> 4326), PMTiles 81.1K
ca_county_boundaries: 58 features (CRS 3857 -> 4326), PMTiles 755K
ca_assembly_districts: 80 features (CRS 4269 -> 4326), PMTiles 965K
ca_senate_districts: 40 features (CRS 4269 -> 4326), PMTiles 815K
ca_cdfw_regions: 7 features (CRS 4326 -> 4326), PMTiles 610K
ca_watershed_boundaries: 140 features (CRS 3857 -> 4326), PMTiles 2.55M
boem_wind_planning: 9833 features (CRS 4326 -> 4326), PMTiles 4.06M
noaa_aquaculture_aoas: 10 features (CRS 6414 -> 4326), PMTiles 7.97K
noaa_ocean_disposal: 2148 features (CRS 4269 -> 4326), PMTiles 816K
ca_ports: 194 features (CRS 4326 -> 4326), PMTiles 75.1K
meow_ecoregions: 232 features (CRS 4326 -> 4326), PMTiles 20.7M
noaa_iea_regions: 1 features (CRS 4326 -> 4326), PMTiles 328K

7 Build consolidated `spatial` + `spatial_attribute` parquet

Code

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)

8 Upload to GCS

Code

# upload PMTiles to calcofi-files-public bucket
cat(glue(
  "Uploading PMTiles to gs://{gcs_bucket}/{gcs_prefix}/\n\n"
))

Uploading PMTiles to gs://calcofi-files-public/_spatial/

Code

sync_to_gcs(
  local_dir = dir_pmtiles,
  gcs_prefix = gcs_prefix,
  bucket = gcs_bucket,
  pattern = "\\.pmtiles$"
)

# A tibble: 17 × 4
   file    action    size reason        
   <chr>   <chr>    <dbl> <chr>         
 1 <rsync> uploaded    NA parallel rsync
 2 <rsync> uploaded    NA parallel rsync
 3 <rsync> uploaded    NA parallel rsync
 4 <rsync> uploaded    NA parallel rsync
 5 <rsync> uploaded    NA parallel rsync
 6 <rsync> uploaded    NA parallel rsync
 7 <rsync> uploaded    NA parallel rsync
 8 <rsync> uploaded    NA parallel rsync
 9 <rsync> uploaded    NA parallel rsync
10 <rsync> uploaded    NA parallel rsync
11 <rsync> uploaded    NA parallel rsync
12 <rsync> uploaded    NA parallel rsync
13 <rsync> uploaded    NA parallel rsync
14 <rsync> uploaded    NA parallel rsync
15 <rsync> uploaded    NA parallel rsync
16 <rsync> uploaded    NA parallel rsync
17 <rsync> uploaded    NA parallel rsync

Code

# 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"
))

Uploading parquets to gs://calcofi-db/ingest/spatial/

Code

sync_to_gcs(
  local_dir = dir_stage,
  sidecar_dir = dir_parquet,
  gcs_prefix = gcs_ingest_prefix,
  bucket = gcs_ingest_bucket
)

# A tibble: 3 × 4
  file    action    size reason        
  <chr>   <chr>    <dbl> <chr>         
1 <rsync> uploaded    NA parallel rsync
2 <rsync> uploaded    NA parallel rsync
3 <rsync> uploaded    NA parallel rsync

9 Write manifest

Code

write_spatial_manifest(parquet_dir = dir_stage, output_dir = dir_parquet)

10 Output summary

Code

cat("### PMTiles\n\n")

10.1 PMTiles

Code

dir_ls(dir_pmtiles, glob = "*.pmtiles") |>
  tibble(path = _) |>
  mutate(
    file = path_file(path),
    size = file_size(path)
  ) |>
  select(file, size) |>
  knitr::kable()

file	size
boem_wind_planning.pmtiles	4.06M
ca_assembly_districts.pmtiles	964.88K
ca_cdfw_regions.pmtiles	610.16K
ca_county_boundaries.pmtiles	754.86K
ca_cowcod_conservation.pmtiles	8.53K
ca_marine_protected_areas.pmtiles	278.9K
ca_maritime_boundaries.pmtiles	119.83K
ca_ports.pmtiles	75.1K
ca_senate_districts.pmtiles	815.27K
ca_swqpa.pmtiles	81.08K
ca_watershed_boundaries.pmtiles	2.55M
meow_ecoregions.pmtiles	20.67M
noaa_aquaculture_aoas.pmtiles	7.97K
noaa_iea_regions.pmtiles	327.56K
noaa_maritime_boundaries.pmtiles	4.68M
noaa_ocean_disposal.pmtiles	816.41K
noaa_onms_sanctuaries.pmtiles	367.03K

Code

cat("\n### Consolidated Parquet\n\n")

10.2 Consolidated Parquet

Code

dir_ls(dir_stage, glob = "*.parquet") |>
  tibble(path = _) |>
  mutate(
    file = path_file(path),
    size = file_size(path)
  ) |>
  select(file, size) |>
  knitr::kable()

file	size
spatial.parquet	53.7M
spatial_attribute.parquet	920.8K

Code

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"
))

10.3 Database Tables

spatial: 13206 features, 19 layers
spatial_attribute: 148461 attribute rows

1 Overview

1.1 Data Flow

1.2 Pre-processing

2 Setup

3 ArcGIS Rest layer(s) to gpkg: CA_watershed-boundaries

4 Clip CA Current LME to US EEZ

5 Helper functions

6 Process datasets

7 Build consolidated spatial + spatial_attribute parquet

8 Upload to GCS

9 Write manifest

10 Output summary

10.1 PMTiles

10.2 Consolidated Parquet

10.3 Database Tables

3 ArcGIS Rest layer(s) to gpkg: `CA_watershed-boundaries`

7 Build consolidated `spatial` + `spatial_attribute` parquet