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. The port follows R's arithmetic, not just its algebra: the same normal-equations solve, the same pivoted QR with lm's 1e-7 rank tolerance, and so the same candidate fits accepted and rejected.
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). This is the one deliberate deviation from R.

How closely it matches R

The port is checked against R three ways, all in the test suite, so none of them needs R to run:

1. R AQEval's own regression fixtures (inst/tinytest/): break locations and confidence intervals, quantified break rows, segment boundaries, date-range statistics and spectral output dimensions (tests/test_regression_r_fixtures.py).
2. A live R run on real AURN data (tests/test_regression_r_break_segments.py). R AQEval 0.6.11, with openair 3.1.0 and strucchange 1.6.0 under R 4.6.1, analysed 20 series: measured and isolated NO2 and NOx from four AURN sites, 2018–2025, averaged to 8 hours, 1 day and 7 days, at h from 0.12 to 0.3. The series R analysed and everything findBreakPoints, testBreakPoints, quantBreakPoints and quantBreakSegments returned are committed as fixtures, including a log of all 4,716 candidate fits in quantBreakSegments' grid searches.
3. openair::timeAverage on the bundled aq.data (tests/test_regression_r_time_average.py).

What is identical to R, and what is only very close:

Things worth knowing behind that table:

• Least-squares results cannot be bit-identical. segmented solves the normal equations, and their rounding depends on the linear-algebra library (BLAS) in use; R gives slightly different digits on different machines too. The decisions that matter — which candidate fits are kept and which wins — came out the same as R's in every case tested, on macOS (Accelerate) and on Linux (OpenBLAS, arm64 and x86-64), with pandas 2.3 and 3.0.
• pandas 3. Up to 0.7.2 the break analysis assumed the date column was stored at nanosecond resolution. pandas 3 no longer defaults to that, and download_aurn_data returns second-resolution dates under it. Two things went wrong. The break-point model behind test_break_points and quant_break_points lost its time trend (every row got the same regressor value), so a different subset of breaks, or none at all, was preferred. And quant_break_segments' time variable came out a few units in the last place away from R's; its grid search starts exactly on data points, so that was enough to change the winning fit, and on the 20 series above it agreed with R on fewer than half. From 0.7.3 the results do not depend on the stored resolution (tested at s, ms, µs and ns) or on the pandas version. Anything run with 0.7.2 and pandas 3 is worth re-running. spectral_frequency had the same assumption, with an effect of about 5×10⁻¹⁰.
• openair 3.1. time_average follows openair 3.1: multiples of a fixed-length unit ("8 hour", "2 day", "7 day") are aligned to the Unix epoch, so 7-day periods start on Thursdays wherever the data begin, and "week" starts on Sundays. Earlier openair releases counted multi-day periods from the first day of the data, as this package did before 0.7.3. R AQEval's tinytest values were recorded that way (break rows 213/317/402 on 2-day aq.data); with openair 3.1.0 R itself returns 214/317/405 for the same call. The tests check both.
• Text is not exact unless parsed exactly. pandas' default CSV float parser is off by one unit in the last place for roughly one value in ten. When comparing with R at this level, read R-exported numbers with pd.read_csv(..., float_precision="round_trip"); the bundled aq.data is loaded that way from 0.7.3 and now equals R's exactly.
• No breaks. As in R, quant_break_segments(data, pollutant, breaks=None) fits a single straight line, so the result of find_break_points can be passed straight in. Leave breaks out to have the breaks found and tested for you.

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: the break-point search is a vectorised implementation of the O(n²) Bai–Perron dynamic program. On the 8-hour BDMA series (~6,000 rows) it runs in well under a second at any h down to 0.05 — including the h=0.12 setting recommended for 8-hour data (R's strucchange takes minutes on the same problem). Memory still grows as O(n²) (~300 MB at n≈6,000), so run the analysis on time-averaged data (8-hour/daily), not raw multi-year hourly series.

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