AQEval Python Package

A Python port of the R package AQEval for air quality time-series analysis.

Overview

AQEval provides tools for the routine detection, characterization, and quantification of discrete changes in air quality time-series, such as identifying the impacts of air quality policy interventions.

The main functions use signal isolation and break-point/segment (BP/S) methods to detect and quantify change events.

Methodology and R parity

This port reproduces the R package's methodology step for step (AQEval v0.6.12; Ropkins et al., JOSS 2026, 10.21105/joss.08839):

1. Find possible points-of-change — find_break_points is a faithful NumPy port of strucchange::breakpoints(y ~ 1, h): the Bai & Perron dynamic-programming segmentation of the intercept-only (mean-shift) model, with the number of breaks selected by BIC and break-date confidence intervals computed with Bai's (1997) asymptotic distribution.
2. Test them — test_break_points fits an OLS model for every subset of the candidate breaks (single shared intercept, per-segment epoch-second regressors, exactly as R's aqe_fitBreakPointsModel); a model is valid only if every coefficient is significant at p < 0.05, and the valid model with the highest adjusted R² wins.
3. Quantify — quant_break_points reports the fitted concentrations either side of each confirmed break; quant_break_segments ports the segmented-package machinery vendored by R AQEval (single Muggeo iteration, it.max = 1, no bootstrap) and grid-searches segmented fits started around the break confidence intervals, quantifying regions-of-change rather than instantaneous steps.
4. Signal isolation — isolate_contribution fits the R default GAM pollutant ~ te(wd,ws) + s(year.day) + s(day.hour) (plus optional background/extra terms). R uses mgcv; this port uses pygam (penalized P-splines with staged GCV lambda selection), the closest maintained pure-Python equivalent. Results agree with mgcv closely but not bit-for-bit (about ±2 µg/m³ on the hourly test fixture); everything else in the package matches the R outputs exactly.

The port is verified two ways:

1. Against the R package's own regression fixtures (inst/tinytest/): break locations and confidence intervals (213/317/402 on 2-day-averaged aq.data), quantified break rows, segment boundaries, date-range statistics and spectral output dimensions all reproduce the R values exactly (tests/test_regression_r_fixtures.py).

2. Head-to-head against a live R AQEval run (R 4.6.1 with strucchange/segmented/mgcv, sourcing the R package code) on the same bundled aq.data:

   Comparison	Result
   Break-points found (2-day, h=0.3; 3-day, h=0.15)	Identical rows, incl. confidence intervals
   Segment boundaries (all rows, all 6 columns)	Identical
   Break-point model trend	r = 1.000000, max diff 8×10⁻⁵ µg/m³
   Break-segment (Muggeo) trend	r = 1.000000, max diff 5×10⁻¹² µg/m³
   Spectral density (13,500 frequencies)	r = 1.000000, max diff 2×10⁻¹⁵
   Isolated signal (3,600 hourly values)	r = 0.992 (mgcv vs pygam)

Conventions carried over from R: break/segment tables hold 1-based row numbers into the input data frame, and report columns use R's dotted names (date.low, c.delta, per.delta, seg.str, ...). As in R, break-point analysis is designed for time-averaged series — resample hourly data first (e.g. time_average(data, "day")); signal isolation runs at the native hourly resolution before averaging.

Installation

# Install from PyPI
pip install aqeval



# Install in development mode
pip install -e .

# Or install dependencies only
pip install -r requirements.txt

Quick Start

from aqeval import (
    load_aq_data,
    isolate_contribution,
    quant_break_points,
    quant_break_segments,
    spectral_frequency,
)
from aqeval.utilities import time_average

# Load sample data
data = load_aq_data()

# Resample to daily averages (openair::timeAverage semantics)
daily = time_average(data, avg_time="day")

# 1. Signal Isolation - Remove seasonal and weather effects
data["dswb_no2"] = isolate_contribution(
    data, 
    "no2", 
    background="bg.no2"
)

# 2. Break-Point Analysis
result = quant_break_points(daily, "no2", h=0.3)
result["plot"].show()  # Interactive Plotly plot

# 3. Break-Segment Analysis
result = quant_break_segments(daily, "no2", h=0.3)
result["plot"].show()

# 4. Spectral Frequency Analysis
result = spectral_frequency(data, "no2")
result["plot"].show()

Main Functions

Signal Processing

• isolate_contribution(): Remove background, seasonal, and weather effects from pollutant time-series using GAM models (pygam; R uses mgcv).
• fit_near_site_model(): Gap-fill one site's series using GAM models built from nearby sites.

Break-Point Detection

• find_break_points(): Detect candidate break-points with the strucchange method (Bai–Perron dynamic programming + BIC + Bai 1997 confidence intervals).
• test_break_points(): Test all break-point subsets; all-coefficients-significant models ranked by adjusted R².
• quant_break_points(): Quantify confirmed break-points with confidence intervals, report and plot.
• quant_break_segments(): Quantify regions-of-change about break-points via Muggeo-style segmented regression seeded from the break confidence intervals.

Analysis Tools

• spectral_frequency(): Perform spectral frequency analysis to identify periodic patterns.
• calc_date_range_stat(): Calculate statistics for specific date ranges.
• calc_rolling_date_range_stat(): Calculate rolling statistics over moving windows.

Site Finding

• find_near_sites(): Find nearby air quality monitoring sites by location.
• find_near_lat_lon(): General-purpose nearest location finder.

Data Access

• load_aq_data(): Load the sample AQEval dataset.
• import_aq_meta(): Import metadata from UK air quality networks.
• download_aurn_data(): Download data from UK air quality networks.
• download_noaa_data(): Download meteorological data from NOAA.

Example Workflow

import pandas as pd
from aqeval import (
    load_aq_data,
    isolate_contribution,
    quant_break_points,
)
from aqeval.utilities import time_average

# 1. Load and prepare data
data = load_aq_data()

# 2. Apply signal isolation (deseasonalization + deweathering + background)
data["isolated_no2"] = isolate_contribution(
    data,
    pollutant="no2",
    background="bg.no2",
    deseason=True,
    deweather=True,
)

# 3. Resample to 14-day averages for smoother analysis
avg_14day = time_average(data, avg_time="14 day")

# 4. Detect and quantify break-points
result = quant_break_points(
    avg_14day, 
    "isolated_no2", 
    h=0.1,
    event={
        "x": "2020-03-23",
        "label": "COVID Lockdown",
        "color": "grey",
    }
)

# 5. View results
print(result["report"])
result["plot"].show()

Tutorial example: BDMA 2018-2023

examples/tutorial_bdma.py replicates the AQEval R tutorial end-to-end with this package:

from aqeval import (download_aurn_data, isolate_contribution,
                    time_average, find_break_points, quant_break_segments)

# 1. AURN data for the BDMA site (R: importAURN(site='bdma', year=2018:2023))
data = download_aurn_data("bdma", 2018, 2023, source="aurn")

# 2. De-weather and de-season NO2, with air temperature as an
#    extra background term (fits no2 ~ s(air_temp)+te(wd,ws)+s(year.day)+s(day.hour))
data["deweatherdeseason"] = isolate_contribution(
    data, "no2", deseason=True, deweather=True, background="air_temp")

# 3. Average to 8-hour resolution (openair::timeAverage equivalent)
data_8h = time_average(data, "8 hour")

# 4. Find and quantify the break segments (h=0.3 for speed; h=0.12 for
#    a full-resolution 8-hour analysis)
break_points = find_break_points(data_8h, "deweatherdeseason", h=0.3)
quant = quant_break_segments(data_8h, "deweatherdeseason", breaks=break_points)
print(quant["report"])

Cross-checked against the identical workflow run in R (each side using its own GAM backend for the isolation step), the two implementations find the same change events:

The isolated hourly series correlate at r = 0.9967 and the fitted segment trends at r = 0.9994 (mean difference 0.06 µg/m³).

Note on performance: like R's strucchange, the break-point search is O(n²) — run it on time-averaged data (8-hour/daily), not raw hourly series. On the 8-hour BDMA series (6,573 rows) h=0.3 takes seconds; smaller h values take longer.

Visualization

All visualization functions return interactive Plotly figures that can be:

• Displayed in Jupyter notebooks
• Saved as HTML files: fig.write_html("plot.html")
• Exported as static images: fig.write_image("plot.png")

Figures Module

The figures module provides a menu of ready-to-use plotting functions. Each function accepts a DataFrame (with a date column and one or more pollutant columns), runs the necessary analysis internally, and returns a single plotly.graph_objects.Figure.

from aqeval.figures import (
    plot_break_points,
    plot_segments_raw,
    plot_segments_isolated,
    plot_raw_vs_isolated_trends,
    plot_raw_vs_isolated_timeseries,
    plot_residuals_timeseries,
    plot_residual_histogram,
    plot_spectral,
    plot_rolling_statistics,
    plot_signal_isolation,
    plot_signal_isolation_trend_comparison,
)

Data preparation

Most figure functions expect hourly data so that signal-isolation (deseasoning / deweathering) can operate at full resolution. Functions that only need daily data (e.g. plot_rolling_statistics) will aggregate internally.

You can use the bundled sample dataset or download live data from a UK air-quality network:

# Option 1 — bundled sample data
from aqeval import load_aq_data
data = load_aq_data()

# Option 2 — download AURN (or SAQN / AQE / WAQN / NI) data
from aqeval.data import download_aurn_data
data = download_aurn_data("LED6", 2018, 2025, source="aurn")
data.columns = data.columns.str.lower()   # normalise column names

Then set the analysis parameters:

POLLUTANT     = "no2"
H             = 0.3     # minimum segment size (fraction of series length)
DESEASON      = True     # remove seasonal patterns
DEWEATHER     = True     # remove weather effects (needs ws / wd columns)
ROLLING_WINDOW = 30      # rolling-mean window in days

1. plot_rolling_statistics — Rolling mean

Daily scatter plot overlaid with a rolling mean line — a quick way to visualise the overall trend and variability before running any formal analysis.

fig = plot_rolling_statistics(data, POLLUTANT, rolling_window=ROLLING_WINDOW)
fig.show()

2. plot_spectral — Spectral frequency analysis

Spectral power-density plot identifying periodic components in the data (e.g. 24-hour diurnal cycle, 7-day weekly pattern, 365-day annual cycle).

fig = plot_spectral(data, POLLUTANT)
fig.show()

3. plot_break_points — Break-point detection & quantification

Detects structural breaks in the time-series using the strucchange method (Bai–Perron dynamic programming with BIC break-count selection) and overlays the fitted break-point model on the observed daily means. Vertical dashed lines indicate detected break-point locations.

fig = plot_break_points(data, POLLUTANT, h=H)
fig.show()

4. plot_segments_raw — Segmented trend analysis (raw data)

Fits segmented linear regression to the raw (observed) daily data, showing the slope within each segment. Dotted lines indicate the confidence interval around each detected break-point.

fig = plot_segments_raw(data, POLLUTANT, h=H)
fig.show()

5. plot_segments_isolated — Segmented trend analysis (deseasoned data)

Same segmented regression, but applied to the deseasoned data. Seasonal patterns are removed first via isolate_contribution, so the trends here reflect underlying emission changes rather than weather/seasonal artefacts.

fig = plot_segments_isolated(data, POLLUTANT, h=H, deseason=DESEASON, deweather=DEWEATHER)
fig.show()

6. plot_raw_vs_isolated_trends — Side-by-side trend comparison

Two-row subplot comparing segmented trends side by side:

• Top: raw data with trend segments and break-point lines
• Bottom: isolated (deseasoned) data with its own trends and break-points

Differences between the two panels reveal where seasonal patterns may mask or create apparent changes.

fig = plot_raw_vs_isolated_trends(data, POLLUTANT, h=H, deweather=DEWEATHER)
fig.show()

7. plot_raw_vs_isolated_timeseries — Overlaid time-series

Overlays the raw and isolated daily mean time-series on a single axis. The gap between the two lines represents the seasonal/weather component that was removed.

fig = plot_raw_vs_isolated_timeseries(data, POLLUTANT, deweather=DEWEATHER)
fig.show()

8. plot_residuals_timeseries — Residuals over time

Colour-coded bar chart of daily residuals (raw − isolated):

• Orange bars: seasonal/weather factors increased concentrations that day
• Blue bars: seasonal/weather factors decreased concentrations

The black line is the rolling mean — if it drifts over time, seasonal patterns may be changing.

fig = plot_residuals_timeseries(data, POLLUTANT, deweather=DEWEATHER, rolling_window=ROLLING_WINDOW)
fig.show()

9. plot_residual_histogram — Residual distribution

Distribution of the residuals (raw − isolated). A symmetric histogram centred near zero indicates the isolation removed a balanced seasonal component. The red dotted line marks the mean.

fig = plot_residual_histogram(data, POLLUTANT, deweather=DEWEATHER, nbins=50)
fig.show()

10. plot_signal_isolation — Signal isolation (two-panel)

Two-panel view of the signal isolation result:

• Top: original daily data (grey scatter) with a rolling mean (blue line)
• Bottom: isolated data (green scatter) with its rolling mean (green line)

Shows at a glance how much seasonal structure has been removed.

fig = plot_signal_isolation(data, POLLUTANT, deseason=DESEASON, deweather=DEWEATHER, rolling_window=ROLLING_WINDOW)
fig.show()

11. plot_signal_isolation_trend_comparison — Isolation trend comparison

The full signal-isolation workflow in one figure: break-point detection and segmented regression on both the original and isolated datasets, displayed in a two-row subplot.

How to interpret the comparison:

• Same break-points in both panels → likely a real emission change
• Break-point only in raw data → may be a weather/seasonal artefact
• Break-point only in isolated data → a real change that was masked by weather
• Different trend magnitudes → weather amplified or dampened the true change

fig = plot_signal_isolation_trend_comparison(data, POLLUTANT, h=H, deseason=DESEASON, deweather=DEWEATHER)
fig.show()

Common parameters

All figure functions accept these optional keyword arguments:

Dependencies

• numpy >= 1.24.0
• pandas >= 2.2.0
• scipy >= 1.10.0
• statsmodels >= 0.14.0
• pygam >= 0.9.1
• plotly >= 5.0.0
• requests >= 2.25.0
• rdata >= 0.8.0
• python-dateutil >= 2.8.0

References

Ropkins, K. and Tate, J. (2021). Early assessment of the impact of the COVID-19 lockdown on air quality in the United Kingdom. Science of The Total Environment, 142374. DOI: 10.1016/j.scitotenv.2020.142374

Ropkins, K., Walker, A., et al. Change Detection of Air Quality Time-Series Using the R Package AQEval. Available at SSRN 4267722.

Ropkins, K., Walker, A., Philips, I., Rushton, C. E., Clark, T. and Tate, J. E. (2026). AQEval: R code for the analysis of discrete change in air quality time-series. Journal of Open Source Software, 11(118), 8839. DOI: 10.21105/joss.08839

Zeileis, A., Kleiber, C., Krämer, W. and Hornik, K. (2003). Testing and dating of structural changes in practice. Computational Statistics & Data Analysis, 44, 109–123.

Muggeo, V. M. (2003). Estimating regression models with unknown break-points. Statistics in Medicine, 22(19), 3055–3071.

License

GPL (>= 3)

Authors

• Original R Package: Karl Ropkins, Anthony Walker, James Tate
• Python Port: Chris Rushton