library(dplyr)
library(tidyr)
library(readr)
library(knitr)
library(DT)
library(lubridate)
library(table.glue)
source("R/utils.R")
source("R/harmonize.R")
source("R/units.R")
source("R/qaqc.R")
source("R/database.R")

Background

Environmental monitoring programs typically submit water quality samples to one or more external laboratories under distinct submission codes. Each code corresponds to a method bundle — for example, a freshwater nutrient panel (E100), an ion chromatography panel (E235), or a seawater metals package (E469S). The same compound can have a different name under each code: what E100 calls “Ammonia, total (as N)” is the same thing E235 calls “Nitrogen - Ammonia as N”.

When results from multiple programs and submission codes are combined into a single dataset, these naming differences produce artificial duplicates, break group-level summaries, and prevent correct joins to reference tables. This vignette walks through a production pipeline for resolving those inconsistencies and producing a clean, analysis-ready dataset.

The pipeline covers:

  1. Analyte name harmonization across lab naming conventions
  2. Dissolved/total fraction parsing
  3. Analysis group classification (Conventional, Nutrients, Metals, etc.)
  4. Freshwater vs. seawater classification by specific conductance
  5. Detection limit flagging and half-DL substitution
  6. Unit standardization (µg/L → mg/L or ng/L)
  7. QAQC: blank screening, field duplicate RPD, holding-time exceedances, dissolved:total ratio check
  8. Field in situ data joined as wide cofactors (sp.cond drives water type; pH, DO, temperature retained for guideline calculations)

1. Data

wq_dat <- read_csv("data/raw/wq_raw.csv",
                   col_types = cols(
                     sample.date = col_date(),
                     result      = col_character(),  # may contain "<" prefix
                     hold.time   = col_character(),  # "T"/"F" must not be parsed as logical
                     dl          = col_double()
                   ))

insitu <- read_csv("data/raw/wq_insitu.csv",
                   col_types = cols(sample.date = col_date()))

In-situ field measurements (specific conductance, pH, DO, temperature, etc.) are recorded once per site visit and stored in a separate file. They must be joined onto the lab rows before water type classification (section 5), because water.type is determined from field specific conductance — not from the lab submission code.

# Left-join in-situ data so sp.cond.uS.cm is available for water type
# classification. Blank sites and duplicates receive NAs (no instrument reading),
# which is handled by the -FB/-TB guard in assign_water_type().
wq_dat <- left_join(wq_dat, insitu,
                    by = c("site.impute" = "site", "sample.date" = "sample.date"))

str(wq_dat)
## spc_tbl_ [1,470 × 29] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ program        : chr [1:1470] "Exposure" "Exposure" "Exposure" "Exposure" ...
##  $ site           : chr [1:1470] "EXP-01" "EXP-01" "EXP-01" "EXP-01" ...
##  $ site.impute    : chr [1:1470] "EXP-01" "EXP-01" "EXP-01" "EXP-01" ...
##  $ sample.id      : chr [1:1470] "EXP-01-0207" "EXP-01-0207" "EXP-01-0207" "EXP-01-0207" ...
##  $ sample.date    : Date[1:1470], format: "2023-02-07" "2023-02-07" ...
##  $ sample.time    : 'hms' num [1:1470] 09:00:00 09:00:00 09:00:00 09:00:00 ...
##   ..- attr(*, "units")= chr "secs"
##  $ matrix         : chr [1:1470] "Water" "Water" "Water" "Water" ...
##  $ sub.matrix     : chr [1:1470] "Freshwater" "Freshwater" "Freshwater" "Freshwater" ...
##  $ analyte        : chr [1:1470] "Aluminum, dissolved" "Aluminum, total" "Antimony, dissolved" "Antimony, total" ...
##  $ test.code      : chr [1:1470] "ICPMS" "ICPMS" "ICPMS" "ICPMS" ...
##  $ lab.id         : chr [1:1470] "Commercial Lab" "Commercial Lab" "Commercial Lab" "Commercial Lab" ...
##  $ result         : chr [1:1470] "17.47" "322.4" "1.518" "2.734" ...
##  $ dl             : num [1:1470] 2 5 0.1 0.2 0.1 0.2 0.5 1 0.02 0.05 ...
##  $ unit           : chr [1:1470] "µg/L" "µg/L" "µg/L" "µg/L" ...
##  $ qualifier      : logi [1:1470] NA NA NA NA NA NA ...
##  $ hold.time      : chr [1:1470] "F" "F" "F" "F" ...
##  $ source         : chr [1:1470] "Contractor" "Contractor" "Contractor" "Contractor" ...
##  $ duplicate      : chr [1:1470] NA NA NA NA ...
##  $ sp.cond.uS.cm  : num [1:1470] 619 619 619 619 619 619 619 619 619 619 ...
##  $ pH             : num [1:1470] 7.2 7.2 7.2 7.2 7.2 7.2 7.2 7.2 7.2 7.2 ...
##  $ sal.ppt        : num [1:1470] 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 ...
##  $ temp.C         : num [1:1470] 8.8 8.8 8.8 8.8 8.8 8.8 8.8 8.8 8.8 8.8 ...
##  $ diss.O2.mg.L   : num [1:1470] 9.7 9.7 9.7 9.7 9.7 9.7 9.7 9.7 9.7 9.7 ...
##  $ diss.O2.percent: num [1:1470] 84 84 84 84 84 84 84 84 84 84 ...
##  $ redox.ORP      : num [1:1470] 227 227 227 227 227 227 227 227 227 227 ...
##  $ turbidity.NTU  : num [1:1470] 8.5 8.5 8.5 8.5 8.5 8.5 8.5 8.5 8.5 8.5 ...
##  $ e.coordinate   : num [1:1470] 497200 497200 497200 497200 497200 ...
##  $ n.coordinate   : num [1:1470] 5448400 5448400 5448400 5448400 5448400 ...
##  $ depth.m        : num [1:1470] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   program = col_character(),
##   ..   site = col_character(),
##   ..   site.impute = col_character(),
##   ..   sample.id = col_character(),
##   ..   sample.date = col_date(format = ""),
##   ..   sample.time = col_time(format = ""),
##   ..   matrix = col_character(),
##   ..   sub.matrix = col_character(),
##   ..   analyte = col_character(),
##   ..   test.code = col_character(),
##   ..   lab.id = col_character(),
##   ..   result = col_character(),
##   ..   dl = col_double(),
##   ..   unit = col_character(),
##   ..   qualifier = col_logical(),
##   ..   hold.time = col_character(),
##   ..   source = col_character(),
##   ..   duplicate = col_character()
##   .. )
##  - attr(*, "problems")=<externalptr>

The problem is visible immediately: different programs submit the same compounds under different analyte names.

2. Analyte Name Harmonization

Each raw analyte name is mapped to a canonical name via a lookup table defined in R/functions/harmonize.R. The lookup covers common variants observed across standard lab submission codes; names not in the table pass through unchanged.

wq_dat <- harmonize_analyte_names(wq_dat)

The full set of harmonized analyte names:

3. Dissolved / Total Fraction Parsing

Metals reported as dissolved or total fractions carry that information as a suffix in the analyte name (e.g. “Copper, total”). The suffix is not part of the analyte identity — it describes the sample preparation — so it is parsed into a separate analysis column and stripped from the name.

wq_dat <- strip_fraction_suffix(wq_dat)

4. Analysis Group Classification

Analytes are assigned to reporting groups that determine which regulatory guidelines apply and how results are presented in the report. The groups are:

Group Examples
Conventional Parameters pH, TSS, DOC, hardness
Nutrients and Biological Indicators ammonia, nitrate, chlorophyll
Major Ions sodium, calcium, chloride
Total Metals copper total, lead total
Dissolved Metals copper dissolved, lead dissolved 
# Phosphorus from nutrient test codes (lower DL, nutrient panel) is mapped to
# "Total Phosphorus" in Nutrients rather than the metals-panel "Phosphorus".
# Adjust these codes to match your lab's submission structure.
phosphorus_nutrient_codes <- c("E372-U", "E372S")

wq_dat <- classify_analysis_group(wq_dat, phosphorus_nutrient_codes)

5. Water Type Classification

Freshwater and seawater samples require different analytical methods (different digestion matrices and detection limits). Classification uses specific conductance when available, with the lab submission code as a cross-check. Blank sites are always classified as freshwater regardless of conductance.

wq_dat <- assign_water_type(wq_dat)

6. Detection Limit Handling

Lab results below the detection limit are reported as character strings prefixed with < (e.g. "<0.001"). Before any arithmetic these must be:

  1. Flagged with a logical dl.f column so censored values can be excluded from parametric statistics.
  2. Substituted with ½ DL as a conventional placeholder for calculations where censored results cannot simply be dropped.
wq_dat <- flag_below_dl(wq_dat)

7. Unit Standardization

Labs report metals in µg/L and nutrients in mg/L. For consistent cross-analyte presentation all results are converted to a common unit per analyte class. Mercury is converted from µg/L to ng/L because its detection threshold is at sub-nanogram concentrations.

Detection limits are scaled by the same conversion factor as results.

wq_dat <- standardize_units(wq_dat)

8. QAQC

8a. Sample hash

A sample.hash uniquely identifies each station–date combination and is used in all subsequent QAQC joining operations.

wq_dat <- wq_dat %>%
  mutate(sample.hash = paste(site.impute, sample.date, sep = "|"))

8b. Analytical Deduplication

When the same analyte appears twice for the same sample (submitted under two test codes, or to two labs), the row with the lower detection limit is retained — the finer threshold is more informative.

dedup_group_vars <- c(
  "site.impute", "sample.date", "sample.id", "water.type", "sample.type",
  "actual.test", "test.code", "duplicate", "analysis", "analyte.impute", "unit.impute"
)

n_before <- nrow(wq_dat)
wq_dat   <- deduplicate_wq(wq_dat, dedup_group_vars)
n_after  <- nrow(wq_dat)

cat(sprintf("Removed %d duplicate rows (%d → %d)\n", n_before - n_after, n_before, n_after))
## Removed 0 duplicate rows (1470 → 1470)

8c. Blank Screening

Field and travel blanks should return results below detection. A detection exceeding 5× DL indicates probable contamination during sampling (field blank) or transport (travel blank) and warrants investigation before the paired environmental samples are used.

blank_results <- screen_blanks(wq_dat)

8d. Field Duplicate RPD

Field duplicates are collected by re-sampling the same location immediately after the primary sample. RPD (relative percent difference) quantifies analytical and sampling precision.

  • RPD = 2 |A - B| / (A + B) × 100
  • Pairs where both results are below DL → reported as <DL
  • Pairs where either result is below 5× DL → NC (not calculable; precision cannot be assessed near the detection threshold)
rpd_table <- compute_rpd(wq_dat)

8e. Holding-Time Exceedances

Some analytes degrade rapidly after collection. Results from samples that exceeded the maximum allowable holding time are flagged; they may still be usable depending on the degree of exceedance and the analyte’s known degradation rate.

wq_dat       <- parse_holding_time(wq_dat)
hold_summary <- holding_time_summary(wq_dat)

8f. Dissolved : Total Metals Ratio

Dissolved metal concentration cannot substantially exceed total metal concentration: dissolved is a sub-fraction of total. A ratio > 1.5 is a red flag for a labelling error, a sample swap, or an analytical matrix effect.

dt_flags <- check_dt_ratio(wq_dat)

9. Week Labelling

Sampling events are grouped into sequential calendar weeks to support week-by-week comparisons in the report.

wq_dat <- label_weeks(wq_dat, n_weeks = 5)

10. Final Dataset

Environmental samples only (blanks and duplicates are presented in the QAQC section). In situ field measurements are retained as wide columns on each lab row — they serve as per-sample cofactors for guideline calculations (e.g. hardness-dependent criteria) and informed the water type classification in section 5. Effluent samples are classified as "Effluent" regardless of conductance.

insitu_wide_cols <- c("pH", "sp.cond.uS.cm", "sal.ppt", "temp.C",
                      "diss.O2.mg.L", "diss.O2.percent", "redox.ORP", "turbidity.NTU")

wq_out <- wq_dat %>%
  filter(is.na(sample.type)) %>%
  select(
    program, sample.id, site.impute,
    sample.date, sample.time, week,
    water.type, actual.test,
    analysis, analyte.impute, test.code,
    result.impute, unit.impute, dl.impute, dl.f,
    source, lab.id, hold.time,
    all_of(insitu_wide_cols),
    e.coordinate, n.coordinate, depth.m
  ) %>%
  mutate(analysis = factor(analysis,
                           levels = analysis_levels[analysis_levels != "Field Measured"])) %>%
  arrange(program, site.impute, sample.date, analysis, analyte.impute)
saveRDS(wq_out, file = "data/processed/wq_harmonized.RDS")

This vignette uses synthetic data generated by data/generate_synthetic_data.R. The pipeline structure reflects work across multiple environmental monitoring programs; all program names, site identifiers, and laboratory references have been replaced with generic labels.