---
title: "Ingest SDMs: NOAA SEFSC GoMex Cetacean & Sea-Turtle Densities"
subtitle: "Monthly hexagon abundance (#/40 km²) → density → suitability [0,100] on the global 0.05° grid"
msens:
target_name: ingest_sdm_gm
workflow_type: ingest
dependency: [build_cell_grid]
output: data/manifests/ingest_sdm_gm.json
dataset: {ds_key: gm, response_type: density, source_authority: NOAA SEFSC, temporal_interval: monthly, native_format: vector}
editor_options:
chunk_output_type: console
---
NOAA SEFSC [Gulf-of-Mexico cetacean & sea-turtle spatial density models](https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.nodc:256800)
(NCEI 0256800) — 19 hexagon shapefiles, one per species/population, each carrying monthly abundance
`{Mon}_n` (individuals per 40 km² hexagon, `-9999` = no data). v8 mirrors the **nc** density ingest
(`ingest_sdm-nc.qmd`): rasterize each (species × month) onto the global 0.05° grid, rescale the
heavy-tailed density to a graded **suitability [0,100]** per species (so months share a scale), plus
an **annual** mean surface. Two dolphins are modelled as split **Oceanic + Shelf** populations whose
densities are **summed** per hexagon before rasterizing.
- **`mdl_key`**: monthly `gm|{sp_id}|{01..12}` + annual `gm|{sp_id}`. `sp_id` = `WORMS:{aphia}` when
`taxa_sci` resolves to a single WoRMS AphiaID, else `GUILD:{taxa_shp}` for species groups.
- **Density → [0,100] is a modeling choice** (v7 never scored gm): per-`sp_id` linear rescale to the
pooled **p99.5** of non-zero density, stored as `val_cap` in `model_gm.csv`. **NOT yet folded into
merge/scoring** (merge treats non-`am` datasets as binary ranges; a scoring change to validate with
`pra_score_delta`).
> **Inputs are on Drive, not in-repo.** Set `GM_SHP_DIR` to the shapefile folder (or hydrate it to
> `~/_big/msens/raw/sdm/gm`). This notebook was authored against the nc pattern but **must be run on
> a host with the Drive inputs** (e.g. `msens1`); it is not executed in the in-repo CI render.
```{r}
#| label: setup
librarian::shelf(dplyr, fs, glue, here, janitor, jsonlite, logger, purrr, readr, readxl,
sf, stringr, terra, tibble, tidyr, worrms, MarineSensitivity/msens, quiet = T)
source(here("libs/paths.R"))
source(here("libs/vars.R"))
sf_use_s2(FALSE)
options(readr.show_col_types = FALSE)
ds_key <- "gm"
dir_dist <- glue("{dir_big_v}/marine-atlas/dist/dataset={ds_key}"); dir_create(dir_dist)
mdl_csv <- glue("{dir_big_v}/marine-atlas/dist/model_{ds_key}.csv")
manifest <- here("data/manifests/ingest_sdm_gm.json"); dir_create(path_dir(manifest))
cog_base <- "https://file.marinesensitivity.org/cog/sdm/raw"
p_cap <- 0.995
NA_VAL <- -9999
# Drive-hosted inputs (override with GM_SHP_DIR); default to a hydrated ~/_big copy if present.
# local default (Ben's Drive); override on any other host by exporting GM_SHP_DIR in the environment
if (!nzchar(Sys.getenv("GM_SHP_DIR")))
Sys.setenv(GM_SHP_DIR = "~/My Drive/projects/msens/data/raw/ncei.noaa.gov - GoMex cetacean & sea turtle SDMs/0256800/2.2/data/0-data/NOAA_SEFSC_Cetacean_SeaTurtle_SDM_shapefiles")
shp_dir <- Sys.getenv("GM_SHP_DIR", unset = path.expand("~/_big/msens/raw/sdm/gm"))
taxa_xls <- if (file_exists(glue("{shp_dir}/spp_gmx.xlsx"))) {
glue("{shp_dir}/spp_gmx.xlsx")
} else glue("{shp_dir}/../spp_gmx.xlsx")
stopifnot("run on a host with the Drive gm inputs (set GM_SHP_DIR)" = dir_exists(shp_dir),
file_exists(cellid_tif))
if (redo_ingest && dir_exists(dir_dist)) { dir_delete(dir_dist); dir_create(dir_dist) }
r_grid <- rast(cellid_tif) # global 0.05° template (values = cell_id)
# write one model surface as (mdl_key, cell_id, val) parquet (V2/zstd). mm = NA -> annual.
pq_path <- function(safe, mm) if (is.na(mm)) glue("{dir_dist}/{ds_key}_{safe}.parquet") else
glue("{dir_dist}/{ds_key}_{safe}_{mm}.parquet")
write_surface <- function(d, sid, mm, safe) {
mk <- if (is.na(mm)) glue("{ds_key}|{sid}") else glue("{ds_key}|{sid}|{mm}")
msens::write_atlas_parquet(
tibble(mdl_key = as.character(mk), cell_id = as.integer(d$cell_id), val = d$val),
pq_path(safe, mm))
}
```
## Models + species crosswalk
```{r}
#| label: models
# one row per shapefile: species (+ Oceanic/Shelf population), monthly abundance columns.
d_taxa <- read_excel(taxa_xls) # cols: taxa_shp, taxa_sci, taxa_doc, taxa_common
m <- tibble(path_shp = dir_ls(shp_dir, glob = "*.shp")) |>
mutate(
base_shp = path_ext_remove(path_file(path_shp)),
taxa_shp = str_replace(base_shp, "(.*)_Monthly_.*", "\\1"),
# the two split-population dolphins: Oceanic_* / Shelf_* -> combine per taxa
population = if_else(str_detect(base_shp, "^(Oceanic|Shelf)_"),
str_replace(base_shp, "^(Oceanic|Shelf)_.*", "\\1"), NA_character_)) |>
left_join(select(d_taxa, taxa_shp, taxa_sci, taxa_common = any_of("taxa_common")), by = "taxa_shp")
# resolve taxa_sci -> single WoRMS AphiaID (WORMS:{id}); species groups fall back to GUILD:{taxa_shp}
resolve_sp_id <- function(sci, shp) {
id <- tryCatch(worrms::wm_name2id(sci), error = function(e) NA_integer_)
if (length(id) == 1 && !is.na(id)) glue("WORMS:{id}") else glue("GUILD:{shp}")
}
m <- m |> mutate(sp_id = map2_chr(taxa_sci, taxa_shp, resolve_sp_id))
log_info("gm shapefiles: {nrow(m)}; sp_id: {n_distinct(m$sp_id)}; split-pop taxa: {sum(!is.na(m$population))}")
knitr::kable(m |> count(sp_id, name = "n_shp") |> arrange(desc(n_shp)) |> head(8))
```
## Rasterize + rescale per species
```{r}
#| label: rasterize
# per hexagon shapefile -> long monthly density (#/km2 = {Mon}_n / 40), summed across populations.
hex_density <- function(paths) {
sf <- bind_rows(lapply(paths, function(p) {
x <- read_sf(p) |> st_transform(4326)
mos <- intersect(glue("{month.abb}_n"), names(x))
x |> select(hexid = HEXID, all_of(mos), geom = geometry)
}))
mos <- intersect(glue("{month.abb}_n"), names(sf))
long <- sf |> st_drop_geometry() |>
pivot_longer(all_of(mos), names_to = "mo_n", values_to = "n") |>
filter(n != NA_VAL) |>
mutate(mm = sprintf("%02d", match(str_remove(mo_n, "_n"), month.abb)),
density = n / 40) # individuals / km2
geoms <- sf |> select(hexid) |> distinct()
list(long = long, geoms = geoms)
}
# rasterize a per-hexagon density onto the 0.05° grid, then read grid cells
cells_from_hex <- function(geoms, dens) { # dens: tibble(hexid, density)
g <- geoms |> left_join(dens, by = "hexid") |> filter(!is.na(density))
if (!nrow(g)) return(tibble(cell_id = integer(), val = numeric()))
r <- terra::rasterize(terra::vect(g), r_grid, field = "density")
msens::cells_from_raster(r, cellid_tif, min_value = 0)
}
sp_ids <- sort(unique(m$sp_id))
reg_rows <- list()
for (sid in sp_ids) {
mi <- m |> filter(sp_id == sid)
safe <- gsub("[^A-Za-z0-9._-]", "_", sid)
if (all(file_exists(pq_path(safe, NA))) && !redo_ingest) next
hd <- hex_density(mi$path_shp) # combine Oceanic+Shelf populations
geoms <- hd$geoms
# sum densities across populations per (hexid, month)
bym <- hd$long |> group_by(hexid, mm) |> summarise(density = sum(density), .groups = "drop")
cap <- as.numeric(quantile(bym$density[bym$density > 0], p_cap, names = FALSE))
if (!is.finite(cap) || cap <= 0) cap <- max(bym$density, na.rm = TRUE)
cells_mo <- list()
for (mm in sort(unique(bym$mm))) {
d <- cells_from_hex(geoms, bym |> filter(mm == !!mm) |> select(hexid, density))
if (!nrow(d)) next
d <- d |> mutate(val = pmin(val / cap * 100, 100)) |> filter(val >= 1)
write_surface(d, sid, mm, safe); cells_mo[[mm]] <- d
}
if (length(cells_mo)) { # annual = mean across months per cell
da <- bind_rows(cells_mo) |> group_by(cell_id) |>
summarise(val = round(mean(val), 2), .groups = "drop")
write_surface(da, sid, NA, safe)
}
reg_rows[[sid]] <- mi |> distinct(sp_id, taxa_sci) |> slice(1) |>
mutate(val_cap = round(cap, 4), months = paste(sort(unique(bym$mm)), collapse = "|"),
populations = paste(sort(unique(na.omit(mi$population))), collapse = "|"))
}
```
## Model registry (`model_gm.csv`)
```{r}
#| label: model-csv
# resumable: a re-render with every surface already on disk leaves reg_rows empty (the rasterize
# loop `next`s), so keep the existing registry rather than rebuilding it from nothing.
if (length(reg_rows) == 0 && file_exists(mdl_csv)) {
model_gm <- read_csv(mdl_csv)
} else {
reg <- bind_rows(reg_rows)
model_gm <- reg |> rowwise() |> do({
r <- .; mos <- strsplit(r$months, "\\|")[[1]]
tibble(
mdl_key = c(glue("{ds_key}|{r$sp_id}|{mos}"), glue("{ds_key}|{r$sp_id}")),
interval = c(mos, "annual"),
sp_id = r$sp_id, scientific_name = r$taxa_sci,
val_cap = r$val_cap, populations = r$populations)
}) |> ungroup()
model_gm <- model_gm |> rowwise() |> filter(file_exists(
glue("{dir_dist}/{ds_key}_{gsub('[^A-Za-z0-9._-]','_',sp_id)}{ifelse(interval=='annual','', paste0('_',interval))}.parquet"))) |>
ungroup()
write_csv(model_gm, mdl_csv)
}
log_info("model_gm.csv: {nrow(model_gm)} mdl_keys ({sum(model_gm$interval=='annual')} annual)")
knitr::kable(model_gm |> count(interval, name = "n"))
```
## Manifest
```{r}
#| label: manifest
h <- msens::hash_parquet(glue("{dir_dist}/*.parquet"))
smry <- arrow::open_dataset(dir_dist) |>
summarise(n_rows = n(), n_mdl = n_distinct(mdl_key),
val_min = round(min(val), 2), val_max = round(max(val), 2)) |> collect()
msens::report_table(smry, caption = "gm density model_cell surface (density → suitability [0,100])")
msens::write_manifest(
manifest, target = "ingest_sdm_gm", content_hash = h,
stats = list(
ds_key = ds_key, n_sp_id = n_distinct(model_gm$sp_id), n_mdl = nrow(model_gm), p_cap = p_cap,
note = "monthly + annual density->[0,100]; Oceanic+Shelf summed; NOT yet in merge/scoring"),
force = msens::force_target("ingest_sdm_gm"))
```
## Compare with AquaMaps
Most of these cetaceans + sea turtles also have an **AquaMaps** model, so the two SDMs can be
contrasted directly. Drag the **swipe** divider left/right: **AquaMaps** (native 0.5° suitability) on
the left vs this **GoMex** surface (0.05°, months averaged to one annual surface) on the right — both
rescaled to 0–100 (spectral). AquaMaps paints broad, coarse suitability across (and beyond) the Gulf;
the SEFSC survey model resolves fine density structure within it. Pick any species that has both
models and drag the round handle on the divider:
```{=html}
<iframe src="compare.html?embed=1" width="100%" height="780" loading="lazy"
style="border:1px solid #2b3a47; border-radius:6px;"
title="AquaMaps vs GoMex SDM swipe compare"></iframe>
```
Full-screen on the [compare page](compare.html){target="_blank"}.