Ingest CalCOFI Bottle Derived Hydrographic Products

Published

2026-10-21

1 Overview

Hydrographic values derived from the CalCOFI Bottle Database rather than measured: spice and the TEOS-10 sigma-theta at every bottle, back to the start of the record (1949). Asked for by Nastassia Patin (spice at each bottle depth, 2026-10-16) and Rasmus Swalethorp (thread “calcofi derived data”, 2026-10-19: “keeping derived values in a separate dataset. Back calculating spice and sigma theta to the beginning using the values in the bottle database”); Ben Best’s reply (2026-10-20): “a separate bottle-derived dataset, back-calculated from the start of the record. From the bottle database’s own temperature and salinity … Spiciness and TEOS-10 sigma theta go in first.” The other bottle-derived values (IntC14, integrated chlorophyll) join this dataset once the DECODR source arrives.

Source: the calcofi_bottle ingest’s temperature (T_degC) and salinity (Salnty) per bottle, with their quality codes (T_qual, S_qual), its depth and its cast’s position — read from that ingest’s staged output when it is current, otherwise from the release (see Input).
The rule is one tested calcofi4db function, bottle_spice_sigma0(): sea pressure from depth and latitude (gsw_p_from_z), absolute salinity, conservative temperature, then gsw_spiciness0 (by ctd_spice(), the rule the CTD series uses) and gsw_sigma0. This notebook only applies it.
Flags first. A bottle whose temperature or salinity the provider flags 8 (suspect) or 9 (missing) gets neither value: a derived value carries no flag, so nothing downstream could undo a flagged input (measurement-bounds skill).
The legacy value stays. The bottle database’s own sigma_theta (computed by DECODR’s algorithm) is untouched in calcofi_bottle; sigma_theta_teos10 is the modern value beside it (Rasmus Swalethorp, 2026-10-16: “keep legacy, add modern, document”). Where the two disagree by more than 0.1 kg m⁻³ the bottle is listed for the provider (Findings).
A known-wrong input is withheld by name. Cruise 2021-05-3322’s salinity column holds sigma-theta (calcofi_bottle Q13), so its bottles get no derived value until the provider answers (Ben Best, 2026-10-20); the notebook asserts none is emitted.

Code

graph LR
  A[calcofi_bottle<br/>temperature + salinity<br/>per bottle, with T_qual / S_qual] --> B[drop 8/9 flags]
  B --> C[bottle_spice_sigma0:<br/>p from depth, SA, CT]
  C --> D[(obs: spiciness0,<br/>sigma_theta_teos10<br/>on the bottle's sample)]

graph LR
  A[calcofi_bottle<br/>temperature + salinity<br/>per bottle, with T_qual / S_qual] --> B[drop 8/9 flags]
  B --> C[bottle_spice_sigma0:<br/>p from depth, SA, CT]
  C --> D[(obs: spiciness0,<br/>sigma_theta_teos10<br/>on the bottle's sample)]

2 Setup

Code

# a development checkout of calcofi4db only when asked (CALCOFI4DB_DIR, as
# build_workflows_index.R does); a pipeline render uses the installed package
if (nzchar(Sys.getenv("CALCOFI4DB_DIR")))
  devtools::load_all(Sys.getenv("CALCOFI4DB_DIR"), quiet = TRUE)
librarian::shelf(
  CalCOFI/calcofi4db, CalCOFI/calcofi4r,
  DBI, dplyr, DT, fs, glue, here, jsonlite, purrr, readr, tibble, tidyr,
  quiet = T)
options(readr.show_col_types = F)
options(DT.options = list(scrollX = TRUE))
stopifnot(
  "calcofi4db >= 4.27.0 carries bottle_spice_sigma0() and add_measurement_source_dataset()" =
    exists("bottle_spice_sigma0", mode = "function") &&
    exists("add_measurement_source_dataset", mode = "function"))

source(here("libs/ingest.R"))
cc <- read_calcofi_meta(here("ingest_calcofi_bottle-derived.qmd"))
provider     <- cc$provider
dataset      <- cc$dataset
dataset_name <- cc$dataset_meta$dataset_name
tables_owned <- cc$tables_owned
ds_key       <- glue("{provider}_{dataset}")
src_key      <- "calcofi_bottle"
dir_label    <- ds_key
dir_meta     <- here(glue("metadata/{provider}/{dataset}"))
dir_parquet  <- here(glue("data/parquet/{dir_label}"))
dir_stage    <- cc_stage_path("parquet", dir_label, create = TRUE)
db_path      <- here(glue("data/wrangling/{dir_label}.duckdb"))
meas_type_csv <- here("metadata/measurement_type.csv")

# the source series, by their calcofi_bottle measurement_type
src_types <- c(temperature = "temperature", salinity = "salinity", legacy = "sigma_theta")
# a legacy-vs-TEOS-10 difference larger than this is listed for the provider; the
# release v2026.10.14 comparison over 199,918 bottles had a median of 0.0037 kg/m3
# and a p95 of 0.0087, so 0.1 is about ten times the p95 and is not a rounding effect
diff_list_kg_m3 <- 0.1

# WITHHELD, by name: cruises whose source T or S is known to be wrong, so nothing is
# derived from them until the provider answers. Ben Best, 2026-10-20: "a value we know is
# wrong never enters a derived product".
# - 2021-05-3322: on every bottle the published Salnty equals the bottle's own sigma_theta
#   (24.19-27.06), so the salinity column holds sigma-theta; calcofi_bottle Q13 (proposed,
#   high) asks the provider. Remove the entry when Q13 is answered and calcofi_bottle
#   re-staged; this dataset's Q03 tracks it.
withheld_cruises <- c("2021-05-3322" = "calcofi_bottle Q13: Salnty holds sigma-theta")

if (overwrite) {
  for (f in c(db_path, paste0(db_path, ".wal"))) if (file_exists(f)) file_delete(f)
}
dir_create(dirname(db_path))
con <- get_duckdb_con(db_path)

3 Input

The preferred input is the calcofi_bottle ingest’s own staged output under cc_stage_dir() — the bytes the next release will be cut from — used only if it is current: its obs and sample row counts must equal what that ingest’s committed manifest.json recorded (its data_hash starts with the row count). Otherwise the notebook reads the promoted release’s obs_env (the temperature, salinity and sigma-theta partitions only) and sample through the resolver helpers (cc_catalog() / cc_release_sources()), never a hand-built path.

Code

dir_src <- cc_stage_path("parquet", src_key)
mf_src  <- read_json(here(glue("data/parquet/{src_key}/manifest.json")))
rows_of <- function(h) sum(vapply(as.list(unlist(h)), function(x) as.numeric(sub(":.*$", "", x)), numeric(1)))
types_in <- paste0("'", src_types, "'", collapse = ", ")

# the staged obs is one file or a hive-partitioned directory; either is read the same way
src_obs_sql <- local({
  if (dir_exists(path(dir_src, "obs")))
    glue("read_parquet('{dir_src}/obs/**/*.parquet', hive_partitioning = true)")
  else glue("read_parquet('{dir_src}/obs.parquet')")
})
stage_current <- local({
  f_smp <- path(dir_src, "sample.parquet")
  if (!dir_exists(dir_src) || !file_exists(f_smp) ||
      !(file_exists(path(dir_src, "obs.parquet")) || dir_exists(path(dir_src, "obs")))) return(FALSE)
  got_obs <- dbGetQuery(con, glue("SELECT COUNT(*) AS n FROM {src_obs_sql}"))$n
  got_smp <- dbGetQuery(con, glue("SELECT COUNT(*) AS n FROM read_parquet('{f_smp}')"))$n
  exp_obs <- rows_of(mf_src$data_hash$obs); exp_smp <- rows_of(mf_src$data_hash$sample)
  ok <- got_obs == exp_obs && got_smp == exp_smp
  cat(glue("stage {dir_src}: obs {got_obs} / {exp_obs}, sample {got_smp} / {exp_smp} -> ",
           "{if (ok) 'CURRENT' else 'NOT current'}"), "\n")
  ok
})

stage /Users/bbest/_big/calcofi/parquet/calcofi_bottle: obs 11135581 / 11135581, sample 931015 / 931015 -> CURRENT

Code

if (stage_current) {
  input_used <- glue("staged calcofi_bottle output ({dir_src}), current against its committed manifest.json")
  dbExecute(con, glue("
    CREATE OR REPLACE VIEW src_obs AS
    SELECT sample_key, measurement_type, measurement_value, measurement_qual,
           grid_key, cruise_key, latitude, longitude, datetime, depth_min_m, depth_max_m
    FROM {src_obs_sql}
    WHERE dataset_key = '{src_key}' AND measurement_type IN ({types_in})"))
  dbExecute(con, glue("
    CREATE OR REPLACE VIEW src_sample AS
    SELECT sample_key, sample_type, site_key FROM read_parquet('{dir_src}/sample.parquet')
    WHERE dataset_key = '{src_key}'"))
} else {
  cat_ <- calcofi4r::cc_catalog("latest")
  input_used <- glue("release {cat_$version} obs_env + sample via cc_release_sources() (the stage was absent or not current)")
  invisible(dbExecute(con, "INSTALL httpfs; LOAD httpfs"))
  src_env <- calcofi4r::cc_release_sources(cat_, "obs_env")
  keep <- grepl(glue("measurement_type=({paste(src_types, collapse = '|')})/"), src_env$urls)
  stopifnot("the release must carry every source partition" = sum(keep) == length(src_types))
  src_env$urls <- src_env$urls[keep]
  # obs_env names the value `value`; the core assembly grain names it measurement_value
  dbExecute(con, paste0("
    CREATE OR REPLACE TABLE src_obs AS
    SELECT sample_key, measurement_type, value AS measurement_value, measurement_qual,
           grid_key, cruise_key, latitude, longitude, datetime, depth_min_m, depth_max_m
    FROM ", calcofi4r::cc_read_parquet_sql(src_env), "
    WHERE dataset_key = '", src_key, "'"))
  dbExecute(con, paste0("
    CREATE OR REPLACE TABLE src_sample AS
    SELECT sample_key, sample_type, site_key FROM ",
    calcofi4r::cc_read_parquet_sql(calcofi4r::cc_release_sources(cat_, "sample")), "
    WHERE dataset_key = '", src_key, "'"))
}

[1] 0

Code

cat(glue("input: {input_used}"), "\n")

input: staged calcofi_bottle output (/Users/bbest/_big/calcofi/parquet/calcofi_bottle), current against its committed manifest.json

4 Bottles

One row per bottle with its temperature, salinity and the legacy sigma-theta, each with its quality code, plus the bottle’s depth and its cast’s position and time (carried on every bottle obs row).

Code

col_of <- function(nm, tp) glue("
  MAX(measurement_value) FILTER (WHERE measurement_type = '{tp}') AS {nm},
  MAX(regexp_replace(CAST(measurement_qual AS VARCHAR), '\\.0+$', ''))
    FILTER (WHERE measurement_type = '{tp}') AS q_{nm},
  COUNT(*) FILTER (WHERE measurement_type = '{tp}') AS n_{nm}")
dbExecute(con, glue("
  CREATE OR REPLACE TABLE btl AS
  SELECT sample_key,
         ANY_VALUE(grid_key) AS grid_key, ANY_VALUE(cruise_key) AS cruise_key,
         ANY_VALUE(latitude) AS latitude, ANY_VALUE(longitude) AS longitude,
         ANY_VALUE(datetime) AS datetime,
         MIN(depth_min_m) AS depth_m, MAX(depth_min_m) AS depth_m_max,
         {paste(imap_chr(src_types, ~ col_of(.y, .x)), collapse = ',\n')}
  FROM src_obs
  GROUP BY sample_key"))

[1] 887667

Code

d_btl_chk <- dbGetQuery(con, "
  SELECT COUNT(*) AS bottles,
         MAX(n_temperature) AS max_t, MAX(n_salinity) AS max_s, MAX(n_legacy) AS max_legacy,
         SUM(CASE WHEN depth_m <> depth_m_max THEN 1 ELSE 0 END) AS n_depth_disagree,
         SUM(CASE WHEN temperature IS NOT NULL AND salinity IS NOT NULL THEN 1 ELSE 0 END) AS with_t_and_s,
         MIN(datetime) FILTER (WHERE temperature IS NOT NULL AND salinity IS NOT NULL) AS first_t_and_s,
         MAX(datetime) FILTER (WHERE temperature IS NOT NULL AND salinity IS NOT NULL) AS last_t_and_s
  FROM btl")
d_btl_chk |> datatable(caption = "Bottles carrying a source series", rownames = FALSE)

Code

stopifnot(
  "a bottle carries one temperature, one salinity and one legacy sigma-theta at most" =
    d_btl_chk$max_t <= 1 && d_btl_chk$max_s <= 1 && d_btl_chk$max_legacy <= 1,
  "a bottle has one depth" = d_btl_chk$n_depth_disagree == 0,
  "every source key is a calcofi_bottle bottle" = dbGetQuery(con, "
    SELECT COUNT(*) AS n FROM btl b LEFT JOIN src_sample s USING (sample_key)
    WHERE s.sample_type IS DISTINCT FROM 'bottle'")$n == 0)

dbGetQuery(con, "
  SELECT 'temperature' AS series, q_temperature AS qual, COUNT(*) AS bottles
  FROM btl WHERE temperature IS NOT NULL GROUP BY ALL
  UNION ALL
  SELECT 'salinity', q_salinity, COUNT(*) FROM btl WHERE salinity IS NOT NULL GROUP BY ALL
  ORDER BY 1, 2") |>
  datatable(caption = "Quality codes on the bottles' temperature and salinity (8 = suspect, 9 = missing are dropped; 6 = ok, taken from the CTD)", rownames = FALSE)

5 Compute the Products

The bottles of a withheld cruise (withheld_cruises in Setup) are set aside before anything is computed: their source salinity is known to be wrong (calcofi_bottle Q13), and a value we know is wrong never enters a derived product (Ben Best, 2026-10-20). They are listed below and in question Q03, and return when Q13 is answered.

Code

d_all <- dbGetQuery(con, "
  SELECT sample_key, cruise_key, temperature, q_temperature, salinity, q_salinity, legacy, q_legacy,
         depth_m, latitude, longitude
  FROM btl WHERE temperature IS NOT NULL AND salinity IS NOT NULL")
is_withheld <- d_all$cruise_key %in% names(withheld_cruises)
d_withheld  <- d_all[is_withheld, ]
d           <- d_all[!is_withheld, ]
d_withheld |>
  count(cruise_key, name = "bottles") |>
  mutate(reason = withheld_cruises[cruise_key],
         salinity_equals_legacy_sigma = vapply(cruise_key, function(ck) with(
           d_withheld[d_withheld$cruise_key == ck, ], sum(abs(salinity - legacy) < 0.01, na.rm = TRUE)),
           numeric(1))) |>
  datatable(caption = "Bottles withheld from the derivation (question Q03)", rownames = FALSE)

Code

stopifnot("every withheld cruise is still in the source (else drop it from the list)" =
            all(names(withheld_cruises) %in% d_all$cruise_key))
t0 <- Sys.time()
r <- bottle_spice_sigma0(d$temperature, d$salinity, d$depth_m, d$longitude, d$latitude,
                         q_temperature = d$q_temperature, q_salinity = d$q_salinity)
d <- bind_cols(d, r)
cat(glue("computed {format(nrow(d), big.mark = ',')} bottles in {format(round(Sys.time() - t0, 1))}"), "\n")

computed 843,828 bottles in 0.9 secs

Code

dbWriteTable(con, "derived", as.data.frame(
  select(d, sample_key, pressure_dbar, spiciness0, sigma_theta_teos10)), overwrite = TRUE)

6 Assert the Flag Rule

A value its provider flags 8 or 9 never enters a derived value; every bottle with a good temperature and salinity at a known depth and position gets both values.

Code

is_flagged <- function(q) sub("\\.0+$", "", as.character(q)) %in% c("8", "9")
flagged <- is_flagged(d$q_temperature) | is_flagged(d$q_salinity)
placed  <- is.finite(d$depth_m) & is.finite(d$latitude) & is.finite(d$longitude)
fa <- tibble(
  withheld_cruise_bottles  = nrow(d_withheld),
  bottles_t_and_s          = nrow(d),
  t_or_s_flagged_8_9       = sum(flagged),
  t_flagged                = sum(is_flagged(d$q_temperature)),
  s_flagged                = sum(is_flagged(d$q_salinity)),
  no_depth_or_position     = sum(!flagged & !placed),
  flagged_with_value       = sum(flagged & (!is.na(d$spiciness0) | !is.na(d$sigma_theta_teos10))),
  good_without_value       = sum(!flagged & placed & (is.na(d$spiciness0) | is.na(d$sigma_theta_teos10))),
  derived                  = sum(!is.na(d$sigma_theta_teos10)))
fa |> datatable(caption = "Flag audit over every bottle with a temperature and a salinity", rownames = FALSE)

Code

stopifnot(
  "a flagged temperature or salinity must yield no derived value" = fa$flagged_with_value == 0,
  "a good bottle at a known depth and position must get both values" = fa$good_without_value == 0,
  "spice and sigma-theta are present together" =
    all(is.na(d$spiciness0) == is.na(d$sigma_theta_teos10)))

7 Register Measurement Types

spiciness0 is reused under this dataset (Ben Best, D8, 2026-10-20): the same quantity by the same rule (ctd_spice()), so a reader filters one key by dataset_key. Its row is defined by calcofi_ctd-derived, which keeps the description, units, derivation and bound; this ingest only adds itself to the row’s _source_datasets (add_measurement_source_dataset()). sigma_theta_teos10 is new (D8): appended with register_measurement_types(), its registry-owned fields set with declare_measurement_fields(), its bound with declare_measurement_bounds(). It is not the bottle’s sigma_theta (DECODR’s legacy algorithm, labelled so in its own row) and it does not join that key’s variable: both come from the same bottles, and a crosswalk unifies series that are not the same physical samples.

No bound is invented or taken from the data. The sigma-theta bound is gsw_sigma0 over the bottle’s registered temperature × salinity bounds at 0 dbar, rounded outward — what the registry already calls physically possible (question Q02). The reused spiciness0 bound (−25 to 18, calcofi_ctd-derived Q07) is checked to cover the same grid.

Code

d_mt0 <- read_measurement_type(meas_type_csv)
src_mt <- function(tp) filter(d_mt0, measurement_type == tp)
t_b <- src_mt("temperature"); s_b <- src_mt("salinity"); sp0 <- src_mt("spiciness0")
stopifnot(nrow(t_b) == 1, nrow(s_b) == 1,
          "spiciness0 is defined by calcofi_ctd-derived and must exist" = nrow(sp0) == 1)

new_type <- tibble(
  measurement_type   = "sigma_theta_teos10",
  description        = "Potential density anomaly (sigma-theta, TEOS-10) per bottle, from the bottle's own temperature and salinity",
  units              = "kg/m3",
  grain              = "obs",
  `_source_column`   = "temperature; salinity; depth",
  `_source_table`    = "obs",
  `_source_datasets` = ds_key,
  `_qual_column`     = NA_character_,
  `_prec_column`     = NA_character_)
d_meas_type <- register_measurement_types(new_type, meas_type_csv)
# derivation, is_canonical, category and the NERC ids are registry-owned: this ingest
# owns the row, so it declares them deliberately (overwrite = TRUE)
d_meas_type <- declare_measurement_fields(
  tibble(
    measurement_type = "sigma_theta_teos10",
    category         = "Physical Oceanography",
    derivation       = paste(
      "calcofi4db::bottle_spice_sigma0(): sea pressure gsw_p_from_z(-depth, latitude), then",
      "gsw_SA_from_SP, gsw_CT_from_t, gsw_sigma0 (TEOS-10, IOC/SCOR/IAPSO 2010) on the bottle's",
      "own temperature (T_degC) and salinity (Salnty); a bottle whose T_qual or S_qual is 8 or 9",
      "has none. The bottle database's own sigma_theta is DECODR's legacy value, kept beside it"),
    is_canonical     = "TRUE",
    # exact match: potential density anomaly {Sigma-0} of the water body by computation from
    # absolute salinity and conservative temperature, UNESCO 2010 (= TEOS-10); checked against
    # the NVS SPARQL endpoint 2026-10-20, not deprecated
    nerc_p01         = "http://vocab.nerc.ac.uk/collection/P01/current/SIGTPR03/",
    units_nerc_p06   = src_mt("sigma_theta")$units_nerc_p06),
  meas_type_csv, categories = read_csv(here("metadata/category.csv"))$category,
  overwrite = TRUE, quiet = TRUE)
d_meas_type <- add_measurement_source_dataset("spiciness0", ds_key, meas_type_csv)

ts_grid <- expand.grid(t = seq(t_b$valid_min, t_b$valid_max, 1),
                       s = seq(s_b$valid_min, s_b$valid_max, 1))
sa_g <- gsw::gsw_SA_from_SP(ts_grid$s, 0, -120, 33)
ct_g <- gsw::gsw_CT_from_t(sa_g, ts_grid$t, 0)
sig_g <- gsw::gsw_sigma0(sa_g, ct_g); sp_g <- gsw::gsw_spiciness0(sa_g, ct_g)
stopifnot(
  "the reused spiciness0 bound must cover the bottle's registered T x S bounds" =
    sp0$valid_min <= min(sp_g) && sp0$valid_max >= max(sp_g))

bounds <- tribble(
  ~measurement_type,    ~valid_min,        ~valid_max,          ~why,
  "sigma_theta_teos10", floor(min(sig_g)), ceiling(max(sig_g)), "gsw_sigma0 over the registered bottle temperature x salinity bounds at 0 dbar, rounded outward (Q02)")
d_meas_type <- declare_measurement_bounds(select(bounds, -why), meas_type_csv, overwrite = TRUE)
bounds |>
  bind_rows(tibble(measurement_type = "spiciness0", valid_min = sp0$valid_min, valid_max = sp0$valid_max,
                   why = "reused: declared by calcofi_ctd-derived (its Q07); covers the bottle T x S grid")) |>
  datatable(caption = "Declared bounds and where each comes from", rownames = FALSE)

Code

stopifnot(
  "this dataset is listed on both of its types" = all(vapply(
    c("spiciness0", "sigma_theta_teos10"), function(tp) ds_key %in% strsplit(
      filter(d_meas_type, measurement_type == tp)$`_source_datasets`, ";")[[1]], logical(1))),
  "spiciness0 keeps calcofi_ctd-derived" = grepl("calcofi_ctd-derived",
    filter(d_meas_type, measurement_type == "spiciness0")$`_source_datasets`, fixed = TRUE))
dbWriteTable(con, "measurement_type", as.data.frame(d_meas_type), overwrite = TRUE)

8 Load Dataset Metadata

Code

this_provider <- provider; this_dataset <- dataset
d_dataset <- ingest_yaml_to_dataset_df(read_ingest_yaml(here())) |>
  filter(provider == this_provider, dataset == this_dataset)
stopifnot(nrow(d_dataset) == 1)
dbWriteTable(con, "dataset", as.data.frame(d_dataset), overwrite = TRUE)

9 Emit Core Tables

Per bottle → obs, on the bottle’s calcofi_bottle sample_key, at the bottle’s depth; position, time, grid and cruise come from the same bottle row, so a derived value lines up with the temperature and salinity it came from and the release climatology covers it with no special case. This dataset emits no sample: the bottle is calcofi_bottle’s, and a second row with the same key would break the global uniqueness assemble_core() asserts.

Code

for (tp in c("spiciness0", "sigma_theta_teos10"))
  append_obs(con, glue("
    SELECT 'env', '{ds_key}', b.sample_key, b.grid_key, b.cruise_key,
           b.latitude, b.longitude, b.datetime, b.depth_m, b.depth_m,
           NULL::VARCHAR, NULL::VARCHAR, '{tp}', x.{tp}, NULL::VARCHAR, NULL::DOUBLE
    FROM btl b JOIN derived x USING (sample_key)
    WHERE x.{tp} IS NOT NULL"))

n_obs <- dbGetQuery(con, "SELECT measurement_type, COUNT(*) AS n FROM obs GROUP BY 1 ORDER BY 1")
n_obs |> datatable(caption = "core projection: obs rows per type", rownames = FALSE)

Code

stopifnot(
  "one obs row per derived bottle and type" = all(n_obs$n == fa$derived),
  "both types emitted" = nrow(n_obs) == 2,
  "no row of a withheld cruise is emitted (calcofi_bottle Q13)" = dbGetQuery(con, glue("
    SELECT COUNT(*) AS n FROM obs
    WHERE cruise_key IN ({paste0(\"'\", names(withheld_cruises), \"'\", collapse = ', ')})
       OR sample_key IN (SELECT sample_key FROM btl
                         WHERE cruise_key IN ({paste0(\"'\", names(withheld_cruises), \"'\", collapse = ', ')}))"))$n == 0)

10 Check Measurement Bounds

Code

b_obs <- check_measurement_bounds(con, "obs", mt = d_meas_type, dataset_key = ds_key)
bounds_datatable(b_obs)

Code

stopifnot(
  "every derived type is bounded" = all(b_obs$status != "undeclared"),
  "no derived value breaks its bound" = all(b_obs$status != "out_of_range"))

11 Findings for the Provider

Values inside the generous bounds are not deleted: a bound catches the impossible. Two lists go to the provider (question Q03): bottles where the TEOS-10 sigma-theta and the bottle database’s own (legacy, DECODR) sigma-theta differ by more than 0.1 kg m⁻³ — a gap that size is not the algorithm, it is one of the inputs (a temperature or salinity the legacy value was not computed from) — and spice below −3 kg m⁻³, which needs a practical salinity near 20, the same class of fault the CTD derived products list (calcofi_ctd-derived Q08).

Cruise 2021-05-3322 is not in these lists because it is withheld (see Compute the Products): at the first render (2026-10-20) it was 918 of the 928 bottles past 0.1 kg m⁻³ and every bottle with spice below −3, because on every bottle the published salinity equals the bottle’s own sigma-theta (calcofi_bottle Q13). Its source code is blank (good, by calcofi_bottle Q01) and the values sit inside the generous bounds, so no flag or bound would have caught it; the withholding is by name, until Q13 is answered.

Code

d_diff <- d |>
  filter(!is.na(sigma_theta_teos10), !is.na(legacy), !is_flagged(q_legacy)) |>
  mutate(diff = sigma_theta_teos10 - legacy)
d_diff_sum <- d_diff |>
  summarize(bottles = n(),
            median_abs = median(abs(diff)), p95_abs = quantile(abs(diff), 0.95),
            p99_abs = quantile(abs(diff), 0.99), max_abs = max(abs(diff)),
            n_gt_0.01 = sum(abs(diff) > 0.01), n_gt_0.1 = sum(abs(diff) > diff_list_kg_m3),
            n_gt_1 = sum(abs(diff) > 1))
d_diff_sum |> mutate(across(ends_with("_abs"), ~ signif(.x, 3))) |>
  datatable(caption = "TEOS-10 minus the bottle database's sigma_theta (kg/m3), over bottles with both and an unflagged legacy value", rownames = FALSE)

Code

d_diff_odd <- d_diff |>
  filter(abs(diff) > diff_list_kg_m3) |>
  left_join(dbGetQuery(con, "SELECT sample_key, datetime FROM btl"), by = "sample_key") |>
  transmute(sample_key, cruise_key, date = as.Date(datetime), depth_m,
            temperature, q_temperature, salinity, q_salinity,
            sigma_theta_legacy = round(legacy, 3), q_legacy,
            sigma_theta_teos10 = round(sigma_theta_teos10, 3), diff = round(diff, 3)) |>
  arrange(desc(abs(diff)))
d_diff_odd |> datatable(caption = glue(
  "Bottles where TEOS-10 and legacy sigma-theta differ by more than {diff_list_kg_m3} kg/m3: ",
  "{nrow(d_diff_odd)} (question Q03)"), rownames = FALSE, filter = "top")

Code

d_spice_odd <- dbGetQuery(con, "
  SELECT b.cruise_key, b.sample_key, b.depth_m, ROUND(x.spiciness0, 2) AS spice,
         ROUND(x.sigma_theta_teos10, 2) AS sigma_theta_teos10
  FROM derived x JOIN btl b USING (sample_key)
  WHERE x.spiciness0 < -3 ORDER BY x.spiciness0")
d_spice_odd |> datatable(caption = glue(
  "Spice below -3 kg/m3: {nrow(d_spice_odd)} bottles (question Q03)"), rownames = FALSE)

12 Relationships

Code

stopifnot(
  "every obs.sample_key is a calcofi_bottle bottle" =
    dbGetQuery(con, "SELECT COUNT(*) AS n FROM obs o LEFT JOIN src_sample s USING (sample_key)
                     WHERE s.sample_type IS DISTINCT FROM 'bottle'")$n == 0,
  "every measurement_type is registered" =
    dbGetQuery(con, "SELECT COUNT(*) AS n FROM obs
                     WHERE measurement_type NOT IN (SELECT measurement_type FROM measurement_type)")$n == 0,
  "one obs row per bottle and type" =
    dbGetQuery(con, "SELECT COUNT(*) AS n FROM (SELECT sample_key, measurement_type FROM obs
                     GROUP BY ALL HAVING COUNT(*) > 1)")$n == 0)

tbls_out <- core_output_tables(con, extra = c("measurement_type", "dataset"))
rels <- core_relationships(tbls_out)
cc_erd(con, tables = c("obs", "measurement_type"), rels = rels,
       colors = list(lightblue = "obs", lightyellow = "measurement_type"))

Table 1: Source → core mapping

source (calcofi_bottle `obs`)	core	transform
`temperature`, `salinity` + `T_qual`, `S_qual`, the bottle’s depth and its cast’s position	`obs.measurement_type = spiciness0`	`bottle_spice_sigma0()` → `ctd_spice()` at `gsw_p_from_z(-depth, latitude)`
the same	`obs.measurement_type = sigma_theta_teos10`	`bottle_spice_sigma0()`: `gsw_sigma0`
`sigma_theta` (legacy, DECODR)	not published here; compared in Findings

13 Data Preview

Code

dbGetQuery(con, "SELECT * FROM obs ORDER BY datetime DESC, sample_key, measurement_type LIMIT 200") |>
  datatable(caption = "obs — the most recent 200 rows", rownames = FALSE, filter = "top")

Code

dbGetQuery(con, "
  SELECT measurement_type, COUNT(*) AS bottles,
         MIN(datetime)::DATE AS first, MAX(datetime)::DATE AS last,
         ROUND(MIN(measurement_value), 3) AS min, ROUND(QUANTILE_CONT(measurement_value, 0.01), 3) AS p01,
         ROUND(MEDIAN(measurement_value), 3) AS median,
         ROUND(QUANTILE_CONT(measurement_value, 0.99), 3) AS p99, ROUND(MAX(measurement_value), 3) AS max
  FROM obs GROUP BY 1 ORDER BY 1") |>
  datatable(caption = "Derived values over every bottle", rownames = FALSE)

14 Validate

validate_for_release() treats every _key column as required. The NULLs below are declared with their reason and count; anything undeclared, or a count that moved, fails the render.

Code

# the expected NULL counts, derived independently of the table being checked (from
# the source bottles that got a value, not from obs)
n_exp <- dbGetQuery(con, "
  SELECT COUNT(*) AS n_obs,
         2 * SUM(CASE WHEN b.grid_key IS NULL THEN 1 ELSE 0 END) AS n_grid_null,
         2 * SUM(CASE WHEN b.cruise_key IS NULL THEN 1 ELSE 0 END) AS n_cruise_null
  FROM btl b JOIN derived x USING (sample_key) WHERE x.sigma_theta_teos10 IS NOT NULL")
# the wrangling tables are not outputs; drop them so the validator sees what ships. The
# src_* inputs are views (stage) or tables (release), and DuckDB refuses DROP VIEW on a
# table, so each is dropped as what it is
for (tb in c("btl", "derived", "src_obs", "src_sample")) {
  kind <- dbGetQuery(con, glue(
    "SELECT table_type FROM information_schema.tables WHERE table_name = '{tb}'"))$table_type
  if (length(kind))
    dbExecute(con, glue("DROP {if (grepl('VIEW', kind[1])) 'VIEW' else 'TABLE'} IF EXISTS {tb}"))
}
results <- validate_for_release(con, checks = "all", strict = FALSE)
nullable <- tribble(
  ~table, ~column,      ~n,                       ~reason,
  "obs",  "taxon_key",  2 * n_exp$n_obs,          "env realm: a derived physical quantity names no organism",
  "obs",  "grid_key",   n_exp$n_grid_null,        "inherited: the calcofi_bottle cast has no grid cell",
  "obs",  "cruise_key", n_exp$n_cruise_null,      "inherited: the calcofi_bottle cast's cruise is not resolved") |>
  filter(n > 0)
found <- tibble(msg = results$errors) |>
  mutate(table  = sub("^Table '([^']+)'.*$", "\\1", msg),
         column = sub("^.*required column '([^']+)'.*$", "\\1", msg),
         n      = suppressWarnings(as.integer(sub("^Table '[^']+' has (\\d+) NULL.*$", "\\1", msg))))
unexplained <- anti_join(found, nullable, by = c("table", "column"))
moved <- inner_join(found, nullable, by = c("table", "column"), suffix = c("", "_decl")) |>
  filter(!is.na(n_decl), n != n_decl)
nullable |> left_join(select(found, table, column, n_found = n), by = c("table", "column")) |>
  datatable(caption = "Declared NULLs vs what validate_for_release() found", rownames = FALSE)

Code

if (nrow(unexplained)) print(unexplained$msg)
stopifnot(
  "a NULL appeared with no declared reason" = nrow(unexplained) == 0,
  "a declared NULL count has changed"       = nrow(moved) == 0)

15 Write Parquet Outputs

Code

d_tbls_rd <- read_csv(glue("{dir_meta}/tbls_redefine.csv"))
d_flds_rd <- read_csv(glue("{dir_meta}/flds_redefine.csv"))
dir_create(dir_parquet)
tbls_out <- core_output_tables(con, extra = c("measurement_type", "dataset"))
parquet_stats <- write_parquet_outputs(
  con = con, output_dir = dir_parquet, tables = tbls_out,
  sort_by = list(obs = c("grid_key", "measurement_type", "depth_min_m")),
  strip_provenance = FALSE)
build_relationships_json(rels = rels, output_dir = dir_parquet, provider = provider, dataset = dataset)

[1] "/Users/bbest/Github/CalCOFI/workflows/data/parquet/calcofi_bottle-derived/relationships.json"

Code

parquet_stats |> mutate(file = basename(path)) |> select(-path) |>
  datatable(caption = "Parquet export statistics")

16 Write Metadata

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,
  set_comments = TRUE, provider = provider, dataset = dataset,
  workflow_url = cc$workflow_url, tables_owned = tables_owned)

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

17 Upload to GCS

Code

if (Sys.getenv("CALCOFI_SKIP_GCS") == "") {
  sync_to_gcs(local_dir = dir_stage, sidecar_dir = dir_parquet,
              gcs_prefix = glue("ingest/{dir_label}"), bucket = "calcofi-db")
} else cat("CALCOFI_SKIP_GCS set — skipping the upload (local verification run)\n")

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

18 Questions for Data Providers

Open: the licence and PI, inherited from the bottle database (Q01); the derived sigma-theta bound (Q02); the bottles where the TEOS-10 and legacy sigma-theta disagree and the spice outliers (Q03); the temperature scale of the early record (Q04).

Code

questions_datatable(
  here(cc$questions_file),
  caption = "Questions for the bottle derived products (ranked)")

19 Cleanup

Code

close_duckdb(con)
cat(glue("input used: {input_used}"), "\n")

input used: staged calcofi_bottle output (/Users/bbest/_big/calcofi/parquet/calcofi_bottle), current against its committed manifest.json

Code

cat(glue("Parquet outputs written to: {dir_stage}; sidecars in {dir_parquet}"), "\n")

Parquet outputs written to: /Users/bbest/_big/calcofi/parquet/calcofi_bottle-derived; sidecars in /Users/bbest/Github/CalCOFI/workflows/data/parquet/calcofi_bottle-derived