Ingest NMFS ESA Critical Habitat → global 0.05° cells

Published

2026-07-13 19:16:59

Rasterize each species’ NMFS ESA critical-habitat polygon(s) onto the global 0.05° cell grid and write one Parquet per species to the marine-atlas release, keyed by mdl_key = "ch_nmfs|{sp_id}" (sp_id = a slug of the scientific name — NMFS critical habitat has no numeric species id).

Method. msens::cells_from_ranges() (exactextractr) marks every cell the (unioned) polygon overlaps, capturing the whole footprint with no land mask (here ocean anyway, critical habitat being marine). The val is the species’ extinction-risk score via msens::compute_er_score() — never hard-coded — from its most-protective ESA status (NMFS:EN = 100, NMFS:TN = 50). This is the shared vector → cells pattern reused by every range dataset (BOTW, FWS, NMFS core areas, turtles).

Source ply.gpkg is derived from the NMFS ESA Critical Habitat geodatabase (see old/ingest_fisheries.noaa.gov_critical-habitat_v7.qmd for that one-time extraction + QA maps). Re-run control: REDO_INGEST=1 rebuilds all models (else the loop resumes, skipping species whose Parquet exists).

1 Design

Figure 1: NMFS ESA critical habitat → extinction-risk value on the global 0.05° grid

2 Setup

Code

librarian::shelf(
  arrow, DBI, dplyr, duckdb, fs, glue, here, jsonlite, logger, readr, sf, terra,
  tibble, MarineSensitivity/msens, quiet = T)
source(here("libs/paths.R"))   # dir_raw, dir_big_v, cellid_tif, ver, s3_atlas
source(here("libs/vars.R"))    # redo_ingest
options(readr.show_col_types = F)
sf_use_s2(FALSE)               # critical-habitat polygons are planar-valid, s2-invalid (self-

Spherical geometry (s2) switched off

Code

                               # crossing loops); planar GEOS unions/valids them (matches ch_fws/rng_fws)

ds_key   <- "ch_nmfs"
ply_geo  <- glue("{dir_raw}/fisheries.noaa.gov/ply.gpkg")   # per-species critical-habitat polygons
dir_dist <- glue("{dir_big_v}/marine-atlas/dist/dataset={ds_key}")
manifest <- here("data/manifests/ingest_ch_nmfs.json")
dir_create(c(dir_dist, path_dir(manifest)))

stopifnot(
  "run build_cell_grid.qmd first (need cellid_tif)" = file_exists(cellid_tif),
  "ply.gpkg missing (derive it, see old/..._v7.qmd)" = file_exists(ply_geo))

3 Read polygons + per-species model crosswalk

One model per species (SCIENAME); its multiple critical-habitat polygons are unioned at rasterize time.

Code

ply <- st_read(ply_geo, quiet = TRUE)
sp  <- ply |>
  st_drop_geometry() |>
  group_by(SCIENAME) |>                       # one model per species (polygons unioned below)
  summarise(
    TAXON = first(TAXON), COMNAME = first(COMNAME),
    # most-protective ESA status among the species' critical-habitat features
    status = if_else(any(grepl("Endangered", LISTSTATUS)), "EN", "TN"), .groups = "drop") |>
  arrange(SCIENAME) |>
  mutate(
    sp_id    = gsub("^_|_$", "", gsub("[^A-Za-z0-9]+", "_", SCIENAME)),
    mdl_key  = mdl_key_raw(ds_key, sp_id),
    er_code  = paste0("NMFS:", status),
    er_score = compute_er_score(er_code))     # val via compute_er_score(), never hard-coded
stopifnot(!any(duplicated(sp$sp_id)))
write_csv(sp, glue("{dir_dist}/../model_{ds_key}.csv"))
nrow(sp)

[1] 38

4 Rasterize each range → one Parquet each

Code

if (redo_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 |> filter(!sp_id %in% done)
log_info("{ds_key}: {length(done)} done; rasterizing {nrow(todo)} of {nrow(sp)} ranges ...")

for (i in seq_len(nrow(todo))) {
  row    <- todo[i, ]
  ply_sp <- ply |> filter(SCIENAME == row$SCIENAME) |> summarise(.groups = "drop")
  d <- cells_from_ranges(ply_sp, cellid_tif, value = row$er_score)   # ER score, not hard-coded
  if (nrow(d) == 0) { log_warn("  {row$sp_id}: no ocean cells covered"); next }
  msens::write_atlas_parquet(
    tibble(mdl_key = row$mdl_key, cell_id = d$cell_id, val = d$val),
    path(dir_dist, row$sp_id, ext = "parquet"))
}

although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar

Code

log_info("{ds_key}: done")

5 Verify

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 = "ch_nmfs model_cell surface")

ch_nmfs model_cell surface
n_models	n_cells	v_min	v_max
38	267897	50	100

6 Manifest (target output)

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_ch_nmfs", content_hash = h,
  stats = list(ds_key = ds_key, ver = ver, n_models = nrow(sp),
               n_cells = smry$n_cells, v_min = smry$v_min, v_max = smry$v_max,
               n_parquet = length(pq), source = path_file(ply_geo)),
  force = msens::force_target("ingest_ch_nmfs"))