Ingest CalCOFI CTD Cast Data

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

1. Web Scraping: Scrapes all CTD .zip download links from calcofi.org

2. Smart Filtering:

   • Downloads all .zip files for archival completeness
   • Only unzips final and preliminary files
   • Skips raw/cast/test files

3. Priority-based Selection:

   • For each cruise, selects final if available, otherwise preliminary
   • Excludes raw/test/prodo cast files

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

5. Standardized Workflow: follows patterns from ingest_calcofi_bottle.qmd — deterministic UUIDs, GCS parquet uploads, metadata sidecars, calcofi4db utilities

2 Setup

# 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)

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

# 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}")
dir_dl <- path_expand(glue("{dir_data}/{provider}/{dataset}/download"))
dir_parquet <- here(glue("data/parquet/{dir_label}"))
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}"))

# overwrite: rebuild the wrangling DB, but keep checkpoint DB + downloads and
# (importantly) keep dir_parquet so write_parquet_outputs can content-hash
# dedup against the prior run — only changed partitions get re-written/uploaded.
# overwrite_all: also delete parquet, checkpoint DB, and RDS intermediates.
if (overwrite) {
  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 (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)
    message("Deleted parquet, checkpoint DB, and RDS intermediates")
  }
}

# restore from checkpoint DB if available (skips expensive read+bind+filter)
if (file_exists(db_checkpoint) && !file_exists(db_path)) {
  file_copy(db_checkpoint, db_path, overwrite = TRUE)
  message(glue("Restored from checkpoint: {db_checkpoint}"))
}
dir_create(c(dir_dl, dirname(db_path), dir_parquet, dir_tmp), recurse = TRUE)

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 RAM
mem_gb <- ps::ps_system_memory()$total / 1024^3 / 2 |> floor()
dbExecute(con, glue("SET memory_limit = '{mem_gb}GB'"))

[1] 0

dbExecute(con, glue("SET temp_directory = '{dir_tmp}'"))

[1] 0

# 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
d_meas_type <- read_csv(
  here("metadata/measurement_type.csv"),
  show_col_types = F
)
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)

3 Check for Resumable State

# detect if parquet outputs are already complete (e.g. prior run failed
# only during GCS upload). If overwrite=TRUE, always rebuild.
parquet_complete <- 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) {
    p <- file.path(dir_parquet, paste0(tbl, ".parquet"))
    d <- file.path(dir_parquet, tbl)
    file_exists(p) || dir_exists(d)
  }, logical(1)))
  if (parquet_ok && !overwrite) {
    parquet_complete <- TRUE
    message(glue(
      "Parquet output already complete ({length(mf$tables)} tables, ",
      "{format(mf$total_rows, big.mark = ',')} rows) — ",
      "skipping computation, resuming at upload"))
  }
}

# if parquet not complete, check for checkpoint (ctd_raw pre-computed)
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
    message(glue(
      "Checkpoint: ctd_raw already loaded ",
      "({format(n_raw, big.mark = ',')} rows) — ",
      "skipping read+bind+filter"))
  }
}

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

4 Scrape CTD Download Links

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

d_zips <- 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 = str_extract(url, "/(\\d{4})/", group = 1) |> as.integer(),
        month = str_extract(url, "-\\d{2}(\\d{2})", group = 1) |> as.integer(),
        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"
        )
      )
    write_csv(d, cache_csv)
    message(glue("Scraped {nrow(d)} zip URLs from {url}"))
    d
  },
  error = function(e) {
    if (file_exists(cache_csv)) {
      message(glue("Website unavailable, using cached URLs from {cache_csv}"))
      read_csv(cache_csv, show_col_types = FALSE)
    } else {
      stop(e)
    }
  }
)

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

d_zips |>
  mutate(
    file_zip = glue("<a href={url}>{file_zip}</a>")
  ) |>
  select(-url) |>
  dt(
    caption = "All zip files to download",
    fname = "ctd_zip_files",
    escape = F
  )

5 Prime Downloads from GCS Source (optional)

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. If that source exposes zips, copy them into dir_dl so the next chunk simply unzips them rather than re-scraping calcofi.org. This is a no-op (falls back to calcofi.org) until the GCS source is populated; set CTD_ZIP_SOURCE="" to disable, or override it to point elsewhere.

# 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))
  # only act if the source actually exposes zips (else fall back to scraping)
  n_src <- tryCatch(
    length(system2(
      "rclone", c("lsf", src, "--include", "*.zip", "--max-depth", "1"),
      stdout = TRUE, stderr = FALSE)),
    error = function(e) 0L)
  if (n_src == 0) {
    message(glue("No zips at {src} — skipping GCS prime (will use calcofi.org)"))
    return(invisible(0L))
  }
  message(glue("Priming {n_src} zip(s) from {src} → {dest_dir}"))
  system2("rclone", c(
    "copy", src, dest_dir,
    "--include", "*.zip", "--max-depth", "1",
    "--transfers", "8", "--checkers", "16"))
  invisible(n_src)
}

prime_zips_from_gcs(ctd_zip_source, dir_dl)

6 Download and Unzip Files

download_and_unzip <- function(url, dest_dir, zip_type, unzip = TRUE) {
  file_zip <- basename(url)
  dest_file <- file.path(dest_dir, file_zip)
  dir_unzip <- file.path(dest_dir, str_remove(file_zip, "\\.zip$"))

  if (file_exists(dest_file)) {
    message(glue("Already exists: {file_zip}"))
  } else {
    message(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)) {
    message(glue("Already unzipped: {file_zip}"))
  } else {
    if (zip_type %in% c("final", "preliminary")) {
      message(glue("Unzipping: {file_zip}"))
      dir_create(dir_unzip, recurse = TRUE)
      unzip(dest_file, exdir = dir_unzip)
    } else if (unzip) {
      message(glue("Skipping unzip (not final/preliminary): {file_zip}"))
    }
  }
}

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

message("Download and extraction complete!")

7 Find and Prioritize CTD Data Files

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(
      str_detect(path_unzip, "Final.*db[_|-]csv")  ~ "final",
      str_detect(path_unzip, "Prelim.*db[_|-]csv") ~ "preliminary",
      # prelim CSVs directly in unzipped dir (no db-csv subfolder)
      str_detect(dir_unzip, "Prelim") &
        path_unzip == paste0(dir_unzip, "/", file_csv) ~ "preliminary",
      # edge case: 2111SR has CSVs only in csvs-plots subfolder
      cruise_key == "2111SR" &
        str_detect(
          path_unzip,
          "Prelim.*csvs-plots.*2111SR.*csv"
        ) ~ "preliminary",
      .default = NA_character_
    ),
    cast_dir = case_when(
      str_detect(file_csv, regex("U\\.csv$", ignore_case = T)) ~ "U",
      str_detect(file_csv, regex("D\\.csv$", ignore_case = T)) ~ "D"
    ),
    priority = case_when(
      data_stage == "final" ~ 1,
      data_stage == "preliminary" ~ 2,
      TRUE ~ 3
    )
  ) |>
  relocate(cruise_key, path_unzip) |>
  arrange(cruise_key, path_unzip) |>
  filter(data_stage %in% c("final", "preliminary"))

# for each cruise, keep only final if available, otherwise preliminary
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)

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

8 Read and Standardize CTD Files

Remove repeat header rows that can occur within files.

d_csv <- d_csv |>
  arrange(basename(path)) |>
  mutate(
    data = map2(path, seq_along(path), \(path, idx) {
      message(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()
      }
      data
    }),
    nrows = map_int(data, ~ if (is.null(.x)) 0 else nrow(.x))
  )

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 = ",")

9 Detect and Correct Column Type Mismatches

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) {
  message("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) {
            message(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
  )

# 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")

10 Format Date-Time Column

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

11 Pseudo-NA Values and Fill Missing Coordinates

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

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

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)
  )

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 = ","
  )

12 Distance Filtering

max_dist_dec_lnst_km <- 10

pts <- d_bind |>
  st_as_sf(coords = c("longitude", "latitude"), remove = F, crs = 4326) |>
  mutate(
    geom_lnst = purrr::map2(lon_lnst, lat_lnst, \(x, y) st_point(c(x, y))) |>
      st_sfc(crs = 4326),
    dist_dec_lnst_km = st_distance(geometry, geom_lnst, by_element = T) |>
      units::set_units(km) |>
      units::drop_units(),
    is_dist_dec_lnst_within_max = dist_dec_lnst_km <= max_dist_dec_lnst_km
  ) |>
  st_join(
    calcofi4r::cc_grid |>
      rename(any_of(c(site_key = "sta_key"))) |>
      select(grid_site = site_key),
    join = st_intersects
  )
st_agr(pts) <- "constant"

12.1 View sample cruise with excess distance points

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

bb_badcr <- st_bbox(pts_badcr) |> 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)
      ) |>
      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
  )

12.2 Filter points

d_pts_cruise_filt_smry <- pts |>
  st_drop_geometry() |>
  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)
  )

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
  )

12.3 View filtered points

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

mapView(bb_cr) +
  mapView(
    pts_filt |>
      group_by(cruise_key) |>
      slice_sample(n = 100) |>
      ungroup() |>
      select(
        cruise_key,
        longitude,
        latitude,
        sta_id,
        line,
        sta,
        dist_dec_lnst_km,
        datetime_start_utc
      ),
    zcol = "cruise_key",
    cex = 5,
    alpha = 0.5
  )