Ingest CalCOFI CTD Cast Data

Published

2026-09-24

1 Overview

This notebook ingests CalCOFI CTD data from https://calcofi.org/data/oceanographic-data/ctd-cast-files/, downloads and unzips final and preliminary CTD files, normalizes into tidy tables (ctd_cast, ctd_measurement, ctd_thin, ctd_summary), and exports to Parquet for the CalCOFI integrated database.

1.1 Key Features

Web Scraping: Scrapes all CTD .zip download links from calcofi.org
Smart Filtering:
- Downloads all .zip files for archival completeness
- Only unzips final and preliminary files
- Skips raw/cast/test files
- Archives are kept on Google Drive; their contents are extracted to local scratch under cc_stage_dir() (see Setup) and never written back
Priority-based Selection:
- For each cruise, selects final if available, otherwise preliminary
- Excludes raw/test/prodo cast files
Tidy Normalization:
- ctd_cast: one row per unique cast (cruise/station/direction)
- ctd_measurement: long-format sensor readings at each depth (supplemental)
- ctd_thin: adaptively-thinned ctd_measurement — single direction, canonical types, ~10 m grid with inflections preserved; the headline CTD table
- ctd_summary: summary stats per station/depth/measurement_type across cast directions
- measurement_type: reference table for measurement codes
Standardized Workflow: follows patterns from ingest_calcofi_bottle.qmd — deterministic UUIDs, GCS parquet uploads, metadata sidecars, calcofi4db utilities
Findings: the report at the end — what is known to be wrong with this data, measured from the parquet that was just written and checked by the same QC rule registry apps/ctd-qaqc runs. It renders even when the heavy path is skipped (see Check for Resumable State), so it can be revised without an hour of recomputation.

2 Setup

Code

# chunk timing hook — prints elapsed time for each chunk
knitr::knit_hooks$set(time_it = function(before, options) {
  if (before) {
    .time_it_t0 <<- Sys.time()
  } else {
    elapsed <- round(difftime(Sys.time(), .time_it_t0, units = "secs"), 1)
    tnow <- format(Sys.time(), "%H:%M:%S")
    message(glue::glue("R chunk {options$label}: {elapsed}s ~ {tnow}"))
  }
})
knitr::opts_chunk$set(time_it = TRUE)

# archives this notebook has to fetch from calcofi.org are 25-125 MB and the
# server delivers ~0.4 MB/s; R's default `timeout` of 60 s aborts every one of
# them at 15-25 MB (run of 2026-08-24: seven consecutive truncations, then the
# render failed on the archives that never arrived). One hour per file.
options(timeout = max(3600, getOption("timeout")))

devtools::load_all(here::here("../calcofi4db"))
devtools::load_all(here::here("../calcofi4r"))
librarian::shelf(
  CalCOFI / calcofi4db,
  CalCOFI / calcofi4r,
  DBI,
  dplyr,
  DT,
  fs,
  glue,
  here,
  httr,
  janitor,
  lubridate,
  mapview,
  plotly,
  purrr,
  ps,
  readr,
  rvest,
  sf,
  stringr,
  tibble,
  tidyr,
  zip,
  quiet = T
)

# common ingest settings (overwrite, dir_data)
source(here("libs/ingest.R"))

# say(): status output to STDOUT so it renders in the notebook.
#
# This notebook previously used message() for ~37 status lines. message() writes
# to stderr, which Quarto does not reliably surface under `code-fold: true` — so
# most of the diagnostics below were being computed and then never shown. The
# repo convention is cat() (see .claude/skills/ingest-new/SKILL.md), but cat() adds no
# trailing newline, so this thin wrapper keeps message()'s line semantics.
say <- function(...) cat(..., "\n", sep = "")

# provider/dataset/metadata read from this file's authoritative YAML block
cc           <- read_calcofi_meta(here("ingest_calcofi_ctd-cast.qmd"))
provider     <- cc$provider
dataset      <- cc$dataset
dataset_name <- cc$dataset_meta$dataset_name
tables_owned <- cc$tables_owned
dir_label <- glue("{provider}_{dataset}")
# ARCHIVES live on Google Drive; EXTRACTED FILES do not. dir_dl holds `.zip` and
# nothing else — Drive is the durable, shareable home for ~370 archives and it
# syncs them fine. Unzipping into it is what Drive cannot survive: the 151
# final/preliminary archives expand to ~124,000 files and ~45 GB, which the Drive
# client never finishes syncing, and mid-sync it evicts files to cloud-only
# placeholders (silently read as 0 rows, see [d_csv]) and mints ` 2.csv` conflict
# copies beside them. That cost release v2026.08.08 ten cruises and needed a
# repair script to undo. Extraction now lands under cc_stage_dir(), outside both
# Drive and the repo, where it is disposable: delete it and the next render
# re-extracts from the archives.
dir_dl <- path_expand(glue("{dir_data}/{provider}/{dataset}/download"))
dir_ext <- cc_stage_path("ctd-cast", "unzip", create = TRUE)
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"))
db_checkpoint <- here(glue("data/wrangling/{dir_label}_checkpoint.duckdb"))
dir_tmp <- here(glue("data/tmp/{dir_label}"))
url <- "https://calcofi.org/data/oceanographic-data/ctd-cast-files/"
dir_meta <- here(glue("metadata/{provider}/{dataset}"))

dir_create(
  c(dir_dl, dir_ext, dirname(db_path), dir_parquet, dir_stage, dir_tmp),
  recurse = TRUE)

# Guard the split rather than merely documenting it: an extracted directory left
# in dir_dl by an older run is not inert. [d_csv] globs the extraction root only,
# so those CSVs no longer reach the pivot, but Drive still tries to sync them —
# which is the problem this split exists to solve. Names the offenders and stops.
dirs_in_dl <- list.dirs(dir_dl, recursive = FALSE, full.names = FALSE)
if (length(dirs_in_dl) > 0)
  stop(glue(
    "{length(dirs_in_dl)} extracted director(ies) found in dir_dl:\n",
    paste0("  ", head(dirs_in_dl, 10), collapse = "\n"),
    if (length(dirs_in_dl) > 10) "\n  ..." else "", "\n",
    "dir_dl holds .zip archives only — extraction belongs in {dir_ext}. ",
    "Remove them with `Rscript scripts/prune_ctd_extracts_from_drive.R`."))

# NOTE: nothing destructive happens here, and no connection is opened. Both wait
# for the resume decision below, which needs the scraped zip list to fingerprint
# the inputs — deleting the wrangling DB before knowing whether the heavy path
# will run would throw away the one table the Gantt appendix reads.

# source archival handled by sync_gd_to_gcs.qmd (rclone): ~18 GB of archives is
# still too large to push through the GD FUSE mount with sync_to_gcs()

# load metadata
# read_measurement_type(), NOT read_csv(): the registry stores empty cells that a
# default read turns back into the literal string "NA", which then survives the
# na.omit() in the measurement-column registry below and becomes a phantom
# quality-flag column name. It also validates the file on the way in, so a
# corrupted write fails here rather than reaching the release. See CLAUDE.md.
d_meas_type <- calcofi4db::read_measurement_type(
  here("metadata/measurement_type.csv")
)
d_flds_rd <- read_csv(glue("{dir_meta}/flds_redefine.csv"), show_col_types = F)
d_tbls_rd <- read_csv(glue("{dir_meta}/tbls_redefine.csv"), show_col_types = F)
d_cruise_corrections <- read_csv(
  glue("{dir_meta}/cruise_key_corrections.csv"), show_col_types = F)

# the `data_stage` vocabulary, MOST ADVANCED FIRST — this ordering IS the
# supersession rule (Final 1m-Binned supersedes Preliminary CTD & Bottle 1m-Binned
# supersedes Preliminary CTD 1m-Binned), so `match()` against it gives the priority
# and no second list can drift out of step with it. Declared here because five
# places downstream need it: the classifier, the per-cruise precedence filter, the
# emit guard, the findings table and the coverage Gantt.
CTD_DATA_STAGES <- c(
  "final", "preliminary_with_bottle", "preliminary_without_bottle")
CTD_STAGE_LABELS <- c(
  final                      = "Final 1m-Binned",
  preliminary_with_bottle    = "Preliminary CTD & Bottle 1m-Binned",
  preliminary_without_bottle = "Preliminary CTD 1m-Binned")

2.1 Archives on Drive, contents on local scratch

	path	holds	synced?
`dir_dl`	`~/My Drive/…/calcofi/ctd-cast/download/`	`.zip` only	yes — Drive → `gs://calcofi-files-public/_sync/…`
`dir_ext`	`~/_big/calcofi/ctd-cast/unzip/` (`cc_stage_dir()`, overridable with `$CALCOFI_STAGE_DIR`)	everything unzipped	no — local, disposable

Roughly 370 archives sync to Drive without trouble. Their contents do not: the 151 final/preliminary archives expand to about 124,000 files and 45 GB, and extracting them into the Drive folder — which is what this notebook did until 2026-08-16 — gave the Drive client a sync job it never finished. The failures were not merely slow, they were silent and they reached the release: Drive evicts a file it has synced to a cloud-only placeholder that still reports its full size but reads as an empty CSV, and it mints 2.csv conflict copies whose names no longer end in a cast-direction letter. Release v2026.08.08 lost ten cruises that way and needed a repair script to walk the tree putting it back.

So the archive and its contents now live on different filesystems. Nothing under dir_ext is precious — delete any part of it and the next render re-extracts that part from dir_dl, which is the copy actually worth keeping and the only one worth syncing. Setup stops the render if an extracted directory reappears in dir_dl.

3 Prime Downloads from the Object Store

The authoritative source for the CTD .zip files is the organization Shared Drive (“CalCOFI Data Folder”), mirrored to gs://calcofi-files-public/_sync/calcofi/ctd-cast/download/ by scripts/sync_gdrive_to_gcs.sh (per-cruise CTDCast/CTDPrelim drops) and scripts/sync_jrw_ctd_to_gcs.sh (Jim Wilkinson’s _CTDFinalDB.zip final archive, 1993-08 onward). Copy anything new into dir_dl so the inventory below can see it — the archive, not its contents (see Archives on Drive, contents on local scratch).

This must run before the inventory, not after. It used to sit after both the scrape and the resume decision, which made a Drive-only zip structurally invisible: d_zips was built purely from calcofi.org, so nothing ever unzipped it and the fingerprint never noticed it. 20-2607SH_CTDPrelim.zip sat on GCS and on disk for over a week that way. Set CTD_ZIP_SOURCE="" to disable, or override it to point elsewhere.

Code

# authoritative zip source (gdrive→gcs); empty string disables priming
ctd_zip_source <- Sys.getenv(
  "CTD_ZIP_SOURCE",
  "gcs-calcofi:calcofi-files-public/_sync/calcofi/ctd-cast/download")

prime_zips_from_gcs <- function(src, dest_dir) {
  # skip if disabled or rclone unavailable
  if (src == "" || Sys.which("rclone") == "")
    return(invisible(0L))
  # system2() does NOT quote its args — it pastes them into a shell command — and
  # dest_dir is under "~/My Drive/…", so an unquoted path split on the space and
  # rclone saw three positional arguments instead of two:
  #   Command copy needs 2 arguments maximum: you provided 3 non flag arguments:
  #   ["gcs-calcofi:…/download" "/Users/bbest/My" "Drive/projects/…/download"]
  # It then failed SILENTLY, because the exit status was never checked: the run
  # printed "Priming N zip(s)", primed nothing, and fell through to scraping
  # calcofi.org. Harmless where the zips are already local; a slow surprise on a
  # fresh machine.
  n_src <- tryCatch(
    length(system2(
      "rclone", c("lsf", shQuote(src), "--include", shQuote("*.zip"),
                  "--max-depth", "1"),
      stdout = TRUE, stderr = FALSE)),
    error = function(e) 0L)
  if (n_src == 0) {
    say(glue("No zips at {src} — skipping GCS prime (will use calcofi.org)"))
    return(invisible(0L))
  }
  say(glue("Priming {n_src} zip(s) from {src} → {dest_dir}"))
  status <- system2("rclone", c(
    "copy", shQuote(src), shQuote(dest_dir),
    "--include", shQuote("*.zip"), "--max-depth", "1",
    "--transfers", "8", "--checkers", "16"))
  if (!identical(status, 0L)) {
    say(glue("GCS prime FAILED (rclone exit {status}) — falling back to ",
             "scraping calcofi.org. This is slower but not fatal."))
    return(invisible(0L))
  }
  invisible(n_src)
}

prime_zips_from_gcs(ctd_zip_source, dir_dl)

Priming 246 zip(s) from gcs-calcofi:calcofi-files-public/_sync/calcofi/ctd-cast/download → /Users/bbest/My Drive/projects/calcofi/data-public/calcofi/ctd-cast/download

4 Build the Zip Inventory

The inventory is the union of two sources, not the web scrape alone:

web — the .zip links published on calcofi.org. One HTTP GET, falling back to the cached list when the site is unreachable.
gcs — zips that reached dir_dl from the Shared Drive but are not on calcofi.org. Two kinds today: a preliminary drop that has not been published yet (20-2607SH_CTDPrelim.zip), and Jim Wilkinson’s _CTDFinalDB.zip final archive, which is not on calcofi.org at all.

This runs on every render, ahead of the resume decision below, because the zip set is one of the two things that decide whether there is anything to recompute.

Code

# a zip's cruise and processing kind are both encoded in its filename, identically
# whichever source it came from. `19-9308NH_CTDFinalDB.zip` → cruise_key 9308NH,
# zip_type final; `20-2607SH_CTDPrelim.zip` → 2607SH, preliminary. Declared once so
# the two sources cannot drift into classifying the same name differently.
add_zip_meta <- function(d) {
  d |>
    mutate(
      cruise_key = str_extract(
        file_zip, "\\d{2}-(\\d{4}[A-Z0-9]{2,4})", group = 1),
      zip_type = case_when(
        str_detect(file_zip, "CTDFinal")         ~ "final",
        str_detect(file_zip, "CTDPrelim")        ~ "preliminary",
        str_detect(file_zip, "CTDCast|CTD_Cast") ~ "raw",
        str_detect(file_zip, "CTDTest")          ~ "test",
        TRUE                                     ~ "unknown"))
}

# --- source 1: calcofi.org ----------------------------------------------------
# read web page, extract all .zip download links; cache to CSV
cache_csv <- file.path(dir_meta, "ctd_zip_urls.csv")

d_zips_web <- tryCatch(
  {
    d <- read_html(url) |>
      html_nodes("a[href$='.zip']") |>
      html_attr("href") |>
      tibble(url = _) |>
      mutate(
        url = if_else(
          str_starts(url, "http"),
          url,
          paste0("https://calcofi.org", url)
        ),
        file_zip = basename(url),
        # year comes from the URL directory (/downloads/2026/), which is
        # authoritative for a published zip; month from the cruise code
        year = str_extract(url, "/(\\d{4})/", group = 1) |> as.integer(),
        month = str_extract(url, "-\\d{2}(\\d{2})", group = 1) |> as.integer()
      ) |>
      add_zip_meta()
    write_csv(d, cache_csv)
    say(glue("Scraped {nrow(d)} zip URLs from {url}"))
    d
  },
  error = function(e) {
    if (file_exists(cache_csv)) {
      say(glue("Website unavailable, using cached URLs from {cache_csv}"))
      read_csv(cache_csv, show_col_types = FALSE)
    } else {
      stop(e)
    }
  }
) |>
  mutate(source = "web")

Scraped 285 zip URLs from https://calcofi.org/data/oceanographic-data/ctd-cast-files/

Code

# --- source 2: zips on disk that calcofi.org does not publish ------------------
# `url` is deliberately NA: there is nothing to download, the file is already in
# dir_dl (put there by prime_zips_from_gcs above, or by hand). download_and_unzip()
# reads that NA as "already local, unzip only".
#
# year/month are parsed from the FILENAME here rather than a URL path. The two
# agree wherever both exist (20-2601RL → /downloads/2026/, 2026-01), but a local
# zip has no URL to read them from. Web rows keep the URL derivation so no
# currently-published value changes.
d_zips_local <- tibble(
  file_zip = list.files(dir_dl, pattern = "\\.zip$")) |>
  anti_join(d_zips_web, by = "file_zip") |>
  add_zip_meta() |>
  mutate(
    url    = NA_character_,
    # "19-" / "20-" gives the century, "YYMM" the rest
    year   = as.integer(str_extract(file_zip, "^(\\d{2})", group = 1)) * 100L +
             as.integer(str_extract(file_zip, "^\\d{2}-(\\d{2})", group = 1)),
    month  = as.integer(str_extract(file_zip, "^\\d{2}-\\d{2}(\\d{2})", group = 1)),
    source = "gcs") |>
  select(url, file_zip, year, month, cruise_key, zip_type, source)

d_zips <- bind_rows(d_zips_web, d_zips_local) |>
  arrange(desc(year), desc(month), file_zip)

stopifnot(!any(d_zips$zip_type == "unknown"))

say(glue("Zip inventory: {nrow(d_zips)} ",
         "({sum(d_zips$source == 'web')} calcofi.org, ",
         "{sum(d_zips$source == 'gcs')} Shared Drive only)"))

Zip inventory: 396 (285 calcofi.org, 111 Shared Drive only)

Code

d_zips |>
  mutate(
    file_zip = if_else(
      is.na(url), file_zip, glue("<a href={url}>{file_zip}</a>"))
  ) |>
  select(-url) |>
  dt(
    caption = "All zip files in the inventory",
    fname = "ctd_zip_files",
    escape = F
  )

5 Check for Resumable State

The heavy path here — download, parse, pivot — is about an hour, and this notebook is re-rendered mostly for reasons that have nothing to do with its inputs: a narrative edit, a new diagnostic, a corrected sentence. Paying an hour to prove that nothing changed is what stops a pipeline notebook being usable as a living document, so the decision is made on evidence rather than on a flag: skip the heavy path when the parquet outputs are complete and the inputs they were built from still fingerprint the same.

The inputs are the published zip URL list scraped above, plus the metadata registries that steer the pivot. calcofi4db::input_fingerprint() hashes them and records the result beside the wrangling database; a missing input hashes as <missing> rather than being skipped, so deleting a corrections file invalidates the outputs instead of reading as “no change”.

When the heavy path is skipped, every chunk from here to the GCS upload is displayed but not evaluated — their code is still the record of how the tables were built, but the tables themselves are not rebuilt, so those sections show no output. The Findings section at the end reads the parquet rather than the wrangling database and therefore always runs; that is what makes this worth doing, and what keeps a skipped render honest rather than merely quiet.

Code

# `overwrite` in libs/ingest.R is TRUE for every ingest and cannot be flipped
# globally — the comment there records that trying it broke a different notebook,
# because the skip path is not uniformly safe. So this notebook decides for itself.
# CTD_FORCE_REBUILD=TRUE (or overwrite_all) runs the heavy path regardless.
fingerprint_path <- here(glue("data/wrangling/{dir_label}_inputs.json"))

# kept OUT of dir_parquet on purpose: sync_to_gcs() mirrors that directory to a
# world-readable bucket, and this is local state about one machine's outputs.
fp <- input_fingerprint(
  # ONLY inputs that can change the OUTPUTS belong here. measurement_qual.csv is
  # deliberately absent: it is a vocabulary the diagnostics decode against, never
  # something the pivot reads, so hashing it would force an hour of rebuild for a
  # documentation edit that cannot alter a single row.
  files = c(
    here("metadata/measurement_type.csv"),
    here("metadata/core_dictionary.csv"),
    file.path(dir_meta, "flds_redefine.csv"),
    file.path(dir_meta, "tbls_redefine.csv"),
    file.path(dir_meta, "cruise_key_corrections.csv"),
    file.path(dir_meta, "flag_overrides.csv"),
    file.path(dir_meta, "metadata_derived.csv")),
  # the engine version is an input too: a change to calcofi4db's projection or
  # parsing would otherwise leave the fingerprint identical and keep outputs that
  # the current code would no longer produce
  #
  # keyed on source+filename, NOT on `url`: a Shared-Drive-only zip has no URL, so
  # hashing `d_zips$url` would let the whole Wilkinson archive (and 2607SH) arrive
  # without invalidating anything — the run would report "inputs unchanged" and
  # skip the rebuild that was the entire point.
  #
  # ...and the archive's SIZE on disk: a provider re-publishing a corrected zip under
  # the same name (20-2607SH_CTDPrelim.zip, 2026-09-14: station numbers >= 100 had
  # been cut to three digits) changes nothing else, so without it the fix would read
  # as "inputs unchanged" and never be ingested (CalCOFI/workflows#104)
  values = c(sort(paste(d_zips$source, d_zips$file_zip,
                        file.size(file.path(dir_dl, d_zips$file_zip)))),
             paste0("calcofi4db ", as.character(packageVersion("calcofi4db")))))
fp_prior <- read_input_fingerprint(fingerprint_path)
fp_same  <- !is.null(fp_prior) && identical(fp_prior$hash, fp$hash)

# parquet is "complete" when every non-supplemental manifest table is on disk
parquet_ok    <- FALSE
manifest_path <- file.path(dir_parquet, "manifest.json")
if (file_exists(manifest_path)) {
  mf       <- jsonlite::read_json(manifest_path)
  expected <- setdiff(mf$tables, unlist(mf$supplemental))
  parquet_ok <- all(vapply(expected, function(tbl) {
    file_exists(file.path(dir_stage, paste0(tbl, ".parquet"))) ||
      dir_exists(file.path(dir_stage, tbl))
  }, logical(1)))
}

force_rebuild <- overwrite_all ||
  isTRUE(as.logical(Sys.getenv("CTD_FORCE_REBUILD", "FALSE")))
parquet_complete <- parquet_ok && fp_same && !force_rebuild

if (parquet_complete) {
  say(glue(
    "Inputs unchanged and parquet complete ({length(mf$tables)} tables, ",
    "{format(mf$total_rows, big.mark = ',')} rows) — ",
    "skipping download/parse/pivot"))
  say(glue("  fingerprint {substr(fp$hash, 1, 12)}, recorded {fp_prior$recorded_at}"))
} else if (force_rebuild) {
  say("Rebuilding: forced (CTD_FORCE_REBUILD / overwrite_all)")
} else if (parquet_ok && !fp_same) {
  say("Rebuilding: inputs changed since the recorded fingerprint —")
  for (f in changed_inputs(fp, fp_prior)) say(glue("  {f}"))
} else {
  say("Rebuilding: no complete parquet output found")
}

Rebuilding: inputs changed since the recorded fingerprint —
  /Users/bbest/Github/CalCOFI/workflows-ctd-new-format/metadata/calcofi/ctd-cast/flag_overrides.csv

Code

# --- only now touch anything -------------------------------------------------
# overwrite keeps the checkpoint DB, the downloads, and (importantly) dir_parquet,
# so write_parquet_outputs() can content-hash dedup against the prior run and only
# changed partitions are re-written and re-uploaded.
if (overwrite && !parquet_complete) {
  if (file_exists(db_path)) file_delete(db_path)
  # clear any stale WAL/tmp from an interrupted run so the restored checkpoint
  # DB opens cleanly (avoids "WAL checkpoint iteration does not match" errors)
  db_wal <- paste0(db_path, ".wal")
  db_tmp <- paste0(db_path, ".tmp")
  if (file_exists(db_wal)) file_delete(db_wal)
  if (dir_exists(db_tmp))  dir_delete(db_tmp)
  if (overwrite_all) {
    if (dir_exists(dir_parquet))    dir_delete(dir_parquet)
    if (dir_exists(dir_stage))      dir_delete(dir_stage)
    if (file_exists(db_checkpoint)) file_delete(db_checkpoint)
    rds_files <- list.files(dir_tmp, pattern = "\\.rds$", full.names = TRUE)
    if (length(rds_files) > 0) file_delete(rds_files)
    say("Deleted parquet, checkpoint DB, and RDS intermediates")
  }
}

# restore from checkpoint DB if available (skips expensive read+bind+filter) — but
# only a checkpoint built from THESE archives: it holds ctd_raw, i.e. the files as
# read, so a replaced archive must not be answered from the old copy of it
fp_zips <- input_fingerprint(
  files  = file.path(dir_meta, "flds_redefine.csv"),
  values = sort(paste(d_zips$file_zip, file.size(file.path(dir_dl, d_zips$file_zip)))))
checkpoint_fp_path <- paste0(db_checkpoint, ".inputs.json")
fp_ckpt <- read_input_fingerprint(checkpoint_fp_path)
if (file_exists(db_checkpoint) &&
    (is.null(fp_ckpt) || !identical(fp_ckpt$hash, fp_zips$hash))) {
  file_delete(db_checkpoint)
  say("Checkpoint discarded: built from a different set of archives")
}
if (!parquet_complete && file_exists(db_checkpoint) && !file_exists(db_path)) {
  file_copy(db_checkpoint, db_path, overwrite = TRUE)
  say(glue("Restored from checkpoint: {db_checkpoint}"))
}

Restored from checkpoint: /Users/bbest/Github/CalCOFI/workflows-ctd-new-format/data/wrangling/calcofi_ctd-cast_checkpoint.duckdb

Code

con <- get_duckdb_con(db_path)
load_duckdb_extension(con, "spatial")
load_duckdb_extension(con, "icu")

# DuckDB's memory limit from what is AVAILABLE, not half of total: this notebook
# shares the machine with R's own ~20 GB of scan frames and whatever else runs (the
# 2026-09-23 run crawled at 3-9x under swap with a fixed half-of-RAM limit). 40 % of
# available, between 2 GB and half of total; DuckDB spills to temp_directory beyond it.
# BUILD_MEMORY_GB overrides it (as libs/duckdb_budget.R does in the sibling repos): on
# 2026-09-24 availability was 5.8 GB at start, so this rule gave DuckDB 2 GB while R itself
# held ~25 GB of scan frames that nothing released (fixed in [rename_fields]).
mem_sys <- ps::ps_system_memory()
mem_gb  <- max(2, floor(min(mem_sys$total / 2, mem_sys$avail * 0.4) / 1024^3))
mem_env <- suppressWarnings(as.numeric(Sys.getenv("BUILD_MEMORY_GB", NA)))
if (!is.na(mem_env) && mem_env > 0) mem_gb <- mem_env
say(glue("DuckDB memory_limit {mem_gb} GB (available {round(mem_sys$avail / 1024^3, 1)} of ",
         "{round(mem_sys$total / 1024^3)} GB; BUILD_MEMORY_GB {if (is.na(mem_env)) 'unset' else mem_env})"))

DuckDB memory_limit 6 GB (available 9 of 24 GB; BUILD_MEMORY_GB 6)

Code

# every connection this notebook opens to the wrangling DB gets the same budget. The
# checkpoint step closes and reopens `con`, and until 2026-09-24 the reopened connection
# silently ran on DuckDB's default (80 % of RAM) instead of this limit: the restarted
# render's R process peaked at 22 GB in [ctd_measurement] with a 6 GB limit "set".
# preserve_insertion_order = false lets the large INSERT ... SELECTs of the pivot stream
# instead of buffering to keep row order, which nothing downstream relies on (every
# table is re-sorted at write_parquet_outputs()).
duck_budget <- function(con) {
  dbExecute(con, glue("SET memory_limit = '{mem_gb}GB'"))
  dbExecute(con, glue("SET temp_directory = '{dir_tmp}'"))
  dbExecute(con, "SET preserve_insertion_order = false")
  invisible(con)
}
duck_budget(con)

# a second, cheaper resume: ctd_raw already pivoted from a prior interrupted run
has_ctd_raw <- FALSE
if (!parquet_complete) {
  has_ctd_raw <- "ctd_raw" %in% DBI::dbListTables(con)
  if (has_ctd_raw) {
    n_raw <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM ctd_raw")$n
    say(glue(
      "Checkpoint: ctd_raw already loaded ",
      "({format(n_raw, big.mark = ',')} rows) — ",
      "skipping read+bind+filter"))
  }
}

Checkpoint: ctd_raw already loaded (8,442,099 rows) — skipping read+bind+filter

Code

# set eval for read+bind+filter chunks
skip_read_bind <- parquet_complete || has_ctd_raw
knitr::opts_chunk$set(eval = !skip_read_bind)

6 Choose Which Archives to Read

Jim Wilkinson’s archive and calcofi.org overlap on 66 cruises, and 130 of the 140 comparable files are identical — 20-0302JD_CTDBTL_001-100D.csv from 20-0302JD_CTDFinalDB.zip is md5-identical to calcofi.org’s copy inside 20-0302JD_CTDFinalQC.zip, and so are the great majority of the rest. But 8 files across 4 cruises are not (1810SR, 1501NH, 1604SH, 1607OS), so this is a preference, not an equivalence.

calcofi.org wins where both publish a final. It is the citable published source, it is what the existing 66 cruises in the release were already built from, and switching them would rewrite years of released data on the strength of an undocumented difference. The Shared Drive archive is read only for the cruises calcofi.org has no final for — which is also what keeps the ~1 hr heavy path from doubling to parse duplicates.

The archive is still mirrored completely to GCS (0.89 GB): it is the basis for any future backfill of calcofi.org, and for answering the question below.

Where the two sources disagree, and why it is a provider question

The jrw_overlap_report chunk below compares every overlapping archive and lists the disagreements rather than asserting they cannot happen — an earlier draft asserted byte-identity and was wrong.

The clearest case is 1810SR, where the Shared Drive copy is 2.1× the size: same 74 casts and the same 82 columns, but 59,726 rows against 29,863, and its timestamps carry seconds (14-Oct-2018 18:00:17) where calcofi.org’s are truncated to the minute (10/14/2018 18:00). On that evidence the Shared Drive copy looks like the fuller product, not a stale one.

That matters mostly for obs_ctd_full — ctd_thin reduces to a ~10 m grid either way — but datetime_start_utc is part of the de-duplication key, so it is not purely cosmetic. Which copy is authoritative for 2003–2019 is the provider’s call, not this notebook’s: see question calcofi_ctd-cast_20. Until it is answered the release keeps what it already had.

Code

# a GCS-only final is redundant when calcofi.org already publishes a final for that
# cruise. Everything else — every preliminary, every raw, and every final for a
# cruise calcofi.org has no final for — is read.
cruises_web_final <- d_zips |>
  filter(source == "web", zip_type == "final") |>
  pull(cruise_key) |>
  unique()

d_zips_skip <- d_zips |>
  filter(source == "gcs", zip_type == "final", cruise_key %in% cruises_web_final)

d_zips_read <- d_zips |>
  anti_join(d_zips_skip, by = "file_zip")

d_zips_new_final <- d_zips_read |>
  filter(source == "gcs", zip_type == "final")

say(glue(
  "Archives to read: {nrow(d_zips_read)} of {nrow(d_zips)} ",
  "({nrow(d_zips_skip)} GCS finals skipped as duplicates of calcofi.org; ",
  "{nrow(d_zips_new_final)} GCS finals are new coverage)"))

d_zips_skip |>
  select(file_zip, cruise_key, year, month) |>
  arrange(year, month) |>
  dt(
    caption = glue(
      "Skipped: {nrow(d_zips_skip)} Shared-Drive final archives whose cruise ",
      "calcofi.org already publishes a final for (byte-identical; mirrored to ",
      "GCS but not parsed)"),
    fname = "ctd_zips_skipped"
  )

7 Download and Unzip Files

Code

# `url` is NA for a Shared-Drive-only zip: there is nothing to fetch, the file is
# already in dest_dir (prime_zips_from_gcs put it there). file_zip is passed
# explicitly rather than derived as basename(url), which would be NA for those.
#
# dest_dir (the archive) and ext_dir (its contents) are separate arguments and
# separate filesystems on purpose — see the dir_dl/dir_ext note in [setup].
download_and_unzip <- function(url, file_zip, dest_dir, ext_dir, zip_type,
                               unzip = TRUE) {
  dest_file <- file.path(dest_dir, file_zip)
  dir_unzip <- file.path(ext_dir, str_remove(file_zip, "\\.zip$"))

  if (file_exists(dest_file)) {
    say(glue("Already exists: {file_zip}"))
  } else if (is.na(url)) {
    # primed source vanished between the prime and here — do not silently proceed
    # with a missing archive
    warning(glue("Expected primed zip is missing: {file_zip}"))
    return(invisible(NULL))
  } else {
    say(glue("Downloading: {file_zip}"))
    tryCatch(
      download.file(url, dest_file, mode = "wb"),
      error = function(e) {
        warning(glue("Failed to download {file_zip}: {e$message}"))
        if (file_exists(dest_file)) {
          file_delete(dest_file)
        }
        return(invisible(NULL))
      }
    )
    if (!file_exists(dest_file)) return(invisible(NULL))
  }

  if (!unzip || !zip_type %in% c("final", "preliminary")) {
    if (unzip) say(glue("Skipping unzip (not final/preliminary): {file_zip}"))
    return(invisible(NULL))
  }

  # "the directory exists" was the whole completeness test, which reads an
  # INTERRUPTED extraction as a finished one — and the extraction root is now
  # disposable local scratch, so a half-written tree is a routine outcome (a
  # cancelled render, a full disk) rather than a rare one. Compare the file count
  # on disk against the archive's own member count instead. That count comes from
  # the zip's central directory (utils::unzip, NOT the masked `zip::unzip` — see
  # [jrw_overlap_report]), so it costs a seek, not an extraction, per archive.
  # `<`, not `!=`: an archive may legitimately be over-populated (an older run of
  # a since-shrunk archive), and re-extracting on that is churn, not repair.
  # an archive REPLACED under the same name (see above) keeps its member count, so
  # the count below would bless the old extraction: key it on the archive's size and
  # mtime instead, recorded beside the extraction
  stamp_file <- file.path(dir_unzip, ".archive_stamp")
  stamp_now  <- paste(file.size(dest_file), as.numeric(file.mtime(dest_file)))
  stale <- dir_exists(dir_unzip) && if (file_exists(stamp_file))
    !identical(read_lines(stamp_file), stamp_now) else
    file.mtime(dest_file) > file.mtime(dir_unzip)
  if (stale) {
    say(glue("Archive changed since it was extracted — re-extracting {file_zip}"))
    dir_delete(dir_unzip)
  }

  n_members <- sum(!str_detect(utils::unzip(dest_file, list = TRUE)$Name, "/$"))
  n_on_disk <- if (dir_exists(dir_unzip)) {
    length(list.files(dir_unzip, recursive = TRUE, all.files = TRUE,
                      no.. = TRUE))
  } else 0L

  if (n_on_disk >= n_members) {
    say(glue("Already unzipped: {file_zip} ({n_on_disk} files)"))
  } else {
    say(glue(
      "Unzipping: {file_zip} ({n_members} files",
      if (n_on_disk > 0) glue(", {n_on_disk} already present") else "", ")"))
    dir_create(dir_unzip, recurse = TRUE)
    # overwrite = TRUE so a partial extraction is completed rather than skipped
    # around; unzip()'s default leaves an existing (possibly truncated) file be.
    unzip(dest_file, exdir = dir_unzip, overwrite = TRUE)
  }
  write_lines(stamp_now, stamp_file)
}

d_zips_read |>
  pwalk(function(url, file_zip, zip_type, ...) {
    download_and_unzip(url, file_zip, dir_dl, dir_ext, zip_type, unzip = TRUE)
  })

say(glue("Download and extraction complete: archives in {dir_dl}, ",
         "extracted files in {dir_ext}"))

# A provider may re-publish a corrected archive under the same name — the 2607
# preliminary file was, on 2026-09-14, after station numbers >= 100 had shipped cut
# to three digits (CalCOFI/workflows#104). The Shared Drive stays the authoritative
# copy (it primes dir_dl above), so this REPORTS a calcofi.org archive whose size
# differs from the local copy rather than replacing it: replacing it here would swap
# an input without the per-series diff a data fix needs. Only archives this notebook
# reads, only an HTTP 200 (calcofi.org lists some archives that 404), and a failed
# HEAD is simply not reported. Fix a hit by updating the Shared Drive copy.
head_size <- function(url) tryCatch({
  r <- curl::curl_fetch_memory(url, handle = curl::new_handle(
    nobody = TRUE, timeout = 30, followlocation = TRUE))
  if (r$status_code != 200) return(NA_real_)
  as.numeric(curl::parse_headers_list(r$headers)[["content-length"]] %||% NA)
}, error = function(e) NA_real_)
d_zip_drift <- d_zips_read |>
  filter(!is.na(url), zip_type %in% c("final", "preliminary"),
         file_exists(file.path(dir_dl, file_zip))) |>
  mutate(local_bytes  = file.size(file.path(dir_dl, file_zip)),
         remote_bytes = map_dbl(url, head_size)) |>
  filter(!is.na(remote_bytes), remote_bytes != local_bytes) |>
  select(file_zip, zip_type, local_bytes, remote_bytes, url)
if (nrow(d_zip_drift) > 0) {
  say(glue("{nrow(d_zip_drift)} read archive(s) differ from calcofi.org — update the ",
           "Shared Drive copy and re-render: {paste(d_zip_drift$file_zip, collapse = ', ')}"))
  dt(d_zip_drift, caption = paste(
    "Read archives whose calcofi.org copy differs in size from the local (Shared Drive)",
    "copy: a re-published file the ingest has NOT taken"), fname = "ctd_zip_drift")
} else say("every read archive matches its calcofi.org size (where calcofi.org answers)")

7.1 Overlap guard

Skipping an archive is a claim that nothing is lost by skipping it, so measure the claim on every overlapping cruise rather than trusting one sample. This reads the zip directory (unzip -l, no extraction) and compares each top-level _CTDBTL_*.csv against calcofi.org’s copy of the same filename by uncompressed size — a few seconds for all 66, where md5 would mean unpacking ~2 GB.

Size equality is a screen, not a proof; what it reliably catches is the case that matters, a Shared-Drive copy holding materially more (or less) than the published one. It does not halt the render: a difference here is a question for the provider, not a defect in this pipeline, and the release keeps what it already had either way.

Code

# JRW ships the db CSVs at the archive TOP LEVEL; calcofi.org buries the same
# filenames under db_csv(s)/ inside the _CTDFinalQC archive. Match on basename —
# several cruises split a cruise across separate runs (013-074, 044-060, ...) so
# the file SETS are not always the same shape and a positional compare would lie.
web_btl <- tibble(
  path = list.files(
    dir_ext, pattern = "_CTDBTL_.*\\.csv$", full.names = TRUE,
    recursive = TRUE)) |>
  filter(str_detect(path, "_CTDFinalQC/db[_-]csvs?/")) |>
  mutate(file_csv = basename(path), web_bytes = as.numeric(file_size(path))) |>
  select(file_csv, web_bytes) |>
  distinct(file_csv, .keep_all = TRUE)

# uncompressed sizes straight out of the zip directory — no extraction.
#
# utils::unzip, NOT the bare name: this notebook shelves the `zip` package, whose
# unzip() masks base R's and takes no `list` argument, so `unzip(list = TRUE)`
# fails with "unused argument (list = TRUE)". Everywhere else in this notebook
# unzip() is called only to EXTRACT, where the two signatures happen to agree —
# which is why the masking had never surfaced before.
jrw_btl <- d_zips_skip$file_zip |>
  set_names() |>
  map_dfr(\(fz) {
    z <- utils::unzip(file.path(dir_dl, fz), list = TRUE)
    tibble(file_csv = basename(z$Name), depth = str_count(z$Name, "/"),
           jrw_bytes = z$Length) |>
      filter(str_detect(file_csv, "_CTDBTL_.*\\.csv$"), depth == 0)
  }, .id = "file_zip") |>
  select(-depth)

d_overlap <- jrw_btl |>
  left_join(web_btl, by = "file_csv") |>
  mutate(
    cruise_key = str_extract(file_zip, "\\d{2}-(\\d{4}[A-Z0-9]{2,4})", group = 1),
    ratio  = round(jrw_bytes / web_bytes, 3),
    status = case_when(
      is.na(web_bytes)        ~ "no matching filename on calcofi.org",
      jrw_bytes == web_bytes  ~ "identical size",
      TRUE                    ~ "DIFFERENT size")) |>
  relocate(cruise_key, file_csv, status, ratio, jrw_bytes, web_bytes)

n_same <- sum(d_overlap$status == "identical size")
n_diff <- sum(d_overlap$status == "DIFFERENT size")
n_none <- sum(d_overlap$status == "no matching filename on calcofi.org")

say(glue(
  "Overlap comparison across {n_distinct(d_overlap$cruise_key)} cruises / ",
  "{nrow(d_overlap)} files: {n_same} identical, {n_diff} different, ",
  "{n_none} with no same-named counterpart (split runs)"))

# Not a hard stop — but a WHOLESALE mismatch would mean the wrong folder was
# synced or the layout changed, which is a defect rather than a data question.
stopifnot(
  "no comparable overlap files at all — check the sync and the archive layout" =
    nrow(d_overlap) > 0,
  "the two sources agree on almost nothing — this is a sync/layout problem, not a data question" =
    n_same >= 0.5 * (n_same + n_diff))

d_overlap |>
  filter(status != "identical size") |>
  arrange(desc(ratio)) |>
  dt(
    caption = glue(
      "Shared-Drive vs calcofi.org finals that do NOT match ({n_diff} different ",
      "size, {n_none} split-run). calcofi.org is used for all of these; see ",
      "question calcofi_ctd-cast_20"),
    fname = "ctd_jrw_overlap_diffs"
  )

8 Find and Prioritize CTD Data Files

CalCOFI’s CTD Cast Files page names six categories across its header row, three of which are 1 m-binned products that supersede one another:

Preliminary CTD 1m-Binned → Preliminary CTD & Bottle 1m-Binned → Final 1m-Binned

but each cruise offers a single _CTDPrelim.zip, replaced in place as processing advances, so the two preliminary tiers are invisible from the zip name. The discriminator is the CSV filename inside the archive — _CTDBTL_ once the bottle merge has happened, plain _CTD_ before it:

archive	CSV	`data_stage`
`_CTDFinalQC.zip` (calcofi.org)	`db_csvs/…_CTDBTL_…csv`	`final`
`_CTDFinalDB.zip` (Shared Drive)	`…_CTDBTL_…csv` at top level	`final`
`_CTDPrelim.zip`	`db-csvs/…_CTDBTL_…csv`	`preliminary_with_bottle`
`_CTDPrelim.zip`	`…_CTD_…csv`	`preliminary_without_bottle`

This matters well beyond labelling. Only the station-corrected salinity and oxygen are canonical — following the provider’s own guidance that “station-corrected data are the best” — and those wait on the bottle merge. A preliminary_without_bottle cruise therefore carries no corrected salinity or oxygen at all, which used to be indistinguishable from a preliminary_with_bottle cruise that has them. Splitting the tier is what lets a consumer tell “not measured” from “not merged yet”.

Code

d_csv <- tibble(
  path = list.files(
    dir_ext,
    pattern = "\\.csv$",
    recursive = TRUE,
    full.names = TRUE
  )
) |>
  mutate(
    file_csv = basename(path),
    path_unzip = str_replace(path, glue("{dir_ext}/"), ""),
    dir_unzip = str_extract(path_unzip, "^[^/]+"),
    cruise_key = str_extract(
      path_unzip,
      "\\d{2}-(\\d{4}[A-Z0-9]{2,4})_.*",
      group = 1
    ),
    data_stage = case_when(
      # final, calcofi.org: db_csv(s)/ inside a *Final* archive — EXCEPT the
      # superseded exports some archives keep beside the canonical file:
      # db-csvs/orig/, orig-csvs/, orig_dbcsvs/, uncorrected/. Every cast in those
      # is also in the top-level file (verified over all 16,576 CSVs on disk,
      # 2026-08-24) — except one: 20-0010NH/db-csvs/orig/ carries cast 001, JRW's
      # "First Cast - switch fluorometer" test dip with a 1992 timestamp and scans
      # to 17,178 m, which the canonical file rightly omits. The old unanchored
      # match classified it `final`; casts 002-074 lost to the shorter canonical
      # path in dedup, cast 001 had no competitor, and v2026.08.14 released it at
      # 14,671 m. NOT "top level only": 20-1104SH's db_csvs/separate_runs/ holds
      # casts 031-036 that exist nowhere else, so subfolders that are not
      # superseded exports stay in and dedup picks the canonical copy.
      str_detect(path_unzip, "Final.*db[_|-]csvs?/") &
        !str_detect(path_unzip, regex("db[_|-]csvs?/(orig|uncorrected)[^/]*/",
                                      ignore_case = TRUE)) ~ "final",
      # final, Shared Drive: JRW's *_CTDFinalDB archives put the db CSVs at the
      # TOP LEVEL. The top-level restriction is load-bearing — these archives also
      # carry orig/, orig_dbcsvs/, SeparateRuns_Fl/ and south-north/ copies of the
      # same casts. Without it both copies classify `final`, and which one survives
      # falls to the `ORDER BY length("_source_file")` tiebreak in dedup_ctd_raw:
      # correct by accident, and only until a path length changes.
      str_detect(dir_unzip, "CTDFinalDB$") &
        path_unzip == paste0(dir_unzip, "/", file_csv) ~ "final",
      # preliminary, bottle-merged: _CTDBTL_ means the bottle merge has run
      str_detect(dir_unzip, "Prelim") &
        str_detect(file_csv, "_CTDBTL_") ~ "preliminary_with_bottle",
      # preliminary, sensor only: _CTD_ (no BTL) — corrected salinity/oxygen and
      # every bottle column are empty in these
      str_detect(dir_unzip, "Prelim") &
        str_detect(file_csv, "_CTD_") ~ "preliminary_without_bottle",
      # NOTE: the hardcoded `cruise_key == "2111SR"` arm that used to sit here is
      # gone, and deliberately. 2111SR puts its db CSVs in csvs-plots/ rather than
      # db-csvs/, so the old *directory*-based rules missed it and it needed a
      # cruise-specific path regex. The rules above key on the FILENAME token and
      # do not care how deep the file sits, so
      # 20-2111SR_CTDPrelim/20-2111SRCTDPrelim/csvs-plots/20-2111SR_CTDBTL_001-075D.csv
      # now classifies as preliminary_with_bottle on its own — which is also
      # accurate than the old arm, which could only say "preliminary". Verified
      # against the archive on disk; a regression would trip the content
      # assertion in [emit_core].
      .default = NA_character_
    ),
    # The direction letter is the LAST character of the name — and "last" has to
    # survive a copy marker. Extraction used to land inside Google Drive, which
    # named a conflict copy `19-9501JD_CTDBTL_001-646D 2.csv`; that ends in
    # " 2.csv", matched neither arm of the old `…D\\.csv$` test, and came through
    # as NA. Not cosmetic: [ctd_measurement] joins on cast_dir and NULL = NULL is
    # never true, and [ctd_thin] needs a downcast. Extraction has since moved to
    # `dir_ext`, off Drive, so nothing mints these any more — the stripping stays
    # because it is two regexes against a filename and the archives themselves
    # are not guaranteed free of a name Drive once produced.
    file_stem = file_csv |>
      str_remove("\\.csv$") |>
      str_remove("\\s*\\(\\d+\\)$") |>          # "name (2).csv"
      str_remove("\\s+\\d+$"),                  # "name 2.csv"
    cast_dir = case_when(
      str_detect(file_stem, regex("U$", ignore_case = T)) ~ "U",
      str_detect(file_stem, regex("D$", ignore_case = T)) ~ "D"
    ),
    # supersession order, straight off CTD_DATA_STAGES so a new stage cannot be
    # added to the vocabulary without also being ranked (an unranked stage would
    # otherwise sort last and be silently discarded by the filter below)
    priority = match(data_stage, CTD_DATA_STAGES)
  ) |>
  relocate(cruise_key, path_unzip) |>
  arrange(cruise_key, path_unzip) |>
  filter(data_stage %in% CTD_DATA_STAGES) |>
  # An archive skipped in [pick_archives] is never unzipped — but this chunk globs
  # the extraction directory, not the inventory, so a copy left behind by an
  # earlier run (before the skip existed, or from a run where it was disabled)
  # would still be found here and would compete with calcofi.org's identical copy.
  # Drop them by their unzip directory so the skip decision holds regardless of
  # what is on disk.
  filter(!dir_unzip %in% str_remove(d_zips_skip$file_zip, "\\.zip$"))

# for each cruise keep only the most-advanced stage available:
# final > preliminary_with_bottle > preliminary_without_bottle
d_priority <- d_csv |>
  group_by(cruise_key) |>
  summarize(
    best_priority = min(priority),
    .groups = "drop"
  )

d_csv <- d_csv |>
  inner_join(d_priority, by = "cruise_key") |>
  filter(priority == best_priority) |>
  select(-best_priority)

# A file whose name does not END in a direction letter is not a db product this
# ingest reads. Four exist — `…084IDnoQC.csv` and `…067Dno001b.csv` pairs — and
# they sit beside the canonical D/U pair inside the same archive. Carrying them
# with `cast_dir = NA` was strictly worse than dropping them: [ctd_measurement]
# joins on cast_dir, `NULL = NULL` is never true, so they became 112,245 cast
# rows that could hold no measurement and still reached `sample` as casts with
# nothing under them. Drop them, and say which.
d_csv_nodir <- d_csv |> filter(is.na(cast_dir))
d_csv       <- d_csv |> filter(!is.na(cast_dir))

if (nrow(d_csv_nodir) > 0)
  d_csv_nodir |>
    select(cruise_key, file_csv, data_stage) |>
    arrange(cruise_key, file_csv) |>
    dt(
      caption = glue(
        "Skipped: {nrow(d_csv_nodir)} file(s) whose name resolves no cast ",
        "direction (alternate renderings beside the canonical D/U pair)"),
      fname = "ctd_files_no_cast_direction")

# Every cruise must keep a downcast. [ctd_thin] takes one direction per physical
# cast and prefers D, so a cruise with no D produces no thin rows and therefore
# no `obs` — and because write_parquet_outputs() DELETES a partition that has
# left the data, and sync_to_gcs() mirrors the deletion, the cruise then leaves
# the release entirely while its casts survive in `sample`. That is precisely how
# v2026.08.08 lost 10 cruises and 874,000 observations with no error anywhere.
cruises_no_downcast <- d_csv |>
  group_by(cruise_key) |>
  summarize(has_downcast = any(cast_dir == "D"), .groups = "drop") |>
  filter(!has_downcast) |>
  pull(cruise_key)

if (length(cruises_no_downcast) > 0)
  stop(glue(
    "{length(cruises_no_downcast)} cruise(s) have no downcast file: ",
    "{paste(cruises_no_downcast, collapse = ', ')}.\n",
    "ctd_thin cannot emit a headline profile for them, so they would reach the ",
    "release as casts with no observations."))

cruises_csv_notzip <- setdiff(
  unique(d_csv$cruise_key),
  unique(d_zips$cruise_key)
) |>
  sort()
stopifnot(length(cruises_csv_notzip) == 0)

cruises_zip_notcsv <- setdiff(
  unique(d_zips$cruise_key),
  unique(d_csv$cruise_key)
) |>
  sort()

d_csv |>
  select(-any_of(c("path", "data", "col_empties", "col_types"))) |>
  relocate(cruise_key, path_unzip) |>
  arrange(cruise_key, path_unzip) |>
  dt(
    caption = "Files to Ingest",
    fname = "ctd_files_to_ingest"
  )

9 Read and Standardize CTD Files

Remove repeat header rows that can occur within files.

Code

d_csv <- d_csv |>
  arrange(basename(path)) |>
  mutate(
    data = map2(path, seq_along(path), \(path, idx) {
      say(glue("Reading {idx}/{nrow(d_csv)}: {basename(path)}"))

      all_lines <- read_lines(path)
      header_line <- all_lines[1]

      # find repeat header rows (excluding row 1)
      repeat_header_rows <- which(all_lines[-1] == header_line)

      if (length(repeat_header_rows) > 0) {
        all_lines <- all_lines[-(repeat_header_rows + 1)]
        tmp_file <- tempfile(fileext = ".csv")
        write_lines(all_lines, tmp_file)
        data <- read_csv(tmp_file, guess_max = Inf, show_col_types = F) |>
          clean_names()
        file_delete(tmp_file)
      } else {
        data <- read_csv(path, guess_max = Inf, show_col_types = F) |>
          clean_names()
      }
      # a flag column the provider left empty throughout a file reads as logical
      # (all NA). Make every flag column numeric so the files bind to one type and
      # COALESCE(provider flag, sensor flag) below never mixes BOOLEAN with DOUBLE.
      #
      # ...and DROP a 2026-layout flag column the file leaves wholly empty. bind_rows
      # broadcasts every column to every scan of the archive, so the 27 new *Q columns
      # of ONE file (the corrected 2607 preliminary, every code empty) became ~30M
      # doubles each and took R past its 24 GB vector limit in
      # [pts_distance_from_lnst] (run of 2026-09-23). An all-empty column carries no
      # flag, and [meas_col_registry] adds any flag column still missing as NULL
      # before the COALESCE, so dropping it loses nothing.
      q_new <- c("temp_ave_q", "salt1_corr_q", "salt2_corr_q", "salt_ave_corr_q",
                 "ox1_cruise_corr_q", "ox1_sta_corr_q", "ox2_cruise_corr_q",
                 "ox2_sta_corr_q", "ox_ave_sta_corr_q", "ox1u_mq", "ox1u_m_cruise_corr_q",
                 "ox1u_m_sta_corr_q", "ox2u_mq", "ox2u_m_cruise_corr_q", "ox2u_m_sta_corr_q",
                 "ox_aveu_m_sta_corr_q", "est_chl_cruise_corr_q", "est_chl_sta_corr_q",
                 "est_no3_cruise_corr_q", "est_no3_sta_corr_q", "batq", "po_t1q", "po_t2q",
                 "dyn_ht_q", "svaq", "ox_sat1q", "ox_sat2q")
      q_empty <- intersect(q_new, names(data))[
        vapply(intersect(q_new, names(data)), \(c) all(is.na(data[[c]])), logical(1))]
      data |>
        mutate(across(where(is.logical) & matches("q$"), as.numeric)) |>
        # the layout marker survives the drop (the cross-check needs it)
        mutate(.src_qc_2026 = "salt_ave_corr_q" %in% names(data)) |>
        select(-any_of(q_empty))
    }),
    nrows = map_int(data, ~ if (is.null(.x)) 0 else nrow(.x)),
    # the file layout: the CTD team's 2026 processing software writes a quality
    # code for every corrected, averaged and derived series (SaltAve_CorrQ,
    # OxAve_StaCorrQ, EstChl_StaCorrQ, …; CalCOFI/workflows#105), the legacy
    # layout flags only the sensors. The cross-check in [combine_sensor_pairs]
    # reads it.
    src_format = map_chr(data, ~ if (isTRUE(.x$.src_qc_2026[1])) "qc_2026" else "legacy"),
    data = map(data, ~ select(.x, -.src_qc_2026))
  )

say(glue("file layouts: {paste(names(table(d_csv$src_format)), table(d_csv$src_format), sep = ' = ', collapse = ', ')}"))
write_rds(d_csv |> transmute(`_source_file` = path_unzip, src_format),
          glue("{dir_tmp}/d_src_format.rds"))

d_csv |>
  arrange(cruise_key, file_csv) |>
  select(cruise_key, file_csv, nrows, src_format) |>
  dt(
    caption = "Number of rows read per file",
    fname = "ctd_rows_per_file"
  ) |>
  formatCurrency("nrows", currency = "", digits = 0, mark = ",")

# A file that reads zero rows is a failure, not an observation — and it is the
# one failure this pipeline could not see.
#
# Extraction used to land inside Google Drive, which evicts a synced file to a
# cloud-only placeholder. `list.files()` still reported the full size (12,254,953
# bytes for 19-9501JD_CTDBTL_001-646D.csv) and `ls -lO` marked it `dataless`, but
# reading it timed out at the filesystem and `read_csv()` returned a **0-row
# tibble with no error**. On the v2026.08.08 run 20 files across 14 cruises read 0
# rows exactly this way. Nothing downstream distinguishes "this cruise has no
# downcast" from "the downcast file was empty".
#
# `dir_ext` is now plain local disk, so that specific cause is gone — but the
# assertion stays, because what it really tests is that the bytes reaching the
# pivot are the bytes in the archive, and a truncated extraction says the same
# thing. The remedy just changed: delete the archive's extraction directory and
# re-render, rather than coaxing a file back out of Drive.
d_csv_empty <- d_csv |> filter(nrows == 0)

if (nrow(d_csv_empty) > 0)
  stop(glue(
    "{nrow(d_csv_empty)} CSV(s) read 0 rows:\n",
    paste0("  ", d_csv_empty$path_unzip, collapse = "\n"), "\n",
    "A CTD db CSV is never legitimately empty — its extraction is incomplete. ",
    "Delete the offending director(ies) under {dir_ext} and re-render; ",
    "[download_and_unzip] re-extracts from the archive in {dir_dl}."))

10 Detect and Correct Column Type Mismatches

Code

d_csv <- d_csv |>
  mutate(
    col_empties = map(data, \(x) {
      tibble(
        col_name = names(x),
        n_empty = map_int(x, \(col) sum(is.na(col)))
      ) |>
        filter(n_empty == nrow(x)) |>
        pull(col_name)
    }),
    col_types = map2(data, col_empties, \(x, y) {
      tibble(
        col_name = names(x),
        col_type = map_chr(x, \(col) class(col)[1])
      ) |>
        filter(!col_name %in% y)
    })
  )

# find most common type for each column across all files
d_types <- d_csv |>
  select(path, col_types) |>
  unnest(col_types) |>
  count(col_name, col_type) |>
  group_by(col_name) |>
  slice_max(n, n = 1, with_ties = FALSE) |>
  ungroup() |>
  select(col_name, expected_type = col_type)

d_mismatches <- d_csv |>
  select(cruise_key, path, col_types) |>
  unnest(col_types) |>
  left_join(d_types, by = "col_name") |>
  filter(col_type != expected_type) |>
  arrange(col_name, path)

if (nrow(d_mismatches) > 0) {
  say("Type mismatches detected - converting columns...")
}

# bind data, converting mismatched columns to expected type
d_bind <- d_csv |>
  mutate(
    data = map2(data, path, \(x, p) {
      x_mismatches <- d_mismatches |>
        filter(path == p)

      if (nrow(x_mismatches) > 0) {
        for (i in 1:nrow(x_mismatches)) {
          col <- x_mismatches$col_name[i]
          expected <- x_mismatches$expected_type[i]
          na_before <- sum(is.na(x[[col]]))

          suppressWarnings({
            x[[col]] <- switch(
              expected,
              "numeric" = as.numeric(x[[col]]),
              "integer" = as.integer(x[[col]]),
              "logical" = as.logical(x[[col]]),
              "character" = as.character(x[[col]]),
              x[[col]]
            )
          })

          na_after <- sum(is.na(x[[col]]))
          na_generated <- na_after - na_before

          if (na_generated > 0) {
            say(glue(
              "  {basename(p)}: {col} ({x_mismatches$col_type[i]} -> {expected}) generated {na_generated} NAs"
            ))
          }

          d_mismatches[
            d_mismatches$path == p & d_mismatches$col_name == col,
            "nas_generated"
          ] <<- na_generated
        }
      }
      x
    })
  ) |>
  unnest(data)

if (nrow(d_mismatches) > 0) {
  d_mismatches |>
    group_by(col_name, expected_type, col_type) |>
    summarize(
      n_files = n(),
      total_nas = sum(nas_generated, na.rm = TRUE),
      files = paste(basename(path), collapse = "; "),
      .groups = "drop"
    ) |>
    arrange(desc(total_nas)) |>
    dt(
      caption = "Type mismatches by column",
      fname = "ctd_type_mismatches_by_column"
    ) |>
    formatCurrency(
      c("n_files", "total_nas"),
      currency = "",
      digits = 0,
      mark = ","
    )
}

# reconcile cruise_key from study where mismatched
d_bind <- d_bind |>
  mutate(
    cruise_key = if_else(
      !is.na(study) &
        cruise_key != study &
        study != "Study",
      study,
      cruise_key
    ),
    `_source_file` = path_unzip
  ) |>
  relocate(`_source_file`, .after = cruise_key) |>
  select(
    -path_unzip,
    -path,
    -file_csv,
    -dir_unzip,
    -priority,
    -project,
    -study,
    -nrows,
    -col_empties,
    -col_types,
    -src_format
  )

# fill missing ord_occ from cast_id
d_bind <- d_bind |>
  mutate(
    ord_occ = if_else(
      is.na(ord_occ) & !is.na(cast_id),
      str_extract(cast_id, "_(\\d{3})", group = 1),
      ord_occ
    )
  )

# normalize raw coordinate column names (CalCOFI CTD CSV ships Lon_Dec/Lat_Dec ->
# clean_names lon_dec/lat_dec) to the canonical dictionary names used downstream
d_bind <- d_bind |>
  rename(longitude = lon_dec, latitude = lat_dec)

# save intermediate for resumability
write_rds(d_bind, glue("{dir_tmp}/d_bind_pre_datetime.rds"), compress = "gz")

11 Format Date-Time Column

Code

d_bind <- read_rds(glue("{dir_tmp}/d_bind_pre_datetime.rds"))

d_bind <- d_bind |>
  mutate(
    date_time_utc = trimws(date_time_utc),
    date_time_format = case_when(
      str_detect(
        date_time_utc,
        "^\\d{1,2}-[A-z]{3}-\\d{4} \\d{1,2}:\\d{1,2}:\\d{1,2}$"
      ) ~ "dmy_hms",
      str_detect(
        date_time_utc,
        "^\\d{1,2}/\\d{1,2}/\\d{4} \\d{1,2}:\\d{1,2}$"
      ) ~ "mdy_hm",
      TRUE ~ "unknown"
    ),
    datetime_start_utc = NA,
    datetime_start_utc = if_else(
      date_time_format == "dmy_hms",
      suppressWarnings(dmy_hms(date_time_utc, tz = "UTC")),
      datetime_start_utc
    ),
    datetime_start_utc = if_else(
      date_time_format == "mdy_hm",
      suppressWarnings(mdy_hm(date_time_utc, tz = "UTC")),
      datetime_start_utc
    )
  ) |>
  relocate(date_time_format, datetime_start_utc, .after = date_time_utc) |>
  arrange(datetime_start_utc, depth)

stopifnot(all(!is.na(d_bind$datetime_start_utc)))

write_rds(d_bind, glue("{dir_tmp}/d_bind_post_datetime.rds"), compress = "gz")

d_bind |>
  group_by(cruise_key, `_source_file`, date_time_format) |>
  summarize(
    n_records = n(),
    datetime_min = min(datetime_start_utc),
    datetime_max = max(datetime_start_utc),
    .groups = "drop"
  ) |>
  arrange(cruise_key, `_source_file`, date_time_format) |>
  pivot_wider(
    names_from = date_time_format,
    values_from = n_records,
    values_fill = NA_real_
  ) |>
  relocate(dmy_hms, mdy_hm, .after = cruise_key) |>
  dt(
    caption = "Date-Time formats detected by cruise_key and source_file",
    fname = "ctd_datetime_formats"
  ) |>
  formatDate(c("datetime_min", "datetime_max"), method = "toLocaleString")

12 Pseudo-NA Values and Fill Missing Coordinates

Code

d_bind <- read_rds(glue("{dir_tmp}/d_bind_post_datetime.rds"))

pseudoNA_values <- c(-9.99e-29, -99)

# --- station numbers cut to three digits: repaired by the cast's own position ---
# The CTD team's preliminary files for 2507SR, 2511SR, 2601RL, 2604SH and 2607SH
# wrote station numbers >= 100 with three digits: 100/110/120 read as 000/010/020
# (`Sta` "000.0", `Sta_ID` "090.0 000.0"). The distance filter then dropped every
# such cast as a position error, so those cruises lost their offshore stations in
# every release through v2026.09.11 (CalCOFI/workflows#104). The provider re-issued
# 2607 on 2026-09-14, and the ONLY difference is `Sta`/`Sta_ID` = station + 100; the
# other four are not re-issued yet (questions.csv Q39). So repair exactly that, and
# only where the cast's recorded GPS position says so: an inshore number (< 30)
# whose cast lies within CTD_TRUNC_KM of station + 100 and farther than that from
# the station as written. Every repair is listed below; [filter_pts] still stops the
# render if a truncated station survives unrepaired.
CTD_TRUNC_KM <- 10   # the distance filter's own cutoff (max_dist_dec_lnst_km)
sta_num  <- suppressWarnings(as.numeric(d_bind$sta))
line_num <- suppressWarnings(as.numeric(d_bind$line))
i_cand <- which(!is.na(sta_num) & !is.na(line_num) & sta_num < 30 & line_num > 0 &
                  !is.na(d_bind$longitude) & !is.na(d_bind$latitude) &
                  d_bind$longitude > -180 & d_bind$latitude > -90)
d_trunc_fix <- tibble(
  .i = i_cand, src = d_bind$`_source_file`[i_cand], cast_id = d_bind$cast_id[i_cand],
  line = line_num[i_cand], sta = sta_num[i_cand],
  lon = d_bind$longitude[i_cand], lat = d_bind$latitude[i_cand]) |>
  group_by(src, cast_id, line, sta) |>
  summarize(lon = median(lon), lat = median(lat), n_scans = n(), .groups = "drop")
if (nrow(d_trunc_fix) > 0) {
  nominal <- function(line, sta) sf_project(
    from = "+proj=calcofi", to = "+proj=longlat +datum=WGS84", pts = cbind(line, sta))
  hav_km <- function(lon1, lat1, lon2, lat2) {
    r <- pi / 180
    a <- sin((lat2 - lat1) * r / 2)^2 +
      cos(lat1 * r) * cos(lat2 * r) * sin((lon2 - lon1) * r / 2)^2
    2 * 6371.0088 * asin(sqrt(a))
  }
  xy0   <- nominal(d_trunc_fix$line, d_trunc_fix$sta)
  xy100 <- nominal(d_trunc_fix$line, d_trunc_fix$sta + 100)
  d_trunc_fix <- d_trunc_fix |>
    mutate(km_as_written = hav_km(lon, lat, xy0[, 1], xy0[, 2]),
           km_plus_100   = hav_km(lon, lat, xy100[, 1], xy100[, 2])) |>
    filter(km_plus_100 <= CTD_TRUNC_KM, km_as_written > CTD_TRUNC_KM)
}
if (nrow(d_trunc_fix) > 0) {
  key_fix <- paste(d_trunc_fix$src, d_trunc_fix$cast_id, d_trunc_fix$line, d_trunc_fix$sta)
  i_fix   <- i_cand[paste(d_bind$`_source_file`[i_cand], d_bind$cast_id[i_cand],
                          line_num[i_cand], sta_num[i_cand]) %in% key_fix]
  sta_new <- sta_num[i_fix] + 100
  d_bind$sta[i_fix] <- if (is.character(d_bind$sta)) sprintf("%05.1f", sta_new) else sta_new
  d_bind$sta_id[i_fix] <- sprintf("%05.1f %05.1f", line_num[i_fix], sta_new)
  say(glue("station numbers repaired (+100, confirmed by position): ",
           "{nrow(d_trunc_fix)} casts, {format(length(i_fix), big.mark = ',')} scans, ",
           "files: {paste(unique(basename(d_trunc_fix$src)), collapse = ', ')}"))
  write_rds(d_trunc_fix, glue("{dir_tmp}/d_station_truncation_repaired.rds"))
  d_trunc_fix |>
    mutate(file = basename(src), sta_repaired = sta + 100) |>
    select(file, cast_id, line, sta_as_written = sta, sta_repaired, n_scans,
           km_as_written, km_plus_100) |>
    dt(caption = paste(
      "Station numbers cut to three digits by the source, repaired to station + 100",
      "where the cast's recorded position confirms it (CalCOFI/workflows#104)"),
      fname = "ctd_station_truncation_repaired")
} else say("no truncated station numbers found")
rm(sta_num, line_num, i_cand)

d_bind <- d_bind |>
  mutate(
    is_lnsta_pseudoNA = ifelse(
      as.numeric(line) == 0 & as.numeric(sta) == 0,
      TRUE,
      FALSE
    ),
    is_lon_pseudoNA = map_lgl(longitude, ~ some(pseudoNA_values, near, .x)),
    is_lat_pseudoNA = map_lgl(latitude, ~ some(pseudoNA_values, near, .x)),
    longitude = if_else(
      is_lon_pseudoNA | is_lat_pseudoNA,
      NA_real_,
      longitude
    ),
    latitude = if_else(
      is_lon_pseudoNA | is_lat_pseudoNA,
      NA_real_,
      latitude
    ),
    lon_lnst = if_else(
      !is.na(line) & !is.na(sta) & !is_lnsta_pseudoNA,
      as.numeric(sf_project(
        from = "+proj=calcofi",
        to = "+proj=longlat +datum=WGS84",
        pts = cbind(x = as.numeric(line), y = as.numeric(sta))
      )[, 1]),
      NA_real_
    ),
    lat_lnst = if_else(
      !is.na(line) & !is.na(sta) & !is_lnsta_pseudoNA,
      as.numeric(sf_project(
        from = "+proj=calcofi",
        to = "+proj=longlat +datum=WGS84",
        pts = cbind(x = as.numeric(line), y = as.numeric(sta))
      )[, 2]),
      NA_real_
    ),
    longitude = if_else(is.na(longitude), lon_lnst, longitude),
    latitude = if_else(is.na(latitude), lat_lnst, latitude)
  )

# A scan with no position AND no CalCOFI line/station to derive one from cannot
# be placed, so it cannot enter `sample`/`obs` (both key on a position). Until
# 2026-08-24 this was a bare hard stop, which held for years because every such
# scan had a line/station. Then calcofi.org's 19-9604JD and 19-9608NH FinalQC
# archives (fetched for the first time that day; JRW's FinalDB copies omit them)
# brought 9,225 scans from the "RATHBURN CORE STN" canyon site: lat/lon -99 and
# the station NAME in the line/sta columns. Those are not CalCOFI grid casts
# (questions.csv Q, "Rathburn"). Report them by file and drop them; the assertion
# below still holds for everything placeable, and the drop is bounded so an
# archive that loses its positions wholesale fails loudly rather than vanishing.
CTD_UNPLACED_MAX_PCT <- 0.2   # 9,225 of 9,135,263 = 0.10% at 2026-08-24; only ever DOWN

d_unplaced <- d_bind |>
  filter(is.na(longitude) | is.na(latitude))
if (nrow(d_unplaced) > 0) {
  d_unplaced |>
    count(cruise_key, `_source_file`, line, sta, name = "n_scans") |>
    arrange(desc(n_scans)) |>
    dt(
      caption = glue(
        "{format(nrow(d_unplaced), big.mark = ',')} scans DROPPED: no position and ",
        "no line/station to derive one from"),
      fname = "ctd_scans_unplaced")
}
stopifnot(
  "unplaceable scans exceed the ratchet: a source lost its positions; investigate before dropping" =
    100 * nrow(d_unplaced) / nrow(d_bind) <= CTD_UNPLACED_MAX_PCT)
d_bind <- d_bind |>
  filter(!is.na(longitude), !is.na(latitude))

stopifnot(sum(is.na(d_bind$longitude)) == 0)
stopifnot(sum(is.na(d_bind$latitude)) == 0)

d_bind |>
  filter(is_lon_pseudoNA | is_lat_pseudoNA) |>
  group_by(cruise_key) |>
  summarize(
    n_lon_pseudoNA = sum(is_lon_pseudoNA),
    n_lat_pseudoNA = sum(is_lat_pseudoNA)
  ) |>
  dt(
    caption = glue(
      "Cruises with pseudo-NAs ({paste(pseudoNA_values, collapse = ',')}) ",
      "set to NA, then filled from line/station coordinates."
    ),
    fname = "ctd_cruises_pseudoNA_lonlat"
  ) |>
  formatCurrency(
    c("n_lon_pseudoNA", "n_lat_pseudoNA"),
    currency = "",
    digits = 0,
    mark = ","
  )

13 Distance Filtering

Code

max_dist_dec_lnst_km <- 10

# THIS USED TO MATERIALIZE EVERY SCAN AS AN sf OBJECT WITH TWO GEOMETRY COLUMNS.
# With 45 more cruises that stopped fitting: `pts` reached the machine's 24 GB
# ceiling and `filter_pts` died with "vector memory limit of 24.0 Gb reached"
# after the distance step alone had run for 74 minutes. Both symptoms came from
# doing point work at the wrong grain, so fix the grain rather than the limit —
# raising mem.maxVSize on a 24 GB machine only buys swapping, and DuckDB is
# already holding half of RAM.
#
# Three changes, none of which alters a single output value:
#
#  1. The nominal line/station coordinate takes ~200 distinct values in the whole
#     dataset — one per CalCOFI grid station — but `purrr::map2(st_point)` built
#     one sfg per SCAN, i.e. ~30M R-level closure calls to produce ~200 distinct
#     points. Build the lookup once and join it on.
#  2. Run per cruise. Every operation here is row-wise (distance by_element) or a
#     point-in-polygon against a static grid, so there is no cross-cruise
#     dependency and batching is exactly equivalent — it just bounds peak memory
#     to one cruise instead of the whole archive.
#  3. Return PLAIN TIBBLES. Nothing downstream wants geometry: `d_filt` already
#     did `st_drop_geometry() |> select(-geom_lnst)`, and the two map previews
#     below rebuild it from lon/lat on a ~100-point sample.

grid_ref <- calcofi4r::cc_grid |>
  rename(any_of(c(site_key = "sta_key"))) |>
  select(grid_site = site_key)

# ~200 rows: the distinct nominal station coordinates, as points, built once
lnst_ref <- d_bind |>
  distinct(lon_lnst, lat_lnst) |>
  filter(!is.na(lon_lnst), !is.na(lat_lnst)) |>
  mutate(
    geom_lnst = purrr::map2(lon_lnst, lat_lnst, \(x, y) st_point(c(x, y))) |>
      st_sfc(crs = 4326))

say(glue("Distance filter: {nrow(lnst_ref)} distinct station coordinates ",
         "across {format(nrow(d_bind), big.mark = ',')} scans"))

# compute the three derived columns for one cruise, and only those columns — the
# batch results are bound and joined back, so d_bind is never duplicated
dist_one_cruise <- function(d) {
  s <- d |>
    select(.row_id, longitude, latitude, lon_lnst, lat_lnst) |>
    left_join(lnst_ref, by = c("lon_lnst", "lat_lnst")) |>
    st_as_sf(coords = c("longitude", "latitude"), remove = FALSE, crs = 4326)

  # a scan with no nominal coordinate gets NA distance, exactly as st_distance
  # against an empty geometry did before
  s$dist_dec_lnst_km <- rep(NA_real_, nrow(s))
  has <- !st_is_empty(s$geom_lnst) & !is.na(s$lon_lnst)
  if (any(has))
    s$dist_dec_lnst_km[has] <- st_distance(
      st_geometry(s)[has], s$geom_lnst[has], by_element = TRUE) |>
      units::set_units(km) |>
      units::drop_units()

  s$is_dist_dec_lnst_within_max <- s$dist_dec_lnst_km <= max_dist_dec_lnst_km
  s$geom_lnst <- NULL
  st_agr(s) <- "constant"

  s |>
    st_join(grid_ref, join = st_intersects) |>
    st_drop_geometry() |>
    select(.row_id, dist_dec_lnst_km, is_dist_dec_lnst_within_max, grid_site)
}

d_bind <- d_bind |> mutate(.row_id = row_number())

# The batch results are three columns plus an id — assign them into d_bind by
# column rather than joining. A `right_join` back would allocate a second copy of
# the whole ~30M-row frame, which on a 24 GB machine is the very thing that
# killed the previous run.
# only the columns dist_one_cruise() reads: group_split() COPIES what it splits, and
# splitting the whole ~110-column frame doubled R's peak (2026-09-23: 24 GB limit)
d_dist <- d_bind |>
  select(.row_id, cruise_key, longitude, latitude, lon_lnst, lat_lnst) |>
  group_split(cruise_key) |>
  purrr::map(dist_one_cruise) |>
  bind_rows() |>
  arrange(.row_id)

stopifnot(
  "distance batches must cover every scan exactly once, in order" =
    identical(d_dist$.row_id, d_bind$.row_id))

d_bind$dist_dec_lnst_km            <- d_dist$dist_dec_lnst_km
d_bind$is_dist_dec_lnst_within_max <- d_dist$is_dist_dec_lnst_within_max
d_bind$grid_site                   <- d_dist$grid_site
rm(d_dist)

pts <- d_bind |> select(-.row_id)

13.1 View sample cruise with excess distance points

Code

badcr_id <- "1507OC"

# `pts` is a plain tibble now (see pts_distance_from_lnst) — geometry is rebuilt
# here, on the sampled subset this map actually draws, rather than being carried
# on all ~30M scans just so two previews can use it
as_pts_sf <- function(d) st_as_sf(
  d, coords = c("longitude", "latitude"), remove = FALSE, crs = 4326)

# the extent from the lon/lat COLUMNS. st_bbox() would have to build geometry for
# every row first, which is the allocation this chunk exists to avoid — and the
# bounding box of a point set is just four numbers.
pts_bbox <- function(d) st_bbox(
  c(xmin = min(d$longitude, na.rm = TRUE), ymin = min(d$latitude, na.rm = TRUE),
    xmax = max(d$longitude, na.rm = TRUE), ymax = max(d$latitude, na.rm = TRUE)),
  crs = st_crs(4326))

pts_badcr <- pts |>
  filter(cruise_key == badcr_id)

bb_badcr <- pts_badcr |> pts_bbox() |> st_as_sfc()

mapView(bb_badcr) +
  mapView(
    pts_badcr |>
      filter(is_dist_dec_lnst_within_max) |>
      slice_sample(n = 1000) |>
      bind_rows(
        pts_badcr |>
          filter(!is_dist_dec_lnst_within_max)
      ) |>
      as_pts_sf() |>
      select(
        cruise_key,
        datetime_start_utc,
        longitude,
        latitude,
        sta_id,
        line,
        sta,
        dist_dec_lnst_km,
        is_dist_dec_lnst_within_max
      ),
    layer.name = glue(
      "Cruise {badcr_id}<br>distance (km)<br>lon/lat to line/station"
    ),
    zcol = "dist_dec_lnst_km",
    cex = 5,
    alpha = 0.5
  )

13.2 Filter points

Code

d_pts_cruise_filt_smry <- pts |>
  group_by(cruise_key) |>
  filter(any(!is_dist_dec_lnst_within_max)) |>
  summarize(
    n_all = n(),
    n_outside_grid = sum(is.na(grid_site)),
    pct_outside_grid = n_outside_grid / n_all,
    n_gt_cutoff = sum(!is_dist_dec_lnst_within_max, na.rm = T),
    avg_dist_gt_cutoff_km = if_else(
      n_gt_cutoff > 0,
      mean(dist_dec_lnst_km[!is_dist_dec_lnst_within_max], na.rm = T),
      NA_real_
    ),
    pct_gt_cutoff = n_gt_cutoff / n_all,
    n_rm = sum(!is_dist_dec_lnst_within_max | is.na(grid_site), na.rm = T),
    pct_rm = n_rm / n_all,
    .groups = "drop"
  ) |>
  filter(n_rm > 0) |>
  arrange(desc(pct_rm))

pts_filt <- pts |>
  filter(
    is_dist_dec_lnst_within_max,
    !is.na(grid_site)
  )

# A station number written with only three digits is not "far from its station",
# it is a different station: 100/110/120 shipped as 000.0/010.0/020.0 in five
# preliminary files (CalCOFI/workflows#104), and the distance filter above dropped
# every such cast as a position error — quietly, one row in a table of cruises.
# [fill_lonlat] repairs the ones the cast's position confirms; this is the backstop:
# any dropped scan that still carries the signature (an inshore number < 30 whose
# cast sits within the cutoff of station + 100) stops the render.
d_trunc <- pts |>
  filter(!is_dist_dec_lnst_within_max, !is.na(sta), !is.na(line),
         suppressWarnings(as.numeric(sta)) < 30) |>
  distinct(cruise_key, `_source_file`, line, sta, longitude, latitude)
if (nrow(d_trunc) > 0) {
  xy100 <- sf_project(
    from = "+proj=calcofi", to = "+proj=longlat +datum=WGS84",
    pts  = cbind(as.numeric(d_trunc$line), as.numeric(d_trunc$sta) + 100))
  hav_km <- function(lon1, lat1, lon2, lat2) {
    r <- pi / 180
    a <- sin((lat2 - lat1) * r / 2)^2 +
      cos(lat1 * r) * cos(lat2 * r) * sin((lon2 - lon1) * r / 2)^2
    2 * 6371.0088 * asin(sqrt(a))
  }
  d_trunc$dist_sta_plus100_km <- hav_km(d_trunc$longitude, d_trunc$latitude,
                                        xy100[, 1], xy100[, 2])
  d_trunc <- filter(d_trunc, dist_sta_plus100_km <= max_dist_dec_lnst_km)
}
if (nrow(d_trunc) > 0)
  stop(glue(
    "{nrow(distinct(d_trunc, `_source_file`, line, sta))} station(s) look TRUNCATED to three ",
    "digits (the cast sits at station + 100): ",
    paste(unique(glue("{basename(d_trunc$`_source_file`)} {d_trunc$line}/{d_trunc$sta}")),
          collapse = "; "),
    ". Ask the CTD team for a corrected file; do not widen the distance filter."))

d_pts_cruise_filt_smry |>
  dt(
    caption = glue(
      "Cruises with rows filtered: outside CalCOFI grid or exceeded ",
      "{max_dist_dec_lnst_km} km cutoff from line/station coordinates."
    ),
    escape = F,
    fname = "ctd_cruises_distance_from_lnst"
  ) |>
  formatCurrency(
    c(
      "n_all",
      "n_gt_cutoff",
      "avg_dist_gt_cutoff_km",
      "n_outside_grid",
      "n_rm"
    ),
    currency = "",
    digits = 0,
    mark = ","
  ) |>
  formatPercentage(
    c("pct_gt_cutoff", "pct_outside_grid", "pct_rm"),
    digits = 2
  )

13.3 View filtered points

Code

bb_cr <- pts_filt |> pts_bbox() |> st_as_sfc()

mapView(bb_cr) +
  mapView(
    # 20 per cruise, not 100: at 100 the widget reached 11.7 MB on ONE line, and pandoc
    # sat on it for 45+ minutes at full CPU (2026-09-24) after every R chunk had finished
    pts_filt |>
      group_by(cruise_key) |>
      slice_sample(n = 20) |>
      ungroup() |>
      as_pts_sf() |>
      select(
        cruise_key,
        longitude,
        latitude,
        sta_id,
        line,
        sta,
        dist_dec_lnst_km,
        datetime_start_utc
      ),
    zcol = "cruise_key",
    cex = 5,
    alpha = 0.5
  )

13.4 Check for duplicates by datetime and depth

Code

d_dupes <- pts_filt |>
  select(datetime_start_utc, depth) |>
  group_by(datetime_start_utc, depth) |>
  summarize(n = n(), .groups = "drop") |>
  filter(n > 1) |>
  arrange(desc(n))

n_distinct_dtime_depth <- pts_filt |>
  select(datetime_start_utc, depth) |>
  n_distinct()

d_dupes |>
  slice(c(1:10, (n() - 9):n())) |>
  dt(
    caption = glue(
      "First and last 10 duplicates by datetime_start_utc, depth ",
      "(nrows = {format(nrow(d_dupes), big.mark = ',')}; ",
      "n_distinct = {format(n_distinct_dtime_depth, big.mark = ',')})"
    ),
    fname = "ctd_duplicates_by_datetime_depth"
  )

14 Rename Fields

Apply field renames from flds_redefine.csv.

Code

# apply column renames from flds_redefine
renames <- d_flds_rd |>
  filter(fld_old != fld_new) |>
  select(fld_old, fld_new)

# `pts_filt` is already a plain tibble and geom_lnst no longer exists — the
# distance chunk drops both rather than carrying geometry across ~30M scans
d_filt <- pts_filt

for (i in seq_len(nrow(renames))) {
  old <- renames$fld_old[i]
  new <- renames$fld_new[i]
  if (old %in% names(d_filt)) {
    d_filt <- d_filt |> rename(!!new := !!old)
  }
}

# drop intermediate columns
d_filt <- d_filt |>
  select(
    -any_of(c(
      "date_time_pst",
      "date_time_utc",
      "date_time_format",
      "is_lon_pseudoNA",
      "is_lat_pseudoNA",
      "is_lnsta_pseudoNA",
      "lon_lnst",
      "lat_lnst",
      "dist_dec_lnst_km",
      "is_dist_dec_lnst_within_max",
      "grid_site"
    ))
  )

# write raw table for pivoting
dbWriteTable(con, "ctd_raw", d_filt, overwrite = TRUE)
# which layout each source file is in, for the provider-average cross-check
dbWriteTable(con, "ctd_src_format", read_rds(glue("{dir_tmp}/d_src_format.rds")),
             overwrite = TRUE)

# The scans live in DuckDB from here on. Release R's copies: d_bind, pts, pts_filt and
# d_filt are each a ~30M-row frame, and keeping them made the R process 35 GB by
# [add_spatial] on 2026-09-24, with DuckDB's own 2 GB on top, until the machine ran out
# of swap. Nothing after this point reads them.
rm(list = intersect(c("d_bind", "pts", "pts_filt", "d_filt", "d_csv"), ls()))
invisible(gc())
say(glue("Loaded {nrow(d_filt)} rows into ctd_raw"))

15 De-duplicate ctd_raw

Some cruises ship the same casts twice in the download archive — a primary copy plus a reorganized copy in a sub-folder (e.g. db_csvs/ alongside db_csvs/south-north/), both classified final by the path regex, so both pass the priority filter (see the duplicates reported above). The copies are identical except for re-rounded coordinates (and occasionally a re-processed measurement value). Collapse to one row per measurement position so the deterministic UUIDs downstream (ctd_cast_uuid, ctd_measurement_uuid, …) are unique.

Code

raw_cols <- dbGetQuery(
  con,
  "SELECT column_name FROM information_schema.columns WHERE table_name = 'ctd_raw'"
)$column_name

# a measurement position is identified by these five columns; duplicate copies
# share all five and differ only in re-rounded coordinates / re-processed values
dedup_key  <- c("cruise_key", "cast_key", "cast_dir", "datetime_start_utc", "depth_m")
stopifnot(
  "ctd_raw missing expected columns for de-duplication" =
    all(c(dedup_key, "_source_file", "latitude", "longitude") %in% raw_cols))

# keep the primary copy: reorganized duplicates sit in deeper source paths, so
# prefer the shortest `_source_file`; coordinates tie-break for full determinism
n_before <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM ctd_raw")$n
dbExecute(con, glue("
  CREATE OR REPLACE TABLE ctd_raw AS
  SELECT * FROM ctd_raw
  QUALIFY ROW_NUMBER() OVER (
    PARTITION BY {paste(dedup_key, collapse = ', ')}
    ORDER BY length(\"_source_file\"), \"_source_file\", latitude, longitude) = 1"))
n_after <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM ctd_raw")$n

cat(glue(
  "ctd_raw de-duplicated: removed {format(n_before - n_after, big.mark = ',')} ",
  "duplicate-position rows ",
  "({round(100 * (n_before - n_after) / max(n_before, 1), 1)}%); ",
  "{format(n_after, big.mark = ',')} rows remain"), "\n")

Code

# must always evaluate: line ~217 may have set eval=FALSE to skip read+bind
# when restoring from checkpoint; this chunk re-enables computation through
# write_parquet. without `eval: true` it would inherit eval=FALSE and be
# skipped, so write_parquet never runs yet downstream chunks expect its output.
if (parquet_complete) {
  knitr::opts_chunk$set(eval = FALSE)
  say("Parquet complete — skipping to upload")
} else {
  knitr::opts_chunk$set(eval = TRUE)
}

16 Cross-Dataset Bridge

Load reference tables (ship, cruise, grid) from the ichthyo ingest, derive ship_key, and validate against known ships/cruises.

Code

# tables declared in calcofi.modifies — loaded as TABLEs (writable)
modifies_tables <- c("ship")  # from calcofi.modifies in YAML frontmatter
load_prior_tables(
  con         = con,
  tables      = modifies_tables,
  parquet_dir = cc_stage_path("parquet", "swfsc_ichthyo")
)

# A tibble: 1 × 3
  table  rows has_geom
  <chr> <dbl> <lgl>   
1 ship     48 FALSE

Code

# cruise + grid as VIEWs (read-only: cruise used for JOIN, grid for assign_grid_key)
load_prior_tables(
  con         = con,
  tables      = c("cruise", "grid"),
  parquet_dir = cc_stage_path("parquet", "swfsc_ichthyo"),
  as_view     = TRUE
)

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

Code

# snapshot PKs of modified tables before modifications
modifies_pks <- list()
for (tbl in modifies_tables) {
  pk_col <- dbGetQuery(con, glue(
    "SELECT column_name FROM information_schema.columns
     WHERE table_name = '{tbl}' ORDER BY ordinal_position LIMIT 1"))$column_name
  modifies_pks[[tbl]] <- list(
    pk_col = pk_col,
    keys   = dbGetQuery(con, glue("SELECT {pk_col} FROM {tbl}"))[[1]])
}

# apply cruise_key corrections for known study-column data quality errors.
# runs unconditionally so it fires for both fresh runs (YYMMKK keys) and
# checkpoint-restored runs (already-converted YYYY-MM-NODC keys).
for (i in seq_len(nrow(d_cruise_corrections))) {
  n_wrong <- dbGetQuery(con, glue(
    "SELECT COUNT(*) AS n FROM ctd_raw
     WHERE cruise_key = '{d_cruise_corrections$cruise_key_raw[i]}'"))$n
  if (n_wrong > 0) {
    dbExecute(con, glue(
      "UPDATE ctd_raw SET cruise_key = '{d_cruise_corrections$cruise_key_correct[i]}'
       WHERE cruise_key = '{d_cruise_corrections$cruise_key_raw[i]}'"))
    say(glue(
      "cruise_key correction: {d_cruise_corrections$cruise_key_raw[i]} → ",
      "{d_cruise_corrections$cruise_key_correct[i]} ({format(n_wrong, big.mark=',')} rows)"))
  }
}

# detect whether cruise_key already in YYYY-MM-NODC format (from checkpoint)
sample_ck <- dbGetQuery(
  con, "SELECT cruise_key FROM ctd_raw LIMIT 1")$cruise_key
cruise_key_already_converted <- grepl("^\\d{4}-\\d{2}-", sample_ck)

if (!cruise_key_already_converted) {
  # 0404NHJD is a combined New Horizon + Jordan Davis cruise; use primary ship NH
  dbExecute(
    con,
    "UPDATE ctd_raw SET cruise_key = LEFT(cruise_key, 6)
     WHERE LENGTH(cruise_key) > 6"
  )

  # derive ship_key from last 2 chars of cruise_key (YYMMKK format from filename)
  dbExecute(con, "ALTER TABLE ctd_raw ADD COLUMN IF NOT EXISTS ship_key VARCHAR")
  dbExecute(con, "UPDATE ctd_raw SET ship_key = RIGHT(cruise_key, 2)")

  # validate ship_key against ship reference table
  ctd_ships <- dbGetQuery(con, "SELECT DISTINCT ship_key FROM ctd_raw")
  ref_ships <- dbGetQuery(con, "SELECT ship_key FROM ship")
  orphan_ships <- setdiff(ctd_ships$ship_key, ref_ships$ship_key)

  if (length(orphan_ships) > 0) {
    say(glue(
      "{length(orphan_ships)} ship_key(s) in CTD not in ship table: ",
      "{paste(orphan_ships, collapse = ', ')}"
    ))

    # run match_ships() with centralized renames for orphan ships
    orphan_tbl <- dbGetQuery(con, glue(
      "SELECT DISTINCT ship_key AS ship_code, NULL AS ship_name
       FROM ctd_raw
       WHERE ship_key IN ({paste(shQuote(orphan_ships, type = 'sh'), collapse = ', ')})")) |>
      mutate(ship_name = as.character(ship_name))

    ship_result <- match_ships(
      unmatched_ships  = orphan_tbl,
      reference_ships  = dbReadTable(con, "ship"),
      ship_renames_csv = here("metadata/ship_renames.csv"),
      fetch_ices       = FALSE)

    # insert interim entries for still-unmatched ships (ship_nodc = "?XX?")
    ensure_interim_ships(con, ship_result)

    dbGetQuery(
      con,
      glue(
        "SELECT ship_key, COUNT(DISTINCT cruise_key) AS n_cruises,
                COUNT(*) AS n_rows
         FROM ctd_raw
         WHERE ship_key IN ({paste(shQuote(orphan_ships, type = 'sh'), collapse = ', ')})
         GROUP BY ship_key
         ORDER BY ship_key"
      )
    ) |>
      dt(
        caption = "Orphan ship_keys (in CTD but not in ship table)",
        fname = "ctd_orphan_ship_keys"
      )
  } else {
    say("All CTD ship_keys found in ship reference table")
  }

  # convert cruise_key from YYMMKK → YYYY-MM-NODC format
  convert_cruise_key_format(con, "ctd_raw", old_key_col = "cruise_key")
  dbExecute(con, "ALTER TABLE ctd_raw DROP COLUMN cruise_key")
  dbExecute(con, "ALTER TABLE ctd_raw RENAME COLUMN cruise_key_new TO cruise_key")
} else {
  say("cruise_key already in YYYY-MM-NODC format (from checkpoint)")
}

cruise_key already in YYYY-MM-NODC format (from checkpoint)

Code

# validate cruise_key against cruise reference table
ctd_cruises <- dbGetQuery(con, "SELECT DISTINCT cruise_key FROM ctd_raw")
ref_cruises <- dbGetQuery(con, "SELECT cruise_key FROM cruise")
orphan_cruises <- setdiff(ctd_cruises$cruise_key, ref_cruises$cruise_key)

if (length(orphan_cruises) > 0) {
  say(glue(
    "{length(orphan_cruises)} cruise_key(s) in CTD not in cruise table: ",
    "{paste(head(orphan_cruises, 10), collapse = ', ')}",
    "{if (length(orphan_cruises) > 10) '...' else ''}"
  ))
  dbGetQuery(
    con,
    glue(
      "SELECT cruise_key, COUNT(*) AS n_rows
       FROM ctd_raw
       WHERE cruise_key IN ({paste(shQuote(orphan_cruises, type = 'sh'), collapse = ', ')})
       GROUP BY cruise_key
       ORDER BY cruise_key"
    )
  ) |>
    dt(
      caption = "Orphan cruise_keys (in CTD but not in cruise table)",
      fname = "ctd_orphan_cruise_keys"
    )
} else {
  say("All CTD cruise_keys found in cruise reference table")
}

24 cruise_key(s) in CTD not in cruise table: 2024-01-33UD, 2026-07-3322, 2018-10-33P4, 2020-07-33P4, 2020-10-33P4, 2025-02-33UD, 2024-04-3322, 2024-08-33P4, 2021-11-33P4, 2016-11-33P4...

17 Save Checkpoint

Save a DuckDB checkpoint after the expensive read+bind+filter+bridge steps so subsequent runs can resume from here.

Code

if (!file_exists(db_checkpoint) || overwrite) {
  close_duckdb(con)
  file_copy(db_path, db_checkpoint, overwrite = TRUE)
  write_input_fingerprint(checkpoint_fp_path, fp_zips)
  con <- get_duckdb_con(db_path)
  duck_budget(con)
  load_duckdb_extension(con, "spatial")
  load_duckdb_extension(con, "icu")
  say(glue("Saved checkpoint: {db_checkpoint}"))
} else {
  say(glue("Checkpoint already exists: {db_checkpoint}"))
}

Saved checkpoint: /Users/bbest/Github/CalCOFI/workflows-ctd-new-format/data/wrangling/calcofi_ctd-cast_checkpoint.duckdb

18 Measurement Column Registry

Resolve which ctd_raw columns are measurements and which are their quality flags, from the shared metadata/measurement_type.csv registry. This split drives everything downstream: measurement columns pivot into ctd_measurement, and whatever is left over is by definition cast-level and becomes ctd_cast. Getting it wrong is how a per-scan column once leaked into ctd_cast and broke its primary key (workflows#53).

ctd_wide was retired here

This section used to also build ctd_wide, a ~1 GB wide-format copy of ctd_raw for ERDDAP. It has no consumer: ERDDAP now serves CTD through EDDTableFromDatabase over DuckDB views plus the netCDF files built below, and the wide file’s whole-file heap read is what OOM’d ERDDAP in the first place (see the serving benchmark). The column-registry logic it carried is load-bearing, so it stays; the table itself is gone.

Code

# identify measurement and quality columns from measurement_type metadata
d_meas_ctd <- d_meas_type |>
  filter(str_detect(`_source_datasets`, "calcofi_ctd-cast"))

# --- every series names the PROVIDER'S OWN flag column (CalCOFI/workflows#105) ---
# The 2026 processing software flags each corrected, averaged and derived series
# (SaltAve_CorrQ, EstChl_StaCorrQ, …) and the registry now points each series at that
# column (libs/build_ctd_measurement_registry.R § 4). A legacy-layout file has no such
# column, so fill it from the sensor flag the series inherited before 2026-09-23:
# COALESCE(provider flag, sensor flag) — the provider's word wherever it has one, and
# never weaker than the old rule where it does not. An average has no sensor to fall
# back on (NA): it is rebuilt from its pair in [combine_sensor_pairs] and keeps the
# provider's flag on it only as the provider's own verdict.
qual_fallback <- c(
  temp_ave_q            = NA,        salt_ave_corr_q       = NA,
  salt1_corr_q          = "salt1q",  salt2_corr_q          = "salt2q",
  ox1_sta_corr_q        = "ox1q",    ox2_sta_corr_q        = "ox2q",
  ox1_cruise_corr_q     = "ox1q",    ox2_cruise_corr_q     = "ox2q",
  ox_ave_sta_corr_q     = NA,
  ox1u_mq               = "ox1q",    ox2u_mq               = "ox2q",
  ox1u_m_sta_corr_q     = "ox1q",    ox2u_m_sta_corr_q     = "ox2q",
  ox1u_m_cruise_corr_q  = "ox1q",    ox2u_m_cruise_corr_q  = "ox2q",
  ox_aveu_m_sta_corr_q  = NA,
  est_chl_sta_corr_q    = "fluor_q", est_chl_cruise_corr_q = "fluor_q",
  est_no3_sta_corr_q    = "isusq",   est_no3_cruise_corr_q = "isusq",
  batq                  = NA,        dyn_ht_q              = NA,
  svaq                  = NA,
  po_t1q                = "temp1q",  po_t2q                = "temp2q",
  ox_sat1q              = "ox1q",    ox_sat2q              = "ox2q")
stopifnot(
  "every provider flag column the registry names has a fallback rule here" =
    all(setdiff(na.omit(d_meas_ctd$`_qual_column`),
                c("", "pr_q", "temp1q", "temp2q", "salt1q", "salt2q", "ox1q", "ox2q",
                  "fluor_q", "isusq", "sig_theta_ts1q", "sig_theta_ts2q", "trans_q",
                  "p_hq", "sparq", "parq")) %in% names(qual_fallback)))
for (q in names(qual_fallback))
  dbExecute(con, glue("ALTER TABLE ctd_raw ADD COLUMN IF NOT EXISTS {q} DOUBLE"))
q_fill <- qual_fallback[!is.na(qual_fallback)]
dbExecute(con, glue("UPDATE ctd_raw SET {paste0(names(q_fill), ' = COALESCE(TRY_CAST(', names(q_fill), ' AS DOUBLE), TRY_CAST(', q_fill, ' AS DOUBLE))', collapse = ', ')}"))

[1] 8442099

Code

cat(glue("provider flag columns: {length(qual_fallback)} ({length(q_fill)} filled from their sensor flag where the file predates them)\n"))

provider flag columns: 27 (20 filled from their sensor flag where the file predates them)

Code

meas_cols <- d_meas_ctd$`_source_column` |> na.omit() |> as.character()
meas_cols <- meas_cols[meas_cols != ""]
qual_cols <- d_meas_ctd$`_qual_column` |> na.omit() |> as.character()
qual_cols <- qual_cols[qual_cols != ""]

raw_cols <- dbGetQuery(
  con,
  "SELECT column_name FROM information_schema.columns WHERE table_name = 'ctd_raw'"
)$column_name

# DIAGNOSTIC: a measurement column in the registry that is absent from ctd_raw is
# a silent no-op in the pivot; a numeric ctd_raw column that is NOT registered
# falls through to cast_cols and can break the ctd_cast PK. Report both.
meas_missing <- setdiff(meas_cols, raw_cols)
cast_leftover <- setdiff(raw_cols, c(meas_cols, qual_cols, "depth_m", "_source_file"))

cat(glue("registered measurement types : {nrow(d_meas_ctd)}\n"))

registered measurement types : 55

Code

cat(glue("  canonical (-> ctd_thin)    : {sum(d_meas_ctd$is_canonical, na.rm = TRUE)}\n"))

canonical (-> ctd_thin)    : 44

Code

cat(glue("  quality-flag columns       : {length(qual_cols)}\n"))

quality-flag columns       : 42

Code

cat(glue("registry columns absent from ctd_raw: {length(meas_missing)}",
         if (length(meas_missing)) glue(" -> {paste(meas_missing, collapse = ', ')}") else "", "\n"))

registry columns absent from ctd_raw: 0

Code

cat(glue("cast-level columns (-> ctd_cast)    : {length(cast_leftover)}\n"))

cast-level columns (-> ctd_cast)    : 14

Code

# the canonical set is what ctd_thin and the thinned netCDF carry, so show it
# explicitly rather than leaving it implicit in a flag column
d_meas_ctd |>
  select(measurement_type, units, is_canonical, source_column = `_source_column`) |>
  arrange(desc(is_canonical), measurement_type) |>
  datatable(
    rownames  = FALSE,
    caption   = "CTD measurement registry — is_canonical drives ctd_thin membership",
    options   = list(pageLength = 10, scrollX = TRUE, dom = "ftip"))

19 Split into Tidy Tables

19.1 ctd_cast — cast-level metadata

Code

# d_meas_ctd, meas_cols, qual_cols come from the Measurement Column Registry above

# use ctd_raw columns (includes ship_key and data_stage added via SQL)
raw_cols <- dbGetQuery(
  con,
  "SELECT column_name FROM information_schema.columns WHERE table_name = 'ctd_raw'"
)$column_name

cast_cols <- setdiff(
  raw_cols,
  c(meas_cols, qual_cols, "depth_m", "_source_file")
)

# extract one row per cast, keyed (cruise_key, cast_key, cast_dir,
# datetime_start_utc). within that key the only cast-level column that varies is
# lat/lon (sub-meter GPS jitter between same-second scans), so order by it for a
# deterministic pick. measurement-like columns (e.g. ox_aveu_m_sta_corr) must be
# registered in measurement_type.csv so they pivot to ctd_measurement and never
# reach cast_cols — otherwise SELECT DISTINCT would emit multiple rows per cast.
dbExecute(
  con,
  glue(
    "
  CREATE OR REPLACE TABLE ctd_cast AS
  SELECT {paste(cast_cols, collapse = ', ')}
  FROM ctd_raw
  QUALIFY ROW_NUMBER() OVER (
    PARTITION BY cruise_key, cast_key, cast_dir, datetime_start_utc
    ORDER BY latitude, longitude) = 1
  ORDER BY datetime_start_utc, cruise_key"
  )
)

[1] 7485734

Code

n_cast <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM ctd_cast")$n
say(glue("ctd_cast: {format(n_cast, big.mark = ',')} rows"))

ctd_cast: 7,485,734 rows

Code

# assign deterministic UUID from composite natural key (md5-based, DuckDB-native)
assign_deterministic_uuids_md5(
  con = con,
  table_name = "ctd_cast",
  id_col = "ctd_cast_uuid",
  key_cols = c("cruise_key", "cast_key", "cast_dir", "datetime_start_utc")
)

# guard: ctd_cast_uuid must be a valid unique PK — fail loudly if a
# measurement-like column ever leaks into cast_cols again (see #53)
stopifnot(
  "ctd_cast_uuid is not unique" =
    dbGetQuery(
      con,
      "SELECT COUNT(*) - COUNT(DISTINCT ctd_cast_uuid) AS n FROM ctd_cast"
    )$n == 0)

19.2 ctd_measurement — long-format sensor readings

Code

# add ctd_cast_uuid to ctd_raw for FK linkage
dbExecute(
  con,
  "
  ALTER TABLE ctd_raw ADD COLUMN IF NOT EXISTS ctd_cast_uuid VARCHAR"
)

[1] 0

Code

dbExecute(
  con,
  "
  UPDATE ctd_raw AS r
  SET ctd_cast_uuid = c.ctd_cast_uuid
  FROM ctd_cast AS c
  WHERE r.cruise_key   = c.cruise_key
    AND r.cast_key     = c.cast_key
    AND r.cast_dir     = c.cast_dir
    AND r.datetime_start_utc = c.datetime_start_utc"
)

[1] 8442099

Code

# --- the two-sensor averages are NOT repaired here -----------------------------
# The source's pre-computed averages (TempAve, SaltAve_Corr, OxAve_StaCorr,
# OxAveuM_StaCorr) are rebuilt from their two sensors by the PROVIDER'S FLAGS in
# § Combine the Sensor Pairs by Their Flags, after the CTD team's accepted flags are
# applied — Rasmus Swalethorp's rule (2026-09-09), calcofi4db::combine_sensor_pair().
# Until 2026-09-11 this spot "repaired" them by the registry's physical bounds alone,
# which let a sensor the source flags 8 (questionable) or 9 (bad) into the average
# whenever its value was physically possible (1,425 temperature and 6,001 salinity
# averages in v2026.09.10) and never saw the accepted flags at all. A physical bound
# is a sentinel guard, not a quality test.

# pivot wide-to-long one measurement type at a time to avoid OOM
dbExecute(con, "DROP TABLE IF EXISTS ctd_measurement")

[1] 0

Code

created <- FALSE
for (i in seq_len(nrow(d_meas_ctd))) {
  row <- d_meas_ctd[i, ]
  src_col <- row$`_source_column`
  meas_type <- row$measurement_type
  qual_col <- row$`_qual_column`
  # a series the ingest BUILDS (oxygen_ml_l_ave_cruise_corr, from its sensor pair in
  # [combine_sensor_pairs]) has no source column: nothing to pivot
  if (is.na(src_col) || src_col == "") next

  # a bare NULL is typed INT32 by DuckDB, so if the FIRST measurement type has no
  # qual column the CREATE TABLE below makes measurement_qual an integer column —
  # and the first later type that does have one fails inserting 'nan'. Type it.
  #
  # The source stores these codes in DOUBLE columns, so a plain cast yields "9.0"
  # rather than "9" and no longer matches the controlled vocabulary in
  # metadata/measurement_qual.csv (the CTD set: 0/1/2/8/9). Strip a trailing ".0" TEXTUALLY
  # rather than casting through INTEGER: a numeric cast would silently round an
  # unexpected "9.5" to 9, whereas this leaves anything non-integral verbatim so it
  # shows up in the audit below instead of being quietly normalized away.
  qual_expr <- if (!is.na(qual_col) && qual_col != "") {
    glue("regexp_replace(CAST({qual_col} AS VARCHAR), '\\.0+$', '')")
  } else {
    "CAST(NULL AS VARCHAR)"
  }

  sql_select <- glue(
    "
    SELECT ctd_cast_uuid, depth_m,
      '{meas_type}' AS measurement_type,
      CAST({src_col} AS DOUBLE) AS measurement_value,
      {qual_expr} AS measurement_qual
    FROM ctd_raw
    WHERE {src_col} IS NOT NULL
      AND NOT isnan(CAST({src_col} AS DOUBLE))
      AND isfinite(CAST({src_col} AS DOUBLE))
    "
  )

  if (!created) {
    dbExecute(con, glue("CREATE OR REPLACE TABLE ctd_measurement AS {sql_select}"))
    created <- TRUE
  } else {
    dbExecute(con, glue("INSERT INTO ctd_measurement {sql_select}"))
  }
  say(glue("  {i}/{nrow(d_meas_ctd)}: {meas_type}"))
}

  1/55: beam_attenuation
  2/55: btl_ammonium
  3/55: btl_chlorophyll_a
  4/55: btl_depth
  5/55: btl_nitrate
  6/55: btl_nitrite
  7/55: btl_phaeopigment
  8/55: btl_phosphate
  9/55: btl_silicate
  10/55: btl_temperature
  11/55: dynamic_height
  12/55: est_chlorophyll_a_cruise_corr
  13/55: est_chlorophyll_a_sta_corr
  14/55: est_nitrate_cruise_corr
  15/55: est_nitrate_sta_corr
  16/55: fluorescence_v
  17/55: isus_v
  18/55: oxygen_btl_ml_l
  19/55: oxygen_btl_umol_kg
  20/55: oxygen_ml_l_1
  21/55: oxygen_ml_l_1_cruise_corr
  22/55: oxygen_ml_l_1_sta_corr
  23/55: oxygen_ml_l_2
  24/55: oxygen_ml_l_2_cruise_corr
  25/55: oxygen_ml_l_2_sta_corr
  27/55: oxygen_ml_l_ave_sta_corr
  28/55: oxygen_saturation_1
  29/55: oxygen_saturation_2
  30/55: oxygen_umol_kg_1
  31/55: oxygen_umol_kg_1_cruise_corr
  32/55: oxygen_umol_kg_1_sta_corr
  33/55: oxygen_umol_kg_2
  34/55: oxygen_umol_kg_2_cruise_corr
  35/55: oxygen_umol_kg_2_sta_corr
  36/55: oxygen_umol_kg_ave_sta_corr
  37/55: par
  38/55: ph
  39/55: potential_temperature_1
  40/55: potential_temperature_2
  41/55: pressure
  42/55: salinity_1
  43/55: salinity_1_corr
  44/55: salinity_2
  45/55: salinity_2_corr
  46/55: salinity_ave_corr
  47/55: salinity_btl
  48/55: sigma_theta_1
  49/55: sigma_theta_2
  50/55: spar
  51/55: specific_volume_anomaly
  52/55: temperature_1
  53/55: temperature_2
  54/55: temperature_ave
  55/55: transmissometer

Code

n_meas_raw <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM ctd_measurement")$n
say(glue("ctd_measurement (pre-sentinel): {format(n_meas_raw, big.mark = ',')} rows"))

ctd_measurement (pre-sentinel): 297,896,021 rows

Code

# --- strip the -99 missing-value sentinel from MEASUREMENTS -------------------
# -99 is this source's documented missing marker. The notebook already knows that:
# `pseudoNA_values <- c(-9.99e-29, -99)` in the fill_lonlat chunk strips it from
# longitude/latitude. It was never applied to the measurement columns, so -99 has
# been flowing through ctd_measurement -> ctd_thin -> obs as though it were a real
# reading — 84,302 rows in release v2026.07.17, including CANONICAL oxygen
# (oxygen_ml_l_ave_sta_corr 953, oxygen_umol_kg_ave_sta_corr 4,294) alongside
# isus_v 40,479, ph 31,493 and spar 6,189.
#
# The NOT isnan / isfinite guard in the pivot above cannot catch it: -99 is a
# perfectly finite double. A -99 mL/L oxygen silently corrupts any mean, minimum
# or anomaly a consumer computes, which is why this is a DELETE and not merely a
# report — a missing value has no business existing as a row in long format.
#
# Matched exactly (= -99, not a tolerance): every observed case is exactly -99.00,
# and an exact test cannot swallow a real reading that merely rounds near it. The
# two genuinely signed types (dynamic_height, specific_volume_anomaly) are
# reported separately below so a real -99 there would be visible rather than
# assumed away.
SENTINEL <- -99

sentinel_tally <- dbGetQuery(con, glue("
  SELECT measurement_type, COUNT(*) AS n,
         COUNT(DISTINCT ctd_cast_uuid) AS n_casts
  FROM ctd_measurement WHERE measurement_value = {SENTINEL}
  GROUP BY 1 ORDER BY n DESC"))

n_sentinel <- sum(sentinel_tally$n)
dbExecute(con, glue(
  "DELETE FROM ctd_measurement WHERE measurement_value = {SENTINEL}"))

[1] 25082834

Code

n_meas <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM ctd_measurement")$n
say(glue("-99 sentinel rows removed: {format(n_sentinel, big.mark = ',')} ",
         "across {nrow(sentinel_tally)} measurement type(s)"))

-99 sentinel rows removed: 25,082,834 across 39 measurement type(s)

Code

say(glue("ctd_measurement: {format(n_meas, big.mark = ',')} rows"))

ctd_measurement: 272,813,187 rows

Code

stopifnot("sentinel delete accounting" = n_meas_raw - n_sentinel == n_meas)

# --- enforce the registry's declared physical bounds -------------------------
# metadata/measurement_type.csv has carried valid_min/valid_max since the CTD
# registry was built, and NOTHING applied them — not here, not at release time.
# So v2026.08.07 shipped ~31k impossible values: ph down to -10, oxygen_ml_l_1 to
# -79.5, salinity_2 up to 46.3, sigma_theta_1 to -50.8, temperature_ave to -47.6.
#
# The -99 DELETE above only catches the sentinel ITSELF. It cannot catch a
# sentinel that has been scaled, offset or averaged upstream, and it says nothing
# about a genuinely wild reading. A declared bound is the general guard.
#
# This was ~40 lines of inline SQL that ran for this dataset alone. It is now
# calcofi4db::check_measurement_bounds() / drop_out_of_bounds(), which every
# ingest calls and release_database.qmd re-checks — because the same defect was
# live in calcofi_mets.sw_ph (494 values at -99, 16.6% of the type) with its
# bounds declared and equally unread.
# --- declare the PAR sensors' physical bounds ---------------------------------
# `par` and `spar` were two of the nine CTD types with no bound, and v2026.09.06
# published what that leaves in at the release grain: `spar` from -3.07e17 to
# 2.01e16 µE/m²/s (31,623 values on ten cruises, 1995–2008 — 6% of the series) and
# `par` with 34,233 values at exactly 9999 — a fill, not a reading — plus 1,065 more
# above 5,000 and eight below -100 (24 cruises, 1994–2025). Photosynthetically
# available radiation at the sea surface cannot exceed the ~2,500 µE/m²/s the sun
# delivers at the top of the atmosphere, and a dark reading sits within a few units
# of zero, so -100…5,000 catches only the impossible (Ben, 2026-09-10: a certain
# bug, not a question). Idempotent, as the bottle's § Declare and Enforce Physical
# Bounds is.
ctd_bounds <- tibble::tribble(
  ~measurement_type, ~valid_min, ~valid_max,
  "par",                   -100,       5000,
  "spar",                  -100,       5000)
d_meas_type <- calcofi4db::declare_measurement_bounds(
  ctd_bounds, here("metadata/measurement_type.csv"))

bounds_pre <- check_measurement_bounds(
  con, "ctd_measurement", mt = d_meas_type)

# report BEFORE deleting: `undeclared` is a finding too, and the seven CTD types
# with no bound are the ones still to take to the provider
bounds_datatable(bounds_pre)

Code

oob_tally <- drop_out_of_bounds(con, "ctd_measurement", mt = d_meas_type)

say(glue("out-of-range rows removed: {format(sum(oob_tally$n_bad), big.mark = ',')} ",
         "across {nrow(oob_tally)} type(s), against ",
         "{sum(bounds_pre$status != 'undeclared')} declared bound(s); ",
         "{sum(bounds_pre$status == 'undeclared')} type(s) still undeclared"))

out-of-range rows removed: 1,808,472 across 38 type(s), against 49 declared bound(s); 5 type(s) still undeclared

Code

# the guard must leave nothing behind: a violation here means the delete and the
# check disagree on a bound
bounds_post <- check_measurement_bounds(
  con, "ctd_measurement", mt = d_meas_type, include_undeclared = FALSE)
stopifnot(
  "ctd_measurement must hold no value outside a declared bound" =
    sum(bounds_post$status == "out_of_range") == 0)

# assign deterministic UUID (md5-based, DuckDB-native — single SQL, no R data transfer)
assign_deterministic_uuids_md5(
  con = con,
  table_name = "ctd_measurement",
  id_col = "ctd_measurement_uuid",
  key_cols = c("ctd_cast_uuid", "depth_m", "measurement_type")
)

# add cruise_key for partitioned parquet output
dbExecute(
  con,
  "ALTER TABLE ctd_measurement ADD COLUMN IF NOT EXISTS cruise_key VARCHAR"
)

[1] 0

Code

dbExecute(
  con,
  "UPDATE ctd_measurement AS m
   SET cruise_key = c.cruise_key
   FROM ctd_cast AS c
   WHERE m.ctd_cast_uuid = c.ctd_cast_uuid"
)

[1] 271004715

Code

# remove orphaned measurements with no parent cast
n_orphan <- dbGetQuery(
  con,
  "SELECT COUNT(*) FROM ctd_measurement WHERE ctd_cast_uuid IS NULL"
)[[1]]
if (n_orphan > 0) {
  dbExecute(con, "DELETE FROM ctd_measurement WHERE ctd_cast_uuid IS NULL")
  cat(glue("{format(n_orphan, big.mark = ',')} orphaned ctd_measurement ",
           "rows removed (NULL ctd_cast_uuid)"), "\n")
}

20 Apply Accepted QC Flags

The CTD team’s PostgreSQL QA/QC database (server/postgis/init/40_ctd.sql, schema ctd) is a working store beside this release, never a dependency this pipeline connects to live — the bridge is a nightly snapshot, gs://calcofi-db/qc/ctd/flag_accepted.parquet (server/scripts/pg_flag_snapshot.R, cron pg_flag_snapshot.sh), of every ctd.flag row a curator has accepted. This chunk downloads that snapshot, joins each accepted flag to the scan it names, and overwrites ctd_measurement.measurement_qual for the match.

Join key. ctd.file.archive is the zip filename and .path is the member path inside it (libs/pg_ctd.R::pg_ctd_discover_files(): archive = paste0(dir_unzip, ".zip"), path = rel) — byte-identical to this notebook’s own `_source_file` (the d_csv chunk’s path_unzip), so archive reconstructs from it as the text before the first / plus .zip. Within one file, ctd.scan.cast_id / .depth are the verbatim source Cast_ID / Depth columns — flds_redefine.csv renames those to this notebook’s cast_key / depth_m, so (archive, _source_file, cast_key, depth_m) is the same identity as (archive, path, cast_id, depth) in the snapshot. ctd.flag.variable is a ctd.scan column name (e.g. salt1_corr), not a measurement_type — mapped through `_source_column` on d_meas_ctd (built above in “Measurement Column Registry”).

Two flag shapes. A scan-level flag carries depth/cast_id (from the snapshot’s LEFT JOIN ctd.scan); a file/cast-level flag (ctd.flag.scan_id IS NULL, the whole file or cast is being marked) does not, and applies to every scan in that file — for variable = '*' that means every measurement type, not one.

Code

flag_url   <- "https://storage.googleapis.com/calcofi-db/qc/ctd/flag_accepted.parquet"
flag_cache_dir <- here("data/cache/ctd-cast")
dir_create(flag_cache_dir)
flag_cache_path <- file.path(flag_cache_dir, "flag_accepted.parquet")

# download-first, skip when unchanged: a HEAD request's Content-Length is
# compared to the cached file's size before spending the GET. A size collision
# with genuinely different content is not a concern here — the snapshot is
# COPY'd fresh each night (server/scripts/pg_flag_snapshot.R), never edited in
# place, so "same size" only happens when the file legitimately did not change.
flag_head <- tryCatch(httr::HEAD(flag_url), error = function(e) NULL)
remote_size <- if (!is.null(flag_head) && httr::status_code(flag_head) == 200)
  suppressWarnings(as.numeric(httr::headers(flag_head)[["content-length"]])) else NA_real_
remote_last_modified <- if (!is.null(flag_head))
  httr::headers(flag_head)[["last-modified"]] else NA_character_

local_size <- if (file_exists(flag_cache_path)) as.numeric(file_size(flag_cache_path)) else NA_real_

need_download <- is.na(local_size) ||
  (!is.na(remote_size) && local_size != remote_size)

if (need_download) {
  stopifnot("could not reach flag_accepted.parquet and no local cache exists" =
              !is.null(flag_head) && httr::status_code(flag_head) == 200)
  utils::download.file(flag_url, flag_cache_path, mode = "wb", quiet = TRUE)
  say(glue("Downloaded flag_accepted.parquet ({round(file_size(flag_cache_path) / 1024, 1)} KB, ",
           "last modified {remote_last_modified %||% 'unknown'})"))
} else {
  say(glue("flag_accepted.parquet unchanged ({round(local_size / 1024, 1)} KB) — using cache"))
}

flag_accepted.parquet unchanged (0.6 KB) — using cache

Code

d_flags <- dbGetQuery(con, glue("SELECT * FROM read_parquet('{flag_cache_path}')"))
n_flags <- nrow(d_flags)

if (n_flags == 0) {

  say(glue("0 accepted flags in snapshot (last modified {remote_last_modified %||% 'unknown'})"))

} else {

  # variable (a ctd.scan column name, e.g. "salt1_corr") -> measurement_type,
  # via the same crosswalk the pivot above used (d_meas_ctd$`_source_column`)
  d_var_xwalk <- d_meas_ctd |>
    select(variable = `_source_column`, measurement_type) |>
    filter(!is.na(variable), variable != "")
  dbWriteTable(con, "_flag_var_xwalk", d_var_xwalk, overwrite = TRUE)
  dbWriteTable(con, "_flags_accepted", d_flags, overwrite = TRUE)

  # the same (archive, _source_file, cast_key, depth_m) identity as ctd.file /
  # ctd.scan (see the section prose above), one row per physical scan
  dbExecute(con, "
    CREATE OR REPLACE TABLE _ctd_raw_keys AS
    SELECT DISTINCT ctd_cast_uuid, depth_m, cast_key, \"_source_file\",
           regexp_extract(\"_source_file\", '^([^/]+)', 1) || '.zip' AS archive
    FROM ctd_raw")
  stopifnot("(archive, _source_file, cast_key, depth_m) is not a unique scan key in ctd_raw" =
    dbGetQuery(con, "
      SELECT COUNT(*) - COUNT(DISTINCT archive || '|' || \"_source_file\" || '|' || cast_key || '|' || depth_m) AS n
      FROM _ctd_raw_keys")$n == 0)

  # scan-level flags: exact (archive, path, cast_id, depth) match
  dbExecute(con, "
    CREATE OR REPLACE TABLE _flags_scan_matched AS
    SELECT f.flag_id, f.qual_code, f.variable, k.ctd_cast_uuid, k.depth_m
    FROM _flags_accepted f
    JOIN _ctd_raw_keys k
      ON k.archive = f.archive AND k.\"_source_file\" = f.path
     AND k.cast_key = f.cast_id AND k.depth_m = f.depth
    WHERE f.depth IS NOT NULL")

  # file/cast-level flags (no scan_id in the ledger — every scan in the file)
  dbExecute(con, "
    CREATE OR REPLACE TABLE _flags_file_matched AS
    SELECT f.flag_id, f.qual_code, f.variable, k.ctd_cast_uuid, k.depth_m
    FROM _flags_accepted f
    JOIN _ctd_raw_keys k
      ON k.archive = f.archive AND k.\"_source_file\" = f.path
    WHERE f.depth IS NULL")

  n_scan_flags <- dbGetQuery(con, "SELECT COUNT(DISTINCT flag_id) AS n FROM _flags_accepted WHERE depth IS NOT NULL")$n
  n_file_flags <- n_flags - n_scan_flags
  n_scan_matched_flags <- dbGetQuery(con, "SELECT COUNT(DISTINCT flag_id) AS n FROM _flags_scan_matched")$n
  n_file_matched_flags <- dbGetQuery(con, "SELECT COUNT(DISTINCT flag_id) AS n FROM _flags_file_matched")$n
  n_unmatched <- (n_scan_flags - n_scan_matched_flags) + (n_file_flags - n_file_matched_flags)

  stopifnot(
    "accepted flag(s) in the snapshot find no scan in ctd_raw — a broken (archive, path, cast_id, depth) bridge" =
      n_unmatched == 0)

  # apply: named-variable rows go through the crosswalk; '*' applies to every
  # measurement_type recorded for that (ctd_cast_uuid, depth_m)
  dbExecute(con, "
    UPDATE ctd_measurement m
    SET measurement_qual = CAST(x.qual_code AS VARCHAR)
    FROM (
      SELECT fm.ctd_cast_uuid, fm.depth_m, w.measurement_type, fm.qual_code
      FROM (SELECT * FROM _flags_scan_matched UNION ALL SELECT * FROM _flags_file_matched) fm
      JOIN _flag_var_xwalk w ON w.variable = fm.variable
      WHERE fm.variable != '*'
      UNION ALL
      SELECT fm.ctd_cast_uuid, fm.depth_m, m2.measurement_type, fm.qual_code
      FROM (SELECT * FROM _flags_scan_matched UNION ALL SELECT * FROM _flags_file_matched) fm
      JOIN ctd_measurement m2
        ON m2.ctd_cast_uuid = fm.ctd_cast_uuid AND m2.depth_m = fm.depth_m
      WHERE fm.variable = '*'
    ) x
    WHERE m.ctd_cast_uuid = x.ctd_cast_uuid AND m.depth_m = x.depth_m
      AND m.measurement_type = x.measurement_type")

  n_scans_affected <- dbGetQuery(con, "
    SELECT COUNT(DISTINCT ctd_cast_uuid || '|' || depth_m) AS n
    FROM (SELECT * FROM _flags_scan_matched UNION ALL SELECT * FROM _flags_file_matched)")$n

  say(glue("accepted CTD QC flags applied: {format(n_flags, big.mark = ',')} flag(s) · ",
           "{format(n_scans_affected, big.mark = ',')} scan(s) affected · 0 unmatched"))

  n_flags_applied <- n_flags  # exported into metadata.json below, for release_database.qmd's n_flags_pending
  dbExecute(con, "DROP TABLE _flag_var_xwalk")
  dbExecute(con, "DROP TABLE _flags_accepted")
  dbExecute(con, "DROP TABLE _ctd_raw_keys")
  dbExecute(con, "DROP TABLE _flags_scan_matched")
  dbExecute(con, "DROP TABLE _flags_file_matched")
}

0 accepted flags in snapshot (last modified Wed, 19 Aug 2026 21:45:55 GMT)

Note

This chunk is shipping unrun this round (2026-09-03): the snapshot is a 600-byte header-only parquet with zero accepted rows (last modified 2026-08-19), so the CTD team has not yet accepted a flag through the ledger. It first has an effect at the next CTD ingest render after the team accepts one. RELEASES.md gets an entry only once n_flags_applied > 0 — a zero-effect chunk shipping is not itself a release-content change.

21 Apply Flag Overrides

The accepted-flag ledger above holds flags a CTD curator has entered against individual scans in the team’s PostgreSQL database. Some flags are known before anyone enters them, most often because the provider’s own cruise notes say a sensor failed. metadata/calcofi/ctd-cast/flag_overrides.csv is where those live: one committed row per (cruise, cast range, direction, series), each with the flag, the reason, the source it rests on and the question that keeps it provisional until the provider flags at source. It is applied after the provider’s flags and the accepted ledger and before any average, so an overridden sensor can never enter temperature_ave or anything built on it: the climatology, the transect shards, measurements.json, the derived dataset.

The first entry (2026-09-24, questions.csv Q37) is cruise 2026-07-3322 casts 1–6. The provider’s notes say the secondary temperature sensor was faulty there, and the data agree: sensor 2 sits a median 11 °C from sensor 1, up to 42 °C, against 0.002 °C on casts 7–72. Every series the provider builds from that temperature is flagged with it: conductivity-derived salinity, sigma-theta, potential temperature, and the sensor-2 oxygens, since oxygen concentration and saturation use the sensor’s own temperature and salinity. The corrected forms of those series are included so the rule still holds when a bottle-merged file for the cruise arrives.

Code

path_overrides <- file.path(dir_meta, "flag_overrides.csv")
d_overrides <- read_csv(path_overrides, show_col_types = FALSE, na = "",
                        col_types = cols(cast_min = "i", cast_max = "i", qual_code = "c",
                                         .default = "c"))
stopifnot(
  "flag_overrides.csv: every row needs a reason, a source and a question" =
    all(nzchar(d_overrides$reason) & nzchar(d_overrides$source) & nzchar(d_overrides$question)),
  "flag_overrides.csv: cast_dir is D, U or *" = all(d_overrides$cast_dir %in% c("D", "U", "*")),
  "flag_overrides.csv: measurement_type must be a registered CTD type" =
    all(d_overrides$measurement_type %in% d_meas_ctd$measurement_type))

dbWriteTable(con, "_flag_overrides", d_overrides, overwrite = TRUE)
dbExecute(con, "
  CREATE OR REPLACE TEMP TABLE _override_hits AS
  SELECT o.cruise_key, o.cast_min, o.cast_max, o.cast_dir, o.measurement_type, o.qual_code,
         c.cast_key, m.ctd_cast_uuid, m.depth_m
  FROM _flag_overrides o
  JOIN ctd_cast c
    ON c.cruise_key = o.cruise_key
   AND (o.cast_dir = '*' OR c.cast_dir = o.cast_dir)
   AND TRY_CAST(regexp_extract(c.cast_key, '_([0-9]{3})', 1) AS INTEGER)
         BETWEEN o.cast_min AND o.cast_max
  JOIN ctd_measurement m
    ON m.ctd_cast_uuid = c.ctd_cast_uuid AND m.measurement_type = o.measurement_type")

[1] 13235

Code

dbExecute(con, "
  UPDATE ctd_measurement m SET measurement_qual = h.qual_code
  FROM _override_hits h
  WHERE m.ctd_cast_uuid = h.ctd_cast_uuid AND m.depth_m = h.depth_m
    AND m.measurement_type = h.measurement_type")

[1] 13235

Code

d_override_report <- d_overrides |>
  left_join(dbGetQuery(con, "
    SELECT cruise_key, cast_min, cast_max, cast_dir, measurement_type,
           COUNT(*) AS n_flagged, COUNT(DISTINCT cast_key) AS n_casts
    FROM _override_hits GROUP BY ALL"),
    by = c("cruise_key", "cast_min", "cast_max", "cast_dir", "measurement_type")) |>
  mutate(across(c(n_flagged, n_casts), \(x) coalesce(as.integer(x), 0L)))
# a series the cruise's file does not carry yet (the corrected forms, before the bottle merge)
# matches nothing, which is fine; an override whose FIRST series matches nothing is a broken
# key, and so is an override row whose cruise has no casts at all in that range
stopifnot("a flag override matches no cast: check cruise_key / cast range / direction" =
  all(d_override_report |> group_by(cruise_key, cast_min, cast_max, cast_dir) |>
        summarize(n = sum(n_flagged), .groups = "drop") |> pull(n) > 0))
dbExecute(con, "DROP TABLE _flag_overrides")

[1] 0

Code

say(glue("flag overrides applied: {nrow(d_overrides)} rule(s), ",
         "{format(sum(d_override_report$n_flagged), big.mark = ',')} value(s) flagged"))

flag overrides applied: 24 rule(s), 13,235 value(s) flagged

Code

datatable(d_override_report |> select(cruise_key, cast_min, cast_max, cast_dir, measurement_type,
                                      qual_code, n_casts, n_flagged, reason, question),
          rownames = FALSE, options = list(pageLength = 12, scrollX = TRUE),
          caption = "Flag overrides (metadata/calcofi/ctd-cast/flag_overrides.csv) and what each flagged")

22 Combine the Sensor Pairs by Their Flags

The canonical CTD temperature, salinity and oxygen are the average of a sensor pair, and the source ships that average pre-computed. It is not trustworthy as shipped: where the second sensor is missing, SaltAve_Corr and TempAve are sometimes exactly half of sensor 1 (1998-07-32NM, 2024-01-33UD: salinity ≈ 17 at 1,000 m beside a sensor reading 34), TempAve reads 35–40 °C at 1–2 m on 2003-04-31JD beside a sensor reading 12, and none of these carries a flag. So every average is rebuilt here from its two sensors, by the rule Rasmus Swalethorp set for this project (2026-09-09), implemented and tested as calcofi4db::combine_sensor_pair():

a sensor flagged 8 (questionable) or 9 (bad) is dropped;
a 1 (use primary) or 2 (use secondary) selects that sensor alone;
otherwise the mean of the two, one alone when the other is missing, and nothing when neither survives.

The flags are each sensor’s own (T1Q/T2Q, S1Q/S2Q, Ox1Q/Ox2Q), which the corrected series inherit from their sensor (libs/build_ctd_measurement_registry.R § 4), after the CTD team’s accepted flags (§ above) have been written onto them. A sensor value removed above as a -99 sentinel or outside its physical bound is absent here, which the rule treats like a 9. A physical bound is only that guard: it never decides that a flagged value is good.

Every average is rebuilt, oxygen included, and an average is written wherever a surviving sensor exists even if the source shipped none (Ben, 2026-09-11: Rasmus is the authority on the rule). The source’s own average survives only on a scan where the source shipped no sensor value at all — there is nothing to rebuild it from.

Code

pair_rule <- tibble::tribble(
  ~ave,                          ~s1,                          ~s2,                          ~raw1,              ~raw2,
  "temperature_ave",             "temperature_1",              "temperature_2",              "temp1",            "temp2",
  "salinity_ave_corr",           "salinity_1_corr",            "salinity_2_corr",            "salt1_corr",       "salt2_corr",
  "oxygen_ml_l_ave_sta_corr",    "oxygen_ml_l_1_sta_corr",     "oxygen_ml_l_2_sta_corr",     "ox1_sta_corr",     "ox2_sta_corr",
  # the source ships no cruise-corrected average; built here (CalCOFI/workflows#106)
  "oxygen_ml_l_ave_cruise_corr", "oxygen_ml_l_1_cruise_corr",  "oxygen_ml_l_2_cruise_corr",  "ox1_cruise_corr",  "ox2_cruise_corr",
  "oxygen_umol_kg_ave_sta_corr", "oxygen_umol_kg_1_sta_corr",  "oxygen_umol_kg_2_sta_corr",  "ox1u_m_sta_corr",  "ox2u_m_sta_corr")

# a sensor series with no flag column would make the rule silently flag-blind
qual_cols <- setNames(d_meas_ctd$`_qual_column`, d_meas_ctd$measurement_type)
stopifnot("every paired sensor series must carry its sensor's flag column" =
            all(!is.na(qual_cols[c(pair_rule$s1, pair_rule$s2)]) &
                  qual_cols[c(pair_rule$s1, pair_rule$s2)] != ""))

combined_sql <- calcofi4db::combine_sensor_pair_sql("p.v1", "p.v2", "p.q1", "p.q2")
q89 <- function(q) glue("regexp_replace(CAST({q} AS VARCHAR), '\\.0+$', '') IN ('8', '9')")
# the same md5 key assign_deterministic_uuids_md5() uses, so a rebuilt row keeps its id
uuid_sql <- "substr(_h, 1, 8) || '-' || substr(_h, 9, 4) || '-' || substr(_h, 13, 4) || '-' ||
             substr(_h, 17, 4) || '-' || substr(_h, 21, 12)"

# the file layout per source file ([rename_fields]); a checkpoint from before
# 2026-09-23 has no such table, which the cross-check reports as 'unknown'
src_format_tbl <- if ("ctd_src_format" %in% dbListTables(con)) "ctd_src_format" else
  "(SELECT CAST(NULL AS VARCHAR) AS _source_file, CAST(NULL AS VARCHAR) AS src_format)"

# the provider's own average, beside ours: agreement per layout. The rule is ours
# (Rasmus, 2026-09-09) and it stands; this table is how a disagreement reaches the
# CTD team (Ben Gire, Kelsey Vogel) as a QC finding instead of being a silent choice.
pair_xcheck <- list()

pair_report <- pmap_dfr(pair_rule, function(ave, s1, s2, raw1, raw2) {
  dbExecute(con, glue("
    CREATE OR REPLACE TEMP TABLE _pair AS
    WITH m AS (
      SELECT ctd_cast_uuid, depth_m,
        any_value(cruise_key) AS cruise_key,
        max(measurement_value) FILTER (WHERE measurement_type = '{s1}')  AS v1,
        max(measurement_qual)  FILTER (WHERE measurement_type = '{s1}')  AS q1,
        max(measurement_value) FILTER (WHERE measurement_type = '{s2}')  AS v2,
        max(measurement_qual)  FILTER (WHERE measurement_type = '{s2}')  AS q2,
        max(measurement_value) FILTER (WHERE measurement_type = '{ave}') AS v_src,
        max(measurement_qual)  FILTER (WHERE measurement_type = '{ave}') AS q_src
      FROM ctd_measurement
      WHERE measurement_type IN ('{s1}', '{s2}', '{ave}')
      GROUP BY ctd_cast_uuid, depth_m),
    r AS (
      -- did the SOURCE ship any sensor value on this scan (before sentinels and bounds)?
      SELECT x.ctd_cast_uuid, x.depth_m,
        bool_or({raw1} IS NOT NULL OR {raw2} IS NOT NULL) AS src_has_sensor,
        any_value(coalesce(f.src_format, 'unknown'))     AS src_format
      FROM ctd_raw x LEFT JOIN {src_format_tbl} f ON f._source_file = x._source_file
      GROUP BY x.ctd_cast_uuid, x.depth_m)
    SELECT m.*, coalesce(r.src_has_sensor, FALSE) AS src_has_sensor,
           coalesce(r.src_format, 'unknown') AS src_format
    FROM m LEFT JOIN r USING (ctd_cast_uuid, depth_m)"))

  dbExecute(con, glue("
    CREATE OR REPLACE TEMP TABLE _pair_new AS
    SELECT p.*,
      CASE WHEN p.src_has_sensor THEN {combined_sql} ELSE p.v_src END AS v_new
    FROM _pair p"))

  r <- dbGetQuery(con, glue("
    SELECT
      COUNT(*) FILTER (WHERE v_src IS NOT NULL)                                   AS n_shipped,
      COUNT(*) FILTER (WHERE v_new IS NOT NULL)                                   AS n_written,
      COUNT(*) FILTER (WHERE v_src IS NOT NULL AND v_new IS NOT NULL
                         AND abs(v_src - v_new) > 1e-4)                           AS n_changed,
      COUNT(*) FILTER (WHERE v_src IS NULL AND v_new IS NOT NULL)                 AS n_filled,
      COUNT(*) FILTER (WHERE v_src IS NOT NULL AND v_new IS NULL)                 AS n_removed,
      COUNT(*) FILTER (WHERE ({q89('q1')} AND v1 IS NOT NULL)
                          OR ({q89('q2')} AND v2 IS NOT NULL))                    AS n_flag_excluded,
      COUNT(*) FILTER (WHERE NOT src_has_sensor AND v_src IS NOT NULL)            AS n_src_only
    FROM _pair_new"))

  pair_xcheck[[ave]] <<- dbGetQuery(con, glue("
    SELECT '{ave}' AS measure, src_format,
      COUNT(*) FILTER (WHERE v_src IS NOT NULL AND v_new IS NOT NULL)             AS n_both,
      COUNT(*) FILTER (WHERE v_src IS NOT NULL AND v_new IS NOT NULL
                         AND abs(v_src - v_new) > 1e-4)                           AS n_disagree,
      max(abs(v_src - v_new)) FILTER (WHERE v_src IS NOT NULL AND v_new IS NOT NULL) AS max_abs_diff,
      COUNT(*) FILTER (WHERE {q89('q_src')})                                      AS n_provider_flag_89,
      COUNT(*) FILTER (WHERE {q89('q_src')} AND v_new IS NOT NULL)                AS n_provider_89_kept_flagged,
      COUNT(*) FILTER (WHERE v_src IS NOT NULL AND v_new IS NULL)                 AS n_provider_only
    FROM _pair_new GROUP BY src_format ORDER BY src_format"))

  dbExecute(con, glue("DELETE FROM ctd_measurement WHERE measurement_type = '{ave}'"))
  dbExecute(con, glue("
    INSERT INTO ctd_measurement BY NAME
    SELECT {uuid_sql} AS ctd_measurement_uuid,
           ctd_cast_uuid, depth_m, '{ave}' AS measurement_type,
           v_new AS measurement_value,
           -- the provider's verdict on ITS average (SaltAve_CorrQ, OxAve_StaCorrQ, …)
           -- rides on ours: a value it flags 8/9 stays flagged for every consumer
           -- that applies qual_ok, even though the pair rule rebuilt it
           q_src AS measurement_qual, cruise_key
    FROM (
      SELECT *, md5(concat_ws('|', COALESCE(CAST(ctd_cast_uuid AS VARCHAR), ''),
                                   COALESCE(CAST(depth_m AS VARCHAR), ''),
                                   '{ave}')) AS _h
      FROM _pair_new WHERE v_new IS NOT NULL)"))

  mutate(r, measure = ave, .before = 1)
})
dbExecute(con, "DROP TABLE IF EXISTS _pair")

[1] 0

Code

dbExecute(con, "DROP TABLE IF EXISTS _pair_new")

[1] 0

Code

say(glue("sensor pairs combined by their flags: ",
         "{format(sum(pair_report$n_written), big.mark = ',')} averages written, ",
         "{format(sum(pair_report$n_changed), big.mark = ',')} changed from the source's, ",
         "{format(sum(pair_report$n_filled), big.mark = ',')} filled where it shipped none, ",
         "{format(sum(pair_report$n_flag_excluded), big.mark = ',')} with a flagged sensor left out"))

sensor pairs combined by their flags: 40,008,341 averages written, 2,745,036 changed from the source's, 13,387,804 filled where it shipped none, 36,836 with a flagged sensor left out

Code

datatable(pair_report, caption = "Averages rebuilt from their sensor pairs by the flag rule",
          rownames = FALSE)

Code

pair_xcheck <- bind_rows(pair_xcheck)
write_csv(pair_xcheck, glue("{dir_tmp}/ctd_pair_xcheck.csv"), na = "")
datatable(pair_xcheck, rownames = FALSE, caption = paste(
  "Cross-check: our flag-rule average against the provider's own (SaltAve_Corr +",
  "SaltAve_CorrQ, OxAve_StaCorr + OxAve_StaCorrQ, …) per file layout. n_disagree > 0 is",
  "a QC finding for the CTD team, never a silent choice; ours is what is written."))

Code

# the rebuilt averages must still sit inside their own bounds, and each scan has one
bounds_pairs <- check_measurement_bounds(
  con, "ctd_measurement", mt = d_meas_type, include_undeclared = FALSE)
stopifnot(
  "a rebuilt average fell outside its declared bound" =
    sum(bounds_pairs$status == "out_of_range" &
          bounds_pairs$measurement_type %in% pair_rule$ave) == 0,
  "one rebuilt average per scan" =
    dbGetQuery(con, glue("
      SELECT COUNT(*) - COUNT(DISTINCT ctd_cast_uuid || '|' || depth_m || '|' || measurement_type) AS n
      FROM ctd_measurement WHERE measurement_type IN ({paste0(\"'\", pair_rule$ave, \"'\", collapse = ', ')})"))$n == 0)

# and the rule itself, asserted on what was written: where a sensor is flagged 8/9 the
# average is the OTHER sensor (or absent) — never a mean that includes the flagged one
n_flag_leak <- sum(pmap_int(pair_rule, function(ave, s1, s2, raw1, raw2) {
  dbGetQuery(con, glue("
    WITH p AS (
      SELECT ctd_cast_uuid, depth_m,
        max(measurement_value) FILTER (WHERE measurement_type = '{s1}')  AS v1,
        max(measurement_qual)  FILTER (WHERE measurement_type = '{s1}')  AS q1,
        max(measurement_value) FILTER (WHERE measurement_type = '{s2}')  AS v2,
        max(measurement_qual)  FILTER (WHERE measurement_type = '{s2}')  AS q2,
        max(measurement_value) FILTER (WHERE measurement_type = '{ave}') AS va
      FROM ctd_measurement WHERE measurement_type IN ('{s1}', '{s2}', '{ave}')
      GROUP BY ctd_cast_uuid, depth_m)
    SELECT COUNT(*) AS n FROM p
    WHERE va IS NOT NULL AND (
      ({q89('q1')} AND v1 IS NOT NULL AND (v2 IS NULL OR {q89('q2')} OR abs(va - v2) > 1e-9)) OR
      ({q89('q2')} AND v2 IS NOT NULL AND (v1 IS NULL OR {q89('q1')} OR abs(va - v1) > 1e-9)))"))$n
}))
stopifnot("an average includes a sensor its provider flags 8/9" = n_flag_leak == 0)

# and the overrides in particular: on every overridden cast the average is the OTHER sensor
# exactly — for Q37, temperature_ave on 2026-07-3322 casts 1-6 equals temperature_1
n_override_leak <- dbGetQuery(con, glue("
  WITH t AS (
    SELECT h.ctd_cast_uuid, h.depth_m
    FROM _override_hits h WHERE h.measurement_type = 'temperature_2'),
  p AS (
    SELECT m.ctd_cast_uuid, m.depth_m,
      max(measurement_value) FILTER (WHERE measurement_type = 'temperature_1')   AS v1,
      max(measurement_qual)  FILTER (WHERE measurement_type = 'temperature_1')   AS q1,
      max(measurement_value) FILTER (WHERE measurement_type = 'temperature_ave') AS va
    FROM ctd_measurement m JOIN t USING (ctd_cast_uuid, depth_m)
    WHERE measurement_type IN ('temperature_1', 'temperature_ave')
    GROUP BY 1, 2)
  SELECT
    COUNT(*) FILTER (WHERE va IS NOT NULL AND (v1 IS NULL OR {q89('q1')} OR abs(va - v1) > 1e-12)) AS n_leak,
    COUNT(*) FILTER (WHERE va IS NOT NULL) AS n_ave,
    max(abs(va - v1)) AS max_abs_diff
  FROM p"))
say(glue("override check: {n_override_leak$n_ave} temperature_ave values on overridden casts, ",
         "max |temperature_ave - temperature_1| = {n_override_leak$max_abs_diff}"))

override check: 1775 temperature_ave values on overridden casts, max |temperature_ave - temperature_1| = 0

Code

stopifnot("an average on an overridden cast includes the overridden sensor" =
            n_override_leak$n_leak == 0)
dbExecute(con, "DROP TABLE IF EXISTS _override_hits")

[1] 0

22.1 ctd_summary — summary stats across cast directions

Code

# one measurement type at a time: every group is (cruise, site, depth, TYPE), so the
# per-type slices are disjoint and their union is exactly the one-statement result. The
# single GROUP BY over ~280M rows needed a hash table larger than the DuckDB budget
# (6 GB, 2026-09-24: "Out of Memory Error ... 5.5 GiB/5.5 GiB used"); it had only ever
# passed because the post-checkpoint connection ran unlimited.
summ_sql <- function(type) glue("
  SELECT
    c.cruise_key,
    c.site_key,
    m.depth_m,
    m.measurement_type,
    AVG(m.measurement_value)    AS avg,
    CASE
      WHEN COUNT(*) = 1 THEN 0
      ELSE COALESCE(STDDEV_SAMP(m.measurement_value), 0)
    END                          AS stddev,
    COUNT(*)                     AS n_obs
  FROM ctd_measurement m
  INNER JOIN ctd_cast c ON m.ctd_cast_uuid = c.ctd_cast_uuid
  WHERE m.measurement_type = '{type}'
    AND NOT isnan(m.measurement_value)
    AND isfinite(m.measurement_value)
  GROUP BY c.cruise_key, c.site_key, m.depth_m, m.measurement_type")
summ_types <- dbGetQuery(con, "SELECT DISTINCT measurement_type FROM ctd_measurement ORDER BY 1")$measurement_type
dbExecute(con, glue("CREATE OR REPLACE TABLE ctd_summary AS {summ_sql(summ_types[1])}"))

[1] 3626441

Code

for (type in summ_types[-1])
  dbExecute(con, glue("INSERT INTO ctd_summary {summ_sql(type)}"))

n_summ <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM ctd_summary")$n
say(glue("ctd_summary: {format(n_summ, big.mark = ',')} rows"))

ctd_summary: 140,809,234 rows

Code

assign_deterministic_uuids_md5(
  con = con,
  table_name = "ctd_summary",
  id_col = "ctd_summary_uuid",
  key_cols = c("cruise_key", "site_key", "depth_m", "measurement_type")
)

22.2 measurement_type — reference table

Code

dbWriteTable(con, "measurement_type", d_meas_type, overwrite = TRUE)
say(glue("measurement_type: {nrow(d_meas_type)} rows"))

measurement_type: 202 rows

Code

# validate: all measurement_types used by this dataset are registered
ctd_types_used <- dbGetQuery(
  con,
  "SELECT DISTINCT measurement_type FROM ctd_measurement"
)$measurement_type

registered <- d_meas_type |>
  filter(str_detect(`_source_datasets`, "calcofi_ctd-cast")) |>
  pull(measurement_type)

unregistered <- setdiff(ctd_types_used, registered)
stopifnot(
  "all measurement_types used by calcofi_ctd-cast must be registered in _source_datasets" = length(
    unregistered
  ) ==
    0
)

22.3 Measurement provenance — instrument and accuracy eras

Which instrument produced a value, to what stated accuracy, and over which years. Recovered from the CalCOFI hydrographic master’s 0-Measurements table into metadata/calcofi/hydro-master/measurement_method.csv; the pipeline had no equivalent. It matters most where a method changed mid-record — water temperature switches from reversing thermometer (through 1993-04-15) to CTD thermistor (from 1993-08-11), which is the boundary this dataset begins at, and phosphate crosses four instruments from 1949.

The registry is one-to-many by design (Temperature has 6 eras), so it is joined by property rather than by exact type: CTD types carry sensor and correction suffixes (temperature_1, salinity_ave_corr, oxygen_ml_l_ave_sta_corr) that the era table knows nothing about, so those are stripped back to the base property before matching.

Code

path_method <- here("metadata/calcofi/hydro-master/measurement_method.csv")

d_method <- read_csv(path_method, show_col_types = F) |>
  filter(!is.na(measurement_type)) |>
  select(property = measurement_type, method, accuracy, year_started, year_ended)

# strip sensor/correction suffixes and the bottle marker to reach the base property:
#   temperature_ave  -> temperature ; oxygen_ml_l_ave_sta_corr -> oxygen_ml_l
#   btl_temperature  -> temperature ; salinity_btl            -> salinity
#   oxygen_btl_ml_l  -> oxygen_ml_l  (the marker is INFIX here, not a prefix/suffix —
#                                     the three positions are why this is not one regex)
base_property <- function(x) {
  x |>
    str_remove("^btl_") |>
    str_remove("_btl$") |>
    str_replace("_btl_", "_") |>
    str_remove("_(cruise|sta)_corr$") |>
    str_remove("_corr$") |>
    str_remove("_ave$") |>
    str_remove("_[12]$")
}

d_prov <- d_meas_ctd |>
  transmute(measurement_type, units, is_canonical,
            property = base_property(measurement_type)) |>
  inner_join(d_method, by = "property", relationship = "many-to-many") |>
  arrange(desc(is_canonical), measurement_type, year_started)

n_typ <- n_distinct(d_prov$measurement_type)
say(glue("CTD types with instrument provenance: {n_typ} of {nrow(d_meas_ctd)} ",
         "({n_distinct(d_prov$property)} distinct properties, ",
         "{nrow(d_prov)} type x era rows)"))

CTD types with instrument provenance: 25 of 55 (9 distinct properties, 38 type x era rows)

Code

say(glue("  unmatched (no era recorded): ",
         "{paste(setdiff(d_meas_ctd$measurement_type, d_prov$measurement_type) |> head(8), collapse = ', ')}",
         if (nrow(d_meas_ctd) - n_typ > 8) " ..." else ""))

  unmatched (no era recorded): beam_attenuation, btl_ammonium, btl_depth, dynamic_height, est_chlorophyll_a_cruise_corr, est_chlorophyll_a_sta_corr, est_nitrate_cruise_corr, est_nitrate_sta_corr ...

Code

d_prov |>
  transmute(measurement_type, canonical = is_canonical, units,
            method, accuracy, from = year_started, to = year_ended) |>
  datatable(
    caption  = paste("CTD measurement provenance — instrument, stated accuracy and",
                     "era, from the CalCOFI hydrographic master's 0-Measurements"),
    rownames = FALSE,
    options  = list(pageLength = 10, scrollX = TRUE, dom = "ftip"))

Documentation only, for now

This table is rendered, not written into metadata.json. Promoting it into the sidecar so consumers get instrument provenance alongside units is a follow-up that belongs in calcofi4db::build_metadata_json(), and 18 of the 35 era rows still need a data manager to confirm their mapping to canonical types (hydro_master_06).

22.4 Drop intermediate table

Code

dbExecute(con, "DROP TABLE IF EXISTS ctd_raw")

[1] 0

23 Measurement Frequency

ctd_measurement is by far the largest table in this dataset — in the prior release the calcofi_ctd-cast parquet totaled 22.16 GB, with ctd_measurement alone at 15.8 GB (233 M rows). Three stacked multipliers drive that size. The queries below quantify each, motivating the adaptively-thinned ctd_thin table built in the next section.

Code

freq_tbl_sizes <- dbGetQuery(con, "
  SELECT 'ctd_cast'         AS table_name, COUNT(*) AS n_rows FROM ctd_cast
  UNION ALL
  SELECT 'ctd_measurement' AS table_name, COUNT(*) AS n_rows FROM ctd_measurement
  UNION ALL
  SELECT 'ctd_summary'     AS table_name, COUNT(*) AS n_rows FROM ctd_summary")

n_meas_rows <- freq_tbl_sizes$n_rows[freq_tbl_sizes$table_name == "ctd_measurement"]
cat(glue("ctd_measurement: {format(n_meas_rows, big.mark = ',')} rows"), "\n")

ctd_measurement: 284,379,888 rows

Code

freq_tbl_sizes |>
  dt(caption = "Normalized CTD table row counts (this run)",
     fname   = "ctd_freq_table_sizes")

23.1 Multiplier 1 — sub-meter depth resolution

CTD profiles record a sample roughly every metre from surface to seafloor (~3,600 m max). Thinning to a 10 m grid is therefore an order-of-magnitude reduction on the depth axis.

Code

# depth interval between consecutive samples within a cast, using `pressure`
# as a representative measurement_type (all types share a cast's depth grid)
freq_depth_int <- dbGetQuery(con, "
  WITH gaps AS (
    SELECT depth_m - LAG(depth_m) OVER (
             PARTITION BY ctd_cast_uuid ORDER BY depth_m) AS d_gap
    FROM ctd_measurement
    WHERE measurement_type = 'pressure')
  SELECT
    ROUND(MIN(d_gap), 2)                 AS min_m,
    ROUND(QUANTILE_CONT(d_gap, 0.25), 2) AS q25_m,
    ROUND(MEDIAN(d_gap), 2)              AS median_m,
    ROUND(QUANTILE_CONT(d_gap, 0.75), 2) AS q75_m,
    ROUND(MAX(d_gap), 2)                 AS max_m
  FROM gaps
  WHERE d_gap IS NOT NULL")

cat(glue(
  "median within-cast depth interval: {freq_depth_int$median_m} m ",
  "(~{round(10 / freq_depth_int$median_m)}x reduction at a 10 m grid)"), "\n")

median within-cast depth interval: 1 m (~10x reduction at a 10 m grid)

Code

freq_depth_int |>
  dt(caption = "Within-cast depth interval (m), representative type `pressure`",
     fname   = "ctd_freq_depth_resolution")

23.2 Multiplier 2 — up- and down-casts both stored

Nearly every physical cast is recorded as a separate down-cast and up-cast. Keeping a single direction roughly halves the rows.

Code

# physical cast = (cruise_key, site_key, cast_base); direction resolved from
# cast_dir, falling back to the cast_key suffix (trailing d/u)
freq_cast_dir <- dbGetQuery(con, "
  WITH phys AS (
    SELECT DISTINCT
      cruise_key, site_key,
      regexp_replace(cast_key, '[duDU]$', '') AS cast_base,
      COALESCE(
        cast_dir,
        CASE
          WHEN RIGHT(cast_key, 1) IN ('d', 'D') THEN 'D'
          WHEN RIGHT(cast_key, 1) IN ('u', 'U') THEN 'U'
        END) AS dir
    FROM ctd_cast)
  SELECT
    COUNT(DISTINCT cruise_key || '|' || site_key || '|' || cast_base)
                                       AS n_physical_casts,
    COUNT(*) FILTER (WHERE dir = 'D')   AS n_with_downcast,
    COUNT(*) FILTER (WHERE dir = 'U')   AS n_with_upcast,
    COUNT(*) FILTER (WHERE dir IS NULL) AS n_unknown_dir
  FROM phys")

cat(glue(
  "{format(freq_cast_dir$n_physical_casts, big.mark = ',')} physical casts; ",
  "{format(freq_cast_dir$n_with_downcast, big.mark = ',')} have a downcast, ",
  "{format(freq_cast_dir$n_with_upcast, big.mark = ',')} an upcast"), "\n")

9,675 physical casts; 9,668 have a downcast, 9,671 an upcast

Code

freq_cast_dir |>
  dt(caption = "Cast direction coverage per physical cast",
     fname   = "ctd_freq_cast_direction")

23.3 Multiplier 3 — redundant measurement-type variants

Many of the measurement types below are variants of the same property (multiple oxygen units and sensor/correction combinations, dual temperature and salinity sensors, derived estimates). ctd_thin keeps the canonical types — the averages, the bottle references, the uncorrected sensors (2026-08-06, so a sensor-only cruise still has salinity and oxygen) and, since 2026-09-10, the bottle-corrected per-sensor series and the chlorophyll / nitrate sensor estimates at both correction levels (Rasmus Swalethorp’s transect-plotter rule) — see the is_canonical flag in metadata/measurement_type.csv, set by libs/build_ctd_measurement_registry.R.

Code

freq_meas_types <- dbGetQuery(con, "
  SELECT
    measurement_type,
    COUNT(*)                      AS n_rows,
    COUNT(DISTINCT ctd_cast_uuid) AS n_casts
  FROM ctd_measurement
  GROUP BY measurement_type
  ORDER BY n_rows DESC")

cat(glue(
  "{nrow(freq_meas_types)} distinct measurement_type values in ctd_measurement"),
  "\n")

55 distinct measurement_type values in ctd_measurement

Code

freq_meas_types |>
  dt(caption = "Rows per measurement_type",
     fname   = "ctd_freq_measurement_types")

24 CTD Thin

ctd_thin is an adaptively-thinned ctd_measurement: one cast direction per physical cast, canonical measurement types only, and a ~10 m depth grid with thermocline / halocline inflections preserved via Ramer–Douglas–Peucker line simplification. It is the headline CTD table; full ctd_measurement is retained as a supplemental output. See the Measurement Frequency section above for the size rationale.

24.1 Verify cast pairing

Code

# physical cast = (cruise_key, site_key, cast_base); direction comes from the
# cast_key suffix (trailing d/u), which is part of the deterministic
# ctd_cast_uuid key and fully consistent per cast. confirm casts are paired
# before thinning by direction.
ctd_pairing <- dbGetQuery(con, "
  WITH phys AS (
    SELECT DISTINCT
      cruise_key, site_key,
      regexp_replace(cast_key, '[duDU]$', '') AS cast_base,
      CASE WHEN RIGHT(cast_key, 1) IN ('d', 'D') THEN 'D'
           WHEN RIGHT(cast_key, 1) IN ('u', 'U') THEN 'U' END AS dir
    FROM ctd_cast)
  SELECT n_dirs, COUNT(*) AS n_physical_casts
  FROM (
    SELECT cruise_key, site_key, cast_base, COUNT(DISTINCT dir) AS n_dirs
    FROM phys GROUP BY 1, 2, 3)
  GROUP BY n_dirs ORDER BY n_dirs")

n_paired <- sum(ctd_pairing$n_physical_casts[ctd_pairing$n_dirs == 2])
n_single <- sum(ctd_pairing$n_physical_casts[ctd_pairing$n_dirs == 1])
cat(glue(
  "{format(n_paired, big.mark = ',')} physical casts have both directions, ",
  "{format(n_single, big.mark = ',')} have one"), "\n")

9,657 physical casts have both directions, 18 have one

Code

# sanity: most casts should be paired — a near-zero n_paired means the pairing
# key is wrong and the single-direction lever would silently keep everything
stopifnot(
  "ctd_thin: cast pairing key looks wrong (no two-direction casts)" =
    n_paired > n_single)

ctd_pairing |>
  dt(caption = "Cast directions per physical cast", fname = "ctd_thin_pairing")

24.2 Choose a single cast direction

Code

# canonical measurement types for ctd_thin (one per property; see is_canonical
# in metadata/measurement_type.csv)
canon_types <- d_meas_type |>
  filter(str_detect(`_source_datasets`, "calcofi_ctd-cast"), is_canonical) |>
  pull(measurement_type)
canon_in <- paste0("'", canon_types, "'", collapse = ", ")
cat(glue("{length(canon_types)} canonical measurement types: ",
         "{paste(canon_types, collapse = ', ')}"), "\n")

44 canonical measurement types: beam_attenuation, btl_ammonium, btl_chlorophyll_a, btl_depth, btl_nitrate, btl_nitrite, btl_phaeopigment, btl_phosphate, btl_silicate, btl_temperature, dynamic_height, est_chlorophyll_a_cruise_corr, est_chlorophyll_a_sta_corr, est_nitrate_cruise_corr, est_nitrate_sta_corr, fluorescence_v, isus_v, oxygen_btl_ml_l, oxygen_btl_umol_kg, oxygen_ml_l_1, oxygen_ml_l_1_cruise_corr, oxygen_ml_l_1_sta_corr, oxygen_ml_l_2, oxygen_ml_l_2_cruise_corr, oxygen_ml_l_2_sta_corr, oxygen_ml_l_ave_cruise_corr, oxygen_ml_l_ave_sta_corr, oxygen_umol_kg_1, oxygen_umol_kg_2, oxygen_umol_kg_ave_sta_corr, par, ph, pressure, salinity_1, salinity_1_corr, salinity_2, salinity_2_corr, salinity_ave_corr, salinity_btl, sigma_theta_1, spar, specific_volume_anomaly, temperature_ave, transmissometer

Code

# choose one direction per physical cast: prefer downcast, fall back to upcast.
# _chosen_dir maps every chosen-direction ctd_cast_uuid to its physical cast
# (phys_cast_id) and direction.
dbExecute(con, "
  CREATE OR REPLACE TEMP TABLE _chosen_dir AS
  WITH resolved AS (
    SELECT DISTINCT
      ctd_cast_uuid, cruise_key, site_key, cast_key,
      regexp_replace(cast_key, '[duDU]$', '') AS cast_base,
      CASE WHEN RIGHT(cast_key, 1) IN ('d', 'D') THEN 'D'
           WHEN RIGHT(cast_key, 1) IN ('u', 'U') THEN 'U' END AS dir
    FROM ctd_cast),
  pick AS (
    SELECT cruise_key, site_key, cast_base, cast_key, dir
    FROM (SELECT DISTINCT cruise_key, site_key, cast_base, cast_key, dir
          FROM resolved)
    QUALIFY ROW_NUMBER() OVER (
      PARTITION BY cruise_key, site_key, cast_base
      ORDER BY CASE dir WHEN 'D' THEN 1 WHEN 'U' THEN 2 ELSE 3 END,
               cast_key) = 1)
  SELECT
    r.ctd_cast_uuid,
    r.cruise_key,
    r.cruise_key || '|' || r.site_key || '|' || r.cast_base AS phys_cast_id,
    p.dir AS cast_dir
  FROM resolved r
  JOIN pick p USING (cruise_key, site_key, cast_base, cast_key)")

[1] 3755547

Code

n_chosen <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM _chosen_dir")$n
cat(glue("_chosen_dir: {format(n_chosen, big.mark = ',')} cast records ",
         "(single direction per physical cast)"), "\n")

_chosen_dir: 3,755,547 cast records (single direction per physical cast)

24.3 10 m depth grid

Code

# depth backbone: every distinct depth a canonical measurement is recorded at,
# per chosen-direction physical cast
dbExecute(con, glue("
  CREATE OR REPLACE TEMP TABLE _profile_depths AS
  SELECT DISTINCT cd.phys_cast_id, m.depth_m
  FROM ctd_measurement m
  JOIN _chosen_dir cd ON m.ctd_cast_uuid = cd.ctd_cast_uuid
  WHERE m.measurement_type IN ({canon_in})"))

[1] 4155030

Code

# 10 m grid: the sample nearest each 10 m node per profile, plus the shallowest
# and deepest sample (RDP needs fixed profile endpoints)
dbExecute(con, "
  CREATE OR REPLACE TEMP TABLE _grid_depths AS
  WITH ranked AS (
    SELECT phys_cast_id, depth_m,
      ROW_NUMBER() OVER (
        PARTITION BY phys_cast_id, ROUND(depth_m / 10.0) * 10.0
        ORDER BY ABS(depth_m - ROUND(depth_m / 10.0) * 10.0), depth_m) AS rn
    FROM _profile_depths),
  endpoints AS (
    SELECT phys_cast_id, MIN(depth_m) AS depth_m FROM _profile_depths GROUP BY 1
    UNION
    SELECT phys_cast_id, MAX(depth_m) AS depth_m FROM _profile_depths GROUP BY 1)
  SELECT phys_cast_id, depth_m FROM ranked WHERE rn = 1
  UNION
  SELECT phys_cast_id, depth_m FROM endpoints")

[1] 430926

Code

n_prof <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM _profile_depths")$n
n_grid <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM _grid_depths")$n
cat(glue("_grid_depths: {format(n_grid, big.mark = ',')} depths on the 10 m grid ",
         "(from {format(n_prof, big.mark = ',')} profile depths)"), "\n")

_grid_depths: 430,926 depths on the 10 m grid (from 4,155,030 profile depths)

24.4 Inflection points (Ramer–Douglas–Peucker)

Code

# Ramer-Douglas-Peucker line simplification: keep depths where a profile deviates
# from the local straight line by more than `eps` (captures thermocline /
# halocline inflections). pure DP — the ~10 m backbone is provided separately by
# _grid_depths, so RDP only adds genuine inflections on top of it.
rdp_keep <- function(x, y, eps) {
  n    <- length(x)
  keep <- rep(FALSE, n)
  if (n == 0L) return(keep)
  keep[c(1L, n)] <- TRUE
  if (n <= 2L) return(keep)
  stack <- list(c(1L, n))
  while (length(stack) > 0L) {
    seg   <- stack[[length(stack)]]
    stack[[length(stack)]] <- NULL
    i <- seg[1L]; j <- seg[2L]
    if (j - i < 2L) next
    idx <- (i + 1L):(j - 1L)
    dx  <- x[j] - x[i]
    dy  <- y[j] - y[i]
    den <- sqrt(dx * dx + dy * dy)
    d <- if (den == 0) {
      abs(y[idx] - y[i])
    } else {
      abs(dy * x[idx] - dx * y[idx] + x[j] * y[i] - y[j] * x[i]) / den
    }
    if (max(d) > eps) {
      k <- idx[which.max(d)]
      keep[k] <- TRUE
      stack[[length(stack) + 1L]] <- c(i, k)
      stack[[length(stack) + 1L]] <- c(k, j)
    }
  }
  keep
}

# per-variable tolerance in measurement units — the primary tuning knob, set so
# the median retained-sample gap lands near the 10 m grid spacing (a looser eps
# keeps fewer inflections; revisit with the spot-check in the verification section)
rdp_eps <- c(temperature_ave = 0.2, salinity_ave_corr = 0.04)

# run RDP per cruise to keep memory flat: pull the two key variables for
# chosen-direction casts, simplify each profile, collect retained depths
cruise_keys <- dbGetQuery(con, "SELECT DISTINCT cruise_key FROM _chosen_dir")$cruise_key

rdp_retained <- purrr::map(cruise_keys, function(ck) {
  d_prof <- dbGetQuery(con, glue("
    SELECT cd.phys_cast_id, m.measurement_type, m.depth_m,
           AVG(m.measurement_value) AS measurement_value
    FROM ctd_measurement m
    JOIN _chosen_dir cd
      ON m.ctd_cast_uuid = cd.ctd_cast_uuid AND cd.cruise_key = '{ck}'
    WHERE m.cruise_key = '{ck}'
      AND m.measurement_type IN ('temperature_ave', 'salinity_ave_corr')
      AND m.measurement_value IS NOT NULL
    GROUP BY cd.phys_cast_id, m.measurement_type, m.depth_m"))
  if (nrow(d_prof) == 0) return(NULL)
  d_prof |>
    arrange(phys_cast_id, measurement_type, depth_m) |>
    group_by(phys_cast_id, measurement_type) |>
    filter(rdp_keep(depth_m, measurement_value,
                    eps = rdp_eps[[measurement_type[1]]])) |>
    ungroup() |>
    distinct(phys_cast_id, depth_m)
}) |>
  purrr::list_rbind() |>
  distinct(phys_cast_id, depth_m)

dbWriteTable(con, "_rdp_retained", rdp_retained,
             temporary = TRUE, overwrite = TRUE)
cat(glue("_rdp_retained: {format(nrow(rdp_retained), big.mark = ',')} ",
         "(profile, depth) pairs flagged as inflections"), "\n")

_rdp_retained: 203,658 (profile, depth) pairs flagged as inflections

24.5 Assemble ctd_thin

Code

# Union three retention reasons, in precedence order: 10 m grid, RDP inflections,
# then bottle-trip depths.
#
# The third is why this section changed. Thinning is designed for DENSE SENSOR
# SCANS: the 10 m backbone and the RDP inflections are both derived from
# temperature/salinity profiles, so a depth is kept because the *sensor* profile
# bends there. Bottle values are not a profile — they are a handful of discrete,
# lab-analysed samples per cast at depths chosen by the watch, and nothing about
# the sensor profile knows where they are. Selecting depths by sensor geometry and
# then subsetting bottle values to those depths silently discards most of them.
#
# Measured against release v2026.07.17: of 133,206 distinct (cast, depth) bottle
# ammonium measurements, only 35,549 (26.7%) sit at a depth this logic retained.
# Roughly half the loss is the cast-direction choice (by design); the rest was
# depth thinning throwing away real lab samples. So flagging btl_ammonium as
# canonical is necessary but NOT sufficient — without this clause it would have
# landed about a quarter of the ammonium record in `obs` and looked like success.
#
# A bottle depth costs nothing to keep (a few rows per cast) and cannot be
# reconstructed by interpolation, so retain every one carrying a canonical value.
canon_btl <- d_meas_ctd |>
  filter(is_canonical, str_starts(measurement_type, "btl_") |
           measurement_type %in% c("salinity_btl", "oxygen_btl_ml_l",
                                   "oxygen_btl_umol_kg")) |>
  pull(measurement_type)
say(glue("bottle-grain canonical types retained whole: ",
         "{if (length(canon_btl)) paste(canon_btl, collapse = ', ') else '(none)'}"))

bottle-grain canonical types retained whole: btl_ammonium, btl_chlorophyll_a, btl_depth, btl_nitrate, btl_nitrite, btl_phaeopigment, btl_phosphate, btl_silicate, btl_temperature, oxygen_btl_ml_l, oxygen_btl_umol_kg, salinity_btl

Code

btl_clause <- if (length(canon_btl)) glue("
  UNION
  SELECT cd.phys_cast_id, m.depth_m, 'bottle' AS retained_reason
  FROM ctd_measurement m
  JOIN _chosen_dir cd ON m.ctd_cast_uuid = cd.ctd_cast_uuid
  WHERE m.measurement_type IN ({paste(sprintf(\"'%s'\", canon_btl), collapse = ', ')})
    AND NOT EXISTS (
      SELECT 1 FROM _grid_depths g
      WHERE g.phys_cast_id = cd.phys_cast_id AND g.depth_m = m.depth_m)
    AND NOT EXISTS (
      SELECT 1 FROM _rdp_retained r
      WHERE r.phys_cast_id = cd.phys_cast_id AND r.depth_m = m.depth_m)") else ""

dbExecute(con, glue("
  CREATE OR REPLACE TEMP TABLE _retained_depths AS
  SELECT phys_cast_id, depth_m, 'grid' AS retained_reason FROM _grid_depths
  UNION
  SELECT r.phys_cast_id, r.depth_m, 'inflection' AS retained_reason
  FROM _rdp_retained r
  WHERE NOT EXISTS (
    SELECT 1 FROM _grid_depths g
    WHERE g.phys_cast_id = r.phys_cast_id AND g.depth_m = r.depth_m)
  {btl_clause}"))

[1] 683896

Code

# ctd_thin: ctd_measurement rows that are a canonical type, on the chosen
# direction, at a retained depth — a pure row subset (values never interpolated).
# the QUALIFY de-dups defensively on (ctd_cast_uuid, depth_m, measurement_type):
# the dedup_ctd_raw chunk above already makes ctd_measurement unique, so this is
# a no-op safeguard that keeps ctd_thin's PK guaranteed-unique on its own.
dbExecute(con, glue("
  CREATE OR REPLACE TABLE ctd_thin AS
  WITH candidates AS (
    SELECT
      m.ctd_cast_uuid,
      m.depth_m,
      m.measurement_type,
      m.measurement_value,
      m.measurement_qual,
      cd.cast_dir,
      rd.retained_reason,
      m.cruise_key
    FROM ctd_measurement m
    JOIN _chosen_dir cd      ON m.ctd_cast_uuid = cd.ctd_cast_uuid
    JOIN _retained_depths rd ON rd.phys_cast_id = cd.phys_cast_id
                            AND rd.depth_m      = m.depth_m
    WHERE m.measurement_type IN ({canon_in}))
  SELECT *
  FROM candidates
  QUALIFY ROW_NUMBER() OVER (
    PARTITION BY ctd_cast_uuid, depth_m, measurement_type
    ORDER BY measurement_qual NULLS LAST, measurement_value) = 1
  ORDER BY cruise_key, ctd_cast_uuid, measurement_type, depth_m"))

[1] 18843115

Code

n_thin <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM ctd_thin")$n
n_meas <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM ctd_measurement")$n
cat(glue(
  "ctd_thin: {format(n_thin, big.mark = ',')} rows ",
  "({round(100 * n_thin / n_meas, 1)}% of ctd_measurement's ",
  "{format(n_meas, big.mark = ',')})"), "\n")

ctd_thin: 18,843,115 rows (6.6% of ctd_measurement's 284,379,888)

Code

# deterministic PK from the natural key, mirroring ctd_measurement
assign_deterministic_uuids_md5(
  con        = con,
  table_name = "ctd_thin",
  id_col     = "ctd_thin_uuid",
  key_cols   = c("ctd_cast_uuid", "depth_m", "measurement_type"))

# drop scratch tables; _chosen_dir is kept for the verification chunk that
# follows and dropped there (downstream metadata helpers enumerate all DB tables)
dbExecute(con, "DROP TABLE IF EXISTS _profile_depths")

[1] 0

Code

dbExecute(con, "DROP TABLE IF EXISTS _grid_depths")

[1] 0

Code

dbExecute(con, "DROP TABLE IF EXISTS _rdp_retained")

[1] 0

Code

dbExecute(con, "DROP TABLE IF EXISTS _retained_depths")

[1] 0

24.6 Verify ctd_thin

Code

# 1. row count + thinning ratio
n_thin <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM ctd_thin")$n
n_meas <- dbGetQuery(con, "SELECT COUNT(*) AS n FROM ctd_measurement")$n
cat(glue("ctd_thin: {format(n_thin, big.mark = ',')} rows = ",
         "{round(100 * n_thin / n_meas, 1)}% of ctd_measurement"), "\n")

ctd_thin: 18,843,115 rows = 6.6% of ctd_measurement

Code

stopifnot(
  "ctd_thin is empty"                          = n_thin > 0,
  "ctd_thin not smaller than ctd_measurement"  = n_thin < n_meas)

# 2. measurement_type is a subset of the canonical list
thin_types <- dbGetQuery(con,
  "SELECT DISTINCT measurement_type FROM ctd_thin")$measurement_type
stopifnot("ctd_thin has non-canonical measurement_type" =
  length(setdiff(thin_types, canon_types)) == 0)

# 3. ctd_thin_uuid is unique
stopifnot("ctd_thin_uuid not unique" = dbGetQuery(con,
  "SELECT COUNT(*) - COUNT(DISTINCT ctd_thin_uuid) AS n FROM ctd_thin")$n == 0)

# 4. faithful subset of ctd_measurement — every ctd_thin row corresponds to a
#    real ctd_measurement row with the same value (thinning selects rows, never
#    interpolates). a NOT EXISTS check, since ctd_measurement has duplicate keys.
stopifnot("ctd_thin rows not a faithful subset of ctd_measurement" =
  dbGetQuery(con, "
    SELECT COUNT(*) AS n
    FROM ctd_thin t
    WHERE NOT EXISTS (
      SELECT 1 FROM ctd_measurement m
      WHERE m.ctd_cast_uuid     = t.ctd_cast_uuid
        AND m.depth_m           = t.depth_m
        AND m.measurement_type  = t.measurement_type
        AND m.measurement_value IS NOT DISTINCT FROM t.measurement_value)")$n == 0)

# 5. retained_reason populated with the expected categories. 'bottle' joins
#    'grid'/'inflection' only when a canonical bottle-grain type exists (see the
#    assemble chunk), so it is required conditionally rather than always.
reason_breakdown <- dbGetQuery(con,
  "SELECT retained_reason, COUNT(*) AS n, COUNT(DISTINCT depth_m) AS n_depths
   FROM ctd_thin GROUP BY 1 ORDER BY n DESC")
expected_reasons <- c("grid", "inflection", if (length(canon_btl)) "bottle")
stopifnot(
  "ctd_thin retained_reason has unexpected values" =
    all(reason_breakdown$retained_reason %in% expected_reasons),
  "ctd_thin retained_reason missing a category" =
    all(expected_reasons %in% reason_breakdown$retained_reason))
# previously computed only for the assertion and never shown
say("retained_reason breakdown:")

retained_reason breakdown:

Code

print(reason_breakdown)

  retained_reason        n n_depths
1            grid 11415516     2077
2      inflection  4356805     1586
3          bottle  3070794      612

Code

# 5b. the point of the 'bottle' reason: every canonical bottle-grain measurement on
#     a chosen-direction cast must survive thinning, not just those that happened to
#     land on a sensor-derived depth. Assert 100% retention rather than trusting it.
if (length(canon_btl)) {
  btl_in <- paste(sprintf("'%s'", canon_btl), collapse = ", ")
  btl_cov <- dbGetQuery(con, glue("
    WITH avail AS (
      SELECT m.ctd_cast_uuid, m.depth_m, m.measurement_type
      FROM ctd_measurement m JOIN _chosen_dir cd USING (ctd_cast_uuid)
      WHERE m.measurement_type IN ({btl_in}))
    SELECT (SELECT COUNT(*) FROM avail) AS available,
           (SELECT COUNT(*) FROM ctd_thin WHERE measurement_type IN ({btl_in})) AS retained"))
  say(glue("bottle-grain canonical: {format(btl_cov$retained, big.mark = ',')} of ",
           "{format(btl_cov$available, big.mark = ',')} retained ",
           "({round(100 * btl_cov$retained / btl_cov$available, 1)}%)"))
  stopifnot("every canonical bottle measurement on the chosen direction must be retained" =
              btl_cov$retained == btl_cov$available)
}

bottle-grain canonical: 1,581,859 of 1,581,859 retained (100%)

Code

# 6. one cast direction per physical cast; no physical cast dropped entirely
coverage <- dbGetQuery(con, glue("
  WITH meas_casts AS (
    SELECT DISTINCT cd.phys_cast_id
    FROM ctd_measurement m
    JOIN _chosen_dir cd ON m.ctd_cast_uuid = cd.ctd_cast_uuid
    WHERE m.measurement_type IN ({canon_in})),
  thin AS (
    SELECT cd.phys_cast_id, COUNT(DISTINCT t.cast_dir) AS n_dir
    FROM ctd_thin t
    JOIN _chosen_dir cd ON t.ctd_cast_uuid = cd.ctd_cast_uuid
    GROUP BY cd.phys_cast_id)
  SELECT
    (SELECT COUNT(*) FROM thin)                  AS n_thin_casts,
    (SELECT COUNT(*) FROM thin WHERE n_dir > 1)  AS n_multi_dir,
    (SELECT COUNT(*) FROM (
       SELECT phys_cast_id FROM meas_casts
       EXCEPT
       SELECT phys_cast_id FROM thin))           AS n_dropped"))
stopifnot(
  "ctd_thin: physical cast with >1 direction" = coverage$n_multi_dir == 0,
  "ctd_thin: physical cast dropped entirely"  = coverage$n_dropped   == 0)
cat(glue("{format(coverage$n_thin_casts, big.mark = ',')} physical casts ",
         "retained (single direction each)"), "\n")

9,675 physical casts retained (single direction each)

Code

# 7. depth-gap distribution between retained samples
gap_stats <- dbGetQuery(con, "
  WITH g AS (
    SELECT t.depth_m - LAG(t.depth_m) OVER (
             PARTITION BY cd.phys_cast_id, t.measurement_type
             ORDER BY t.depth_m) AS gap_m
    FROM ctd_thin t
    JOIN _chosen_dir cd ON t.ctd_cast_uuid = cd.ctd_cast_uuid)
  SELECT
    ROUND(MEDIAN(gap_m), 2)              AS median_gap_m,
    ROUND(QUANTILE_CONT(gap_m, 0.95), 2) AS p95_gap_m,
    ROUND(MAX(gap_m), 2)                 AS max_gap_m
  FROM g WHERE gap_m IS NOT NULL")
cat(glue("depth gaps between retained samples: median {gap_stats$median_gap_m} m, ",
         "p95 {gap_stats$p95_gap_m} m, max {gap_stats$max_gap_m} m"), "\n")

depth gaps between retained samples: median 7 m, p95 11 m, max 2983 m

Code

# 8. every ctd_thin type is registered as canonical in measurement_type
stopifnot("ctd_thin types not registered as canonical in measurement_type" =
  length(dbGetQuery(con, "
    SELECT DISTINCT t.measurement_type
    FROM ctd_thin t
    LEFT JOIN measurement_type mt ON t.measurement_type = mt.measurement_type
    WHERE mt.measurement_type IS NULL OR NOT mt.is_canonical")$measurement_type) == 0)

# done with _chosen_dir — drop before downstream metadata helpers enumerate tables
dbExecute(con, "DROP TABLE IF EXISTS _chosen_dir")

[1] 0

Code

reason_breakdown |>
  dt(caption = "ctd_thin rows by retained_reason", fname = "ctd_thin_reason")

24.7 Spot-check thinned profiles

Overlay the full ctd_measurement profile (line) against the ctd_thin samples (points) for the six temperature profiles with the most retained inflection points — confirming thermocline curvature is tracked while flat segments stay sparse.

Code

# a ctd_cast_uuid is a single depth sample, not a profile — group by the
# physical cast (cruise_key + cast_key). pick six temperature profiles with the
# most retained inflection points (the hardest cases to thin faithfully).
cast_map <- "(SELECT DISTINCT ctd_cast_uuid, cast_key FROM ctd_cast)"

spot <- dbGetQuery(con, glue("
  WITH thin AS (
    SELECT t.cruise_key, cm.cast_key, t.retained_reason
    FROM ctd_thin t
    JOIN {cast_map} cm ON t.ctd_cast_uuid = cm.ctd_cast_uuid
    WHERE t.measurement_type = 'temperature_ave')
  SELECT cruise_key, cast_key
  FROM thin
  GROUP BY cruise_key, cast_key
  HAVING COUNT(*) >= 20
  ORDER BY COUNT(*) FILTER (WHERE retained_reason = 'inflection') DESC,
           cruise_key, cast_key
  LIMIT 6"))
spot$cast <- paste("cast", seq_len(nrow(spot)))
spot_in <- if (nrow(spot) > 0) {
  paste(sprintf("('%s', '%s')", spot$cruise_key, spot$cast_key), collapse = ", ")
} else "('', '')"

# full ctd_measurement profile (de-duplicated) and the thinned samples
d_full <- dbGetQuery(con, glue("
  SELECT t.cruise_key, cm.cast_key, t.depth_m,
         AVG(t.measurement_value) AS measurement_value
  FROM ctd_measurement t
  JOIN {cast_map} cm ON t.ctd_cast_uuid = cm.ctd_cast_uuid
  WHERE t.measurement_type = 'temperature_ave'
    AND (t.cruise_key, cm.cast_key) IN ({spot_in})
  GROUP BY t.cruise_key, cm.cast_key, t.depth_m
  ORDER BY t.cruise_key, cm.cast_key, t.depth_m"))
d_thin <- dbGetQuery(con, glue("
  SELECT t.cruise_key, cm.cast_key, t.depth_m, t.measurement_value,
         t.retained_reason
  FROM ctd_thin t
  JOIN {cast_map} cm ON t.ctd_cast_uuid = cm.ctd_cast_uuid
  WHERE t.measurement_type = 'temperature_ave'
    AND (t.cruise_key, cm.cast_key) IN ({spot_in})
  ORDER BY t.cruise_key, cm.cast_key, t.depth_m"))

# short facet labels keyed on cruise_key + cast_key
prof_key  <- function(d) paste(d$cruise_key, d$cast_key)
cast_lab  <- setNames(spot$cast, prof_key(spot))
d_full$cast <- cast_lab[prof_key(d_full)]
d_thin$cast <- cast_lab[prof_key(d_thin)]

ggplot2::ggplot() +
  ggplot2::geom_path(
    data = d_full, ggplot2::aes(measurement_value, depth_m),
    color = "grey70", linewidth = 0.3) +
  ggplot2::geom_point(
    data = d_thin,
    ggplot2::aes(measurement_value, depth_m, color = retained_reason),
    size = 1.1) +
  ggplot2::scale_y_reverse() +
  ggplot2::scale_color_manual(
    values = c(grid = "#0077cc", inflection = "#ff6600")) +
  ggplot2::facet_wrap(~ cast, scales = "free", nrow = 2) +
  ggplot2::labs(
    x = "temperature_ave (degC)", y = "depth (m)", color = "retained",
    title = "Full ctd_measurement profile (line) vs ctd_thin samples (points)") +
  ggplot2::theme_minimal(base_size = 9) +
  ggplot2::theme(legend.position = "bottom")

25 Schema Diagram

Code

# define PK/FK relationships for visualization and relationships.json
ctd_rels <- list(
  primary_keys = list(
    ctd_cast        = "ctd_cast_uuid",
    ctd_measurement = "ctd_measurement_uuid",
    ctd_thin        = "ctd_thin_uuid",
    ctd_summary     = "ctd_summary_uuid",
    measurement_type = "measurement_type"),
  foreign_keys = list(
    list(table = "ctd_measurement", column = "ctd_cast_uuid", ref_table = "ctd_cast", ref_column = "ctd_cast_uuid"),
    list(table = "ctd_thin", column = "ctd_cast_uuid", ref_table = "ctd_cast", ref_column = "ctd_cast_uuid"),
    list(table = "ctd_measurement", column = "measurement_type", ref_table = "measurement_type", ref_column = "measurement_type"),
    list(table = "ctd_thin", column = "measurement_type", ref_table = "measurement_type", ref_column = "measurement_type"),
    list(table = "ctd_summary", column = "measurement_type", ref_table = "measurement_type", ref_column = "measurement_type")))

cc_erd(con, rels = ctd_rels)

26 Add Spatial

Code

add_point_geom(con, "ctd_cast", lon_col = "longitude", lat_col = "latitude")

# grid already loaded in Cross-Dataset Bridge section
assign_grid_key(con, "ctd_cast")

   status       n
1 in_grid 7485734

27 Derived View

Code

dbExecute(
  con,
  "
  CREATE OR REPLACE VIEW ctd_cast_derived AS
  SELECT *,
    datetime_start_utc AT TIME ZONE 'UTC' AT TIME ZONE 'US/Pacific' AS datetime_pst,
    EXTRACT(YEAR FROM datetime_start_utc)  AS year,
    EXTRACT(MONTH FROM datetime_start_utc) AS month,
    EXTRACT(DOY FROM datetime_start_utc)   AS julian_day
  FROM ctd_cast"
)

[1] 0

28 Load Dataset Metadata

Code

# dataset registry built from authoritative ingest_*.qmd YAML (was dataset.csv)
d_dataset <- ingest_yaml_to_dataset_df(read_ingest_yaml(here()))
dbWriteTable(con, "dataset", d_dataset, overwrite = TRUE)
say(glue("dataset: {nrow(d_dataset)} datasets registered"))

dataset: 16 datasets registered

29 Data Quality Diagnostics

Range and sentinel audit over ctd_measurement. This exists because the v2026.07.17 release shipped values that cannot be real — salinity_ave_corr spanning −45 to 1016 PSU, oxygen maxima around 1e9–1e10 — alongside -99 sentinels that were ingested as genuine measurements. Nothing in the pipeline flagged either, so this section makes both visible at ingest time.

Empty on a fast re-render — see Findings instead

This section audits the intermediate ctd_measurement table, which only exists when the heavy path ran. On a re-render with unchanged inputs its chunks are shown but not evaluated, so they have no output (see Check for Resumable State).

That is not a gap in coverage: the equivalent checks over the published obs and obs_ctd_full are in Findings, which reads the parquet and therefore always runs. This section is the pre-write audit; that one is the as-shipped report.

-99 is fixed at the source; the extremes are not

-99 is handled. The pivot above now deletes it, because the notebook already treated -99 as this source’s missing marker for longitude/latitude — the rule had simply never been extended to the measurement columns. The audit below should therefore report zero -99 rows; if it does not, that fix has regressed.

The extremes are still only reported. They do not stopifnot(): until the providers confirm whether values like 2.1e9 mL/L oxygen are bad scans or a units error (question 02), dropping them would be guessing. The plausible ranges below are a first pass and belong in metadata/measurement_type.csv as valid_min/valid_max once agreed — which would also let the CF netCDF writer emit them as real valid_min/valid_max variable attributes.

Code

# audit trail for the -99 deletion in the pivot above
if (nrow(sentinel_tally)) {
  sentinel_tally |>
    mutate(pct_of_removed = round(100 * n / sum(n), 1)) |>
    datatable(
      caption  = glue("-99 sentinel rows removed from ctd_measurement ",
                      "({format(n_sentinel, big.mark = ',')} total)"),
      rownames = FALSE,
      options  = list(pageLength = 10, scrollX = TRUE, dom = "tip"))
} else {
  say("no -99 sentinel values found in ctd_measurement")
}

Code

# Plausible physical ranges now come from the REGISTRY, not a tribble local to this
# chunk — `metadata/measurement_type.csv` columns valid_min / valid_max, populated by
# libs/build_ctd_measurement_registry.R. Same values as the inline table they replace;
# declaring them once means they are reviewable in a diff, emittable by the CF netCDF
# writer as real valid_min/valid_max variable attributes, and available to any QC rule
# rather than only to this audit.
#
# Still deliberately generous, and still report-only: they catch impossible values,
# they are not the agreed oceanographic ranges, and nothing below drops a row.
plaus <- d_meas_type |>
  filter(!is.na(valid_min) | !is.na(valid_max)) |>
  select(measurement_type, lo = valid_min, hi = valid_max)

diag_sql <- glue("
  SELECT measurement_type,
         COUNT(*)                                           AS n,
         COUNT(*) FILTER (WHERE measurement_value = -99)     AS n_neg99,
         COUNT(*) FILTER (WHERE isnan(measurement_value)
                             OR NOT isfinite(measurement_value)) AS n_nonfinite,
         MIN(measurement_value)                             AS v_min,
         MAX(measurement_value)                             AS v_max
  FROM ctd_measurement GROUP BY 1")
d_diag <- dbGetQuery(con, diag_sql) |>
  left_join(plaus, by = "measurement_type") |>
  mutate(
    pct_neg99   = round(100 * n_neg99 / n, 2),
    # count out-of-range only where a plausible range is declared, and exclude the
    # -99 sentinel so it is reported once (as a sentinel) rather than twice
    has_range   = !is.na(lo),
    out_of_range = NA_integer_)

for (i in which(d_diag$has_range)) {
  r <- d_diag[i, ]
  d_diag$out_of_range[i] <- dbGetQuery(con, glue("
    SELECT COUNT(*) AS n FROM ctd_measurement
    WHERE measurement_type = '{r$measurement_type}'
      AND measurement_value <> -99
      AND (measurement_value < {r$lo} OR measurement_value > {r$hi})"))$n
}

say(glue("types audited                     : {nrow(d_diag)}"))

types audited                     : 55

Code

say(glue("  with a declared plausible range : {sum(d_diag$has_range)}"))

  with a declared plausible range : 50

Code

say(glue("  WITHOUT a range (unchecked)     : {sum(!d_diag$has_range)}"))

  WITHOUT a range (unchecked)     : 5

Code

say(glue("-99 sentinel values (all types)   : ",
         "{format(sum(d_diag$n_neg99), big.mark = ',')}"))

-99 sentinel values (all types)   : 0

Code

say(glue("out-of-range values (ranged types): ",
         "{format(sum(d_diag$out_of_range, na.rm = TRUE), big.mark = ',')}"))

out-of-range values (ranged types): 0

Code

d_diag |>
  transmute(
    measurement_type, n,
    `-99` = n_neg99, `%-99` = pct_neg99, `non-finite` = n_nonfinite,
    min = signif(v_min, 4), max = signif(v_max, 4),
    range = if_else(has_range, glue("{lo}..{hi}"), glue("—")),
    `out of range` = out_of_range) |>
  arrange(desc(`-99` + coalesce(`out of range`, 0L)), measurement_type) |>
  datatable(
    caption  = paste("ctd_measurement range + sentinel audit —",
                     "rows with -99 or out-of-range values sort first"),
    rownames = FALSE,
    options  = list(pageLength = 12, scrollX = TRUE, dom = "ftip")) |>
  formatStyle("-99",          color = styleInterval(0, c("inherit", "#d03b3b"))) |>
  formatStyle("out of range", color = styleInterval(0, c("inherit", "#d03b3b")))

29.1 Quality flags — decoded against the controlled vocabulary

measurement_qual was a verbatim pass-through for years: carried through every table, never interpreted, with no vocabulary anywhere in the repo. There is one now — metadata/measurement_qual.csv, recovered from the CalCOFI hydrographic master’s own field documentation (6 = data OK but taken from CTD, 8 = value is suspect, 9 = missing data) and confirmed to apply to CTD, whose flags use the same codes.

Decoded and reported — deliberately not acted on

Nothing here drops or rewrites a flagged value. Three questions are open with the data providers (metadata/calcofi/hydro-master/questions.csv, two of them blockers): what the undocumented single-digit codes mean, whether a NULL flag means “good” or “never assessed”, and whether S_qual’s 253 distinct values in the bottle master indicate a bitmask or a corrupted column. Acting on a flag whose semantics are unconfirmed would be guessing with released data.

The headline variables carry no flags at all, and that is the more important gap: source flags attach to the component sensors (Temp1Q, Salt1Q, Ox1Q, Ox2Q), while the canonical types are the averages — so quality information is lost in the mean. Deciding that propagation rule (worst-of? both-must-pass?) is ctd-qaqc work, not something to invent here.

Code

# qual_code is forced to character: every value in the vocabulary is numeric, so a
# type-guessing read makes it a double, and joining that to the VARCHAR
# measurement_qual is a hard error rather than a silent miss.
d_qual_vocab <- read_csv(
  here("metadata/measurement_qual.csv"),
  col_types = cols(qual_code = col_character(), .default = col_guess())) |>
  select(measurement_qual = qual_code, label, description, is_documented)

d_qual_obs <- dbGetQuery(con, "
  SELECT measurement_qual, COUNT(*) AS n,
         COUNT(DISTINCT measurement_type) AS n_types
  FROM ctd_measurement
  WHERE measurement_qual IS NOT NULL
  GROUP BY 1 ORDER BY n DESC") |>
  left_join(d_qual_vocab, by = "measurement_qual") |>
  mutate(
    label         = coalesce(label, "UNRECOGNIZED"),
    is_documented = coalesce(is_documented, FALSE))

n_flagged   <- sum(d_qual_obs$n)
n_unrecog   <- sum(d_qual_obs$n[d_qual_obs$label == "UNRECOGNIZED"])
n_canon_flg <- dbGetQuery(con, glue("
  SELECT COUNT(DISTINCT measurement_type) AS n FROM ctd_measurement
  WHERE measurement_qual IS NOT NULL
    AND measurement_type IN ({canon_in})"))$n

say(glue("flagged measurements        : {format(n_flagged, big.mark = ',')}"))

flagged measurements        : 365,496

Code

say(glue("  codes seen                : {nrow(d_qual_obs)}"))

  codes seen                : 5

Code

say(glue("  NOT in the vocabulary     : {format(n_unrecog, big.mark = ',')}"))

  NOT in the vocabulary     : 1,092

Code

say(glue("canonical types carrying any flag: {n_canon_flg} of {length(canon_types)}"))

canonical types carrying any flag: 24 of 44

Code

# a code outside the vocabulary is a finding, not noise — surface it loudly rather
# than letting a left_join quietly render it as NA in the table below
if (n_unrecog > 0) {
  say(glue("  unrecognized code(s): ",
           "{paste(d_qual_obs$measurement_qual[d_qual_obs$label == 'UNRECOGNIZED'], collapse = ', ')}"))
}

  unrecognized code(s): nan

Code

d_qual_obs |>
  transmute(code = measurement_qual, label, description,
            documented = is_documented, n, `types` = n_types) |>
  datatable(
    caption  = paste("ctd_measurement quality flags decoded against",
                     "metadata/measurement_qual.csv"),
    rownames = FALSE,
    options  = list(pageLength = 10, scrollX = TRUE, dom = "tip")) |>
  formatStyle("documented",
              color = styleEqual(c(TRUE, FALSE), c("inherit", "#d03b3b")))

30 Validate and Enforce Types

Code

validate_for_release(con)

$passed
[1] TRUE

$checks
# A tibble: 11 × 4
   check      table            status message       
   <chr>      <chr>            <chr>  <glue>        
 1 row_counts cruise           pass   691 rows      
 2 row_counts ctd_cast         pass   7485734 rows  
 3 row_counts ctd_cast_derived pass   7485734 rows  
 4 row_counts ctd_measurement  pass   284379888 rows
 5 row_counts ctd_src_format   pass   377 rows      
 6 row_counts ctd_summary      pass   140809234 rows
 7 row_counts ctd_thin         pass   18843115 rows 
 8 row_counts dataset          pass   16 rows       
 9 row_counts grid             pass   218 rows      
10 row_counts measurement_type pass   202 rows      
11 row_counts ship             pass   48 rows       

$errors
character(0)

$warnings
character(0)

Code

enforce_column_types(con, d_flds_rd = d_flds_rd)

# A tibble: 0 × 5
# ℹ 5 variables: table <chr>, column <chr>, from_type <chr>, to_type <chr>,
#   success <lgl>

31 Preview Tables

One datatable() per table, each in its own chunk — preview_tables() in a loop renders only the first DT widget (see .claude/skills/ingest-new/SKILL.md), which is why several of these previews were previously invisible. Every preview is horizontally scrollable (scrollX) and column-capped: ctd_cast is ~40 columns wide and overflowed the page as a static table.

Code

# shared preview: LIMIT in SQL (never collect a 216M-row table), scrollX so wide
# tables scroll inside their own box instead of overflowing the article, and
# geometry dropped (the R DuckDB driver cannot marshal GEOMETRY).
preview_tbl <- function(tbl, n = 100, cols = NULL) {
  all_cols <- dbGetQuery(con, glue(
    "SELECT column_name, data_type FROM information_schema.columns
      WHERE table_name = '{tbl}'"))
  keep <- all_cols$column_name[all_cols$data_type != "GEOMETRY"]
  if (!is.null(cols)) keep <- intersect(keep, cols)
  n_row <- dbGetQuery(con, glue("SELECT COUNT(*) AS n FROM {tbl}"))$n
  d <- dbGetQuery(con, glue(
    "SELECT {paste(keep, collapse = ', ')} FROM {tbl} LIMIT {n}"))
  datatable(
    d, rownames = FALSE,
    caption = glue("{tbl} — {format(n_row, big.mark = ',')} rows x ",
                   "{nrow(all_cols)} cols (showing first {nrow(d)} rows",
                   if (!is.null(cols)) ", selected cols" else "", ")"),
    options = list(pageLength = 5, scrollX = TRUE, dom = "tip"))
}

Code

# ~40 columns; show the identity + position core and let scrollX handle the rest
preview_tbl("ctd_cast", cols = c(
  "ctd_cast_uuid", "cruise_key", "site_key", "cast_key", "cast_dir",
  "datetime_start_utc", "latitude", "longitude", "ship_key", "data_stage",
  "ord_occ"))

Code

preview_tbl("ctd_measurement")

Code

preview_tbl("ctd_thin")

Code

preview_tbl("ctd_summary")

Code

preview_tbl("measurement_type")

32 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.

Two grains ship from here. The default obs carries the thinned ctd_thin series (~5.5M rows) at realm = 'env', which is what every consumer reads. The full-resolution scan set becomes obs_ctd_full (~216M rows), a supplemental output that is hosted and catalogued but excluded from the default table list — cc_get_db(supplemental = TRUE) opts into it. Building it was previously the release’s job, reaching back into this ingest’s ctd_measurement parquet; it belongs here, with the dataset that produces it.

sample is one row per physical cast (cast_key), not per scan: ctd_cast is per-scan, so the arm deduplicates to the cast before minting the event.

Code

ds_key <- "calcofi_ctd-cast"

# This projection lives here, in the notebook that owns the dataset, rather than in
# a switch(dataset_key, …) arm inside calcofi4db. No taxa: CTD is the env realm.
#
# 1. sample — one row per PHYSICAL cast. ctd_cast is per-scan, so the QUALIFY is
#    load-bearing: without it `sample` would get one row per scan (see the
#    scan-grain gotcha in CLAUDE.md).
#    16 columns, not 15: `data_stage` is append_sample()'s optional trailing
#    column, so the final/preliminary distinction the source insists on ("for
#    non-publication use") reaches consumers instead of dying here (question
#    calcofi_ctd-cast_14).
append_sample(con, glue("
  SELECT * FROM (
    SELECT {ns_key(ds_key, 'cast', 'cast_key')} AS sample_key, 'cast' AS sample_type,
           NULL::VARCHAR AS parent_sample_key,
           {ns_key(ds_key, 'cast', 'cast_key')} AS root_sample_key,
           '{ds_key}' AS dataset_key, grid_key, site_key, cruise_key,
           TRY_CAST(ord_occ AS INTEGER) AS order_occ, latitude, longitude,
           CAST(datetime_start_utc AS TIMESTAMP) AS datetime,
           NULL::DOUBLE AS depth_min_m, NULL::DOUBLE AS depth_max_m,
           NULL::VARCHAR AS tow_type, data_stage
    FROM ctd_cast
  ) q QUALIFY row_number() OVER (PARTITION BY sample_key ORDER BY datetime) = 1"))

# 2. obs — the THINNED scan set is the headline (ctd_thin, not ctd_measurement);
#    full resolution goes to the supplemental obs_ctd_full below.
append_obs(con, glue("
  SELECT 'env', '{ds_key}',
         {ns_key(ds_key, 'cast', 'cc.cast_key')},
         cc.grid_key, cc.cruise_key, cc.latitude, cc.longitude,
         CAST(cc.datetime_start_utc AS TIMESTAMP), t.depth_m, t.depth_m,
         NULL::VARCHAR, NULL::VARCHAR, t.measurement_type, t.measurement_value,
         t.measurement_qual, NULL::DOUBLE
  FROM ctd_thin t JOIN ctd_cast cc ON t.ctd_cast_uuid = cc.ctd_cast_uuid"))

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

core projection — sample=19330 obs=18843115

Code

n_obs <- dbGetQuery(con, "SELECT COUNT(*) FROM obs")[[1]]
n_exp <- dbGetQuery(con,
  "SELECT COUNT(*) FROM ctd_thin t JOIN ctd_cast cc ON t.ctd_cast_uuid = cc.ctd_cast_uuid")[[1]]
n_cast <- dbGetQuery(con, "SELECT COUNT(DISTINCT cast_key) FROM ctd_cast")[[1]]
stopifnot(
  "obs must be one row per thinned scan" = n_obs == n_exp,
  "sample must be one row per PHYSICAL cast, not per scan" =
    dbGetQuery(con, "SELECT COUNT(*) FROM sample")[[1]] == n_cast,
  "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,
  # ctd_cast is filtered to CTD_DATA_STAGES upstream, so a NULL here means the
  # 16th positional column silently fell off the arm rather than that a cast has
  # no stage
  "every sample row must carry a known data_stage" =
    dbGetQuery(con, glue(
      "SELECT COUNT(*) FROM sample
       WHERE data_stage NOT IN ({stages_sql}) OR data_stage IS NULL",
      stages_sql = paste0("'", CTD_DATA_STAGES, "'", collapse = ", ")))[[1]] == 0,
  # the whole point of ctd_thin: obs must not silently become the full scan set
  "obs must be fed by ctd_thin, not ctd_measurement" =
    n_obs < dbGetQuery(con, "SELECT COUNT(*) FROM ctd_measurement")[[1]])
cat(glue("obs parity: {format(n_obs, big.mark=',')} thinned rows across ",
         "{format(n_cast, big.mark=',')} casts"), "\n")

obs parity: 18,843,115 thinned rows across 19,330 casts

Code

# --- data_stage is classified by FILENAME; verify it against CONTENT ----------
# The `_CTDBTL_` / `_CTD_` token is a convention, not a guarantee, and 2111SR is
# classified by a hardcoded cruise rule that no token supports at all. What
# actually distinguishes the tiers is whether the bottle merge has run, so test
# that directly: a `preliminary_without_bottle` cast must have NO bottle-derived
# or bottle-corrected observations, and the other two tiers must have some.
#
# This is what catches a re-organized or mis-named archive — a case the filename
# rule cannot see, and which would otherwise reach consumers as a cruise silently
# missing its salinity and oxygen.
bottle_dependent <- c(
  "salinity_ave_corr", "oxygen_ml_l_ave_sta_corr", "oxygen_umol_kg_ave_sta_corr",
  "salinity_btl", "oxygen_btl_ml_l", "btl_depth")

d_stage_content <- dbGetQuery(con, glue("
  SELECT s.data_stage,
         COUNT(DISTINCT s.sample_key)                             AS n_casts,
         COUNT(*) FILTER (WHERE o.measurement_type IN ({btl_sql})) AS n_bottle_obs
  FROM sample s LEFT JOIN obs o USING (sample_key)
  GROUP BY s.data_stage",
  btl_sql = paste0("'", bottle_dependent, "'", collapse = ", ")))

stopifnot(
  "a preliminary_without_bottle cast must carry NO bottle-merged measurement" =
    all(d_stage_content$n_bottle_obs[
      d_stage_content$data_stage == "preliminary_without_bottle"] == 0),
  "final and preliminary_with_bottle casts must carry bottle-merged values" =
    all(d_stage_content$n_bottle_obs[
      d_stage_content$data_stage %in%
        c("final", "preliminary_with_bottle")] > 0))

d_stage_content |>
  mutate(stage_label = CTD_STAGE_LABELS[data_stage]) |>
  arrange(match(data_stage, CTD_DATA_STAGES)) |>
  dt(
    caption = paste(
      "`data_stage` verified against content: bottle-merged observations are",
      "present in the two merged tiers and absent from the sensor-only tier"),
    fname = "ctd_data_stage_content"
  )

Code

# --- obs_ctd_full: supplemental full-resolution scans -----------------------
# heavy (~216M rows); set BUILD_OBS_CTD_FULL=FALSE for a fast structural render.
build_obs_ctd_full <- as.logical(Sys.getenv("BUILD_OBS_CTD_FULL", "TRUE"))
if (build_obs_ctd_full) {
  n_full <- append_obs(con, obs_tbl = "obs_ctd_full", select_sql = "
    SELECT 'env' realm, 'calcofi_ctd-cast' dataset_key,
           'calcofi_ctd-cast:cast:' || CAST(cc.cast_key AS VARCHAR) sample_key,
           cc.grid_key, cc.cruise_key, cc.latitude, cc.longitude,
           CAST(cc.datetime_start_utc AS TIMESTAMP) datetime,
           m.depth_m depth_min_m, m.depth_m depth_max_m,
           NULL::VARCHAR taxon_key, NULL::VARCHAR life_stage,
           m.measurement_type, m.measurement_value, m.measurement_qual,
           NULL::DOUBLE measurement_prec
    FROM ctd_measurement m
    JOIN ctd_cast cc ON m.ctd_cast_uuid = cc.ctd_cast_uuid")
  stopifnot(
    "obs_ctd_full must cover every scan" = n_full ==
      dbGetQuery(con, "SELECT COUNT(*) FROM ctd_measurement m
                       JOIN ctd_cast cc ON m.ctd_cast_uuid = cc.ctd_cast_uuid")[[1]])
  cat(glue("obs_ctd_full: {format(n_full, big.mark=',')} rows"), "\n")

  # Verify the bounds guard reached the SUPPLEMENTAL table too. It should — this
  # is built from the same `ctd_measurement` the guard cleaned — but "should"
  # is exactly what went unchecked: v2026.08.07 published an obs_ctd_full whose
  # `ph` ran to -2.98 (5,963 values below the declared floor) while the staged
  # ingest output for the same partition was clean. The released bytes and the
  # ingest's own bytes disagreed and nothing compared them.
  #
  # Assert rather than report: if this table ever stops deriving from the guarded
  # table, that is a silent regression across 216M published rows.
  b_full <- check_measurement_bounds(con, "obs_ctd_full", mt = d_meas_type)
  bounds_datatable(b_full)
  cat(glue("obs_ctd_full bounds: {sum(b_full$status == 'ok')} ok, ",
           "{sum(b_full$status == 'undeclared')} undeclared ",
           "(the derived _sta_corr/_cruise_corr oxygen types — Q24)"), "\n")
  stopifnot(
    "obs_ctd_full must inherit the bounds guard applied to ctd_measurement" =
      sum(b_full$status == "out_of_range") == 0)
} else {
  cat("obs_ctd_full skipped (BUILD_OBS_CTD_FULL=FALSE)\n")
}

obs_ctd_full: 284,379,888 rows 
obs_ctd_full bounds: 50 ok, 5 undeclared (the derived _sta_corr/_cruise_corr oxygen types — Q24)

32.1 No cruise may leave `obs` and keep its casts

The guards above are upstream and specific — an empty read, a name with no direction letter, a cruise with no downcast. This one is downstream and general: whatever the cause, a cruise that reaches this point with casts and no observations must not be written, because write_parquet_outputs() deletes the partition that has left the data and sync_to_gcs() mirrors the deletion, so the cruise disappears from the release while its casts remain in sample. Every foreign key still resolves — FK validation runs child → parent, and a parent with no children breaks nothing — which is why v2026.08.08 published the loss.

The correct count for this dataset is zero, so it is asserted rather than ratcheted. sample legitimately carries about half its cast rows without observations (one row per cast per direction, against one direction kept in obs), which is why the grain is the cruise and not the sample.

Code

# `d_cruise_cov`, never `cov` — that shadows stats::cov for the rest of the
# render, and a short name is exactly how the release chunk clobbered its own
# `d_cov` and died after a 50-minute freeze.
d_cruise_cov <- check_cruise_coverage(con, max_orphan_cruises = 0L)

d_cruise_cov |>
  dt(
    caption = "Cruise coverage: cruises carrying samples with no observations",
    fname   = "ctd_cruise_coverage")

33 Write Parquet

Code

# collect mismatches for manifest
mismatches <- list(
  ships             = collect_ship_mismatches(con, "ctd_cast"),
  measurement_types = collect_measurement_type_mismatches(
    con, here("metadata/measurement_type.csv")),
  cruise_keys       = collect_cruise_key_mismatches(con, "ctd_cast"))

# core shards + the supplemental outputs this dataset owns. ctd_cast/ctd_thin/
# ctd_measurement/ctd_summary are superseded by sample/obs/obs_ctd_full, and
# ctd_wide is retired (ERDDAP reads DuckDB views + netCDF now).
tbls_out <- core_output_tables(
  con, extra = c("obs_ctd_full", "measurement_type"))
parquet_stats <- write_parquet_outputs(
  con          = con,
  output_dir   = dir_parquet,
  tables       = tbls_out,
  partition_by = list(
    obs          = "cruise_key",
    obs_ctd_full = "cruise_key"),
  # Sort key MUST match `core_sort` in release_database.qmd, so the shard the
  # release reads is already clustered the way it re-exports it. Two parts
  # matter and both were measured, not guessed:
  #   grid_key FIRST co-locates every column that is a function of the cast —
  #     sample_key, hex_id, latitude, longitude, datetime — which is most of the
  #     row width. measurement_type-first clusters one 54-value column and
  #     scatters the rest: 4.61 GB vs 1.22 GB for obs_ctd_full.
  #   datetime LAST is not decoration. (grid_key, depth_min_m, measurement_type)
  #     leaves large tie groups, and rows inside a tie land in arbitrary order,
  #     scattering lat/lon/datetime again. Adding the tiebreak made a partition
  #     27.55 -> 20.20 MB (CTD) and 23.22 -> 16.95 MB (mets), i.e. ~27% below
  #     what the release's own re-export produced. Do not drop it.
  sort_by = list(
    obs          = c("measurement_type", "depth_min_m"),
    obs_ctd_full = c("grid_key", "depth_min_m", "measurement_type", "datetime"),
    sample       = "hilbert:longitude,latitude"),
  strip_provenance = FALSE,
  mismatches       = mismatches,
  supplemental     = c("obs_ctd_full")
)

parquet_stats |>
  dt(fname = "ctd_parquet_stats")

34 Write Metadata JSON

Code

metadata_path <- 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"),
                           glue("{dir_meta}/metadata_derived.csv")),
  output_dir           = dir_parquet,
  tables               = tbls_out,
  provider             = provider,
  dataset              = dataset,
  workflow_url         = cc$workflow_url,
  tables_owned         = tables_owned
)

metadata.json documentation gaps (4 tables, 74 columns) — these render blank in cc_describe_table() / cc_db_catalog():
columns with no description_md: 30    measurement_type.valid_min, measurement_type.valid_max, measurement_type.valid_depth_min_m, measurement_type.valid_depth_max_m, measurement_type.derivation, measurement_type.grain, measurement_type.denominator, measurement_type.category, measurement_type.variable, measurement_type.nerc_p01, measurement_type.units_nerc_p06, obs_ctd_full.obs_id (+18 more)
measurement columns with no units: 13    measurement_type.valid_min, measurement_type.valid_max, measurement_type.valid_depth_min_m, measurement_type.valid_depth_max_m, measurement_type.derivation, measurement_type.denominator, measurement_type.category, measurement_type.variable, measurement_type.nerc_p01, measurement_type.units_nerc_p06, obs_ctd_full.depth_min_m, obs_ctd_full.depth_max_m (+1 more)
  backfill via metadata/{provider}/{dataset}/flds_redefine.csv, then re-run

Code

# write relationships.json sidecar with PKs/FKs
build_relationships_json(
  rels       = core_relationships(tbls_out),
  output_dir = dir_parquet,
  provider   = provider,
  dataset    = dataset
)

[1] "/Users/bbest/Github/CalCOFI/workflows-ctd-new-format/data/parquet/calcofi_ctd-cast/relationships.json"

Code

# n_flags_applied (from "Apply Accepted QC Flags" above) is not a
# build_metadata_json() field, so it is patched onto the written sidecar
# rather than widening that helper's signature for one dataset's one number.
# release_database.qmd's qc_flags_pending chunk reads it to report how many of
# the snapshot's accepted flags this release does NOT yet reflect.
d_metadata_json <- jsonlite::read_json(metadata_path)
d_metadata_json$n_flags_applied <- if (exists("n_flags_applied")) n_flags_applied else 0L
jsonlite::write_json(d_metadata_json, metadata_path, auto_unbox = TRUE, pretty = TRUE, null = "null")

Both sidecars are read back from disk rather than from the objects the chunk above returned, and both viewers always evaluate. They describe files that exist in either mode, so gating them on the heavy path would have silently dropped two of this notebook’s published outputs from a fast render — which is exactly the kind of quiet degradation a skip path must not introduce.

Code

path_meta_json <- file.path(dir_parquet, "metadata.json")
if (file_exists(path_meta_json)) {
  listviewer::jsonedit(
    jsonlite::fromJSON(path_meta_json, simplifyVector = FALSE), mode = "view")
} else {
  say(glue("metadata.json not found at {path_meta_json}"))
}

Code

path_rels_json <- file.path(dir_parquet, "relationships.json")
if (file_exists(path_rels_json)) {
  listviewer::jsonedit(
    jsonlite::fromJSON(path_rels_json, simplifyVector = FALSE), mode = "view")
} else {
  say(glue("relationships.json not found at {path_rels_json}"))
}

35 Export Modified Dependency Deltas

Code

# export _new delta sidecars for modified dependency tables (from calcofi.modifies)
# these are NOT in the manifest — picked up by build_release_table_registry()
for (tbl in modifies_tables) {
  pk_col   <- modifies_pks[[tbl]]$pk_col
  keys_old <- modifies_pks[[tbl]]$keys
  keys_new <- dbGetQuery(con, glue("SELECT {pk_col} FROM {tbl}"))[[1]]
  additions <- setdiff(keys_new, keys_old)

  if (length(additions) > 0) {
    pq_path <- file.path(dir_stage, paste0(tbl, "_new.parquet"))
    vals <- paste(shQuote(additions, "sh"), collapse = ", ")
    export_parquet(con,
      glue("SELECT * FROM {tbl} WHERE {pk_col} IN ({vals})"), pq_path)
    say(glue("{length(additions)} new {tbl} row(s) -> {tbl}_new.parquet"))
  } else {
    say(glue("No new {tbl} rows — {tbl}_new not exported"))
  }
}

No new ship rows — ship_new not exported

36 Upload to GCS

Code

# always re-enable eval for upload and cleanup
knitr::opts_chunk$set(eval = TRUE)

Code

# Record the fingerprint as soon as the OUTPUTS are complete — after parquet,
# metadata and the delta sidecars, and BEFORE the upload.
#
# It used to be written last, after the GCS sync and the appendices, to avoid
# claiming "these inputs produced those outputs" for a run that died mid-pivot.
# That reasoning is right; the placement was too conservative. The outputs are
# finished HERE, and the upload is a separate concern that can fail for reasons
# having nothing to do with them — as it did: a 2 h 45 m run wrote every table
# correctly, then lost one object of a 3.2 GB transfer to a slow link, and threw
# away the whole rebuild because the fingerprint had not been reached yet. The
# identical upload succeeded by hand a minute later.
#
# Writing it here is still honest, because `parquet_ok` in [check_resume] does
# not take this file's word for anything: it independently verifies that every
# non-supplemental table named in the manifest is actually on disk. And the sync
# is no longer gated on the rebuild, so a fingerprint written before a failed
# upload cannot leave the mirror stale.
if (!parquet_complete) {
  write_input_fingerprint(fingerprint_path, fp)
  say(glue("Recorded input fingerprint {substr(fp$hash, 1, 12)} → {fingerprint_path}"))
}

Recorded input fingerprint 3ca4e41aa853 → /Users/bbest/Github/CalCOFI/workflows-ctd-new-format/data/wrangling/calcofi_ctd-cast_inputs.json

Code

# ALWAYS sync, even when the parquet was unchanged.
#
# This used to be gated on `!parquet_complete`, which coupled two independent
# facts: "the outputs were rebuilt" and "GCS matches the outputs". A run that
# rebuilt everything and then lost the upload to a network blip left GCS stale,
# and the *next* run — seeing complete parquet and an unchanged fingerprint —
# skipped the sync too, so the mirror stayed stale indefinitely with nothing
# reporting it.
#
# rsync is the right shape for this: it compares before it transfers, so a sync
# over an unchanged directory is a listing and no bytes. Paying that on every
# render is far cheaper than a silently stale mirror, and it makes the upload
# idempotent — safe to re-run after any failure.
sync_to_gcs(
  local_dir    = dir_stage,
  sidecar_dir  = dir_parquet,
  gcs_prefix   = glue("ingest/{dir_label}"),
  bucket       = "calcofi-db",
  delete_stale = TRUE
)

# A tibble: 25,629 × 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
 6 <rsync> uploaded    NA parallel rsync
 7 <rsync> uploaded    NA parallel rsync
 8 <rsync> uploaded    NA parallel rsync
 9 <rsync> uploaded    NA parallel rsync
10 <rsync> uploaded    NA parallel rsync
# ℹ 25,619 more rows

37 Findings

What this ingest just published, and the condition it is in. The sections above are the build; this is the report — the place a reader should be able to learn what is known to be wrong with CalCOFI CTD data without opening an app, a rule file or a release.

Three things make it worth having here rather than only in ctd-qaqc:

it reads the parquet that was just written, not the wrangling database, so it is checking what shipped rather than what was in memory before shipping;
it runs the same rule registry (metadata/qc_rules/) through the same engine (calcofi4db::qc_run_all()) that the app runs, so a number here and a number there cannot disagree;
it renders on a fast re-run (see Check for Resumable State), so the narrative can be revised without an hour of recomputation.

It sits after the GCS upload deliberately. Everything above is the ingest’s actual product and is already delivered by the time this runs; a broken diagnostic must not be able to discard it — the lesson the Gantt appendix records.

Code

# in-memory DuckDB with VIEWS over this run's parquet. Views cost nothing until
# queried, and reading the outputs rather than the wrangling tables is what makes
# this section render identically on a full run and on a fast re-run.
fcon <- get_duckdb_con(":memory:")
dbExecute(fcon, "SET threads TO 4")

[1] 0

Code

pq_view <- function(tbl, path) {
  if (!file_exists(path) && !dir_exists(path)) return(FALSE)
  src <- if (dir_exists(path))
    glue("read_parquet('{path}/**/*.parquet', hive_partitioning = true)") else
      glue("read_parquet('{path}')")
  dbExecute(fcon, glue("CREATE OR REPLACE VIEW {tbl} AS SELECT * FROM {src}"))
  TRUE
}
has_obs  <- pq_view("obs",          file.path(dir_stage, "obs"))
has_smp  <- pq_view("sample",       file.path(dir_stage, "sample.parquet"))
has_full <- pq_view("obs_ctd_full", file.path(dir_stage, "obs_ctd_full"))
stopifnot("no obs/sample parquet to report on" = has_obs && has_smp)

# Reference inputs the rules join against — the quality vocabulary, the harmonic
# climatology and station bottom depths mined from the Access master, and a
# seafloor depth per cast from the GEBCO crop that apps/ctd-viz already commits.
# None of it is part of the release. A missing input is left as a MISSING TABLE so
# its rules ERROR rather than returning zero rows and reading as clean.
gebco_tif <- normalizePath(
  file.path(here(), "..", "apps/ctd-viz/data/gebco_calcofi.tif"), mustWork = FALSE)
say("QC reference data staged:")

QC reference data staged:

Code

qc_ref <- qc_stage_reference(fcon, here(), gebco_tif = gebco_tif)

  measurement_type       202 rows (116 with a declared range)
  measurement_qual       15 codes
  climatology_harmonic   5,880 rows
  station                75 rows
  standard_depth         14 rows
  station_class          7 rows
  mld_sigma              620 rows
  nutclinedepth          585 rows
  sample_seafloor        19,330 casts (19,330 with a depth)

Code

f_head <- dbGetQuery(fcon, "
  SELECT COUNT(*)                                       AS n_obs,
         COUNT(DISTINCT sample_key)                     AS n_casts,
         COUNT(DISTINCT cruise_key)                     AS n_cruises,
         COUNT(DISTINCT measurement_type)               AS n_types,
         CAST(MIN(datetime) AS DATE)                    AS d0,
         CAST(MAX(datetime) AS DATE)                    AS d1,
         COUNT(*) FILTER (WHERE measurement_value = -99) AS n_neg99
  FROM obs")
n_smp  <- dbGetQuery(fcon, "SELECT COUNT(*) AS n FROM sample")$n
f_full <- if (has_full) dbGetQuery(fcon,
  "SELECT COUNT(*) AS n, COUNT(*) FILTER (WHERE measurement_value = -99) AS n_neg99
   FROM obs_ctd_full") else NULL

say(glue("obs           : {format(f_head$n_obs, big.mark = ',')} rows, ",
         "{f_head$n_types} measurement types, ",
         "{format(f_head$n_casts, big.mark = ',')} casts, ",
         "{f_head$n_cruises} cruises, {f_head$d0} to {f_head$d1}"))

obs           : 18,843,115 rows, 44 measurement types, 9,674 casts, 134 cruises, 1993-08-11 to 2026-07-13

Code

say(glue("sample        : {format(n_smp, big.mark = ',')} cast events"))

sample        : 19,330 cast events

Code

if (!is.null(f_full))
  say(glue("obs_ctd_full  : {format(f_full$n, big.mark = ',')} rows (supplemental)"))

obs_ctd_full  : 284,379,888 rows (supplemental)

37.1 The `-99` sentinel — fixed at the source, and now guarded

-99 is this source’s documented missing marker. The notebook had always known that — pseudoNA_values <- c(-9.99e-29, -99) strips it from longitude and latitude — but the rule was never extended to the measurement columns, so -99 flowed through ctd_measurement → ctd_thin → obs as though it were a reading. Release v2026.07.17 shipped 84,302 of them, including canonical oxygen (oxygen_ml_l_ave_sta_corr 953, oxygen_umol_kg_ave_sta_corr 4,294) alongside isus_v 40,479, ph 31,493 and spar 6,189.

The NOT isnan / isfinite guard in the pivot could never catch it: -99 is a perfectly finite double. A -99 mL/L oxygen silently corrupts any mean, minimum or anomaly a consumer computes, which is why the pivot now deletes the rows rather than reporting them — a missing value has no business existing as a row in long format.

The count below is an assertion, not a statistic: it is the regression guard for that fix, and the ctd_sentinel_neg99 rule is the same test run from the registry.

Code

say(glue("-99 in obs          : {f_head$n_neg99}   (was 84,302 in v2026.07.17)"))

-99 in obs          : 0   (was 84,302 in v2026.07.17)

Code

if (!is.null(f_full))
  say(glue("-99 in obs_ctd_full : {f_full$n_neg99}"))

-99 in obs_ctd_full : 0

Code

stopifnot(
  "the -99 sentinel is back in obs — the pivot fix has regressed" =
    f_head$n_neg99 == 0,
  "the -99 sentinel is back in obs_ctd_full" =
    is.null(f_full) || f_full$n_neg99 == 0)

37.2 `btl_*` promoted to canonical — and the depth-thinning trap

The CTD cast files ship the bottle values taken on the same cast: nutrients, chlorophyll, and the lab references (BTL_Temp, SaltB, OxB) that sensor calibration is judged against. Eleven of those twelve types were registered but flagged is_canonical = FALSE, so they were excluded from ctd_thin and therefore from obs — only btl_ammonium was canonical, which read like a stray edit rather than a decision. Making the rest canonical is what enables the bottle-vs-sensor calibration rules below; without the reference side, those rules have nothing to compare and (correctly) report skip, not pass.

Flagging them canonical was necessary but not sufficient

Thinning is designed for dense sensor scans: the ~10 m backbone and the RDP inflections are both derived from the temperature/salinity profile, so a depth is kept because the sensor profile bends there. Bottle values are not a profile — they are a handful of discrete lab samples at depths chosen by the watch, and nothing about the sensor profile knows where they are.

Selecting depths by sensor geometry and then subsetting bottle values to those depths silently discards most of them. Measured against v2026.07.17: of 133,206 distinct (cast, depth) bottle ammonium measurements, only 35,549 (26.7%) sat at a retained depth. So promoting the types alone would have landed about a quarter of the bottle record in obs and looked like a success.

ctd_thin therefore keeps a third retention reason — bottle — that retains every depth carrying a canonical bottle-grain value. The table below verifies it: within the cast direction obs covers, retention must be 100%.

Code

# bottle-grain types carry the marker as a prefix (btl_ammonium), a suffix
# (salinity_btl) or an infix (oxygen_btl_ml_l) — three positions, hence two LIKEs
btl_where <- "(measurement_type LIKE 'btl\\_%' ESCAPE '\\'
               OR measurement_type LIKE '%\\_btl%' ESCAPE '\\')"

f_btl <- dbGetQuery(fcon, glue("
  SELECT COUNT(*) FILTER (WHERE {btl_where}) AS n_btl, COUNT(*) AS n FROM obs"))
say(glue("bottle-grain rows in obs: {format(f_btl$n_btl, big.mark = ',')} of ",
         "{format(f_btl$n, big.mark = ',')} ",
         "({round(100 * f_btl$n_btl / f_btl$n, 1)}%)"))

bottle-grain rows in obs: 1,581,859 of 18,843,115 (8.4%)

Code

# retention WITHIN the direction obs keeps. The other direction is not loss: the
# source duplicates every bottle value onto both cast records (question 05), so
# comparing against the full scan set would report a 50% "loss" that is really the
# duplicate copy being dropped.
if (has_full) {
  f_ret <- dbGetQuery(fcon, glue("
    WITH chosen AS (SELECT DISTINCT sample_key FROM obs),
    f AS (SELECT DISTINCT sample_key, depth_min_m, measurement_type
          FROM obs_ctd_full
          WHERE {btl_where} AND sample_key IN (SELECT sample_key FROM chosen)),
    t AS (SELECT DISTINCT sample_key, depth_min_m, measurement_type
          FROM obs WHERE {btl_where})
    SELECT COUNT(*) AS n_avail,
           COUNT(*) FILTER (WHERE t.sample_key IS NOT NULL) AS n_kept
    FROM f LEFT JOIN t USING (sample_key, depth_min_m, measurement_type)"))
  say(glue("bottle values retained by thinning: ",
           "{format(f_ret$n_kept, big.mark = ',')} of ",
           "{format(f_ret$n_avail, big.mark = ',')} on the chosen direction ",
           "({round(100 * f_ret$n_kept / f_ret$n_avail, 1)}%)"))
  stopifnot(
    "depth thinning is dropping bottle values again (retained_reason = 'bottle')" =
      f_ret$n_kept == f_ret$n_avail)
}

bottle values retained by thinning: 1,566,483 of 1,566,483 on the chosen direction (100%)

Code

dbGetQuery(fcon, glue("
  SELECT measurement_type, COUNT(*) AS n,
         COUNT(DISTINCT sample_key) AS n_casts,
         ROUND(MIN(measurement_value), 3) AS v_min,
         ROUND(MAX(measurement_value), 3) AS v_max
  FROM obs WHERE {btl_where} GROUP BY 1 ORDER BY n DESC")) |>
  datatable(
    caption  = paste("Bottle-grain measurement types now carried in the default",
                     "obs — the reference side of every calibration check"),
    rownames = FALSE,
    options  = list(pageLength = 12, scrollX = TRUE, dom = "tip"))

37.3 One direction per physical cast: half of `sample` carries no `obs`

sample is one row per cast_key, and cast_key carries the direction suffix (…d / …u) — so a station occupation that recorded both a down- and an upcast contributes two sample rows. obs carries ctd_thin, which by design keeps one direction per physical cast. The consequence is arithmetic, but it is not stated anywhere a consumer would look: roughly half of the CTD rows in the release’s sample have no observations attached to them in the default tier.

That is not data loss — the other direction is in the supplemental obs_ctd_full — but a consumer counting casts from sample and a consumer counting casts from obs will get answers that differ by about a factor of two, and neither is wrong. Which direction should be authoritative is question calcofi_ctd-cast_06.

Code

f_dir <- dbGetQuery(fcon, "
  SELECT COUNT(*) AS n_sample,
         COUNT(*) FILTER (
           WHERE sample_key IN (SELECT DISTINCT sample_key FROM obs)) AS n_with_obs
  FROM sample")
say(glue("sample rows                 : {format(f_dir$n_sample, big.mark = ',')}"))

sample rows                 : 19,330

Code

say(glue("  with any obs (default tier): {format(f_dir$n_with_obs, big.mark = ',')} ",
         "({round(100 * f_dir$n_with_obs / f_dir$n_sample, 1)}%)"))

  with any obs (default tier): 9,674 (50%)

Code

say(glue("  no obs (other direction)   : ",
         "{format(f_dir$n_sample - f_dir$n_with_obs, big.mark = ',')}"))

  no obs (other direction)   : 9,656

37.4 Quality codes, as published

measurement_qual was a verbatim pass-through for years — carried through every table, never interpreted, with no vocabulary anywhere in the repo. There is one now: metadata/measurement_qual.csv, recovered from the CalCOFI hydrographic master’s own field documentation (6 = data OK but taken from CTD, 8 = value is suspect, 9 = missing data) and confirmed to apply to CTD, whose flags use the same codes. The codes were also being stored as "9.0" by a double→string cast, so none of them matched the vocabulary; stripping that textually (not via an INTEGER cast, which would round an unexpected "9.5") turned 3,860 spurious “unrecognized code” hits into zero.

Two things below are worth a reader’s attention, and both are open with the providers rather than acted on here:

The headline variables carry no flags at all. Source flags attach to the component sensors (Temp1Q, Salt1Q, Ox1Q, Ox2Q) while the canonical types are the averages, so the quality information is lost in the mean. A consumer cannot currently tell whether a headline CTD value is trustworthy (question calcofi_ctd-cast_09).
Every row flagged 9 = “missing data” carries a value — including physically impossible ones. Either the flag or the value is wrong, and nothing in the source says which (question calcofi_ctd-cast_13).

Code

f_qual <- dbGetQuery(fcon, "
  SELECT o.measurement_qual AS code, q.label, q.description,
         COUNT(*) AS n,
         COUNT(o.measurement_value)              AS n_with_value,
         COUNT(DISTINCT o.measurement_type)      AS n_types,
         ROUND(MIN(o.measurement_value), 3)      AS v_min,
         ROUND(MAX(o.measurement_value), 3)      AS v_max
  FROM obs o LEFT JOIN measurement_qual q ON o.measurement_qual = q.qual_code
  WHERE o.measurement_qual IS NOT NULL
  GROUP BY 1, 2, 3 ORDER BY n DESC")

n_typed <- dbGetQuery(fcon, "
  SELECT COUNT(DISTINCT measurement_type) AS n_all,
         COUNT(DISTINCT measurement_type) FILTER (
           WHERE measurement_qual IS NOT NULL) AS n_flagged FROM obs")

say(glue("flagged obs rows       : {format(sum(f_qual$n), big.mark = ',')} of ",
         "{format(f_head$n_obs, big.mark = ',')} ",
         "({round(100 * sum(f_qual$n) / f_head$n_obs, 3)}%)"))

flagged obs rows       : 23,786 of 18,843,115 (0.126%)

Code

say(glue("measurement types with any flag: {n_typed$n_flagged} of {n_typed$n_all}"))

measurement types with any flag: 24 of 44

Code

n_missing_flag <- sum(f_qual$n_with_value[f_qual$code == "9"])
if (length(n_missing_flag) && n_missing_flag > 0)
  say(glue("rows flagged 9 = 'missing data' that nonetheless carry a value: ",
           "{format(n_missing_flag, big.mark = ',')}"))

rows flagged 9 = 'missing data' that nonetheless carry a value: 8,172

Code

f_qual |>
  mutate(label = coalesce(label, "UNRECOGNIZED")) |>
  datatable(
    caption  = paste("Quality codes in the published obs, decoded against",
                     "metadata/measurement_qual.csv"),
    rownames = FALSE,
    options  = list(pageLength = 10, scrollX = TRUE, dom = "tip")) |>
  formatStyle("label",
              color = styleEqual("UNRECOGNIZED", "#d03b3b", default = "inherit"))

37.5 What the source documentation says about these numbers

Each source zip ships its own documentation (CTD-CSV-Format.pdf, CTD Data Files.pdf, FilenamingAndAscHeaderNotes.txt, and a PreliminaryDataReadme.txt on preliminary cruises). Combed 2026-08-01; the full account is in the QA/QC protocol. Four points change how a consumer should read what this ingest publishes.

The canonical types are the providers’ own recommendation, not our preference. Every property is published three times — SBE-processed, _CruiseCorr (one regression per cruise, from 4-second averages against ~1400 bottles) and _StaCorr (a regression per cast, from 1 m bin averages against that cast’s ~20 bottles). The source states station-corrected data are the best, particularly for nitrate, which is why *_sta_corr are the types flagged canonical here. Salinity is corrected separately, from bottle salinities below 350 m, applied at all depths.

A missing *_sta_corr is often correct. Station correction needs a 500 m cast with roughly ten or more bottles. On winter and spring cruises, lines north of 76.7 may be occupied with as few as 0–12 bottles and receive cruise-corrected values only. Counting nulls as defects will mis-read a whole class of casts.

Values are 1 m bin averages, not raw scans. obs_ctd_full is full published resolution. That is also why the loop-edit rule’s 1 m threshold is a noise floor rather than a tolerance.

data_stage now reaches the release

The source is explicit that preliminary data are for non-publication use, and specific about what moves: “Temperature and salinities may change very little but oxygen, nitrate and chlorophyll data may change significantly after post-cruise calibrations.”

Until v2026.08 the released sample carried no such column, so a consumer could not tell a final cruise from a preliminary one and the caveat above could not reach them — question calcofi_ctd-cast_14. Core sample now has data_stage, supplied by this ingest as the optional trailing column of append_sample() (calcofi4db 3.4.0). The other 15 ingests, none of which draw the distinction, leave it NULL rather than having to change.

It carries three values, not two, matching the provider’s own published progression — final ← preliminary_with_bottle ← preliminary_without_bottle. The split matters most for exactly the quantities the caveat names: oxygen and chlorophyll are bottle-corrected, so a preliminary_without_bottle cruise has not merely provisional values for them but none at all. Collapsing both preliminary tiers into one label made that indistinguishable from a cruise whose sensors failed.

Code

# how much of what we publish is preliminary? Read from the wrangling DB when the
# heavy path ran; otherwise say so rather than guessing.
if ("ctd_cast" %in% DBI::dbListTables(con)) {
  dbGetQuery(con, "
    SELECT data_stage, COUNT(DISTINCT cruise_key) AS n_cruises,
           COUNT(DISTINCT cast_key) AS n_casts
    FROM ctd_cast GROUP BY 1") |>
    # ordered by processing stage, not by size: the point of the table is where a
    # cruise sits in the progression
    arrange(match(data_stage, CTD_DATA_STAGES)) |>
    mutate(stage_label = CTD_STAGE_LABELS[data_stage], .after = data_stage) |>
    datatable(
      caption = paste("Processing stage in this ingest — carried into the",
                      "release on `sample.data_stage`"),
      rownames = FALSE, options = list(dom = "t"))
} else {
  say(paste("data_stage not available on a fast re-render (ctd_cast is not",
            "materialized); it is charted in the Gantt appendix on a full run."))
}

37.6 QC rules over this ingest’s output

The rules are data, not code: metadata/qc_rules/rules.csv plus one sql/*.sql file per rule, versioned with the pipeline that produces the data they check, so a data manager can review one in a diff. They are executed here by calcofi4db::qc_run_all() — the same engine apps/ctd-qaqc uses, moved into the package once it had two callers so the app and this notebook cannot report different numbers for the same rule.

Every rule targets obs / sample, never a per-dataset table. That is what lets the identical registry run against a release, against this parquet, and (next) against an uploaded file.

skip is not pass

A rule whose input is absent returns zero rows, which is indistinguishable from “the data is clean”. The engine therefore reports skip with a reason, never a silent pass — this is not hypothetical: the three bottle-vs-sensor rules did exactly that against v2026.07.30, which predates the btl_* promotion above. Rules scoped to a single cruise (they read the full-resolution obs_ctd_full) skip here too, and are run separately below.

Code

dir_rules <- here("metadata/qc_rules")
rules_all <- qc_read_rules(dir_rules, active_only = FALSE)
rules_act <- filter(rules_all, active)
rules_now <- filter(rules_act, coalesce(scope, "all") == "all")

f_types <- qc_present_types(fcon, "calcofi_ctd-cast")
qc_res  <- qc_run_all(fcon, rules_now, limit = 500L, present_types = f_types)
qc_sum  <- qc_summarize(qc_res, rules_now)

say(glue("rules in registry : {nrow(rules_all)} ",
         "({nrow(rules_act)} active, {nrow(rules_all) - nrow(rules_act)} parked)"))

rules in registry : 22 (19 active, 3 parked)

Code

say(glue("run here          : {nrow(rules_now)} release-wide ",
         "(+{nrow(rules_act) - nrow(rules_now)} cruise-scoped, below)"))

run here          : 14 release-wide (+5 cruise-scoped, below)

Code

for (s in c("FAIL", "ERROR", "flag", "skip", "pass")) {
  n <- sum(qc_sum$status == s)
  if (n > 0) say(glue("  {format(s, width = 6)}: {n}"))
}

  flag  : 9
  pass  : 5

Code

say(glue("total findings    : ",
         "{format(sum(qc_sum$n, na.rm = TRUE), big.mark = ',')}"))

total findings    : 29,388

Code

qc_sum |>
  transmute(rule_key, status, findings = n, severity, type = rule_type, target,
            description, `sec` = elapsed_s, note) |>
  datatable(
    caption  = paste("QC rule results over this ingest's parquet —",
                     "registry metadata/qc_rules/, engine calcofi4db::qc_run_all()"),
    rownames = FALSE,
    options  = list(pageLength = 15, scrollX = TRUE, dom = "tip")) |>
  formatStyle(
    "status",
    color = styleEqual(c("pass", "flag", "FAIL", "ERROR", "skip"),
                       c("#1a7f37", "#b06000", "#d03b3b", "#d03b3b", "#6e7781")),
    fontWeight = styleEqual(c("FAIL", "ERROR"), c("bold", "bold"),
                            default = "normal"))

37.6.1 Values outside their declared range

Ranges live in metadata/measurement_type.csv (valid_min / valid_max), moved there from an inline tribble so they are reviewable in a diff, emittable by the CF netCDF writer as real variable attributes, and usable by any rule rather than only by this notebook. They are deliberately generous: they catch values that cannot be real, they do not police oceanography. Nothing below is dropped.

Code

f_range <- dbGetQuery(fcon, "
  SELECT o.measurement_type,
         COUNT(*)                                    AS n_bad,
         COUNT(DISTINCT o.cruise_key)                AS n_cruises,
         ROUND(MIN(o.measurement_value), 3)          AS v_min,
         ROUND(MAX(o.measurement_value), 3)          AS v_max,
         m.valid_min, m.valid_max
  FROM obs o JOIN measurement_type m USING (measurement_type)
  WHERE m.valid_min IS NOT NULL
    AND (o.measurement_value < m.valid_min OR o.measurement_value > m.valid_max)
  GROUP BY 1, 6, 7 ORDER BY n_bad DESC")

f_range |>
  left_join(
    dbGetQuery(fcon, "SELECT measurement_type, COUNT(*) AS n_total FROM obs GROUP BY 1"),
    by = "measurement_type") |>
  transmute(measurement_type, n_bad, n_total,
            `% bad` = round(100 * n_bad / n_total, 2),
            n_cruises, v_min, v_max,
            range = glue("{valid_min}..{valid_max}")) |>
  datatable(
    caption  = "Published values outside their declared plausible range",
    rownames = FALSE,
    options  = list(pageLength = 15, scrollX = TRUE, dom = "tip"))

37.6.2 Impossible temperatures

The single finding most worth a provider’s time. temperature_ave reaches 60.4 °C in the published product, and this is not a sensor warming on deck: most of the values above the 40 °C ceiling are at depth, they span more than two decades and more than one ship, and where a bottle thermometer sampled the same cast at the same depth it disagrees by 30–40 °C. Three independent lines — the range rule, the depth distribution, and the bottle reference — say the same thing, which is why this is stated as a finding rather than a suspicion.

Code

f_t <- dbGetQuery(fcon, "
  SELECT COUNT(*) FILTER (WHERE measurement_value > 40)                       AS n_hot,
         COUNT(*) FILTER (WHERE measurement_value > 40 AND depth_min_m > 2)   AS n_hot_deep,
         COUNT(*) FILTER (WHERE measurement_value < -2)                       AS n_cold,
         COUNT(DISTINCT cruise_key) FILTER (
           WHERE measurement_value > 40 OR measurement_value < -2)            AS n_cruises
  FROM obs WHERE measurement_type = 'temperature_ave'")
say(glue("temperature_ave > 40 degC : {f_t$n_hot} ",
         "({f_t$n_hot_deep} of them deeper than 2 m — not a deck warm-up)"))

temperature_ave > 40 degC : 0 (0 of them deeper than 2 m — not a deck warm-up)

Code

say(glue("temperature_ave < -2 degC : {f_t$n_cold}"))

temperature_ave < -2 degC : 0

Code

say(glue("cruises affected          : {f_t$n_cruises}"))

cruises affected          : 0

Code

dbGetQuery(fcon, "
  SELECT cruise_key,
         COUNT(*)                             AS n,
         COUNT(DISTINCT sample_key)           AS n_casts,
         ROUND(MIN(depth_min_m), 1)           AS depth_min,
         ROUND(MAX(depth_min_m), 1)           AS depth_max,
         ROUND(MIN(measurement_value), 2)     AS t_min,
         ROUND(MAX(measurement_value), 2)     AS t_max
  FROM obs
  WHERE measurement_type = 'temperature_ave'
    AND (measurement_value > 40 OR measurement_value < -2)
  GROUP BY 1 ORDER BY n DESC") |>
  datatable(
    caption  = paste("Physically impossible temperature_ave by cruise —",
                     "the table to send with question calcofi_ctd-cast_02"),
    rownames = FALSE,
    options  = list(pageLength = 15, scrollX = TRUE, dom = "tip"))

Code

# the corroboration: where a bottle thermometer sampled the same cast at the same
# depth, what did IT read? This is only possible because btl_temperature is now
# canonical (above) — before that the reference side was not in obs at all.
f_tbtl <- dbGetQuery(fcon, "
  WITH hot AS (
    SELECT sample_key, cruise_key, depth_min_m, measurement_value AS v_sensor
    FROM obs
    WHERE measurement_type = 'temperature_ave' AND measurement_value > 40)
  SELECT h.cruise_key,
         COUNT(*)                              AS n_pairs,
         ROUND(MIN(h.depth_min_m), 1)          AS depth_min,
         ROUND(MAX(h.depth_min_m), 1)          AS depth_max,
         ROUND(MIN(b.measurement_value), 2)    AS bottle_min,
         ROUND(MAX(b.measurement_value), 2)    AS bottle_max,
         ROUND(MIN(h.v_sensor), 2)             AS sensor_min,
         ROUND(MAX(h.v_sensor), 2)             AS sensor_max
  FROM hot h
  JOIN obs b ON b.sample_key = h.sample_key AND b.depth_min_m = h.depth_min_m
            AND b.measurement_type = 'btl_temperature'
  GROUP BY 1 ORDER BY n_pairs DESC")

if (nrow(f_tbtl)) {
  datatable(
    f_tbtl,
    caption  = paste("Bottle thermometer vs sensor on the same cast and depth,",
                     "where the sensor reads above 40 degC"),
    rownames = FALSE,
    options  = list(pageLength = 10, scrollX = TRUE, dom = "tip"))
} else {
  say("no bottle temperature at the same cast+depth as any >40 degC sensor value")
}

no bottle temperature at the same cast+depth as any >40 degC sensor value

37.6.3 Negative pH, clustered by cruise

Not previously reported, and the largest single range violation in the release: ph accounts for the great majority of out-of-range values, and the bad values are not scattered — they cluster tightly within a cruise, often thousands of rows sitting in a narrow band around −2.5 to −3. Scattered impossible values look like bad scans; a whole cruise pinned at −2.7 looks like a column that was never converted from sensor voltage to pH units. The dataset has two explicit voltage columns (isus_v, fluorescence_v) but no ph_v twin, so a raw-voltage ph would have nowhere else to go. That is a hypothesis for the providers, not a conclusion (question calcofi_ctd-cast_12).

Code

dbGetQuery(fcon, "
  SELECT cruise_key,
         COUNT(*)                                                        AS n_bad,
         ROUND(MIN(measurement_value), 2)                                AS v_min,
         ROUND(MAX(measurement_value), 2)                                AS v_max,
         COUNT(*) FILTER (WHERE measurement_value < 0)                   AS n_negative,
         ROUND(MAX(measurement_value) - MIN(measurement_value), 2)       AS spread
  FROM obs
  WHERE measurement_type = 'ph' AND (measurement_value < 6 OR measurement_value > 9)
  GROUP BY 1 ORDER BY n_bad DESC") |>
  datatable(
    caption  = paste("Out-of-range pH by cruise — a narrow `spread` on thousands",
                     "of rows is the signature of an unconverted column"),
    rownames = FALSE,
    options  = list(pageLength = 10, scrollX = TRUE, dom = "tip"))

37.6.4 Full-resolution profile rules, one cruise at a time

Spikes, loop edits and up/down disagreement are invisible in obs: thinning keeps one direction and roughly one sample per 10 m, which is exactly the structure these rules examine. They read the supplemental obs_ctd_full instead, which is hive-partitioned by cruise_key so a scoped query prunes to one cruise rather than scanning 212M rows.

The cruise below is chosen by evidence, not by recency: it is the one with the most out-of-range canonical values, so the profile rules run where there is most reason to look. This is a sample, not a sweep — the remaining cruises are the app’s job.

Code

rules_cruise <- filter(rules_act, coalesce(scope, "all") == "cruise")

if (has_full && nrow(rules_cruise)) {
  worst_cruise <- dbGetQuery(fcon, "
    SELECT o.cruise_key, COUNT(*) AS n
    FROM obs o JOIN measurement_type m USING (measurement_type)
    WHERE m.valid_min IS NOT NULL
      AND (o.measurement_value < m.valid_min OR o.measurement_value > m.valid_max)
    GROUP BY 1 ORDER BY n DESC LIMIT 1")$cruise_key

  # Since the two-sensor average repair and the bounds guard landed, there is
  # routinely NO out-of-range cruise — the query returns zero rows and
  # `$cruise_key` is character(0). That is the good outcome, and it used to
  # abort the render: character(0) reached qc_run_rule()'s scope guard, where
  # `&&` on a zero-length operand yields NA and `if (NA)` stops the chunk.
  #
  # "Nothing is out of range" is not a reason to skip the profile rules. They
  # test spike, gradient, density inversion and stuck-sensor shapes, none of
  # which a range check can see — so fall back to the cruise with the most
  # full-resolution scans, i.e. the most to look at.
  scope_reason <- "most out-of-range values"
  if (!length(worst_cruise) || is.na(worst_cruise[1])) {
    worst_cruise <- dbGetQuery(fcon, "
      SELECT cruise_key, COUNT(*) AS n FROM obs_ctd_full
      GROUP BY 1 ORDER BY n DESC LIMIT 1")$cruise_key
    scope_reason <- "no cruise has out-of-range values; most full-resolution scans"
  }
  say(glue("cruise scoped to: {worst_cruise} ({scope_reason})"))

  qc_res_cr <- qc_run_all(
    fcon, rules_cruise, limit = 500L, present_types = f_types,
    scope_values = list(cruise_key = worst_cruise))
  qc_sum_cr <- qc_summarize(qc_res_cr, rules_cruise)
} else {
  worst_cruise <- NA_character_
  qc_sum_cr    <- NULL
  say("obs_ctd_full not present — the profile rules are SKIPPED, not passed")
}

cruise scoped to: 2022-04-3322 (no cruise has out-of-range values; most full-resolution scans)

Code

if (!is.null(qc_sum_cr)) {
  qc_sum_cr |>
    transmute(rule_key, status, findings = n, severity, type = rule_type,
              description, `sec` = elapsed_s, note) |>
    datatable(
      caption  = glue("Profile rules over obs_ctd_full for cruise {worst_cruise}"),
      rownames = FALSE,
      options  = list(pageLength = 10, scrollX = TRUE, dom = "t")) |>
    formatStyle(
      "status",
      color = styleEqual(c("pass", "flag", "FAIL", "ERROR", "skip"),
                         c("#1a7f37", "#b06000", "#d03b3b", "#d03b3b", "#6e7781")))
} else {
  say("no cruise-scoped rule results this render")
}

37.7 Questions for Data Providers

Open questions on this dataset, tracked in metadata/calcofi/ctd-cast/questions.csv and surfaced here so they travel with the workflow rather than living in someone’s inbox. Everything above that could not be resolved from the data alone is written down here, with the measured evidence that raised it, so that asking is a matter of sending a link rather than reconstructing an investigation.

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 CTD data providers (ranked)")

Code

close_duckdb(fcon)

38 Cleanup

Code

close_duckdb(con)

# remove checkpoint after successful completion (fresh next run)
if (file_exists(db_checkpoint)) {
  file_delete(db_checkpoint)
  say(glue("Removed checkpoint: {db_checkpoint}"))
}

Removed checkpoint: /Users/bbest/Github/CalCOFI/workflows-ctd-new-format/data/wrangling/calcofi_ctd-cast_checkpoint.duckdb

Code

# the fingerprint is recorded in [record_fingerprint], immediately after the
# outputs are complete and before the GCS upload — see the reasoning there

39 Appendix: Interactive Gantt Chart

This appendix cannot fail the ingest

error: true is deliberate. Everything above this point — the parquet outputs and the GCS sync — is the ingest’s actual product, and it is all complete by the time this renders. On 2026-07-30 a stale path in this cosmetic chart aborted the whole target after a successful 15-minute run, so targets recorded a failure and release_database never started, over a chart. A decorative appendix should degrade to a visible error message, not discard finished work.

Code

# On a fast re-render (inputs unchanged) the heavy path never ran, so ctd_cast was
# never built in the wrangling DB. Gate the chart on its presence rather than
# letting a cosmetic appendix throw — the ingest's real product is already done.
con_pq <- get_duckdb_con(db_path)
load_duckdb_extension(con_pq, "spatial")
has_ctd_cast <- "ctd_cast" %in% DBI::dbListTables(con_pq)
if (!has_ctd_cast) {
  close_duckdb(con_pq)
  say(paste0(
    "ctd_cast not materialized (inputs unchanged, heavy path skipped) — no Gantt ",
    "chart this render. Re-run with CTD_FORCE_REBUILD=TRUE to rebuild it."))
}

Code

# Read ctd_cast from the WRANGLING DB, not from parquet.
#
# This chunk used to read '{dir_parquet}/ctd_cast.parquet'. That file stopped
# being written at the consolidated-core cut-over — Write Parquet now exports
# sample / obs / obs_ctd_full / measurement_type, with ctd_cast, ctd_thin,
# ctd_measurement and ctd_summary superseded. The breakage went unnoticed only
# because this target was excluded from `_targets.R` for the whole period since,
# so re-enabling it surfaced a latent failure: the ingest completed all of its
# real work (parquet written, GCS synced) and then died in a cosmetic appendix.
#
# sample.parquet is the natural modern substitute but has no `data_stage`, and
# the final-vs-preliminary split is the point of the chart. ctd_cast is still a
# real table in the wrangling DB at render time, so read it there — the connection
# is opened by the guard chunk above.
cruise_spans <- tbl(con_pq, "ctd_cast") |>
  group_by(cruise_key, data_stage) |>
  summarize(
    beg_date = min(datetime_start_utc, na.rm = TRUE) |> as.Date(),
    end_date = max(datetime_start_utc, na.rm = TRUE) |> as.Date(),
    .groups = "drop"
  ) |>
  collect() |>
  filter(!is.na(beg_date), !is.na(end_date))

# split multi-year cruises into separate rows
cruise_spans <- cruise_spans |>
  rowwise() |>
  mutate(
    years_spanned = list(year(beg_date):year(end_date))
  ) |>
  unnest(years_spanned) |>
  mutate(
    year = years_spanned,
    adj_beg_date = if_else(
      years_spanned == year(beg_date),
      beg_date,
      as.Date(paste0(years_spanned, "-01-01"))
    ),
    adj_end_date = if_else(
      years_spanned == year(end_date),
      end_date,
      as.Date(paste0(years_spanned, "-12-31"))
    ),
    begin_jday = yday(adj_beg_date),
    end_jday = yday(adj_end_date),
    stage_label = data_stage,
    hover_text = glue(
      "Cruise: {cruise_key}<br>",
      "Stage: {data_stage}<br>",
      "Begin: {format(beg_date, '%Y-%m-%d')}<br>",
      "End: {format(end_date, '%Y-%m-%d')}<br>",
      "Duration: {as.numeric(difftime(end_date, beg_date, units = 'days'))} days"
    )
  ) |>
  ungroup() |>
  arrange(adj_beg_date)

# green = done, warm = in progress; the two preliminary tiers share a hue so the
# eye reads "not final yet" first and the tier second
colors <- c(
  final                  = "#23d355ff",
  preliminary_with_bottle = "#A23B72",
  preliminary_without_bottle        = "#E8862E")

p <- plot_ly()

for (i in 1:nrow(cruise_spans)) {
  row <- cruise_spans[i, ]
  p <- p |>
    add_trace(
      type = "scatter",
      mode = "lines",
      x = c(row$begin_jday, row$end_jday),
      y = c(row$year, row$year),
      line = list(
        color = colors[row$stage_label],
        width = 8
      ),
      text = row$hover_text,
      hoverinfo = "text",
      showlegend = FALSE,
      legendgroup = row$stage_label
    )
}

# one legend entry per stage actually present, in processing order — hand-written
# traces drifted from `colors` when the third tier was added
for (stg in CTD_DATA_STAGES) {
  if (!stg %in% cruise_spans$stage_label) next
  p <- p |>
    add_trace(
      type = "scatter",
      mode = "lines",
      x = c(NA, NA),
      y = c(NA, NA),
      line = list(color = colors[[stg]], width = 8),
      name = CTD_STAGE_LABELS[[stg]],
      legendgroup = stg,
      showlegend = TRUE
    )
}

p <- p |>
  layout(
    title = "CalCOFI CTD Cruise Timeline",
    xaxis = list(
      title = "Day of Year",
      range = c(1, 365),
      dtick = 30,
      tickangle = 0
    ),
    yaxis = list(
      title = "Year",
      autorange = "reversed",
      dtick = 1
    ),
    hovermode = "closest",
    plot_bgcolor = "#f8f9fa",
    paper_bgcolor = "white",
    legend = list(
      x = 1.02,
      y = 1,
      xanchor = "left",
      yanchor = "top"
    )
  )

p

Code

close_duckdb(con_pq)

Session Info

Code

devtools::session_info()

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.6.1 (2026-06-24)
 os       macOS Sequoia 15.7.1
 system   aarch64, darwin23
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Europe/Rome
 date     2026-09-24
 pandoc   3.8.3 @ /opt/homebrew/bin/ (via rmarkdown)
 quarto   1.8.25 @ /usr/local/bin/quarto

─ Packages ───────────────────────────────────────────────────────────────────
 ! package            * version    date (UTC) lib source
   abind                1.4-8      2024-09-12 [1] CRAN (R 4.6.0)
   arrow                25.0.1     2026-08-23 [1] CRAN (R 4.6.1)
   assertthat           0.2.1      2019-03-21 [1] CRAN (R 4.6.0)
   backports            1.5.1      2026-04-03 [1] CRAN (R 4.6.0)
   base64enc            0.1-6      2026-02-02 [1] CRAN (R 4.6.0)
   bit                  4.6.0      2025-03-06 [1] CRAN (R 4.6.0)
   bit64                4.8.6      2026-09-01 [1] CRAN (R 4.6.1)
   blob                 1.3.0      2026-01-14 [1] CRAN (R 4.6.0)
   brio                 1.1.5      2024-04-24 [1] CRAN (R 4.6.0)
   broom                1.0.13     2026-05-14 [1] CRAN (R 4.6.0)
   bslib                0.12.0     2026-08-04 [1] CRAN (R 4.6.1)
   cachem               1.1.0      2024-05-16 [1] CRAN (R 4.6.0)
 P calcofi4db         * 4.16.1     2026-09-23 [?] load_all()
 P calcofi4r          * 1.24.2     2026-09-10 [?] load_all()
   chromote             0.5.1      2025-04-24 [1] CRAN (R 4.6.0)
   class                7.3-24     2026-08-03 [1] CRAN (R 4.6.1)
   classInt             0.4-11     2025-01-08 [1] CRAN (R 4.6.0)
   cli                  3.6.6      2026-04-09 [1] CRAN (R 4.6.0)
   codetools            0.2-20     2024-03-31 [1] CRAN (R 4.6.1)
   crayon               1.5.3      2024-06-20 [1] CRAN (R 4.6.0)
   crosstalk            1.2.2      2025-08-26 [1] CRAN (R 4.6.0)
   curl                 8.0.0      2026-08-25 [1] CRAN (R 4.6.1)
   data.table           1.18.6.1   2026-08-24 [1] CRAN (R 4.6.1)
   DBI                * 1.3.0      2026-02-25 [1] CRAN (R 4.6.0)
   dbplyr               2.6.0      2026-06-17 [1] CRAN (R 4.6.0)
   desc                 1.4.3      2023-12-10 [1] CRAN (R 4.6.0)
   devtools             2.5.2      2026-04-30 [1] CRAN (R 4.6.0)
   DiagrammeR           1.0.12     2026-04-27 [1] CRAN (R 4.6.0)
   DiagrammeRsvg        0.1        2016-02-04 [1] CRAN (R 4.6.0)
   digest               0.6.39     2025-11-19 [1] CRAN (R 4.6.0)
   dm                   1.1.2      2026-05-17 [1] CRAN (R 4.6.0)
   dplyr              * 1.2.1      2026-04-03 [1] CRAN (R 4.6.0)
   DT                 * 0.34.0     2025-09-02 [1] CRAN (R 4.6.0)
   duckdb               1.5.5      2026-07-25 [1] CRAN (R 4.6.1)
   dygraphs             1.1.1.6    2018-07-11 [1] CRAN (R 4.6.0)
   e1071                1.7-17     2025-12-18 [1] CRAN (R 4.6.0)
   ellipsis             0.3.3      2026-04-04 [1] CRAN (R 4.6.0)
   evaluate             1.0.5      2025-08-27 [1] CRAN (R 4.6.0)
   farver               2.1.2      2024-05-13 [1] CRAN (R 4.6.0)
   fastmap              1.2.0      2024-05-15 [1] CRAN (R 4.6.0)
   fs                 * 2.1.0      2026-04-18 [1] CRAN (R 4.6.0)
   fuzzyjoin            0.1.8      2026-02-20 [1] CRAN (R 4.6.0)
   gargle               1.6.1      2026-01-29 [1] CRAN (R 4.6.0)
   generics             0.1.4      2025-05-09 [1] CRAN (R 4.6.0)
   geojsonsf            2.0.5      2025-11-26 [1] CRAN (R 4.6.0)
   ggplot2            * 4.0.3      2026-04-22 [1] CRAN (R 4.6.0)
   glue               * 1.8.1      2026-04-17 [1] CRAN (R 4.6.0)
   googledrive          2.1.2      2025-09-10 [1] CRAN (R 4.6.0)
   gsw                  1.2-0      2024-08-19 [1] CRAN (R 4.6.0)
   gtable               0.3.6      2024-10-25 [1] CRAN (R 4.6.0)
   here               * 1.0.2      2025-09-15 [1] CRAN (R 4.6.0)
   highcharter          0.9.5      2026-04-22 [1] CRAN (R 4.6.0)
   hms                  1.1.4      2025-10-17 [1] CRAN (R 4.6.0)
   htmltools            0.5.9      2025-12-04 [1] CRAN (R 4.6.0)
   htmlwidgets          1.6.4      2023-12-06 [1] CRAN (R 4.6.0)
   httpuv               1.6.17     2026-03-18 [1] CRAN (R 4.6.0)
   httr               * 1.4.9      2026-09-01 [1] CRAN (R 4.6.1)
   httr2                1.3.0      2026-07-13 [1] CRAN (R 4.6.1)
   igraph               2.3.3      2026-06-26 [1] CRAN (R 4.6.1)
   isoband              0.3.0      2025-12-07 [1] CRAN (R 4.6.0)
   janitor            * 2.2.1      2024-12-22 [1] CRAN (R 4.6.0)
   jquerylib            0.1.4      2021-04-26 [1] CRAN (R 4.6.0)
   jsonlite             2.0.0      2025-03-27 [1] CRAN (R 4.6.0)
   KernSmooth           2.23-27    2026-08-12 [1] CRAN (R 4.6.1)
   knitr                1.52       2026-09-06 [1] CRAN (R 4.6.1)
   labeling             0.4.3      2023-08-29 [1] CRAN (R 4.6.0)
   later                1.4.8      2026-03-05 [1] CRAN (R 4.6.0)
   lattice              0.23-1     2026-08-12 [1] CRAN (R 4.6.1)
   leafem               0.2.5      2025-08-28 [1] CRAN (R 4.6.0)
   leaflet              2.2.3      2025-09-04 [1] CRAN (R 4.6.0)
   librarian            1.8.1      2021-07-12 [1] CRAN (R 4.6.0)
   lifecycle            1.0.5      2026-01-08 [1] CRAN (R 4.6.0)
   listviewer           4.0.0      2023-09-30 [1] CRAN (R 4.6.0)
   lubridate          * 1.9.5      2026-02-04 [1] CRAN (R 4.6.0)
   magrittr             2.0.5      2026-04-04 [1] CRAN (R 4.6.0)
   mapgl                0.5.0.9000 2026-09-01 [1] Github (bbest/mapgl@484e869)
   mapview            * 2.11.4     2025-09-08 [1] CRAN (R 4.6.0)
   markdown             2.0        2025-03-23 [1] CRAN (R 4.6.0)
   Matrix               1.7-6      2026-07-25 [1] CRAN (R 4.6.1)
   memoise              2.0.1      2021-11-26 [1] CRAN (R 4.6.0)
   mgcv                 1.9-4      2025-11-07 [1] CRAN (R 4.6.1)
   mime                 0.13       2025-03-17 [1] CRAN (R 4.6.0)
   nlme                 3.1-171    2026-09-01 [1] CRAN (R 4.6.1)
   otel                 0.2.0      2025-08-29 [1] CRAN (R 4.6.0)
   pillar               1.11.1     2025-09-17 [1] CRAN (R 4.6.0)
   pkgbuild             1.4.8      2025-05-26 [1] CRAN (R 4.6.0)
   pkgconfig            2.0.3      2019-09-22 [1] CRAN (R 4.6.0)
   pkgload              1.5.3      2026-06-15 [1] CRAN (R 4.6.0)
   plotly             * 4.12.1     2026-07-22 [1] CRAN (R 4.6.1)
   png                  0.1-9      2026-03-15 [1] CRAN (R 4.6.0)
   processx             3.9.0      2026-04-22 [1] CRAN (R 4.6.0)
   promises             1.5.0      2025-11-01 [1] CRAN (R 4.6.0)
   proxy                0.4-29     2025-12-29 [1] CRAN (R 4.6.0)
   ps                 * 1.9.3      2026-04-20 [1] CRAN (R 4.6.0)
   purrr              * 1.2.2      2026-04-10 [1] CRAN (R 4.6.0)
   quantmod             0.4.29     2026-06-28 [1] CRAN (R 4.6.1)
   R6                   2.6.1      2025-02-15 [1] CRAN (R 4.6.0)
   rappdirs             0.3.4      2026-01-17 [1] CRAN (R 4.6.0)
   raster               3.6-32     2025-03-28 [1] CRAN (R 4.6.0)
   RColorBrewer         1.1-3      2022-04-03 [1] CRAN (R 4.6.0)
   Rcpp                 1.1.2      2026-07-05 [1] CRAN (R 4.6.1)
   readr              * 2.2.0      2026-02-19 [1] CRAN (R 4.6.0)
   rlang                1.3.0      2026-07-05 [1] CRAN (R 4.6.1)
   rlist                0.4.6.2    2021-09-03 [1] CRAN (R 4.6.0)
   rmarkdown            2.32       2026-09-01 [1] CRAN (R 4.6.1)
   rnaturalearth        1.2.0      2026-01-19 [1] CRAN (R 4.6.0)
   rnaturalearthhires   1.0.0.9000 2026-08-29 [1] Github (ropensci/rnaturalearthhires@e4736f6)
   RPostgres            1.4.10     2026-02-16 [1] CRAN (R 4.6.0)
   rprojroot            2.1.1      2025-08-26 [1] CRAN (R 4.6.0)
   rstudioapi           0.19.0     2026-06-11 [1] CRAN (R 4.6.0)
   rvest              * 1.0.5      2025-08-29 [1] CRAN (R 4.6.0)
   S7                   0.2.2      2026-04-22 [1] CRAN (R 4.6.0)
   sass                 0.4.10     2025-04-11 [1] CRAN (R 4.6.0)
   satellite            1.0.6      2025-08-21 [1] CRAN (R 4.6.0)
   scales               1.4.0      2025-04-24 [1] CRAN (R 4.6.0)
   selectr              0.6-0      2026-06-23 [1] CRAN (R 4.6.1)
   sessioninfo          1.2.4      2026-06-04 [1] CRAN (R 4.6.0)
   sf                 * 1.1-2      2026-07-23 [1] CRAN (R 4.6.1)
   shiny                1.14.0     2026-06-21 [1] CRAN (R 4.6.0)
   shinyWidgets         0.9.1      2026-03-09 [1] CRAN (R 4.6.0)
   snakecase            0.11.1     2023-08-27 [1] CRAN (R 4.6.0)
   sp                   2.2-3      2026-07-19 [1] CRAN (R 4.6.1)
   stars                0.7-3      2026-07-20 [1] CRAN (R 4.6.1)
   stringi              1.8.9      2026-08-04 [1] CRAN (R 4.6.1)
   stringr            * 1.6.0      2025-11-04 [1] CRAN (R 4.6.0)
   terra                1.9-46     2026-08-22 [1] CRAN (R 4.6.1)
   testthat           * 3.3.2      2026-01-11 [1] CRAN (R 4.6.0)
   tibble             * 3.3.1      2026-01-11 [1] CRAN (R 4.6.0)
   tidyr              * 1.3.2      2025-12-19 [1] CRAN (R 4.6.0)
   tidyselect           1.2.1      2024-03-11 [1] CRAN (R 4.6.0)
   timechange           0.4.0      2026-01-29 [1] CRAN (R 4.6.0)
   TTR                  0.24.4     2023-11-28 [1] CRAN (R 4.6.0)
   tzdb                 0.5.0      2025-03-15 [1] CRAN (R 4.6.0)
   units                1.0-1      2026-03-11 [1] CRAN (R 4.6.0)
   usethis              3.2.1      2025-09-06 [1] CRAN (R 4.6.0)
   utf8                 1.2.6      2025-06-08 [1] CRAN (R 4.6.0)
   uuid                 1.2-2      2026-01-23 [1] CRAN (R 4.6.0)
   V8                   8.2.0      2026-04-21 [1] CRAN (R 4.6.0)
   vctrs                0.7.3      2026-04-11 [1] CRAN (R 4.6.0)
   viridisLite          0.4.3      2026-02-04 [1] CRAN (R 4.6.0)
   visNetwork           2.1.4      2025-09-04 [1] CRAN (R 4.6.0)
   vroom                1.7.1      2026-03-31 [1] CRAN (R 4.6.0)
   websocket            1.4.4      2025-04-10 [1] CRAN (R 4.6.0)
   withr                3.0.3      2026-06-19 [1] CRAN (R 4.6.0)
   xfun                 0.60       2026-07-09 [1] CRAN (R 4.6.1)
   xml2                 1.6.0      2026-06-22 [1] CRAN (R 4.6.1)
   xtable               1.8-8      2026-02-22 [1] CRAN (R 4.6.0)
   xts                  0.14.2     2026-02-28 [1] CRAN (R 4.6.0)
   yaml                 2.3.12     2025-12-10 [1] CRAN (R 4.6.0)
   zip                * 3.0.2      2026-08-04 [1] CRAN (R 4.6.1)
   zoo                  1.9-0      2026-07-31 [1] CRAN (R 4.6.1)

 [1] /Library/Frameworks/R.framework/Versions/4.6/Resources/library

 * ── Packages attached to the search path.
 P ── Loaded and on-disk path mismatch.

──────────────────────────────────────────────────────────────────────────────

1 Overview

1.1 Key Features

2 Setup

2.1 Archives on Drive, contents on local scratch

3 Prime Downloads from the Object Store

4 Build the Zip Inventory

5 Check for Resumable State

6 Choose Which Archives to Read

7 Download and Unzip Files

7.1 Overlap guard

8 Find and Prioritize CTD Data Files

9 Read and Standardize CTD Files

10 Detect and Correct Column Type Mismatches

11 Format Date-Time Column

12 Pseudo-NA Values and Fill Missing Coordinates

13 Distance Filtering

13.1 View sample cruise with excess distance points

13.2 Filter points

13.3 View filtered points

13.4 Check for duplicates by datetime and depth

14 Rename Fields

15 De-duplicate ctd_raw

16 Cross-Dataset Bridge

17 Save Checkpoint

18 Measurement Column Registry

19 Split into Tidy Tables

19.1 ctd_cast — cast-level metadata

19.2 ctd_measurement — long-format sensor readings

20 Apply Accepted QC Flags

21 Apply Flag Overrides

22 Combine the Sensor Pairs by Their Flags

22.1 ctd_summary — summary stats across cast directions

22.2 measurement_type — reference table

22.3 Measurement provenance — instrument and accuracy eras

22.4 Drop intermediate table

23 Measurement Frequency

23.1 Multiplier 1 — sub-meter depth resolution

23.2 Multiplier 2 — up- and down-casts both stored

23.3 Multiplier 3 — redundant measurement-type variants

24 CTD Thin

24.1 Verify cast pairing

24.2 Choose a single cast direction

24.3 10 m depth grid

24.4 Inflection points (Ramer–Douglas–Peucker)

24.5 Assemble ctd_thin

24.6 Verify ctd_thin

24.7 Spot-check thinned profiles

25 Schema Diagram

26 Add Spatial

27 Derived View

28 Load Dataset Metadata

29 Data Quality Diagnostics

29.1 Quality flags — decoded against the controlled vocabulary

30 Validate and Enforce Types

31 Preview Tables

32 Emit Core Tables

32.1 No cruise may leave obs and keep its casts

33 Write Parquet

34 Write Metadata JSON

35 Export Modified Dependency Deltas

36 Upload to GCS

37 Findings

37.1 The -99 sentinel — fixed at the source, and now guarded

37.2 btl_* promoted to canonical — and the depth-thinning trap

37.3 One direction per physical cast: half of sample carries no obs

37.4 Quality codes, as published

37.5 What the source documentation says about these numbers

37.6 QC rules over this ingest’s output

37.6.1 Values outside their declared range

37.6.2 Impossible temperatures

37.6.3 Negative pH, clustered by cruise

37.6.4 Full-resolution profile rules, one cruise at a time

37.7 Questions for Data Providers

38 Cleanup

39 Appendix: Interactive Gantt Chart

32.1 No cruise may leave `obs` and keep its casts

37.1 The `-99` sentinel — fixed at the source, and now guarded

37.2 `btl_*` promoted to canonical — and the depth-thinning trap

37.3 One direction per physical cast: half of `sample` carries no `obs`