TwQualityControl

Import, standardize, QA/QC, and summarize water-temperature time series for stream monitoring programs. This R package was developed to support data processing for the PSF Stream Temperature Database, a centralized repository of stream temperature records across British Columbia managed by the Pacific Salmon Foundation.

• Package: TwQualityControl
• License: MIT
• Version: 0.0.0.9000

---

Installation

This package is not on CRAN. Install from GitHub:

# install.packages("remotes")
remotes::install_github("salmonwatersheds/tw_quality_control")

Note: R build tools may be required (Rtools on Windows, Xcode CLI tools on macOS). Dependencies installed automatically: dplyr, lubridate, zoo, stats, graphics, grDevices.

---

Quick Start

A complete copy-paste example showing the full pipeline:

library(TwQualityControl)

# --- Generate some sample data ---
set.seed(1)
n <- 200
dt0 <- as.POSIXct("2024-07-01 00:00:00", tz = "America/Vancouver")
raw <- data.frame(
  datetime = dt0 + seq(0, by = 3600, length.out = n),
  tw_value = 15 + sin(seq_len(n) / 6) + rnorm(n, 0, 0.2)
)
# Inject problems for QC to detect:
raw$tw_value[20]    <- raw$tw_value[19] + 4  # spike
raw$tw_value[30:37] <- 12                    # flatline
raw$tw_value[100]   <- NA                    # missing value

# --- Step 1: Format ---
fmt <- format_tw_timeseries(raw, verbose_print = FALSE)

# --- Step 2: QA/QC ---
qa <- qaqc_tw_timeseries(fmt$timeseries,
                          sampling_interval_mins = 60,
                          verbose_print = FALSE)

# --- Step 3: Visualize flags ---
plot_qaqc_flags(qa)

# --- Step 4: Summarize ---
daily   <- daily_timeseries_summary(qa)
monthly <- monthly_timeseries_summary(qa)
yearly  <- yearly_timeseries_summary(qa)

---

Workflow

The package provides six exported functions used in sequence. Raw data is formatted and standardized, passed through automated QA/QC checks, and then summarized at daily, monthly, and yearly time scales for analysis and database upload.

---

1. format_tw_timeseries()

Ingests raw data and returns a standardized list ready for QA/QC. Because stream temperature data arrives from many different sources (government agencies, community groups, research programs), column names, datetime formats, and timezones vary widely. This function auto-detects these properties when not specified, normalizing everything into a consistent format.

Key arguments:

Argument	Description
dat	Input data.frame
colname_datetime	Datetime column name (auto-detected if NULL)
colname_tw_value	Temperature column name (auto-detected if NULL)
datetime_format	Format string, e.g. "%Y-%m-%d %H:%M:%S" (guessed if NULL)
input_tz	Olson timezone, e.g. "America/Vancouver" (guessed if NULL)
timestep_mins	Sampling interval in minutes (estimated if NULL)
from_year, to_year	Expected year range for validation
verbose_print	Print diagnostics and generate plots

Returns a named list with two elements:

• station -- one-row data.frame with start_year, end_year, timezone, sampling_interval_mins
• timeseries -- data.frame with datetime (POSIXct) and tw_value (numeric)

Important: All timestamps are converted to UTC before return, regardless of the original timezone. This is by design -- the database stores everything in UTC.

fmt <- format_tw_timeseries(
  dat = raw,
  colname_datetime = "datetime_local",
  colname_tw_value = "water_temperature",
  input_tz         = "America/Vancouver",
  verbose_print    = TRUE
)

fmt$station
#>   start_year end_year       timezone sampling_interval_mins
#> 1       2024     2024 America/Vancouver                   15

head(fmt$timeseries)
#>              datetime tw_value
#> 1 2024-07-01 15:00:00     12.1
#> 2 2024-07-01 15:15:00     12.3
#> ...

---

2. qaqc_tw_timeseries()

Appends seven QA/QC flag columns to the timeseries. Each flag is one of "P" (pass), "S" (suspect), "F" (fail), or "X" (not evaluated). Automated flagging catches common data quality issues -- sensor drift, logger malfunctions, air exposure events -- that would otherwise bias summary statistics and thermal metric calculations.

Key arguments:

Argument	Description
timeseries	data.frame from format_tw_timeseries()$timeseries
sampling_interval_mins	Sampling interval; estimated from data if NULL
fill_missing_timestamps	Insert rows for expected-but-missing timestamps (TRUE/FALSE)
evaluation_window_size_mins	Rolling window size for SD-based checks (default 300 = 5 hours)
low_thres_suspect, low_thres_fail	Low temperature thresholds (default 0, -0.2)
high_thres_suspect, high_thres_fail	High temperature thresholds (default 30, 35)
spike_suspect, spike_fail	Spike step-change thresholds (default 2, 3 °C)
flatline_suspect, flatline_fail	Flatline rolling-SD thresholds (default 0.02, 0 °C)
changerate_suspect, changerate_fail	Change-rate rolling-SD thresholds (default 2, 3 °C)
verbose_print	Print flag counts and generate the QC plot

Returns the input data.frame with seven additional flag columns.

qa <- qaqc_tw_timeseries(
  timeseries = fmt$timeseries,
  sampling_interval_mins = 60,
  verbose_print = TRUE
)

How spike, flatline, and change rate checks work

All three checks operate on the raw temperature values (tw_value) and use the evaluation_window_size_mins parameter to determine the rolling window size k (number of data points in the window). For example, with evaluation_window_size_mins = 300 (5 hours) and sampling_interval_mins = 60 (hourly data), k = 300 / 60 = 5 points.

Spike detection (spike_suspect = 2, spike_fail = 3)

Spikes are detected by computing the absolute step-change between each point and its immediate neighbours. For each observation at index i, the function computes:

• dx_prev = |tw_value[i] - tw_value[i-1]| (step from previous point)
• dx_next = |tw_value[i+1] - tw_value[i]| (step to next point)

If either step-change >= spike_fail (default 3 °C), the row is flagged "F". If either >= spike_suspect (default 2 °C), it is flagged "S". Because both the forward and backward differences are checked, a single-point spike flags 3 rows: the point before, the spike itself, and the point after.

Example: Suppose five consecutive hourly readings are [14.0, 14.2, 18.5, 14.1, 14.3]. The point at index 3 (18.5 °C) has dx_prev = |18.5 - 14.2| = 4.3 and dx_next = |14.1 - 18.5| = 4.4. Both exceed the fail threshold of 3 °C. Additionally, index 2 sees dx_next = 4.3 >= 3 and index 4 sees dx_prev = 4.4 >= 3, so all three rows (indices 2, 3, and 4) receive an "F" flag.

Flatline detection (flatline_suspect = 0.02, flatline_fail = 0)

Flatlines indicate periods where the sensor is reporting little or no variation, which can indicate a malfunctioning logger, a sensor out of water, or frozen conditions. The check computes the trailing rolling population SD over the last k points (the window defined by evaluation_window_size_mins).

• If rolling_sd <= flatline_fail (default 0 °C, i.e., zero variance), the row is flagged "F".
• If rolling_sd <= flatline_suspect (default 0.02 °C), the row is flagged "S".

Example: With hourly data and a 5-hour window (k = 5), suppose the last 5 readings are [12.0, 12.0, 12.0, 12.0, 12.0]. The rolling SD is 0.0, which is <= flatline_fail (0), so the row is flagged "F". If the readings were [12.00, 12.01, 12.00, 12.01, 12.00], the rolling SD is ~0.005, which is <= flatline_suspect (0.02), so the row is flagged "S".

Change rate detection (changerate_suspect = 2, changerate_fail = 3)

Change rate flags periods of unusually high variability, which can indicate sensor interference, air exposure, or data corruption. The check computes rolling population SD over both a trailing window (last k points) and a leading window (next k points).

• If either the trailing or leading rolling SD >= changerate_fail (default 3 °C), the row is flagged "F".
• If either >= changerate_suspect (default 2 °C), the row is flagged "S".

Example: With hourly data and a 5-hour window (k = 5), suppose the last 5 readings are [10.0, 15.0, 8.0, 16.0, 9.0]. The rolling SD of these values is ~3.3 °C, which exceeds changerate_fail (3 °C), so the row is flagged "F". Using both trailing and leading windows ensures that erratic behaviour is caught regardless of whether it occurs just before or just after the current observation.

Note on rolling SD: The rolling standard deviation used for flatline and change rate is the population SD (sqrt(E[x^2] - E[x]^2), divides by k), not the sample SD that R's sd() returns (divides by k-1). The population SD is always slightly smaller: population_sd = sample_sd * sqrt((k-1)/k). For the default 5-hour window this difference is negligible with sub-hourly data (k >= 20, < 3%) but more pronounced with hourly data (k = 5, ~11%). Keep this in mind when comparing against sd() computed externally or when tuning threshold values.

---

3. plot_qaqc_flags()

Visualizes the flagged timeseries with shaded regions for suspect (grey) and fail (yellow) flags. This provides a quick visual check to confirm that automated flags are targeting real data quality issues and not flagging legitimate temperature variability. Can be called independently to re-visualize without re-running QC.

plot_qaqc_flags(qa)

You can also pass a subset of flag columns:

plot_qaqc_flags(qa, flag_cols = c("spike", "flatline"))

---

4. daily_timeseries_summary()

Computes daily statistics from a QC'd timeseries. Use the flag parameters to control which records are included. Missing dates within the record are filled with NA so that rolling-window calculations (such as the 7-day average daily mean) are computed correctly without bridging data gaps.

Output columns: date, year, month, day, julian_day, minimum, maximum, mean, weekly_adm (7-day rolling average daily mean), median, range, standard_deviation, count_records

daily <- daily_timeseries_summary(
  timeseries = qa,
  time_dup  = c("P", "S"),   # include pass and suspect
  missing   = c("P", "S"),
  spike     = c("P", "S"),   # exclude fail flags for spikes
  verbose_print = TRUE
)

---

5. monthly_timeseries_summary()

Aggregates the timeseries to monthly statistics, including the day of the month when monthly minimum and maximum temperatures were recorded. Monthly summaries are useful for identifying seasonal patterns and comparing temperature regimes across years.

Output columns: year, month, minimum, day_of_month_minimum, maximum, day_of_month_maximum, mean, median, range, standard_deviation, count_records

monthly <- monthly_timeseries_summary(timeseries = qa)

---

6. yearly_timeseries_summary()

Computes annual thermal ecology metrics from a QC'd timeseries. These metrics are widely used in fisheries science to characterize thermal habitat suitability for salmon and other cold-water species -- for example, MWAT and MWMT are standard indicators of summer thermal stress.

Output columns:

Column	Description
year	Year
minimum, maximum, mean, median	Annual temperature statistics
day_minimum, day_maximum	Julian day of annual min/max
range, standard_deviation	Annual range and SD
count_records	Number of records used
max_weekly_avg_temp	Max 7-day rolling average of daily mean (MWAT)
max_weekly_max_temp	Max 7-day rolling average of daily max (MWMT)
max_weekly_min_temp	Max 7-day rolling average of daily min
aug_mean_daily_mean	Mean of August daily means
seven_day_equivalent_constant	7DEC: median of MWAT and MWMT
days_above_15C	Days where daily max exceeds 15 C
days_above_19C	Days where daily max exceeds 19 C
daily_mean_95th_pct	95th percentile of daily mean temperatures
max_daily_range, min_daily_range	Extremes of daily temperature range

The weekly_window_method parameter controls how the 7-day rolling window is computed: "trailing" (default, current day + 6 prior) or "centered" (3 days either side).

yearly <- yearly_timeseries_summary(
  timeseries = qa,
  weekly_window_method = "trailing",
  spike = c("P", "S"),       # exclude spike-flagged records
  flatline = c("P", "S"),
  verbose_print = TRUE
)

---

QA/QC Flag System

Flag values

Flag	Meaning	When assigned
P	Pass	Value passes the check
S	Suspect	Value exceeds a suspect threshold but not fail
F	Fail	Value exceeds a fail threshold
X	Not evaluated	Prerequisite missing (e.g., tw_value is NA)

Flag columns

Column	Check	Suspect trigger	Fail trigger
time_dup	Duplicate timestamps	--	Duplicate detected
missing	Missing temperature	--	tw_value is NA
low_thres	Below low threshold	Below low_thres_suspect (0 °C)	Below low_thres_fail (-0.2 °C)
high_thres	Above high threshold	Above high_thres_suspect (30 °C)	Above high_thres_fail (35 °C)
spike	Step-change to either neighbour (flags the point and its neighbours)	Step >= spike_suspect (2 °C)	Step >= spike_fail (3 °C)
flatline	Low rolling SD	SD <= flatline_suspect (0.02 °C)	SD <= flatline_fail (0 °C)
changerate	High rolling SD (trailing or leading)	SD >= changerate_suspect (2 °C)	SD >= changerate_fail (3 °C)

Using flags in summaries

The daily, monthly, and yearly summary functions accept the same seven flag parameters. Each parameter takes a character vector of flags to include. Rows with flags not in the vector are excluded from the summary.

# Strict filtering: only use records that passed all checks
yearly_strict <- yearly_timeseries_summary(
  timeseries = qa,
  time_dup   = c("P"),
  missing    = c("P"),
  low_thres  = c("P"),
  high_thres = c("P"),
  spike      = c("P"),
  flatline   = c("P"),
  changerate = c("P")
)

# Lenient filtering: include everything except fail
yearly_lenient <- yearly_timeseries_summary(
  timeseries = qa,
  time_dup   = c("P", "S"),
  missing    = c("P", "S"),
  low_thres  = c("P", "S"),
  high_thres = c("P", "S"),
  spike      = c("P", "S"),
  flatline   = c("P", "S"),
  changerate = c("P", "S")
)

---

Notes

Timezone handling: format_tw_timeseries() converts all timestamps to UTC before returning. If input_tz is not specified, the function attempts to guess whether the data is in local time or UTC using three heuristics (summer diurnal peak timing, November DST fall-back, March DST spring-forward). When local time is detected, it is assumed to be "America/Vancouver". For stations outside the Pacific Northwest, always supply input_tz explicitly.

Data without QC flags: The summary functions can accept data without QC flag columns (e.g., a raw data.frame with just datetime and tw_value). Missing flag columns are automatically created with all "P" values, so no records are filtered out.

Timestamp rounding: If qaqc_tw_timeseries() detects that >50% of expected timestamps appear missing, it rounds timestamps to the nearest minute before retrying. A diagnostic message is printed when this occurs. Sub-minute precision is not expected in this data.

---

Acknowledgements

Developed by Matthew Bayly (M.J. Bayly Analytics Ltd.) in collaboration with the Pacific Salmon Foundation.