Ingest CalCOFI Lobster Phyllosoma

Published

2026-08-14

1 Overview

Source: EDI knb-lter-cce.188.4 — spiny lobster (Panulirus interruptus) phyllosoma larvae counts by stage, CalCOFI net tows 1951-2008 (PI Koslow). Provider calcofi.

2 Setup

Code

devtools::load_all(here::here("../calcofi4db"))
devtools::load_all(here::here("../calcofi4r"))
librarian::shelf(
  CalCOFI/calcofi4db, CalCOFI/calcofi4r,
  DBI, dplyr, DT, fs, glue, here, janitor, jsonlite, knitr, lubridate, purrr,
  readr, sf, stringr, tibble, tidyr, units, quiet = T)
options(readr.show_col_types = F); options(DT.options = list(scrollX = TRUE))
source(here("libs/ingest.R"))
cc <- read_calcofi_meta(here("ingest_calcofi_phyllosoma.qmd"))
provider <- cc$provider; dataset <- cc$dataset; tables_owned <- cc$tables_owned
dir_label <- glue("{provider}_{dataset}")
dir_parquet <- here(glue("data/parquet/{dir_label}"))
dir_stage <- cc_stage_path("parquet", dir_label, create = TRUE)
db_path <- here(glue("data/wrangling/{dir_label}.duckdb"))
if (overwrite && file_exists(db_path)) file_delete(db_path)
dir_create(dirname(db_path)); con <- get_duckdb_con(db_path); load_duckdb_extension(con, "spatial")
meas_type_csv <- here("metadata/measurement_type.csv"); d_meas_type <- read_measurement_type(meas_type_csv)

3 Download from EDI

Code

dir_dl <- here(glue("data/cache/{dir_label}")); dir_create(dir_dl)
phyllo_csv <- file.path(dir_dl, "phyllosoma.csv")
edi_url <- "https://pasta.lternet.edu/package/data/eml/knb-lter-cce/188/4/2e66465fa17f78cdc15a2bfc9dce652d"
# Gate on file existence, NOT on `overwrite`. libs/ingest.R states the contract
# plainly -- "downloads are always skipped if files exist regardless of this flag"
# -- but `if (overwrite || ...)` inverted it, and since overwrite defaults to TRUE
# this re-downloaded from EDI on EVERY run. Three pipeline restarts inside ~40
# minutes on 2026-07-30 got us rate-limited: pasta.lternet.edu stayed up (200 at
# the host) while this data object returned 403, and the failed write left the
# cache empty, so the ingest could no longer fall back to a local copy either.
# `overwrite_all` is the flag that means "redo intermediates", so force-refresh
# belongs there.
if (overwrite_all || !file_exists(phyllo_csv)) download.file(edi_url, phyllo_csv, quiet = TRUE)
d_raw <- read_csv(phyllo_csv, col_types = cols(.default = "c")) |> clean_names()
cat(glue("Downloaded {nrow(d_raw)} phyllosoma tows"), "\n")

Downloaded 1859 phyllosoma tows

Code

sync_to_gcs(local_dir = dir_dl, gcs_prefix = glue("archive/{provider}/{dataset}"),
            bucket = "calcofi-files-public", exclude = c(".DS_Store"))

# A tibble: 0 × 4
# ℹ 4 variables: file <chr>, action <chr>, size <dbl>, reason <chr>

4 Build Tow + Measurement Tables

Code

d <- d_raw |>
  mutate(tow_id = row_number(), .before = 1) |>
  mutate(
    date      = suppressWarnings(as.Date(tow_collection_day)),
    latitude  = suppressWarnings(as.numeric(latitude_o)),
    longitude = suppressWarnings(as.numeric(longitude_o)),
    line      = suppressWarnings(as.numeric(station_line)),
    station   = suppressWarnings(as.numeric(station_number)),
    site_key  = if_else(is.na(line)|is.na(station), NA_character_, sprintf("%05.1f %05.1f", line, station)),
    datetime_start_utc = as_datetime(date))

d_tow <- d |>
  transmute(
    tow_id, cruise_orig = year_month_of_tow, ship_name = ship, sorting_lab,
    date, datetime_start_utc, line, station, site_key, latitude, longitude,
    max_tow_depth_m = suppressWarnings(as.numeric(max_tow_depth_m)),
    volume_filtered = suppressWarnings(as.numeric(volume_water_filtered_ml_1000m3)),
    aliquot_pct = suppressWarnings(as.numeric(aliquot_percent)),
    aliquot_adjustment = aliquot_adjustment_value, study_flag)
dbWriteTable(con, "phyllosoma_tow", d_tow, overwrite = TRUE)

stage_cols <- c("total_phyllosoma", paste0("stage_", 1:11))
new_names  <- c("total_phyllosoma", paste0("phyllosoma_stage_", 1:11))
d_meas <- d |>
  select(tow_id, all_of(stage_cols)) |>
  rename_with(~ new_names, all_of(stage_cols)) |>
  pivot_longer(all_of(new_names), names_to = "measurement_type", values_to = "measurement_value") |>
  mutate(measurement_value = suppressWarnings(as.double(measurement_value)),
         measurement_qual = NA_character_) |>
  filter(!is.na(measurement_value)) |>
  mutate(phyllosoma_measurement_id = row_number(), .before = 1)
dbWriteTable(con, "phyllosoma_measurement", d_meas, overwrite = TRUE)
cat(glue("phyllosoma_tow {nrow(d_tow)}, phyllosoma_measurement {nrow(d_meas)}"), "\n")

phyllosoma_tow 1859, phyllosoma_measurement 22308

5 Resolve Keys + Spatial

Code

load_prior_tables(con, parquet_dir = cc_stage_path("parquet", "swfsc_ichthyo"),
                  tables = c("ship","cruise","grid"), geom_tables = c("grid"), as_view = TRUE)

# A tibble: 3 × 3
  table   rows has_geom
  <chr>  <dbl> <lgl>   
1 cruise   691 FALSE   
2 grid     218 TRUE    
3 ship      48 FALSE

Code

# source "Ship" is a NODC ship code (e.g. 31PT), not a name -> join on ship_nodc
d_ship <- dbGetQuery(con, "SELECT ship_key, ship_nodc FROM ship")
zt <- dbGetQuery(con, "SELECT tow_id, ship_name AS ship_nodc, date FROM phyllosoma_tow") |>
  left_join(d_ship, by = "ship_nodc") |>
  mutate(cruise_key = if_else(is.na(ship_nodc)|is.na(date), NA_character_,
                              as.character(glue("{format(date,'%Y-%m')}-{ship_nodc}"))))
valid_ck <- dbGetQuery(con, "SELECT DISTINCT cruise_key FROM cruise")$cruise_key
zt <- zt |> mutate(cruise_key = if_else(cruise_key %in% valid_ck, cruise_key, NA_character_))
dbWriteTable(con, "zk", zt |> select(tow_id, ship_key, cruise_key), overwrite = TRUE)
dbExecute(con, "ALTER TABLE phyllosoma_tow ADD COLUMN IF NOT EXISTS ship_key VARCHAR")

[1] 0

Code

dbExecute(con, "ALTER TABLE phyllosoma_tow ADD COLUMN IF NOT EXISTS cruise_key VARCHAR")

[1] 0

Code

dbExecute(con, "UPDATE phyllosoma_tow t SET ship_key=z.ship_key, cruise_key=z.cruise_key FROM zk z WHERE t.tow_id=z.tow_id")

[1] 1859

Code

dbExecute(con, "DROP TABLE zk")

[1] 0

Code

cat(glue("ship {sum(!is.na(zt$ship_key))}/{nrow(zt)}; cruise_key {sum(!is.na(zt$cruise_key))}/{nrow(zt)}"), "\n")

ship 1816/1859; cruise_key 1634/1859

Code

add_point_geom(con, "phyllosoma_tow", lon_col = "longitude", lat_col = "latitude")
assign_grid_key(con, "phyllosoma_tow") |> datatable(caption = "Grid assignment")

6 Measurement Types + Finalize

7 Emit Core Tables

Project this dataset into the shared consolidated core model (design_env-bio-consolidation.md). These core tables are this ingest’s output: release_database.qmd concatenates the per-dataset shards rather than re-deriving the core, so there is exactly one projection to keep correct.

The source is one column per stage plus a total. Only total_phyllosoma is an occurrence headline; the eleven phyllosoma_stage_N columns are a stage distribution within that occurrence, so they become obs_attribute rows (measurement_type = 'stage', bin_value = N) rather than eleven separate observations. metadata/measurement_taxon.csv supplies the split and the taxon_key for Panulirus interruptus.

Code

ds_key <- "calcofi_phyllosoma"
# filter the crosswalk to THIS dataset -- an unfiltered read leaks other datasets'
# taxa into this shard (the retired emit_core_tables() wrapper filtered internally)
mt_taxon <- read_csv(here("metadata/measurement_taxon.csv"),
                     col_types = cols(worms_id = "i", itis_id = "i",
                                      bin_value = "d", .default = "c")) |>
  filter(dataset_key == ds_key)
tx_over  <- read_csv(here("metadata/taxon_override.csv"), show_col_types = FALSE)

# This projection lives here, in the notebook that owns the dataset, not in a
# switch(dataset_key, ...) arm inside calcofi4db. The reusable SHAPES stay in the
# package (sample_arm_self / compat_measurement_sql), so this is a declaration.

# cross-reference: resolve each taxon against BOTH authorities (cached in
# metadata/taxon_xref.csv, so a re-run costs no API calls). This fills the
# `worms_id` COLUMN on itis:-keyed taxa without touching their key — a consumer
# joining on worms_id used to match ZERO rows for every seabird and marine
# mammal — backfills `itis_id` the other way, replaces an id its authority has
# deprecated so the key is always an accepted id, and fetches the real
# `taxonomic_status` with the date it was checked. Must precede the lineage
# fetch, which should ask about the accepted id, not the deprecated one.
ensure_taxon_xref(con, mt_taxon, tx_over,
                  cache_csv = here("metadata/taxon_xref.csv"))

# lineage: fetch each taxon's WoRMS/ITIS classification (cached in
# metadata/taxon_lineage.csv, so a re-run costs no API calls) and stage it as the
# `taxon` hierarchy build_taxon_reference() reads. Without it a crosswalk- or
# vocabulary-resolved taxon reaches the release with a key and a name and NOTHING
# else — no rank, no parent_taxon_key, no classification — so hierarchy rollups
# ("all Decapoda") silently match nothing and no error is raised anywhere.
ensure_taxon_lineage(con, mt_taxon, tx_over,
                     cache_csv = here("metadata/taxon_lineage.csv"))
n_taxon    <- build_taxon_reference(con, mt_taxon, tx_over)
n_ds_taxon <- build_dataset_taxon(con,   mt_taxon, tx_over)
# stage the crosswalk WITH its derived taxon_key — never dbWriteTable() the raw
# CSV, which has no taxon_key column at all (and a 'worms:' || worms_id string
# built inline would mis-key any ITIS-resolved taxon)
ensure_measurement_taxon(con, mt_taxon, dataset_key = ds_key)

append_sample(con, sample_arm_self(
  ds_key, "phyllosoma_tow", "tow_id", "tow", site_expr = "site_key",
  depth_min = "0::DOUBLE", depth_max = "max_tow_depth_m"))

# obs — the occurrence headline is ONLY total_phyllosoma (target='obs'); the
# INNER join to the crosswalk is what excludes the per-stage columns
append_obs(con, glue("
  SELECT 'bio', '{ds_key}', {ns_key(ds_key, 'tow', 'tw.tow_id')},
         tw.grid_key, tw.cruise_key, tw.latitude, tw.longitude,
         CAST(tw.datetime_start_utc AS TIMESTAMP), 0::DOUBLE, tw.max_tow_depth_m,
         mx.taxon_key, mx.life_stage, mx.measurement_type, m.measurement_value,
         m.measurement_qual, NULL::DOUBLE
  FROM phyllosoma_measurement m JOIN phyllosoma_tow tw USING (tow_id)
  JOIN _measurement_taxon mx ON mx.dataset_key = '{ds_key}'
                            AND mx.raw_measurement_type = m.measurement_type
                            AND mx.target = 'obs'"))

# obs_attribute — stage frequency: phyllosoma_stage_N -> ('stage', bin_value=N,
# count). Zero-count stages are dropped: a stage column of 0 within a tow that
# did catch phyllosoma is an absence already implied by the headline.
append_obs_attribute(con, glue("
  SELECT '{ds_key}', {ns_key(ds_key, 'tow', 'tw.tow_id')},
         mx.taxon_key, mx.life_stage, mx.measurement_type, mx.bin_value,
         NULL::VARCHAR, CAST(m.measurement_value AS INTEGER), NULL::VARCHAR
  FROM phyllosoma_measurement m JOIN phyllosoma_tow tw USING (tow_id)
  JOIN _measurement_taxon mx ON mx.dataset_key = '{ds_key}'
                            AND mx.raw_measurement_type = m.measurement_type
                            AND mx.target = 'attribute'
  WHERE m.measurement_value > 0"))

core <- list(
  sample        = dbGetQuery(con, "SELECT COUNT(*) FROM sample")[[1]],
  obs           = dbGetQuery(con, "SELECT COUNT(*) FROM obs")[[1]],
  obs_attribute = dbGetQuery(con, "SELECT COUNT(*) FROM obs_attribute")[[1]],
  taxon         = n_taxon,
  dataset_taxon = n_ds_taxon)
cat(glue(
  "core projection — sample={core$sample %||% 0} obs={core$obs %||% 0} ",
  "obs_attribute={core$obs_attribute %||% 0} ",
  "taxon={core$taxon %||% 0} dataset_taxon={core$dataset_taxon %||% 0}\n"))

core projection — sample=1859 obs=1859 obs_attribute=369 taxon=13 dataset_taxon=1

Code

# the headline must be the total only — one row per tow — and the
# stage bins must not leak into obs
n_obs <- core$obs
n_exp <- dbGetQuery(con, "
  SELECT COUNT(*) FROM phyllosoma_measurement m JOIN phyllosoma_tow tw USING (tow_id)
  WHERE m.measurement_type = 'total_phyllosoma'")[[1]]
n_stage_in_obs <- dbGetQuery(con,
  "SELECT COUNT(*) FROM obs WHERE measurement_type LIKE 'phyllosoma_stage%'")[[1]]
d_bins <- dbGetQuery(con, "
  SELECT MIN(bin_value) AS lo, MAX(bin_value) AS hi, COUNT(DISTINCT bin_value) AS n
  FROM obs_attribute WHERE measurement_type = 'stage'")
stopifnot(
  "obs must carry the total only, one row per tow" = n_obs == n_exp,
  "stage columns must go to obs_attribute, not obs"             = n_stage_in_obs == 0,
  "phyllosoma obs must all carry a taxon_key" =
    dbGetQuery(con, "SELECT COUNT(*) FROM obs WHERE taxon_key IS NULL")[[1]] == 0,
  "every obs.taxon_key must RESOLVE in taxon (not merely be non-NULL)" =
    dbGetQuery(con, "SELECT COUNT(*) FROM obs o LEFT JOIN taxon t USING (taxon_key)
                     WHERE t.taxon_key IS NULL")[[1]] == 0,
  "stage bins must span 1-11"  = d_bins$lo == 1 && d_bins$hi == 11,
  "every obs.sample_key must resolve in sample" =
    dbGetQuery(con, "SELECT COUNT(*) FROM obs o LEFT JOIN sample s USING (sample_key)
                     WHERE s.sample_key IS NULL")[[1]] == 0)
cat(glue("obs parity: {format(n_obs, big.mark = ',')} headline rows; ",
         "{format(core$obs_attribute, big.mark = ',')} stage bins ",
         "across {d_bins$n} stages ({d_bins$lo}-{d_bins$hi})"), "\n")

obs parity: 1,859 headline rows; 369 stage bins across 7 stages (1-11)

Code

# serve the retired per-dataset table names as VIEWs over the core, so
# in-notebook consumers and ad-hoc queries keep working against the old names
# (exact for every column the core models)
invisible(dbExecute(con, "DROP TABLE IF EXISTS phyllosoma_measurement_src"))
invisible(dbExecute(con, "ALTER TABLE phyllosoma_measurement RENAME TO phyllosoma_measurement_src"))
invisible(dbExecute(con, glue(
  "CREATE OR REPLACE VIEW phyllosoma_measurement AS
   {compat_measurement_sql(ds_key, 'tow', 'tow_id', 'phyllosoma_measurement_id')}")))
cat(glue("compat view phyllosoma_measurement over obs: ",
         "{dbGetQuery(con, 'SELECT COUNT(*) FROM phyllosoma_measurement')[[1]]} rows"), "\n")

compat view phyllosoma_measurement over obs: 1859 rows

Code

phyllo_types <- tibble(
  measurement_type = new_names,
  description = c("Total phyllosoma larvae count", paste0("Phyllosoma stage ", 1:11, " count")),
  units = "count", `_source_column` = new_names, `_source_table` = "phyllosoma_measurement",
  `_source_datasets` = "calcofi_phyllosoma", `_qual_column` = NA_character_, `_prec_column` = NA_character_)
new_types <- phyllo_types |> filter(!measurement_type %in% d_meas_type$measurement_type)
if (nrow(new_types) > 0) { d_meas_type <- bind_rows(d_meas_type, new_types); write_csv(d_meas_type, meas_type_csv, na = "") }
dbWriteTable(con, "measurement_type", d_meas_type, overwrite = TRUE)

d_dataset <- ingest_yaml_to_dataset_df(read_ingest_yaml(here()))
dbWriteTable(con, "dataset", d_dataset, overwrite = TRUE)

phyllo_rels <- list(
  primary_keys = list(phyllosoma_tow = "tow_id", phyllosoma_measurement = "phyllosoma_measurement_id",
                      measurement_type = "measurement_type"),
  foreign_keys = list(
    list(table="phyllosoma_measurement", column="tow_id", ref_table="phyllosoma_tow", ref_column="tow_id"),
    list(table="phyllosoma_measurement", column="measurement_type", ref_table="measurement_type", ref_column="measurement_type")))
# the SOURCE shape, documenting the wrangling above; the published tables are
# the consolidated core, so relationships.json comes from core_relationships()
cc_erd(con, tables = c("phyllosoma_tow","phyllosoma_measurement","measurement_type","dataset"), rels = phyllo_rels,
       colors = list(lightblue = c("phyllosoma_tow","phyllosoma_measurement"), lightyellow = "measurement_type", white = "dataset"))

Code

results <- validate_for_release(con, checks = "all", strict = FALSE)
cat("Validation:", ifelse(results$passed, "PASSED", "FAILED"), "\n")

Validation: FAILED

Code

dir_create(dir_parquet)
tbls_out <- core_output_tables(con, extra = c("measurement_type", "dataset"))
write_parquet_outputs(con = con, output_dir = dir_parquet,
  tables = tbls_out,
  sort_by = list(obs = c("grid_key", "measurement_type")), strip_provenance = FALSE)

# A tibble: 7 × 5
  table             rows file_size path                              partitioned
  <chr>            <dbl>     <dbl> <chr>                             <lgl>      
1 sample            1859     26453 /Users/bbest/_big/calcofi/parque… FALSE      
2 obs               1859     35569 /Users/bbest/_big/calcofi/parque… FALSE      
3 obs_attribute      369      3287 /Users/bbest/_big/calcofi/parque… FALSE      
4 taxon               13      4055 /Users/bbest/_big/calcofi/parque… FALSE      
5 dataset_taxon        1      1324 /Users/bbest/_big/calcofi/parque… FALSE      
6 measurement_type   200     12063 measurement_type.parquet          FALSE      
7 dataset             16     11099 dataset.parquet                   FALSE

Code

build_relationships_json(rels = core_relationships(tbls_out), output_dir = dir_parquet,
  provider = provider, dataset = dataset)

[1] "/Users/bbest/Github/CalCOFI/workflows/data/parquet/calcofi_phyllosoma/relationships.json"

Code

d_tbls_rd <- read_csv(here("metadata/calcofi/phyllosoma/tbls_redefine.csv"))
d_flds_rd <- read_csv(here("metadata/calcofi/phyllosoma/flds_redefine.csv"))
build_metadata_json(con = con, d_tbls_rd = d_tbls_rd, d_flds_rd = d_flds_rd,
  metadata_derived_csv = c(here("metadata/core_dictionary.csv"),
                           here("metadata/calcofi/phyllosoma/metadata_derived.csv")),
  output_dir = dir_parquet, tables = tbls_out,
  set_comments = TRUE, provider = provider, dataset = dataset,
  workflow_url = cc$workflow_url, tables_owned = tables_owned)

metadata.json documentation gaps (7 tables, 102 columns) — these render blank in cc_describe_table() / cc_db_catalog():
tables with no description_md: 2    measurement_type, dataset
columns with no description_md: 32    dataset.provider, dataset.dataset, dataset.dataset_name, dataset.dataset_name_short, dataset.category, dataset.color, dataset.description, dataset.citation_main, dataset.citation_others, dataset.link_calcofi_org, dataset.link_data_source, dataset.link_others (+20 more)
measurement columns with no units: 26    dataset.dataset_name_short, dataset.category, dataset.color, dataset.citation_main, dataset.citation_others, dataset.link_calcofi_org, dataset.link_data_source, dataset.link_others, dataset.tables, dataset.coverage_temporal, dataset.coverage_spatial, dataset.license (+14 more)
  backfill via metadata/{provider}/{dataset}/flds_redefine.csv, then re-run

[1] "/Users/bbest/Github/CalCOFI/workflows/data/parquet/calcofi_phyllosoma/metadata.json"

Code

sync_to_gcs(local_dir = dir_stage, sidecar_dir = dir_parquet, gcs_prefix = glue("ingest/{dir_label}"), bucket = "calcofi-db")

# A tibble: 5 × 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

8 Questions for Data Providers

Code

# one validated read + render for every ingest: the vocabulary and the column
# order live in calcofi4db, not in 16 hand-written factor() calls
questions_datatable(
  here(cc$questions_file),
  caption = "Questions for the CalCOFI phyllosoma data providers (ranked)")

Code

close_duckdb(con); cat(glue("Parquet outputs written to: {dir_parquet}"), "\n")

Parquet outputs written to: /Users/bbest/Github/CalCOFI/workflows/data/parquet/calcofi_phyllosoma