Ingest SDMs: NCCOS Atlantic & Pacific Seabird Densities

Seasonal seabird density (#/km²) → suitability [0,100] on the global 0.05° grid

NCCOS at-sea seabird density models, Atlantic (Winship et al. 2018) + Pacific (Leirness et al. 2021), as seasonal COGs (n_per_km2, 2-km, oblique-mercator). v8 rasterizes each (species × season) onto the global 0.05° cell grid and rescales the heavy-tailed density to a graded suitability [0,100] per species (so seasons/regions share a scale), plus an annual mean surface for scoring.

• Source COGs (in-repo): data/sdm/raw/nc_atl_birds_dens/*.tif (140) + nc_pac_birds_dens/*.tif (135); band 1 = n_per_km2 (bands 2/3 = CV/CI). Also public at file.marinesensitivity.org/cog/sdm/raw/{ds}/{sp_code}_{season}.tif (served natively via titiler).
• mdl_key: seasonal nc|{sp_id}|{season} + annual nc|{sp_id} (sp_id = GBIF:…, or GUILD:{sp_code} for the multi-species guilds that lack a single taxon id). Regions (atl/pac) are spatially disjoint and unioned per sp_id.

Density → [0,100] is a modeling choice, not a v7-faithful recipe (v7 never scored nc). Here: per-sp_id linear rescale to the pooled p99.5 of non-zero density (caps outliers), stored as val_cap in model_nc.csv for reversibility. NOT yet folded into merge/scoring — see the note at the end (merge treats non-am datasets as binary ranges; folding density in is a scoring change to validate with pra_score_delta).

1 Setup

librarian::shelf(arrow, dplyr, fs, glue, here, jsonlite, logger, purrr, readr, sf, terra, tibble,
                 tidyr, MarineSensitivity/msens, quiet = T)
source(here("libs/paths.R"))
source(here("libs/vars.R"))
options(readr.show_col_types = FALSE)
`%||%` <- function(a, b) if (is.null(a) || length(a) == 0 || all(is.na(a))) b else a

ds_key   <- "nc"
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_nc.json"); dir_create(path_dir(manifest))
cog_base <- "https://file.marinesensitivity.org/cog/sdm/raw"
p_cap    <- 0.995                                        # density percentile used as the 100 point
stopifnot(file_exists(cellid_tif))

if (redo_ingest && dir_exists(dir_dist)) { dir_delete(dir_dist); dir_create(dir_dist) }

# write one model surface as (mdl_key, cell_id, val) parquet. ssn = NA -> annual key/file.
pq_path <- function(safe, ssn) if (is.na(ssn)) glue("{dir_dist}/{ds_key}_{safe}.parquet") else
                                              glue("{dir_dist}/{ds_key}_{safe}_{ssn}.parquet")
write_parquet_row <- function(d, sid, ssn, safe) {
  mk <- if (is.na(ssn)) glue("{ds_key}|{sid}") else glue("{ds_key}|{sid}|{ssn}")
  msens::write_atlas_parquet(                              # V2/zstd/80MB row groups (Phase A)
    tibble(mdl_key = as.character(mk), cell_id = as.integer(d$cell_id), val = d$val),
    pq_path(safe, ssn))
}

2 Models + species crosswalk

# one row per (region, sp_code, season) at band 1 (density)
m <- read_csv(here("data/nc_models.csv")) |>
  filter(var == "n_per_km2") |>
  transmute(region_ds = ds_key, region, sp_code, season, cog_url, bidx)

# sp_code -> sp_id / names, unioned across the two regions' crosswalks; guilds get GUILD:{sp_code}
read_xw <- function(f) {
  if (!file_exists(f)) return(tibble())
  read_csv(f) |> transmute(sp_code, sp_id = as.character(sp_id), sp_common, sp_scientific)
}
xw <- bind_rows(
  read_xw(here("data/sdm/derived/nc_atl_birds_dens/mw_spp.csv")),
  read_xw(here("data/sdm/derived/nc_pac_birds_dens/m_spp.csv")),
  read_xw(here("data/sdm/derived/nc_pac_birds_dens/m_spp-multiple.csv")),
  read_xw(here("data/sdm/derived/nc_atl_birds_dens/mw_spp-multiple.csv"))) |>
  distinct(sp_code, .keep_all = TRUE)

m <- m |> left_join(xw, by = "sp_code") |>
  mutate(
    is_guild = is.na(sp_id) | grepl("/", sp_scientific %||% ""),
    sp_id    = ifelse(is.na(sp_id), glue("GUILD:{sp_code}"), sp_id),
    scientific_name = sp_scientific %||% sp_code)
log_info("nc models: {nrow(m)} (species×season×region); {n_distinct(m$sp_id)} sp_id; guilds: {sum(m$is_guild & !duplicated(m$sp_id))}")
knitr::kable(m |> count(region, season, name = "n") |> pivot_wider(names_from = season, values_from = n))

region	fall	spring	summer	winter
Atlantic	39	39	35	27
Pacific	33	42	36	24

3 Rasterize + rescale per species

# local source cog for a model row (reproject omerc -> 4326, band 1)
read_cog <- function(region_ds, sp_code, season) {
  f <- here(glue("data/sdm/raw/{region_ds}/{sp_code}_{season}.tif"))
  if (!file_exists(f)) return(NULL)
  r <- rast(f, lyrs = 1)
  terra::project(r, "EPSG:4326", method = "bilinear")
}

# process per sp_id: pool all its (season×region) densities for a shared p99.5 cap, then write one
# seasonal surface per season (regions unioned) + an annual mean. resumable (skip existing parquet).
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)
  done <- all(file_exists(glue("{dir_dist}/{ds_key}_{safe}_{unique(mi$season)}.parquet"))) &&
          file_exists(glue("{dir_dist}/{ds_key}_{safe}.parquet"))
  if (done && !redo_ingest) next

  # read every source raster for this sp_id, pool non-zero densities -> cap
  rs   <- pmap(list(mi$region_ds, mi$sp_code, mi$season), read_cog)
  keep <- !vapply(rs, is.null, logical(1)); mi <- mi[keep, ]; rs <- rs[keep]
  if (!length(rs)) { log_warn("no source rasters for {sid}"); next }
  vals <- unlist(lapply(rs, function(r) { v <- values(r, mat = FALSE); v[!is.na(v) & v > 0] }))
  cap  <- as.numeric(quantile(vals, p_cap, names = FALSE))
  if (!is.finite(cap) || cap <= 0) cap <- max(vals, na.rm = TRUE)

  # per season: union regions' cells, rescale density->[0,100] by the shared cap
  cells_season <- list()
  for (ssn in unique(mi$season)) {
    idx <- which(mi$season == ssn)
    d <- map_dfr(idx, function(j) {
      r <- rs[[j]]; r <- terra::clamp(r / cap * 100, upper = 100, values = TRUE)
      msens::cells_from_raster(r, cellid_tif, min_value = 1)
    })
    if (!nrow(d)) next
    d <- d |> group_by(cell_id) |> summarise(val = max(val), .groups = "drop")   # union regions (disjoint)
    write_parquet_row(d, sid, ssn, safe)
    cells_season[[ssn]] <- d
  }
  # annual = mean across the available seasons (per cell)
  if (length(cells_season)) {
    da <- bind_rows(cells_season, .id = "season") |>
      group_by(cell_id) |> summarise(val = round(mean(val), 2), .groups = "drop")
    write_parquet_row(da, sid, NA, safe)
  }
  reg_rows[[sid]] <- mi |> distinct(sp_id, sp_code, scientific_name, is_guild) |>
    slice(1) |> mutate(val_cap = round(cap, 4),
                       seasons = paste(sort(unique(mi$season)), collapse = "|"),
                       regions = paste(sort(unique(mi$region)), collapse = "|"))
}

4 Model registry (model_nc.csv)

reg <- bind_rows(reg_rows)
# expand to one mdl_key per (seasonal + annual) surface actually written
model_nc <- reg |> rowwise() |> do({
  r <- .; ssns <- strsplit(r$seasons, "\\|")[[1]]; safe <- gsub("[^A-Za-z0-9._-]","_",r$sp_id)
  tibble(
    mdl_key = c(glue("{ds_key}|{r$sp_id}|{ssns}"), glue("{ds_key}|{r$sp_id}")),
    interval = c(ssns, "annual"),
    sp_id = r$sp_id, sp_code = r$sp_code, scientific_name = r$scientific_name,
    is_guild = r$is_guild, val_cap = r$val_cap, regions = r$regions,
    cog_url = c(glue("{cog_base}/nc_{ifelse(r$regions=='Atlantic','atl','pac')}_birds_dens/{r$sp_code}_{ssns}.tif"), NA))
}) |> ungroup() |>
  filter(file_exists(glue("{dir_dist}/{ds_key}_{gsub('[^A-Za-z0-9._-]','_',sp_id)}{ifelse(interval=='annual','', paste0('_',interval))}.parquet")))

write_csv(model_nc, mdl_csv)
log_info("model_nc.csv: {nrow(model_nc)} mdl_keys ({sum(model_nc$interval=='annual')} annual)")
knitr::kable(model_nc |> count(interval, name = "n"))

interval	n
annual	79
fall	61
spring	67
summer	62
winter	44

5 Manifest

# content-addressed, deterministic manifest (Phase A) — no wall-clock, so an identical re-ingest
# produces byte-identical bytes and does NOT invalidate downstream targets.
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 = "nc density model_cell surface (density → suitability [0,100])")

nc density model_cell surface (density → suitability [0,100])

n_rows	n_mdl	val_min	val_max
12069856	313	1	100

msens::write_manifest(
  manifest, target = "ingest_sdm_nc", content_hash = h,
  stats = list(
    ds_key = ds_key, n_sp_id = n_distinct(model_nc$sp_id), n_mdl = nrow(model_nc), p_cap = p_cap,
    note = "seasonal + annual density->[0,100]; NOT yet folded into merge/scoring"),
  force = msens::force_target("ingest_sdm_nc"))

References

Leirness, J. B., Josh Adams, Lisa T. Ballance, Michael Coyne, Jonathan J. Felis, Trevor Joyce, David M. Pereksta, Arliss J. Winship, Christopher FG Jeffrey, and D. Ainley. 2021. “Modeling at-Sea Density of Marine Birds to Support Renewable Energy Planning on the Pacific Outer Continental Shelf of the Contiguous United States.” https://repository.library.noaa.gov/view/noaa/49073.

Winship, Arliss J., Brian Patrick Kinlan, Timothy Paul White, Jeffrey Leirness, and John Christensen. 2018. “Modeling at-Sea Density of Marine Birds to Support Atlantic Marine Renewable Energy Planning.”