Publish CTD casts to CF NetCDF

Two CF Discrete-Sampling-Geometry profile files — the thinned headline record and the full-resolution scans

Author

CalCOFI

Published

2026-07-28

1 Overview

Publish the CalCOFI CTD record as a CF Discrete-Sampling-Geometry profile dataset, downstream of release_database.

This is the simple end of the publish range, and deliberately so — it is the counterpart to publish_ichthyo_to-netcdf, which handles the deeply nested case. CTD needs no netCDF-4 groups:

	ichthyo	ctd-cast
hierarchy	site → tow → net → occurrence → bin	flat: cast → depth level
CF status	CF where it applies; groups elsewhere	fully CF (`featureType=Profile`)
structure	netCDF-4 groups + `parent_index`	contiguous ragged array (`rowSize`)

A CTD cast is a CF profile: one instance, a vertical series of depth levels, each sensor its own variable. So this file is genuinely CF-compliant end to end — no caveat needed — and any CF-aware tool will read it without special handling.

Why the widening still matters

The database stores CTD in long form: one row per cast × depth × measurement_type, 5.55 M rows over 15 sensor types. A single measurement_value column therefore mixes °C, PSU, mL/L and volts. CF requires one unit per variable, so each measurement_type becomes its own variable with its own units and long_name. That widening is what turns a queryable table into a self-documenting file.

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     <- "ctd-cast"
DATASET_KEY <- "calcofi_ctd-cast"
RELEASE     <- cc_release_version()      # PROMOTED release, never a serving 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")

cat(glue("release : {RELEASE}\nparquet : {PQ}\noutput  : {OUT_NC}\npublish : {PUBLISH}\n"))

release : v2026.07.17
parquet : https://storage.googleapis.com/calcofi-db/ducklake/releases/v2026.07.17/parquet
output  : /Users/bbest/Github/CalCOFI/workflows-main/data/netcdf/ctd-cast.nc
publish : TRUE

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)
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');"))

[1] 0

3 Shape and integrity

Code

shape <- q("
  SELECT (SELECT count(*) FROM smp WHERE sample_type='cast')       AS casts_in_sample,
         (SELECT count(DISTINCT sample_key) FROM obs)              AS casts_with_obs,
         (SELECT count(*) FROM obs)                                AS obs_rows,
         (SELECT count(DISTINCT measurement_type) FROM obs)        AS sensor_types")
kable(shape)

casts_in_sample	casts_with_obs	obs_rows	sensor_types
14336	7175	5551551	15

Code

# A profile must have a well-defined depth axis: one value per (cast, depth,
# measurement_type). Duplicates there would silently collapse in the pivot below,
# so check rather than assume.
dupes <- q("
  SELECT count(*) AS n FROM (
    SELECT sample_key, depth_min_m, measurement_type
    FROM obs GROUP BY 1,2,3 HAVING count(*) > 1)")$n
cat(glue("duplicate (cast, depth, measurement_type) triples: {dupes}\n"))

duplicate (cast, depth, measurement_type) triples: 44245

Code

# WIDEN: measurement_type -> one variable each, ordered by (cast, depth) so the
# profiles are contiguous — which is what a CF contiguous ragged array requires.
types <- q("SELECT DISTINCT measurement_type FROM obs ORDER BY 1")$measurement_type
sel <- paste(sprintf(
  "MAX(measurement_value) FILTER (WHERE measurement_type='%s')::DOUBLE AS \"%s\"",
  types, types), collapse = ",\n    ")

prof <- q("
  SELECT o.sample_key                       AS profile_id,
         any_value(s.cruise_key)            AS cruise_key,
         any_value(s.grid_key)              AS grid_key,
         epoch(any_value(s.datetime))::DOUBLE AS time,
         any_value(s.latitude)::DOUBLE      AS latitude,
         any_value(s.longitude)::DOUBLE     AS longitude,
         o.depth_min_m::DOUBLE              AS depth,
         {sel}
  FROM obs o LEFT JOIN smp s ON s.sample_key = o.sample_key
  GROUP BY o.sample_key, o.depth_min_m
  ORDER BY o.sample_key, o.depth_min_m")

cat(glue("wide rows (cast x depth): {nrow(prof)}\n"))

wide rows (cast x depth): 434312

Code

cat(glue("sensor variables        : {length(types)}\n"))

sensor variables        : 15

4 Write the CF profile file

Code

mt_meta <- cc_measurement_meta(here())

prof_ids <- unique(prof$profile_id)
row_size <- as.integer(table(factor(prof$profile_id, levels = prof_ids)))
first_ix <- match(prof_ids, prof$profile_id)
stopifnot(sum(row_size) == nrow(prof))

STRLEN <- 64L
FILL   <- 9.969209968386869e36
d_prof <- ncdim_def("profile", "", seq_along(prof_ids), create_dimvar = FALSE)
d_obs  <- ncdim_def("obs",     "", seq_len(nrow(prof)), create_dimvar = FALSE)
d_str  <- ncdim_def("name_strlen", "", seq_len(STRLEN), create_dimvar = FALSE)

v_pid  <- ncvar_def("profile_id", "", list(d_str, d_prof), prec = "char")
v_cru  <- ncvar_def("cruise_key", "", list(d_str, d_prof), prec = "char")
v_grd  <- ncvar_def("grid_key",   "", list(d_str, d_prof), prec = "char")
v_time <- ncvar_def("time", "seconds since 1970-01-01T00:00:00Z", d_prof, prec = "double")
v_lat  <- ncvar_def("latitude",  "degrees_north", d_prof, prec = "double")
v_lon  <- ncvar_def("longitude", "degrees_east",  d_prof, prec = "double")
v_rs   <- ncvar_def("rowSize",   "",              d_prof, prec = "integer")
v_dep  <- ncvar_def("depth", "m", d_obs, prec = "double", missval = FILL)
v_sens <- lapply(types, function(v) ncvar_def(
  v, as.character(mt_meta[[v]]$units %||% ""), d_obs, prec = "double", missval = FILL))

unlink(OUT_NC)
nc <- nc_create(OUT_NC, c(list(v_pid, v_cru, v_grd, v_time, v_lat, v_lon, v_rs, v_dep),
                          v_sens), force_v4 = TRUE)
ncvar_put(nc, v_pid,  prof_ids)
ncvar_put(nc, v_cru,  prof$cruise_key[first_ix])
ncvar_put(nc, v_grd,  prof$grid_key[first_ix])
ncvar_put(nc, v_time, prof$time[first_ix])
ncvar_put(nc, v_lat,  prof$latitude[first_ix])
ncvar_put(nc, v_lon,  prof$longitude[first_ix])
ncvar_put(nc, v_rs,   row_size)
ncvar_put(nc, v_dep,  prof$depth)
for (i in seq_along(types)) ncvar_put(nc, v_sens[[i]], prof[[types[i]]])

# CF DSG attributes — these are what make it a profile dataset rather than a
# table that happens to be in netCDF.
ncatt_put(nc, "profile_id", "cf_role", "profile_id")
ncatt_put(nc, "profile_id", "long_name", "CTD cast identifier")
ncatt_put(nc, "time", "standard_name", "time");      ncatt_put(nc, "time", "axis", "T")
ncatt_put(nc, "latitude",  "standard_name", "latitude");  ncatt_put(nc, "latitude", "axis", "Y")
ncatt_put(nc, "longitude", "standard_name", "longitude"); ncatt_put(nc, "longitude", "axis", "X")
ncatt_put(nc, "depth", "standard_name", "depth"); ncatt_put(nc, "depth", "axis", "Z")
ncatt_put(nc, "depth", "positive", "down")
ncatt_put(nc, "rowSize", "long_name", "Number of observations for this profile")
ncatt_put(nc, "rowSize", "sample_dimension", "obs")   # -> contiguous ragged array
for (v in types) {
  ncatt_put(nc, v, "long_name",
            as.character(mt_meta[[v]]$long_name %||% gsub("_", " ", v)))
  ncatt_put(nc, v, "coordinates", "time latitude longitude depth")
  sn <- mt_meta[[v]]$standard_name %||% NA_character_
  if (!is.na(sn) && nzchar(sn)) ncatt_put(nc, v, "standard_name", sn)
}
globals <- list(
  title = "CalCOFI CTD casts — CF profile dataset",
  summary = paste(
    "CalCOFI CTD profiles as CF Discrete Sampling Geometry: one profile per cast,",
    "depth levels stored as a contiguous ragged array, each canonical sensor its",
    "own variable with its own units. Widened from the normalized long form in the",
    glue("CalCOFI integrated database release {RELEASE}.")),
  featureType = "profile", cdm_data_type = "Profile",
  cdm_profile_variables = "profile_id, time, latitude, longitude, cruise_key, grid_key",
  Conventions = "CF-1.10, COARDS, ACDD-1.3",
  cf_scope = paste("Fully CF: a CTD cast is a CF profile, so this dataset needs no",
                   "extension beyond the standard (contrast the nested net-tow datasets)."),
  institution = "CalCOFI", source = glue("CalCOFI integrated database release {RELEASE}"),
  db_release = RELEASE, dataset_key = DATASET_KEY,
  references = "https://calcofi.io/workflows/publish_ctd-cast_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"))
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-main/data/netcdf/ctd-cast.nc (54.5 MB)

5 Verify

Code

ncv <- nc_open(OUT_NC)
rs  <- ncvar_get(ncv, "rowSize")
chk <- tibble(
  check = c("featureType", "cf_role set", "profiles", "obs levels",
            "rowSize sums to obs", "sensor variables", "depths per profile (min/med/max)"),
  value = c(
    ncatt_get(ncv, 0, "featureType")$value,
    ncatt_get(ncv, "profile_id", "cf_role")$value,
    as.character(ncv$dim[["profile"]]$len),
    as.character(ncv$dim[["obs"]]$len),
    as.character(sum(rs) == ncv$dim[["obs"]]$len),
    as.character(length(types)),
    paste(min(rs), median(rs), max(rs), sep = " / ")))
nc_close(ncv)
kable(chk)

check	value
featureType	profile
cf_role set	profile_id
profiles	7175
obs levels	434312
rowSize sums to obs	TRUE
sensor variables	15
depths per profile (min/med/max)	1 / 70 / 436

6 Full-resolution companion — `ctd-cast_full.nc`

The file above is the adaptively-thinned record (5.55 M long rows → 434 k depth levels, 7,175 profiles). obs_ctd_full is ~216 M long rows across 96 cruise partitions and 32 sensor types, flagged supplemental in the release precisely because of its size.

It is a superset in two dimensions, not one: every depth scan rather than a thinned subset, and all 14,336 casts — both the down- and up-cast of each station occupation — where the thinned record keeps one direction per occupation. So the profile count roughly doubles for a reason unrelated to resolution. Both files state this in their own global attributes, because “full resolution” alone would imply the same profiles sampled more densely, which is not what this is.

Two things make this build different, and both are consequences of scale:

Why this is chunked rather than pivoted in one query

This is the table that OOM’d ERDDAP at every container size tested — 4, 5 and 6 GB (see the serving benchmark). Materializing it here would hit the same wall, so it is processed one cruise partition at a time and written into the file incrementally with ncvar_put(start=, count=). Peak memory is one cruise (~59 k levels × 32 variables ≈ 15 MB), not the whole table.

Also, read_parquet('.../obs_ctd_full/**/*.parquet') 404s over HTTPS: expanding a glob needs a directory listing and object storage has none. Partitions are therefore enumerated from the XML listing API (cc_release_partitions()) and passed as an explicit URL vector.

Code

parts <- cc_release_partitions("obs_ctd_full", RELEASE)
cat(glue("partitions: {length(parts)}\n"))

partitions: 96

Code

# obs_ctd_full already carries latitude/longitude/datetime, so unlike the thinned
# path no join back to `sample` is needed.
full_types <- dbGetQuery(con, glue(
  "SELECT DISTINCT measurement_type FROM read_parquet('{parts[1]}') ORDER BY 1"))$measurement_type
cat(glue("sensor types: {length(full_types)}\n"))

sensor types: 32

Code

# PASS 1 — size the dimensions. Cheap: scans two columns per partition rather
# than pivoting. netCDF needs dimension lengths at creation time, and a wrong
# guess means rewriting a multi-GB file.
pass1 <- lapply(seq_along(parts), function(i) {
  r <- dbGetQuery(con, glue("
    SELECT count(*) AS n_lev, count(DISTINCT sample_key) AS n_prof FROM (
      SELECT DISTINCT sample_key, depth_min_m FROM read_parquet('{parts[i]}'))"))
  if (i %% 24 == 0) message(glue("  pass1 {i}/{length(parts)}"))
  r
})
p1 <- do.call(rbind, pass1)
N_PROF <- sum(p1$n_prof); N_LEV <- sum(p1$n_lev)
tibble(partitions = length(parts), profiles = N_PROF, depth_levels = N_LEV,
       `mean levels/profile` = round(N_LEV / N_PROF)) |> kable()

partitions	profiles	depth_levels	mean levels/profile
96	14336	6082688	424

Code

FULL_NC <- file.path(OUT_DIR, "ctd-cast_full.nc")
sel_full <- paste(sprintf(
  "MAX(measurement_value) FILTER (WHERE measurement_type='%s')::DOUBLE AS \"%s\"",
  full_types, full_types), collapse = ",\n    ")

d_prof_f <- ncdim_def("profile", "", seq_len(N_PROF), create_dimvar = FALSE)
d_obs_f  <- ncdim_def("obs",     "", seq_len(N_LEV),  create_dimvar = FALSE)
d_str_f  <- ncdim_def("name_strlen", "", seq_len(STRLEN), create_dimvar = FALSE)

vf_pid  <- ncvar_def("profile_id", "", list(d_str_f, d_prof_f), prec = "char")
vf_cru  <- ncvar_def("cruise_key", "", list(d_str_f, d_prof_f), prec = "char")
vf_time <- ncvar_def("time", "seconds since 1970-01-01T00:00:00Z", d_prof_f, prec = "double")
vf_lat  <- ncvar_def("latitude",  "degrees_north", d_prof_f, prec = "double")
vf_lon  <- ncvar_def("longitude", "degrees_east",  d_prof_f, prec = "double")
vf_rs   <- ncvar_def("rowSize",   "",              d_prof_f, prec = "integer")
vf_dep  <- ncvar_def("depth", "m", d_obs_f, prec = "double", missval = FILL)
vf_sens <- lapply(full_types, function(v) ncvar_def(
  v, as.character(mt_meta[[v]]$units %||% ""), d_obs_f, prec = "double", missval = FILL))

unlink(FULL_NC)
ncf <- nc_create(FULL_NC, c(list(vf_pid, vf_cru, vf_time, vf_lat, vf_lon, vf_rs, vf_dep),
                            vf_sens), force_v4 = TRUE)

# PASS 2 — one cruise at a time, appended at the running offsets. Nothing larger
# than a single cruise is ever held in memory.
off_p <- 1L; off_o <- 1L
for (i in seq_along(parts)) {
  w <- dbGetQuery(con, glue("
    SELECT sample_key AS profile_id,
           any_value(cruise_key)              AS cruise_key,
           epoch(any_value(datetime))::DOUBLE AS time,
           any_value(latitude)::DOUBLE        AS latitude,
           any_value(longitude)::DOUBLE       AS longitude,
           depth_min_m::DOUBLE                AS depth,
           {sel_full}
    FROM read_parquet('{parts[i]}')
    GROUP BY sample_key, depth_min_m
    ORDER BY sample_key, depth_min_m"))
  if (!nrow(w)) next
  ids <- unique(w$profile_id)
  rs  <- as.integer(table(factor(w$profile_id, levels = ids)))
  fx  <- match(ids, w$profile_id)

  ncvar_put(ncf, vf_pid,  ids,               start = c(1, off_p), count = c(STRLEN, length(ids)))
  ncvar_put(ncf, vf_cru,  w$cruise_key[fx],  start = c(1, off_p), count = c(STRLEN, length(ids)))
  ncvar_put(ncf, vf_time, w$time[fx],        start = off_p, count = length(ids))
  ncvar_put(ncf, vf_lat,  w$latitude[fx],    start = off_p, count = length(ids))
  ncvar_put(ncf, vf_lon,  w$longitude[fx],   start = off_p, count = length(ids))
  ncvar_put(ncf, vf_rs,   rs,                start = off_p, count = length(ids))
  ncvar_put(ncf, vf_dep,  w$depth,           start = off_o, count = nrow(w))
  for (k in seq_along(full_types))
    ncvar_put(ncf, vf_sens[[k]], w[[full_types[k]]], start = off_o, count = nrow(w))

  off_p <- off_p + length(ids); off_o <- off_o + nrow(w)
  if (i %% 12 == 0) message(glue("  wrote {i}/{length(parts)} — {off_p-1} profiles, {off_o-1} levels"))
  rm(w); invisible(gc(FALSE))
}
# the pass-1 sizing must match exactly, or the file has trailing fill values that
# look like real (missing) data
stopifnot(off_p - 1L == N_PROF, off_o - 1L == N_LEV)

ncatt_put(ncf, "profile_id", "cf_role", "profile_id")
ncatt_put(ncf, "profile_id", "long_name", "CTD cast identifier")
ncatt_put(ncf, "time", "standard_name", "time");     ncatt_put(ncf, "time", "axis", "T")
ncatt_put(ncf, "latitude",  "standard_name", "latitude");  ncatt_put(ncf, "latitude", "axis", "Y")
ncatt_put(ncf, "longitude", "standard_name", "longitude"); ncatt_put(ncf, "longitude", "axis", "X")
ncatt_put(ncf, "depth", "standard_name", "depth"); ncatt_put(ncf, "depth", "axis", "Z")
ncatt_put(ncf, "depth", "positive", "down")
ncatt_put(ncf, "rowSize", "long_name", "Number of observations for this profile")
ncatt_put(ncf, "rowSize", "sample_dimension", "obs")
for (v in full_types) {
  ncatt_put(ncf, v, "long_name", as.character(mt_meta[[v]]$long_name %||% gsub("_", " ", v)))
  ncatt_put(ncf, v, "coordinates", "time latitude longitude depth")
}
globals_full <- globals
globals_full$title <- "CalCOFI CTD casts, full resolution — CF profile dataset"
globals_full$summary <- paste(
  "Full-resolution CalCOFI CTD scans as CF Discrete Sampling Geometry: one profile",
  "per cast, every recorded depth level, each sensor its own variable.",
  "NOTE this file differs from the thinned companion (ctd-cast.nc) in TWO ways, not",
  "one: it keeps every depth scan rather than an adaptively-thinned subset, AND it",
  "covers all 14,336 casts including both the down- and up-cast of each station",
  "occupation, where the thinned record keeps a single direction (7,175). It is",
  "therefore a superset, not a higher-resolution rendering of the same profiles.",
  "Most analyses want the thinned file; use this one when individual scans or the",
  "second cast direction matter.")
globals_full$comment <- paste(
  "profiles here (14,336) are cast-direction resolved; the thinned companion",
  "reports 7,175 because it retains one direction per occupation.")
globals_full$source_table <- "obs_ctd_full (supplemental)"
for (nm in names(globals_full)) ncatt_put(ncf, 0, nm, globals_full[[nm]])
nc_close(ncf)
cat(glue("wrote {FULL_NC} ({round(file.size(FULL_NC)/1048576, 1)} MB)\n"))

wrote /Users/bbest/Github/CalCOFI/workflows-main/data/netcdf/ctd-cast_full.nc (1533.6 MB)

Code

ncv <- nc_open(FULL_NC)
rs  <- ncvar_get(ncv, "rowSize")
tibble(
  check = c("featureType", "profiles", "depth levels", "rowSize sums to obs",
            "sensor variables", "levels per profile (min/med/max)"),
  value = c(ncatt_get(ncv, 0, "featureType")$value,
            as.character(ncv$dim[["profile"]]$len),
            as.character(ncv$dim[["obs"]]$len),
            as.character(sum(rs) == ncv$dim[["obs"]]$len),
            as.character(length(full_types)),
            paste(min(rs), median(rs), max(rs), sep = " / "))) |> kable()

check	value
featureType	profile
profiles	14336
depth levels	6082688
rowSize sums to obs	TRUE
sensor variables	32
levels per profile (min/med/max)	1 / 515 / 3630

Code

nc_close(ncv)

7 Publish

Code

# Two products, published as SEPARATE datasets: they are the same casts at
# different resolutions, and a user choosing between them should see two entries
# with their own sizes and provenance rather than one entry with a hidden variant.
publish_one <- function(nc_path, ds, src_tables, scope) {
  plan <- cc_netcdf_plan(nc_path, ds, RELEASE)
  man  <- cc_netcdf_manifest(plan, ds, RELEASE, RELEASE, file.size(nc_path),
                             source_tables = src_tables, cf_scope = scope)
  cat(glue("{ds}: sha256 {substr(plan$sha256,1,12)}… identical_to={plan$identical_to} upload={plan$upload}\n"))
  res <- cc_netcdf_publish(nc_path, ds, RELEASE, plan, man)
  cat(glue("  -> {res$url}\n"))
  invisible(res)
}
publish_one(OUT_NC,  DATASET,           c("sample", "obs"),  globals$cf_scope)

ctd-cast: sha256 176139ea93fb… identical_to=NA upload=TRUE

-> https://storage.calcofi.io/calcofi-files-public/netcdf/ctd-cast/v2026.07.17/ctd-cast.nc

Code

publish_one(FULL_NC, "ctd-cast_full",   "obs_ctd_full",      globals_full$cf_scope)

ctd-cast_full: sha256 609dc0a54982… identical_to=NA upload=TRUE

-> https://storage.calcofi.io/calcofi-files-public/netcdf/ctd-cast_full/v2026.07.17/ctd-cast_full.nc

Code

dbDisconnect(con, shutdown = TRUE)

1 Overview

2 Setup

3 Shape and integrity

4 Write the CF profile file

5 Verify

6 Full-resolution companion — ctd-cast_full.nc

7 Publish

6 Full-resolution companion — `ctd-cast_full.nc`