Building a Full Air Quality QC Pipeline with pandas and PyOD

This recipe assembles the four quality stages an air-quality network needs — range validation, QC flagging on a 0–3 scale, drift correction, and PyOD anomaly scoring — into one stateless class that turns raw PM2.5 telemetry into an analysis-ready DataFrame. The ordering matters as much as the algorithms: range checks and flags run on the raw signal so faults are caught on their own terms, drift correction touches only the readings that survived, and the anomaly detector trains solely on rows already marked valid. The result is a single frame carrying the corrected value, a QC flag, and an anomaly score per observation, ready for spatial storage and export. It draws on the method choices explained in the anomaly detection methods overview.

Why Stage Order Governs Correctness

Each stage depends on the invariants the previous one establishes, and getting the order wrong corrupts the output silently rather than loudly. If you correct drift before flagging missing data, a communication dropout enters the rolling baseline as a run of implied-zero readings and permanently biases the correction offset. If you fit the PyOD detector on the full frame including known-bad rows, the model learns that failures are normal and goes blind to the next one. And if you range-check after drift correction, a genuine out-of-range spike may be pulled back inside the plausible envelope by the baseline subtraction and escape the flag it deserved. The pipeline therefore enforces a strict sequence: raw range check → QC flags → drift correction on valid rows → anomaly scoring on clean rows.

The QC flag column is the connective tissue. Every stage both reads and writes it: the range check sets 2 for hard violations and 3 for missing, drift correction masks flags 2 and 3 out of its rolling window, and the anomaly scorer trains only where the flag is 0. A single integer column, not scattered boolean masks, keeps that contract auditable.

Production-Ready Implementation

The AirQualityQC class below is self-contained. Instantiate it with the physical bounds for your sensor, call process on a raw PM2.5 frame, and receive a new frame with pm25_corrected, qc_flag, and anomaly_score columns. It never mutates its input.

# python 3.11 · pandas==2.2.2 · numpy==1.26.4 · scikit-learn==1.4.2 · pyod==1.1.3
from __future__ import annotations
from dataclasses import dataclass
import pandas as pd
import numpy as np
from pyod.models.ecod import ECOD


@dataclass
class AirQualityQC:
    """
    End-to-end QC for PM2.5 telemetry: range check -> QC flags (0-3)
    -> drift correction -> PyOD anomaly scoring.

    Flags follow EPA AQS convention: 0 valid, 1 questionable, 2 invalid,
    3 missing. The pipeline is stateless and deterministic given pinned
    library versions.
    """
    value_col: str = "pm25"
    time_col: str = "timestamp"
    hard_min: float = -5.0        # small negative tolerance for sensor noise
    hard_max: float = 1000.0      # µg/m³ — above this the optical sensor saturates
    soft_min: float = 0.0
    soft_max: float = 500.0
    drift_window: str = "12h"
    contamination: float = 0.01

    def process(self, df: pd.DataFrame) -> pd.DataFrame:
        df = df.copy()
        df[self.time_col] = pd.to_datetime(df[self.time_col], utc=True)
        df = df.sort_values(self.time_col).reset_index(drop=True)
        df = self._range_and_flags(df)
        df = self._correct_drift(df)
        df = self._score_anomalies(df)
        df["cal_version"] = "qc-v1"
        df["processed_at"] = pd.Timestamp.now(tz="UTC")
        return df

    # --- Stage 1 + 2: range check and QC flags on the RAW signal ---
    def _range_and_flags(self, df: pd.DataFrame) -> pd.DataFrame:
        v = df[self.value_col]
        flags = pd.Series(0, index=df.index, dtype=int)
        flags[v.isna()] = 3                                   # missing
        hard = (v < self.hard_min) | (v > self.hard_max)
        flags[hard & (flags == 0)] = 2                        # invalid
        soft = (v < self.soft_min) | (v > self.soft_max)
        flags[soft & (flags == 0)] = 1                        # questionable
        df["qc_flag"] = flags
        return df

    # --- Stage 3: drift correction, masking non-valid rows ---
    def _correct_drift(self, df: pd.DataFrame) -> pd.DataFrame:
        # Only readings that passed hard checks feed the rolling baseline.
        usable = df[self.value_col].where(df["qc_flag"].isin([0, 1]))
        baseline = (
            usable.set_axis(df[self.time_col])
            .rolling(self.drift_window, min_periods=6)
            .mean()
            .to_numpy()
        )
        anchor = np.nanmin(baseline) if np.isfinite(baseline).any() else 0.0
        drift_offset = baseline - anchor
        corrected = df[self.value_col].to_numpy() - drift_offset
        df["pm25_corrected"] = np.where(df["qc_flag"] == 2, np.nan, corrected)
        return df

    # --- Stage 4: PyOD scoring, trained only on clean rows ---
    def _score_anomalies(self, df: pd.DataFrame) -> pd.DataFrame:
        clean = df.loc[df["qc_flag"] == 0, ["pm25_corrected"]].dropna()
        df["anomaly_score"] = np.nan
        if len(clean) >= 50:
            det = ECOD(contamination=self.contamination)
            det.fit(clean.to_numpy())
            scorable = df["pm25_corrected"].notna()
            X = df.loc[scorable, ["pm25_corrected"]].to_numpy()
            df.loc[scorable, "anomaly_score"] = det.decision_function(X)
        return df

A few design choices are load-bearing. Drift correction anchors to the minimum rolling baseline rather than the first value, so a startup transient does not define the zero point — this mirrors the approach in correcting temperature sensor drift using rolling averages, adapted for PM2.5’s faster baseline shifts. The detector is ECOD because it is parameter-free and edge-friendly; swapping in IForest from pyod.models.iforest is a one-line change when you feed it several channels and want joint-implausibility scoring, as weighed in Isolation Forest vs Z-Score for sensor anomaly detection.

Running it 

qc = AirQualityQC(drift_window="12h", contamination=0.02)
out = qc.process(raw_pm25_df)   # columns: timestamp, pm25, latitude, longitude
print(out[["timestamp", "pm25", "pm25_corrected", "qc_flag", "anomaly_score"]].head())

The output frame is the input contract for downstream spatial work. Because it carries latitude, longitude, qc_flag, and anomaly_score per row, it drops straight into exporting QC-flagged data to GeoJSON for QGIS for map-based review.

Parameter Tuning Guide

The defaults suit a mid-range optical PM2.5 sensor. Adjust per hardware class and reporting cadence.

Sensor / context drift_window contamination hard_max (µg/m³) Notes
Low-cost optical PM2.5 (1 min) 6h–12h 0.02 1000 Faster baseline shift from lens contamination
Reference-grade PM2.5 (1 h) 24h–48h 0.005 2000 Stable optics; wider window avoids over-correction
PM10 (optical) 12h 0.02 2000 Coarse fraction; raise hard_max
Indoor PM2.5 6h 0.03 500 Cooking spikes are real events, not faults — gate them
Wildfire-season deployment 12h 0.05 5000 Expect legitimate extremes; loosen bounds and contamination

Rule of thumb: set drift_window to one to two diurnal cycles and contamination to the fault fraction you actually observe in maintenance logs. Do not raise contamination to silence alarms — that hides events; instead pass scores through the spatial coherence gate.

Verification and Testing

The pipeline’s correctness rests on stage ordering and flag semantics, so the tests target exactly those invariants: a hard-out-of-range value must be flag 2 with a null corrected value, a missing value must be flag 3, and drift must be measurably reduced on the valid rows.

import numpy as np
import pandas as pd


def _synthetic_pm25(hours: int = 48) -> pd.DataFrame:
    ts = pd.date_range("2026-01-01", periods=hours * 60, freq="1min", tz="UTC")
    base = 12 + 5 * np.sin(np.linspace(0, hours / 12 * np.pi, len(ts)))
    drift = np.linspace(0, 8, len(ts))          # +8 µg/m³ lens-contamination drift
    pm25 = base + drift + np.random.default_rng(1).normal(0, 1, len(ts))
    pm25[100] = 5000.0                            # hard-out-of-range spike
    pm25[200] = np.nan                            # missing
    return pd.DataFrame({"timestamp": ts, "pm25": pm25,
                         "latitude": 47.6, "longitude": -122.3})


def test_flags_and_drift_reduction():
    out = AirQualityQC(drift_window="12h").process(_synthetic_pm25())

    assert out.loc[100, "qc_flag"] == 2                    # hard violation
    assert np.isnan(out.loc[100, "pm25_corrected"])        # invalid -> null
    assert out.loc[200, "qc_flag"] == 3                    # missing

    valid = out[out["qc_flag"].isin([0, 1])]
    # Corrected series should have less end-to-end trend than the raw series
    raw_trend = valid["pm25"].tail(500).mean() - valid["pm25"].head(500).mean()
    cor_trend = valid["pm25_corrected"].tail(500).mean() - valid["pm25_corrected"].head(500).mean()
    assert abs(cor_trend) < abs(raw_trend), "drift correction should reduce trend"


def test_scores_only_on_valid_rows():
    out = AirQualityQC().process(_synthetic_pm25())
    # The invalid row carries a flag but no anomaly score
    assert np.isnan(out.loc[100, "anomaly_score"])
    assert out["anomaly_score"].notna().sum() > 0

For production promotion, extend these with a determinism check: run process twice on the same input and assert the qc_flag and anomaly_score columns are identical, guaranteeing reproducibility for later re-analysis.

Gotchas

Correcting drift before flagging outages. A gap flagged 3 must be masked out of the rolling baseline. Feed it in as implied data and the baseline sags, over-correcting every reading after the outage. The pipeline masks flags 2 and 3 before computing the baseline for exactly this reason.

Training PyOD on the whole frame. Including flag 1–3 rows in fit teaches the detector that questionable and invalid readings are normal. Always fit on qc_flag == 0 only, then score the rest.

Anchoring drift to the first value. If the sensor powers on mid-spike, the first rolling value is a bad zero point. Anchor to the minimum stable baseline instead, as the class does with np.nanmin.

Treating scores as a fault verdict. anomaly_score marks statistical unusualness, not a broken sensor. An indoor cooking spike or a wildfire plume scores high but is a real event. Route scores through the spatial coherence gate before alerting.