Metadata-Version: 2.4
Name: aqeval
Version: 0.7.1
Summary: Air Quality Evaluation – Python port of the R package AQEval
Author: Chris Rushton
Author-email: Karl Ropkins <k.ropkins@its.leeds.ac.uk>
Maintainer: Chris Rushton
License-Expression: GPL-3.0-or-later
Project-URL: Homepage, https://github.com/chris-r-uol/aqeval_python
Project-URL: Repository, https://github.com/chris-r-uol/aqeval_python
Project-URL: R Original, https://github.com/karlropkins/AQEval
Project-URL: R Documentation, https://karlropkins.github.io/AQEval/
Keywords: air-quality,time-series,break-point,signal-isolation,atmospheric-science
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Atmospheric Science
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24.0
Requires-Dist: pandas>=2.2.0
Requires-Dist: scipy>=1.10.0
Requires-Dist: statsmodels>=0.14.0
Requires-Dist: pygam>=0.10.1
Requires-Dist: plotly>=5.0.0
Requires-Dist: requests>=2.25.0
Requires-Dist: rdata>=0.8.0
Requires-Dist: python-dateutil>=2.8.0
Provides-Extra: dev
Requires-Dist: pytest>=7.0.0; extra == "dev"
Requires-Dist: pytest-cov>=4.0.0; extra == "dev"
Requires-Dist: black>=23.0.0; extra == "dev"
Requires-Dist: flake8>=6.0.0; extra == "dev"
Requires-Dist: mypy>=1.0.0; extra == "dev"
Requires-Dist: build; extra == "dev"
Requires-Dist: twine; extra == "dev"
Dynamic: license-file

# AQEval Python Package

[![CI](https://github.com/chris-r-uol/aqeval_python/actions/workflows/ci.yml/badge.svg)](https://github.com/chris-r-uol/aqeval_python/actions/workflows/ci.yml)

A Python port of the R package [AQEval](https://github.com/karlropkins/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](https://doi.org/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](https://pygam.readthedocs.io) (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

```bash
# Install from PyPI
pip install aqeval



# Install in development mode
pip install -e .

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

## Quick Start

```python
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

```python
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`](examples/tutorial_bdma.py) replicates the
[AQEval R tutorial](https://github.com/chris-r-uol/AQEval_Tutorial)
end-to-end with this package:

```python
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:

| Quantity | R (mgcv) | Python (pygam) |
|---|---|---|
| Break 1 (lower, bpt, upper; 1-based 8h rows) | 2116, 2189, 2240 | 2116, 2189, 2239 |
| Break 2 | 4557, 4701, 4824 | 4557, 4701, 4825 |
| Change event 1 | 41.60 → 34.17 (−17.9%), Nov 2019–Mar 2020 | 41.54 → 34.28 (−17.5%), Dec 2019–Mar 2020 |
| Change event 2 | 37.91 → 33.16 (−12.5%), Apr–May 2022 | 37.82 → 33.22 (−12.2%), Apr–May 2022 |

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`.

```python
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:

```python
# 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:

```python
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.

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

| Parameter | Default | Description |
|-----------|---------|-------------|
| `rolling_window` | `30` | Rolling-mean window in days |
| `title` | auto | Custom plot title |
| `height` | `400` | Figure height in pixels |

---

### 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).

```python
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.

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

| Parameter | Default | Description |
|-----------|---------|-------------|
| `h` | `0.15` | Minimum segment size as a fraction of data length |
| `title` | auto | Custom plot title |
| `height` | `500` | Figure height in pixels |

---

### 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.

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

| Parameter | Default | Description |
|-----------|---------|-------------|
| `h` | `0.15` | Minimum segment size fraction |
| `show_break_ci` | `True` | Show break-point CI lines |

---

### 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.

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

| Parameter | Default | Description |
|-----------|---------|-------------|
| `deseason` | `True` | Remove seasonal patterns |
| `deweather` | `False` | Remove weather effects (needs `ws`/`wd` columns) |

---

### 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.

```python
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.

```python
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.

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

| Parameter | Default | Description |
|-----------|---------|-------------|
| `rolling_window` | `30` | Rolling-mean window in days |

---

### 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.

```python
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.

```python
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

```python
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:

| Parameter | Default | Description |
|-----------|---------|-------------|
| `title` | auto-generated | Custom plot title string |
| `height` | varies (400–700) | Figure height in pixels |

## 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](https://doi.org/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](https://ssrn.com/abstract=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](https://doi.org/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
