Publish ichthyoplankton to CF NetCDF

One self-documenting file preserving the cruise → site → tow → net → occurrence → size/stage-bin hierarchy

Author

CalCOFI

Published

2026-07-28

1 Overview

Publish the SWFSC ichthyoplankton dataset as a single self-documenting netCDF-4 file, downstream of release_database in the pipeline.

Two things distinguish this from a table dump, and they are the whole reason the step exists:

1 · The hierarchy is preserved, not flattened. Ichthyoplankton data nests: a site occupation holds tows, a tow holds nets, a net yields taxon occurrences, and each occurrence resolves into length- and stage-frequency bins. Effort (volume filtered, standard haul factor) is measured per net. Flatten that into one table and each net’s volume_sampled repeats onto every one of its size bins — 76,512 real values become 369,978 repeated ones, and any naive SUM() over the flat table silently over-counts effort by ~5×. netCDF-4 groups store each level once and link children to parents by explicit index.

2 · Variables are widened for CF. The database is normalized — one row per observation with measurement_type naming the quantity and measurement_value holding the number — so a single column mixes quantities with different units. CF requires one unit and one standard_name per variable, so each measurement_type becomes its own variable carrying units from metadata/measurement_type.csv.

CF scope — stated honestly

CF conventions cover profile and trajectory geometries cleanly. There is no CF standard for a tow → net → taxon → size-bin hierarchy. This file is therefore CF where CF applies (coordinates, units, standard names, time) and netCDF-4 convention elsewhere (groups, explicit parent_index links). The global attribute cf_scope says exactly this inside the file, so a consumer is never misled into expecting full CF compliance.

1.1 Structure

Code

graph TD
  S[site<br/>61,104] --> T[tow<br/>75,506]
  T --> N["net<br/>76,512<br/><i>+ effort widened to variables</i>"]
  N --> O[occurrence<br/>459,286<br/>taxon × life_stage]
  O --> L[length_bin<br/>241,871]
  O --> G[stage_bin<br/>128,107]

2 Setup

Code

librarian::shelf(DBI, duckdb, dplyr, glue, ncdf4, readr, jsonlite, digest,
                 knitr, tibble, here, quiet = TRUE)
here <- here::here
source(here("libs/publish_netcdf.R"))

DATASET   <- "ichthyo"                 # file/URL identity (no provider prefix)
DATASET_KEY <- "swfsc_ichthyo"         # provenance stamp inside the database
RELEASE   <- cc_release_version()      # the PROMOTED release, never a local tree
PQ        <- cc_release_parquet(RELEASE)
OUT_DIR   <- here("data/netcdf"); dir.create(OUT_DIR, recursive = TRUE, showWarnings = FALSE)
OUT_NC    <- file.path(OUT_DIR, glue("{DATASET}.nc"))
PUBLISH   <- identical(Sys.getenv("CALCOFI_PUBLISH"), "true")  # opt-in upload

cat(glue("release : {RELEASE}\n"))

release : v2026.07.17

Code

cat(glue("parquet : {PQ}\n"))

parquet : https://storage.googleapis.com/calcofi-db/ducklake/releases/v2026.07.17/parquet

Code

cat(glue("output  : {OUT_NC}\n"))

output  : /Users/bbest/Github/CalCOFI/workflows/data/netcdf/ichthyo.nc

Code

cat(glue("publish : {PUBLISH} (set CALCOFI_PUBLISH=true to upload)\n"))

publish : TRUE (set CALCOFI_PUBLISH=true to upload)

Code

con <- dbConnect(duckdb())
for (s in c("INSTALL httpfs", "LOAD httpfs", "SET memory_limit='6GB'",
            "SET enable_progress_bar=false")) try(dbExecute(con, s), silent = TRUE)

# Read the release over HTTPS rather than any local copy: the first NetCDF we
# published was built from the ERDDAP serving tree and silently shipped a
# month-old snapshot. Reading the promoted release makes that failure impossible.
q <- function(sql, ...) dbGetQuery(con, glue(sql, ..., .envir = parent.frame()))
dbExecute(con, glue("
  CREATE VIEW smp  AS SELECT * FROM read_parquet('{PQ}/sample.parquet')
                      WHERE dataset_key = '{DATASET_KEY}';
  CREATE VIEW obs  AS SELECT * FROM read_parquet('{PQ}/obs/dataset_key={DATASET_KEY}/data_0.parquet');
  CREATE VIEW oatt AS SELECT * FROM read_parquet('{PQ}/obs_attribute.parquet')
                      WHERE dataset_key = '{DATASET_KEY}';
  CREATE VIEW smea AS SELECT * FROM read_parquet('{PQ}/sample_measurement.parquet')
                      WHERE dataset_key = '{DATASET_KEY}';
  CREATE VIEW tax  AS SELECT * FROM read_parquet('{PQ}/taxon.parquet');
  CREATE VIEW cru  AS SELECT * FROM read_parquet('{PQ}/cruise.parquet');"))

[1] 0

3 Data integrity — checked before anything is written

The release is the source of truth, but it is not assumed to be clean. These checks run before the file is built so a structural problem surfaces here rather than as a silently wrong hierarchy inside a published product.

Code

# obs has no unique NATURAL key: (sample_key, taxon_key, life_stage) repeats for a
# handful of rows with DIFFERENT abundance values, so joining children on that
# composite would fan out. obs_id is the only unique identifier — children are
# therefore linked through it.
dup_obs <- q("
  SELECT sample_key, taxon_key, life_stage, count(*) AS n,
         string_agg(measurement_value::VARCHAR, ' | ') AS values
  FROM obs GROUP BY 1,2,3 HAVING count(*) > 1 ORDER BY n DESC")

# obs_attribute carries no obs_id, so it must be matched on the composite. Rows
# whose combination has no abundance record are ORPHANS — they are kept, with a
# missing parent index, rather than dropped, because silently discarding measured
# size bins would be worse than admitting they are unparented.
orphan_n <- q("
  SELECT count(*) AS n FROM oatt a
  LEFT JOIN (SELECT DISTINCT sample_key, taxon_key, life_stage FROM obs) o
    ON  o.sample_key = a.sample_key AND o.taxon_key = a.taxon_key
    AND coalesce(o.life_stage,'') = coalesce(a.life_stage,'')
  WHERE o.sample_key IS NULL")$n

cat(glue("duplicate (sample,taxon,life_stage) keys in obs : {nrow(dup_obs)}\n"))

duplicate (sample,taxon,life_stage) keys in obs : 3

Code

cat(glue("obs_attribute rows with no parent occurrence     : {orphan_n}\n"))

obs_attribute rows with no parent occurrence     : 3186

Code

if (nrow(dup_obs)) kable(dup_obs, caption = "Duplicate occurrence keys (release data)")

Duplicate occurrence keys (release data)
sample_key	taxon_key	life_stage	n	values
swfsc_ichthyo:net:1ae84df2-02eb-ef11-a09e-2cea7fa0979c	worms:127292	larva	2	1.0 \| 6.0
swfsc_ichthyo:net:cc2c47f2-02eb-ef11-a09e-2cea7fa0979c	worms:125561	larva	2	7.0 \| 2.0
swfsc_ichthyo:net:eb2b47f2-02eb-ef11-a09e-2cea7fa0979c	worms:125561	larva	2	11.0 \| 1.0

Known release-data anomalies, carried not hidden

As of v2026.07.17 the ichthyo slice contains 3 duplicate occurrence keys (same net, taxon and life stage, different abundance) and 3,186 orphan attribute rows. Neither is introduced here. They are handled explicitly: duplicates keep distinct obs_ids and both rows are published; orphan bins are published with parent_index = -1 and counted in the global attribute n_orphan_attributes. Nothing is dropped, and nothing is silently merged.

4 Build the levels

Code

# ---- site / tow / net: the sample adjacency, one row per level ----------------
site <- q("SELECT sample_key, cruise_key, grid_key, latitude, longitude, datetime
           FROM smp WHERE sample_type='site' ORDER BY sample_key")
tow  <- q("SELECT sample_key, parent_sample_key, tow_type, depth_min_m, depth_max_m,
                  latitude, longitude, datetime
           FROM smp WHERE sample_type='tow' ORDER BY sample_key")
net  <- q("SELECT sample_key, parent_sample_key FROM smp WHERE sample_type='net'
           ORDER BY sample_key")

# ---- effort: widen sample_measurement into NET variables ---------------------
# This is the CF widening at the event level: each measurement_type becomes its
# own column with its own units, one value per net, stored ONCE.
eff_types <- q("SELECT DISTINCT measurement_type FROM smea ORDER BY 1")$measurement_type
eff_sel <- paste(sprintf(
  "MAX(measurement_value) FILTER (WHERE measurement_type='%s') AS \"%s\"",
  eff_types, eff_types), collapse = ",\n    ")
net_eff <- q("SELECT sample_key, {eff_sel} FROM smea GROUP BY sample_key")
net <- net |> left_join(net_eff, by = "sample_key")

# ---- occurrence: one row per obs_id ------------------------------------------
occ <- q("SELECT o.obs_id, o.sample_key, o.taxon_key, o.life_stage,
                 o.measurement_value AS abundance, o.measurement_qual,
                 t.scientific_name, t.worms_id
          FROM obs o LEFT JOIN tax t USING (taxon_key)
          ORDER BY o.sample_key, o.taxon_key, o.life_stage, o.obs_id")

# ---- bins: split by measurement_type so each has ONE unit --------------------
# body_length bins are millimetres; stage bins are an ordinal code. Keeping them
# in one variable would reintroduce exactly the mixed-units problem that makes the
# normalized form un-CF-able.
bin_sql <- "
  SELECT a.obs_attribute_id, a.sample_key, a.taxon_key, a.life_stage,
         a.bin_value, a.bin_label, a.count
  FROM oatt a WHERE a.measurement_type = '{mt}'
  ORDER BY a.sample_key, a.taxon_key, a.life_stage, a.bin_value"
len_bin <- q(bin_sql, mt = "body_length")
stg_bin <- q(bin_sql, mt = "stage")

tibble::tibble(
  level = c("site","tow","net","occurrence","length_bin","stage_bin"),
  rows  = c(nrow(site), nrow(tow), nrow(net), nrow(occ), nrow(len_bin), nrow(stg_bin))
) |> kable()

level	rows
site	61104
tow	75506
net	76512
occurrence	459286
length_bin	241871
stage_bin	128107

Code

# Parent links as 1-based indices into the parent level. match() on the parent's
# ORDERED key vector is what turns "these rows share a value" into an explicit,
# checkable pointer.
tow$parent_index <- match(tow$parent_sample_key, site$sample_key)
net$parent_index <- match(net$parent_sample_key, tow$sample_key)
occ$parent_index <- match(occ$sample_key,        net$sample_key)

# occurrence key for attaching bins; duplicates resolve to the FIRST obs_id, which
# is why the duplicate count is reported above rather than hidden
occ_key <- paste(occ$sample_key, occ$taxon_key, occ$life_stage, sep = "|")
bin_parent <- function(b) {
  idx <- match(paste(b$sample_key, b$taxon_key, b$life_stage, sep = "|"), occ_key)
  idx[is.na(idx)] <- -1L          # orphan: explicit missing, never dropped
  as.integer(idx)
}
len_bin$parent_index <- bin_parent(len_bin)
stg_bin$parent_index <- bin_parent(stg_bin)

stopifnot(!anyNA(tow$parent_index), !anyNA(net$parent_index), !anyNA(occ$parent_index))
cat(glue("orphan length bins: {sum(len_bin$parent_index == -1)}\n"))

orphan length bins: 40

Code

cat(glue("orphan stage bins : {sum(stg_bin$parent_index == -1)}\n"))

orphan stage bins : 3146

5 Write the netCDF-4 file

Code

mt_meta <- cc_measurement_meta(here())

# ncdf4 exposes no group API, but a slash-separated variable name creates a REAL
# netCDF-4 group (verified with ncdump and h5dump, not just by reading it back).
d_site <- ncdim_def("site_n", "", seq_len(nrow(site)),    create_dimvar = FALSE)
d_tow  <- ncdim_def("tow_n",  "", seq_len(nrow(tow)),     create_dimvar = FALSE)
d_net  <- ncdim_def("net_n",  "", seq_len(nrow(net)),     create_dimvar = FALSE)
d_occ  <- ncdim_def("occ_n",  "", seq_len(nrow(occ)),     create_dimvar = FALSE)
d_len  <- ncdim_def("len_n",  "", seq_len(nrow(len_bin)), create_dimvar = FALSE)
d_stg  <- ncdim_def("stg_n",  "", seq_len(nrow(stg_bin)), create_dimvar = FALSE)

drop_key <- function(df, ...) df |> select(-any_of(c("sample_key", "parent_sample_key",
                                                     "parent_index", ...)))
v_site <- nc_level_vars("site",       drop_key(site), d_site, var_meta = mt_meta)
v_tow  <- nc_level_vars("tow",        drop_key(tow),  d_tow,  d_site, tow$parent_index, mt_meta)
v_net  <- nc_level_vars("net",        drop_key(net),  d_net,  d_tow,  net$parent_index, mt_meta)
v_occ  <- nc_level_vars("occurrence", drop_key(occ),  d_occ,  d_net,  occ$parent_index, mt_meta)
v_len  <- nc_level_vars("length_bin", drop_key(len_bin, "obs_attribute_id","taxon_key","life_stage"),
                        d_len, d_occ, len_bin$parent_index, mt_meta)
v_stg  <- nc_level_vars("stage_bin",  drop_key(stg_bin, "obs_attribute_id","taxon_key","life_stage"),
                        d_stg, d_occ, stg_bin$parent_index, mt_meta)

unlink(OUT_NC)
nc <- nc_create(OUT_NC, c(unname(v_site), unname(v_tow), unname(v_net),
                          unname(v_occ), unname(v_len), unname(v_stg)),
                force_v4 = TRUE)
nc_level_put(nc, "site",       drop_key(site), v_site, var_meta = mt_meta)
nc_level_put(nc, "tow",        drop_key(tow),  v_tow,  tow$parent_index, mt_meta, "site")
nc_level_put(nc, "net",        drop_key(net),  v_net,  net$parent_index, mt_meta, "tow")
nc_level_put(nc, "occurrence", drop_key(occ),  v_occ,  occ$parent_index, mt_meta, "net")
nc_level_put(nc, "length_bin", drop_key(len_bin, "obs_attribute_id","taxon_key","life_stage"),
             v_len, len_bin$parent_index, mt_meta, "occurrence")
nc_level_put(nc, "stage_bin",  drop_key(stg_bin, "obs_attribute_id","taxon_key","life_stage"),
             v_stg, stg_bin$parent_index, mt_meta, "occurrence")

CF_SCOPE <- paste(
  "CF-1.10 where CF applies (coordinates, units, standard names, time).",
  "The site/tow/net/occurrence/bin hierarchy uses netCDF-4 groups with explicit",
  "parent_index links; CF defines no feature type for this nesting.")
globals <- list(
  title = "CalCOFI ichthyoplankton (SWFSC) — hierarchical CF NetCDF",
  summary = paste(
    "Ichthyoplankton occurrences with the sampling hierarchy preserved: site ->",
    "tow -> net -> taxon occurrence -> length/stage frequency bin. Event-level",
    "effort (volume_sampled, std_haul_factor, prop_sorted, plankton biomass) is",
    "stored ONCE per net. Children reference parents by 1-based parent_index, so",
    "summing an effort variable over a child group would double-count."),
  Conventions = "CF-1.10, ACDD-1.3", cf_scope = CF_SCOPE,
  institution = "CalCOFI", source = glue("CalCOFI integrated database release {RELEASE}"),
  db_release = RELEASE, dataset_key = DATASET_KEY,
  references = "https://calcofi.io/workflows/publish_ichthyo_to-netcdf.html",
  license = "CC-BY 4.0", creator_name = "CalCOFI", creator_url = "https://calcofi.io",
  date_created = format(Sys.time(), "%Y-%m-%dT%H:%M:%SZ", tz = "UTC"),
  n_duplicate_occurrence_keys = nrow(dup_obs),
  n_orphan_attributes = as.integer(sum(len_bin$parent_index == -1) +
                                   sum(stg_bin$parent_index == -1)),
  anomaly_note = paste(
    "Counts above are anomalies present in the source release, carried through",
    "deliberately: duplicate occurrence keys retain distinct obs_id and both rows",
    "are published; orphan bins carry parent_index = -1. Nothing was dropped."))
for (nm in names(globals)) ncatt_put(nc, 0, nm, globals[[nm]])
nc_close(nc)

cat(glue("wrote {OUT_NC} ({round(file.size(OUT_NC)/1048576, 1)} MB)\n"))

wrote /Users/bbest/Github/CalCOFI/workflows/data/netcdf/ichthyo.nc (170.7 MB)

6 Verify the written file

Round-trip checks: a published file that merely exists is not evidence it is correct. These re-open the file and confirm the hierarchy survived.

Code

ncv <- nc_open(OUT_NC)
grp <- unique(sub("/.*$", "", names(ncv$var)))
chk <- tibble::tibble(
  check = c("groups present", "net count", "occurrence count",
            "effort stored once per net", "parent_index in range"),
  value = c(
    paste(grp, collapse = ", "),
    as.character(ncv$dim[["net_n"]]$len),
    as.character(ncv$dim[["occ_n"]]$len),
    as.character(length(ncvar_get(ncv, "net/volume_sampled"))),
    { pi_ <- ncvar_get(ncv, "occurrence/parent_index")
      as.character(all(pi_ >= 1 & pi_ <= ncv$dim[["net_n"]]$len)) }))
nc_close(ncv)
kable(chk)

check	value
groups present	site, tow, net, occurrence, length_bin, stage_bin
net count	76512
occurrence count	459286
effort stored once per net	76512
parent_index in range	TRUE

Code

# The point of the structure, made concrete: effort summed at the right level vs
# the flat-table answer a naive join would give.
eff_once <- sum(net$volume_sampled, na.rm = TRUE)
eff_flat <- sum(net$volume_sampled[occ$parent_index], na.rm = TRUE)
tibble::tibble(
  `sum(volume_sampled) over nets (correct)` = round(eff_once),
  `... if repeated onto occurrences (flat)` = round(eff_flat),
  `inflation factor` = round(eff_flat / eff_once, 1)) |> kable()

sum(volume_sampled) over nets (correct)	… if repeated onto occurrences (flat)	inflation factor
21054883	181489941	8.6

7 Publish

Versioned by DB release, with bytes written once: if this build is byte-identical to an earlier release, the release folder still gets a manifest and index but the file is not re-uploaded, and storage.calcofi.io redirects the release-scoped .nc URL to the canonical object.

Code

plan <- cc_netcdf_plan(OUT_NC, DATASET, RELEASE)
man  <- cc_netcdf_manifest(plan, DATASET, RELEASE, RELEASE, file.size(OUT_NC),
                           source_tables = c("sample","obs","obs_attribute",
                                             "sample_measurement","taxon"),
                           cf_scope = CF_SCOPE)
cat(glue("sha256      : {substr(plan$sha256, 1, 16)}...\n"))

sha256      : 8caf1fd7430919e2...

Code

cat(glue("identical_to: {plan$identical_to}\n"))

identical_to: NA

Code

cat(glue("upload      : {plan$upload}\n"))

upload      : TRUE

Code

res <- cc_netcdf_publish(OUT_NC, DATASET, RELEASE, plan, man)
cat(glue("published   : {res$url}\n"))

published   : https://storage.calcofi.io/calcofi-files-public/netcdf/ichthyo/v2026.07.17/ichthyo.nc

Code

dbDisconnect(con, shutdown = TRUE)