Ingest BirdLife BOTW → global 0.05° cells (whole range)

Published

2026-07-13 17:47:57

Rasterize each bird species’ BirdLife BOTW range onto the global 0.05° cell grid and write one Parquet per species, keyed by mdl_key = "bl|{sisid}".

Whole range, no land mask. A bird’s global home range is largely terrestrial, so v8 captures the entire distribution (land + ocean) — a departure from prior versions that masked to the ocean. msens::cells_from_ranges() (exactextractr) is fast and tolerant of BOTW’s messy geometry (self-touching / MULTISURFACE loops), using extant polygons (presence ∈ {1,2,3}). The val is the IUCN extinction-risk score via msens::compute_er_score("IUCN:{code}") (never hard-coded; CR=50, EN=25, VU=5, NT=2, LC/DD=1) — the Red List category normalized to that vocabulary. How marine a range is (pct_marine) is a derived metric computed at merge, where the cell table’s is_ocean joins cleanly for every dataset — not a filter here.

Parallelized with parallel::mclapply (fork); resumable (skip species whose Parquet exists). REDO_INGEST=1 rebuilds all.

1 Design

Figure 1: BirdLife BOTW ranges → whole-range extinction-risk value on the global 0.05° grid

2 Setup

Code

librarian::shelf(
  arrow, DBI, dplyr, duckdb, fs, glue, here, jsonlite, logger, readr, RSQLite, sf,
  terra, tibble, MarineSensitivity/msens, quiet = T)
source(here("libs/paths.R"))   # dir_raw, dir_big_v, cellid_tif, ver
source(here("libs/vars.R"))    # redo_ingest
options(readr.show_col_types = F)
# BOTW stores some ranges as curved MULTISURFACE, which GEOS's st_make_valid cannot
# handle (it errored/hung on ~159 wide-ranging birds). Linearize curves -> MULTIPOLYGON
# on read so cells_from_ranges can rasterize them.
Sys.setenv(OGR_STROKE_CURVE = "TRUE")

ds_key    <- "bl"
# prefer a fast LOCAL copy off Google Drive — per-species queries throttle badly on Drive
botw_gpkg <- {
  f_local <- path.expand("~/_big/msens/raw/BOTW_2024_2.gpkg")
  if (file.exists(f_local)) f_local else glue("{dir_raw}/birdlife.org/BOTW_GPKG_2024_2/BOTW_2024_2.gpkg")
}
dir_dist  <- glue("{dir_big_v}/marine-atlas/dist/dataset={ds_key}")
manifest  <- here("data/manifests/ingest_birdlife_botw.json")
n_workers <- max(1L, parallel::detectCores())
dir_create(c(dir_dist, path_dir(manifest)))

stopifnot(
  "run build_cell_grid.qmd first (need cellid_tif)" = file_exists(cellid_tif),
  "BOTW gpkg not found"                             = file_exists(botw_gpkg))

# ensure sisid is indexed so the ~11k per-species queries are fast (0.03s vs a 10.9s
# full scan of the 17,379-row all_species table — 300x)
if (file.access(botw_gpkg, mode = 2) == 0) {  # writable copy
  con_gpkg <- dbConnect(RSQLite::SQLite(), botw_gpkg)
  try(dbExecute(con_gpkg, "CREATE INDEX IF NOT EXISTS idx_all_species_sisid ON all_species(sisid)"),
      silent = TRUE)
  dbDisconnect(con_gpkg)
}

3 Species list + model crosswalk

One model per sisid (extant polygons unioned at rasterize time); skip extinct (EX). IUCN category from the checklist is carried as merge metadata.

Code

a_attr <- st_read(botw_gpkg, quiet = TRUE,
  query = "SELECT DISTINCT sisid, sci_name FROM all_species WHERE presence IN (1,2,3)")
chk <- st_read(botw_gpkg, quiet = TRUE,
  query = "SELECT ScientificName AS sci_name, FamilyName AS family, CommonName AS common,
             IUCN_RedList_Category_2024 AS iucn
             FROM main_BL_HBW_Checklist_V9") |> st_drop_geometry()

# normalize the IUCN Red List category to compute_er_score()'s vocabulary
# (CR/EN/VU/NT/LC/DD): "CR (PE)"/"CR (PEW)" -> CR; drop extinct (EX/EW); blank -> LC.
# Value = the IUCN extinction-risk score (never hard-coded); merge may add MBTA etc.
sp <- a_attr |> st_drop_geometry() |>
  left_join(chk, by = "sci_name") |>
  mutate(iucn = toupper(trimws(iucn))) |>
  filter(is.na(iucn) | !iucn %in% c("EX", "EW")) |>
  mutate(
    code = case_when(
      is.na(iucn) | iucn == ""            ~ "LC",
      grepl("^CR", iucn)                  ~ "CR",
      iucn %in% c("EN","VU","NT","LC","DD") ~ iucn,
      TRUE                                ~ "LC"),
    er_code  = paste0("IUCN:", code),
    er_score = compute_er_score(er_code),
    mdl_key  = mdl_key_raw(ds_key, sisid)) |>
  arrange(sci_name)
stopifnot(!any(duplicated(sp$sisid)))
write_csv(sp, glue("{dir_dist}/../model_{ds_key}.csv"))
nrow(sp)

[1] 10995

4 Rasterize each range → one Parquet each

Code

if (redo_ingest && dir_exists(dir_dist)) { dir_delete(dir_dist); dir_create(dir_dist) }
done <- as.integer(path_ext_remove(path_file(dir_ls(dir_dist, glob = "*.parquet"))))
todo <- sp |> filter(!sisid %in% done)
log_info("{ds_key}: {length(done)} done; rasterizing {nrow(todo)} of {nrow(sp)} ranges on {n_workers} cores ...")

one <- function(i) {
  sisid  <- todo$sisid[i]
  out_pq <- fs::path(dir_dist, sisid, ext = "parquet")
  p <- tryCatch(sf::st_read(botw_gpkg, quiet = TRUE, query = sprintf(
    "SELECT sisid, geom FROM all_species WHERE sisid = %d AND presence IN (1,2,3)", sisid)),
    error = function(e) NULL)
  if (is.null(p) || nrow(p) == 0) return("empty")
  d <- tryCatch(cells_from_ranges(p, cellid_tif, value = todo$er_score[i]),
                error = function(e) NULL)
  if (is.null(d) || nrow(d) == 0) return("empty")
  msens::write_atlas_parquet(
    tibble::tibble(mdl_key = todo$mdl_key[i], cell_id = d$cell_id, val = d$val), out_pq)
  "done"
}
res <- unlist(parallel::mclapply(seq_len(nrow(todo)),
  function(i) tryCatch(one(i), error = function(e) paste0("error: ", conditionMessage(e))),
  mc.cores = n_workers, mc.preschedule = FALSE))
print(table(res))

res
 done 
10995

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

bl model_cell surface
n_models	n_cells	v_max
10995	1693299770	50

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

--- title: "Ingest BirdLife BOTW → global 0.05° cells (whole range)" msens: target_name: ingest_birdlife_botw workflow_type: ingest dependency: [build_cell_grid] output: data/manifests/ingest_birdlife_botw.json dataset: {ds_key: bl, response_type: range, source_authority: BirdLife, temporal_interval: static, native_format: vector} editor_options: chunk_output_type: console --- Rasterize each bird species' **BirdLife BOTW** range onto the **global 0.05° cell grid** and write one Parquet per species, keyed by `mdl_key = "bl|{sisid}"`. **Whole range, no land mask.** A bird's global home range is largely terrestrial, so v8 captures the **entire** distribution (land + ocean) — a departure from prior versions that masked to the ocean. `msens::cells_from_ranges()` (exactextractr) is fast and tolerant of BOTW's messy geometry (self-touching / MULTISURFACE loops), using extant polygons (`presence ∈ {1,2,3}`). The **val is the IUCN extinction-risk score** via `msens::compute_er_score("IUCN:{code}")` (never hard-coded; CR=50, EN=25, VU=5, NT=2, LC/DD=1) — the Red List category normalized to that vocabulary. How marine a range is (`pct_marine`) is a derived metric computed at **merge**, where the `cell` table's `is_ocean` joins cleanly for every dataset — not a filter here. Parallelized with `parallel::mclapply` (fork); resumable (skip species whose Parquet exists). `REDO_INGEST=1` rebuilds all. ## Design ```{mermaid} %%| label: fig-design %%| fig-cap: "BirdLife BOTW ranges → whole-range extinction-risk value on the global 0.05° grid" flowchart LR src["BOTW gpkg all_species (presence∈1:3)"] --> chk["join BL/HBW checklist → IUCN Red List category"] chk --> er["val = compute_er_score('IUCN:code') CR50 EN25 VU5 NT2 LC/DD1"] er --> ee["cells_from_ranges() exact_extract vs 0.05° cellid COG"] ee --> pq[("dist/dataset=bl/{sisid}.parquet (mdl_key, cell_id, val) · V2/zstd")] pq --> mf["hash_parquet → content-addressed manifest"] ``` ## Setup ```{r} #| label: setup librarian::shelf( arrow, DBI, dplyr, duckdb, fs, glue, here, jsonlite, logger, readr, RSQLite, sf, terra, tibble, MarineSensitivity/msens, quiet = T) source(here("libs/paths.R")) # dir_raw, dir_big_v, cellid_tif, ver source(here("libs/vars.R")) # redo_ingest options(readr.show_col_types = F) # BOTW stores some ranges as curved MULTISURFACE, which GEOS's st_make_valid cannot # handle (it errored/hung on ~159 wide-ranging birds). Linearize curves -> MULTIPOLYGON # on read so cells_from_ranges can rasterize them. Sys.setenv(OGR_STROKE_CURVE = "TRUE") ds_key <- "bl" # prefer a fast LOCAL copy off Google Drive — per-species queries throttle badly on Drive botw_gpkg <- { f_local <- path.expand("~/_big/msens/raw/BOTW_2024_2.gpkg") if (file.exists(f_local)) f_local else glue("{dir_raw}/birdlife.org/BOTW_GPKG_2024_2/BOTW_2024_2.gpkg") } dir_dist <- glue("{dir_big_v}/marine-atlas/dist/dataset={ds_key}") manifest <- here("data/manifests/ingest_birdlife_botw.json") n_workers <- max(1L, parallel::detectCores()) dir_create(c(dir_dist, path_dir(manifest))) stopifnot( "run build_cell_grid.qmd first (need cellid_tif)" = file_exists(cellid_tif), "BOTW gpkg not found" = file_exists(botw_gpkg)) # ensure sisid is indexed so the ~11k per-species queries are fast (0.03s vs a 10.9s # full scan of the 17,379-row all_species table — 300x) if (file.access(botw_gpkg, mode = 2) == 0) { # writable copy con_gpkg <- dbConnect(RSQLite::SQLite(), botw_gpkg) try(dbExecute(con_gpkg, "CREATE INDEX IF NOT EXISTS idx_all_species_sisid ON all_species(sisid)"), silent = TRUE) dbDisconnect(con_gpkg) } ``` ## Species list + model crosswalk One model per `sisid` (extant polygons unioned at rasterize time); skip extinct (`EX`). IUCN category from the checklist is carried as merge metadata. ```{r} #| label: model_csv a_attr <- st_read(botw_gpkg, quiet = TRUE, query = "SELECT DISTINCT sisid, sci_name FROM all_species WHERE presence IN (1,2,3)") chk <- st_read(botw_gpkg, quiet = TRUE, query = "SELECT ScientificName AS sci_name, FamilyName AS family, CommonName AS common, IUCN_RedList_Category_2024 AS iucn FROM main_BL_HBW_Checklist_V9") |> st_drop_geometry() # normalize the IUCN Red List category to compute_er_score()'s vocabulary # (CR/EN/VU/NT/LC/DD): "CR (PE)"/"CR (PEW)" -> CR; drop extinct (EX/EW); blank -> LC. # Value = the IUCN extinction-risk score (never hard-coded); merge may add MBTA etc. sp <- a_attr |> st_drop_geometry() |> left_join(chk, by = "sci_name") |> mutate(iucn = toupper(trimws(iucn))) |> filter(is.na(iucn) | !iucn %in% c("EX", "EW")) |> mutate( code = case_when( is.na(iucn) | iucn == "" ~ "LC", grepl("^CR", iucn) ~ "CR", iucn %in% c("EN","VU","NT","LC","DD") ~ iucn, TRUE ~ "LC"), er_code = paste0("IUCN:", code), er_score = compute_er_score(er_code), mdl_key = mdl_key_raw(ds_key, sisid)) |> arrange(sci_name) stopifnot(!any(duplicated(sp$sisid))) write_csv(sp, glue("{dir_dist}/../model_{ds_key}.csv")) nrow(sp) ``` ## Rasterize each range → one Parquet each ```{r} #| label: interpolate if (redo_ingest && dir_exists(dir_dist)) { dir_delete(dir_dist); dir_create(dir_dist) } done <- as.integer(path_ext_remove(path_file(dir_ls(dir_dist, glob = "*.parquet")))) todo <- sp |> filter(!sisid %in% done) log_info("{ds_key}: {length(done)} done; rasterizing {nrow(todo)} of {nrow(sp)} ranges on {n_workers} cores ...") one <- function(i) { sisid <- todo$sisid[i] out_pq <- fs::path(dir_dist, sisid, ext = "parquet") p <- tryCatch(sf::st_read(botw_gpkg, quiet = TRUE, query = sprintf( "SELECT sisid, geom FROM all_species WHERE sisid = %d AND presence IN (1,2,3)", sisid)), error = function(e) NULL) if (is.null(p) || nrow(p) == 0) return("empty") d <- tryCatch(cells_from_ranges(p, cellid_tif, value = todo$er_score[i]), error = function(e) NULL) if (is.null(d) || nrow(d) == 0) return("empty") msens::write_atlas_parquet( tibble::tibble(mdl_key = todo$mdl_key[i], cell_id = d$cell_id, val = d$val), out_pq) "done" } res <- unlist(parallel::mclapply(seq_len(nrow(todo)), function(i) tryCatch(one(i), error = function(e) paste0("error: ", conditionMessage(e))), mc.cores = n_workers, mc.preschedule = FALSE)) print(table(res)) ``` ## Verify ```{r} #| label: verify pq <- dir_ls(dir_dist, glob = "*.parquet") con <- dbConnect(duckdb()) smry <- dbGetQuery(con, glue( "SELECT count(DISTINCT mdl_key) n_models, count(*) n_cells, 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 = "bl model_cell surface") ``` ## Manifest (target output) ```{r} #| label: manifest # 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_birdlife_botw", content_hash = h, stats = list(ds_key = ds_key, ver = ver, n_species = nrow(sp), n_models = smry$n_models, n_cells = smry$n_cells, v_max = smry$v_max, n_parquet = length(pq), source = path_file(botw_gpkg)), force = msens::force_target("ingest_birdlife_botw")) ```