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
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
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 chunkknitr::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)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 blockcc <-read_calcofi_meta(here("ingest_calcofi_ctd-cast.qmd"))provider <- cc$providerdataset <- cc$datasetdataset_name <- cc$dataset_meta$dataset_nametables_owned <- cc$tables_owneddir_label <-glue("{provider}_{dataset}")dir_dl <-path_expand(glue("{dir_data}/{provider}/{dataset}/download"))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, dirname(db_path), dir_parquet, dir_stage, dir_tmp), recurse =TRUE)# 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)# the ctd-cast download/ dir is 62GB+ — too large for inline# sync_to_gcs() through the GD FUSE mount# 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")
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.
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 primingctd_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 unavailableif (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 CSVcache_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 codeyear =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 272 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 restyear =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: 383 (272 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, "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.values =c(sort(paste(d_zips$source, 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 diskparquet_ok <-FALSEmanifest_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_rebuildif (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}"))} elseif (force_rebuild) {say("Rebuilding: forced (CTD_FORCE_REBUILD / overwrite_all)")} elseif (parquet_ok &&!fp_same) {say("Rebuilding: inputs changed since the recorded fingerprint —")for (f inchanged_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/metadata/measurement_type.csv
<values>
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)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}"))}con <-get_duckdb_con(db_path)load_duckdb_extension(con, "spatial")load_duckdb_extension(con, "icu")# limit memory to half of system RAM so DuckDB spills to disk instead of consuming all RAMmem_gb <- ps::ps_system_memory()$total /1024^3/2|>floor()dbExecute(con, glue("SET memory_limit = '{mem_gb}GB'"))
# a second, cheaper resume: ctd_raw already pivoted from a prior interrupted runhas_ctd_raw <-FALSEif (!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")$nsay(glue("Checkpoint: ctd_raw already loaded ","({format(n_raw, big.mark = ',')} rows) — ","skipping read+bind+filter")) }}# set eval for read+bind+filter chunksskip_read_bind <- parquet_complete || has_ctd_rawknitr::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.
ImportantWhere 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)"))
Archives to read: 293 of 383 (90 GCS finals skipped as duplicates of calcofi.org; 21 GCS finals are new coverage)
Code
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.download_and_unzip <-function(url, file_zip, dest_dir, zip_type, unzip =TRUE) { dest_file <-file.path(dest_dir, file_zip) dir_unzip <-file.path(dest_dir, str_remove(file_zip, "\\.zip$"))if (file_exists(dest_file)) {say(glue("Already exists: {file_zip}")) } elseif (is.na(url)) {# primed source vanished between the prime and here — do not silently proceed# with a missing archivewarning(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 &&dir_exists(dir_unzip)) {say(glue("Already unzipped: {file_zip}")) } else {if (zip_type %in%c("final", "preliminary")) {say(glue("Unzipping: {file_zip}"))dir_create(dir_unzip, recurse =TRUE)unzip(dest_file, exdir = dir_unzip) } elseif (unzip) {say(glue("Skipping unzip (not final/preliminary): {file_zip}")) } }}d_zips_read |>pwalk(function(url, file_zip, zip_type, ...) {download_and_unzip(url, file_zip, dir_dl, zip_type, unzip =TRUE) })
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_dl, 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)"))
Overlap comparison across 90 cruises / 188 files: 148 identical, 10 different, 30 with no same-named counterpart (split runs)
Code
# 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:
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_dl,pattern ="\\.csv$",recursive =TRUE,full.names =TRUE )) |>mutate(file_csv =basename(path),path_unzip =str_replace(path, glue("{dir_dl}/"), ""),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* archivestr_detect(path_unzip, "Final.*db[_|-]csv") ~"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 runstr_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 thesestr_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. `dir_dl` sits inside Google Drive, which names 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. Strip the marker, then read the# letter off the stem.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 download 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_bottled_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")
Code
# 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.
d_csv |>arrange(cruise_key, file_csv) |>select(cruise_key, file_csv, nrows) |>dt(caption ="Number of rows read per file",fname ="ctd_rows_per_file" ) |>formatCurrency("nrows", currency ="", digits =0, mark =",")
Code
# A file that reads zero rows is a failure, not an observation — and it is the# one failure this pipeline could not see.## `dir_dl` sits inside Google Drive, which evicts a synced file to a cloud-only# placeholder. `list.files()` still reports the full size (12,254,953 bytes for# 19-9501JD_CTDBTL_001-646D.csv) and `ls -lO` marks it `dataless`, but reading it# times out at the filesystem and `read_csv()` returns a **0-row tibble with no# error**. On the v2026.08.08 run 20 files across 14 cruises read 0 rows exactly# this way, and the only copy holding data was the `… 2.csv` conflict copy Drive# had made beside each one. Nothing downstream distinguishes "this cruise has no# downcast" from "the downcast file was empty".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. Check for a Google Drive ","cloud-only placeholder (`ls -lO` reports `dataless`) and materialize it ","with `cat file > /dev/null`, or move the download directory off Drive."))
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 filesd_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...")}
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 oncelnst_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"))
Distance filter: 213 distinct station coordinates across 9,646,676 scans
Code
# 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 duplicateddist_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 <-NULLst_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.d_dist <- d_bind |>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_kmd_bind$is_dist_dec_lnst_within_max <- d_dist$is_dist_dec_lnst_within_maxd_bind$grid_site <- d_dist$grid_siterm(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 itas_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 )