Statistical Analysis of Air Quality Data in R and Python: Tools, Packages and Working Code

Spreadsheets are where air quality analysis goes to become unreproducible. A year of hourly data for eight pollutants is over 70,000 cells; a rolling 8-hour maximum computed with a dragged formula is a defect waiting to be discovered; and when next year's data arrives, the whole thing has to be rebuilt by hand. Scripted analysis in R or Python solves all three problems at once: the method is written down, the result can be regenerated in seconds, and every step is reviewable.

This guide covers the practical toolkit — which packages do what, working code for the operations that come up on every project, and the reproducibility practices that make an analysis defensible rather than merely correct.

Statistical analysis of air quality data in R and Python uses purpose-built packages — openair in R and pandas with scipy and scikit-learn in Python — to validate monitoring data, aggregate it to regulatory averaging periods, compute compliance and AQI statistics, run non-parametric trend tests, produce diagnostic plots such as polar plots and calendar heatmaps, and generate reproducible reports directly from the raw data.

R or Python? The Honest Answer

Both are excellent. The choice depends on what surrounds the analysis.

	R	Python
Domain packages	openair — purpose-built for air quality, unmatched	No direct equivalent; assemble from general libraries
Statistics	Deeper and more mature, especially non-parametric methods	Strong and sufficient via scipy and statsmodels
Data handling	dplyr and data.table, both excellent	pandas, excellent
Machine learning	Good (tidymodels, ranger)	Best in class — scikit-learn, XGBoost, PyTorch
Geospatial	sf, terra — very strong	geopandas, rasterio — very strong
Reporting	R Markdown and Quarto, outstanding	Jupyter and Quarto, very good
Production pipelines	Workable	Better — deployment, APIs, scheduling
Learning curve for analysts	Gentler for statistics-first users	Gentler for programming-first users

Practical recommendation. For descriptive analysis, diagnostic plots and trend testing, R with openair will get you there faster than anything else — many operations that take fifty lines of Python are one function call. For machine-learning forecasts, production data pipelines, and integration with web applications, Python is the better home. Many teams use both, and there is no reason not to.

The R Toolkit

openair

David Carslaw's openair package is the reason many air quality analysts use R at all. It implements, as single functions, analyses that would otherwise be substantial projects:

Function	What it does
polarPlot()	Concentration as a bivariate function of wind speed and direction — the best single source-diagnostic plot available
windRose()	Wind rose, and pollutionRose() for concentration-weighted roses
timeVariation()	Diurnal, weekday and monthly profiles with confidence intervals, in one call
calendarPlot()	Daily values as a calendar heatmap, ideal for episode identification
TheilSen()	Theil-Sen trend with Mann-Kendall significance, deseasonalised
smoothTrend()	Non-parametric smoothed trend with uncertainty bands
timeAverage()	Correct aggregation with a data-capture threshold enforced
selectByDate()	Convenient temporal subsetting
percentileRose()	Concentration percentiles by wind direction
trajCluster()	HYSPLIT back-trajectory clustering
corPlot()	Correlation matrix across pollutants

Supporting packages: dplyr and lubridate for data manipulation and dates, trend and Kendall for additional trend tests, zoo for rolling windows, ggplot2 for custom figures, sf for spatial work, and openairmaps for interactive mapping.

A worked R example

The following covers the operations that appear on nearly every project.

library(openair)
library(dplyr)
library(lubridate)

# ---- 1. Load and structure -------------------------------------------
aq <- read.csv("caaqms_station.csv", stringsAsFactors = FALSE) %>%
  mutate(date = ymd_hms(timestamp, tz = "Asia/Kolkata")) %>%
  rename(ws = wind_speed, wd = wind_dir, pm25 = PM2.5, pm10 = PM10)

# ---- 2. Validation ---------------------------------------------------
aq <- aq %>%
  mutate(
    pm25 = ifelse(pm25 < 0 | pm25 > 1000, NA, pm25),
    pm10 = ifelse(pm10 < 0 | pm10 > 2000, NA, pm10),
    # PM2.5 can never exceed PM10 - flag both as suspect
    ratio_bad = !is.na(pm25) & !is.na(pm10) & pm25 > pm10,
    pm25 = ifelse(ratio_bad, NA, pm25),
    pm10 = ifelse(ratio_bad, NA, pm10),
    # stalled anemometer: zero speed with a frozen direction
    wd = ifelse(ws < 0.5, NA, wd)
  )

# ---- 3. Data capture -------------------------------------------------
capture <- aq %>%
  mutate(month = floor_date(date, "month")) %>%
  group_by(month) %>%
  summarise(
    n_possible = n(),
    pm25_valid = sum(!is.na(pm25)),
    pm25_capture = round(100 * pm25_valid / n_possible, 1)
  )
print(capture)

# ---- 4. Aggregation with an enforced capture threshold ---------------
daily <- timeAverage(aq, avg.time = "day", data.thresh = 75)

# ---- 5. NAAQS exceedance statistics ----------------------------------
exceed <- daily %>%
  summarise(
    valid_days      = sum(!is.na(pm25)),
    annual_mean     = round(mean(pm25, na.rm = TRUE), 1),
    exceed_days     = sum(pm25 > 60, na.rm = TRUE),
    exceed_pct      = round(100 * exceed_days / valid_days, 1),
    p98             = round(quantile(pm25, 0.98, na.rm = TRUE), 1)
  )
print(exceed)

# ---- 6. Diagnostic plots ---------------------------------------------
timeVariation(aq, pollutant = "pm25", normalise = FALSE)
polarPlot(aq, pollutant = "pm25", statistic = "mean")
calendarPlot(daily, pollutant = "pm25", year = 2025)
pollutionRose(aq, pollutant = "pm25")

# ---- 7. Trend test ---------------------------------------------------
TheilSen(aq, pollutant = "pm25", deseason = TRUE, autocor = TRUE)

Seven blocks, and you have validated data, a capture statement, correct aggregation, compliance statistics, four diagnostic plots and a trend test with autocorrelation handled. That is the case for openair in a nutshell.

The Python Toolkit

Package	Role
pandas	Data structures, time indexing, resampling, rolling windows
numpy	Numerical operations
scipy.stats	Statistical tests, distributions
statsmodels	Regression, time series, seasonal decomposition
pymannkendall	Mann-Kendall variants including seasonal and pre-whitened
scikit-learn	Regression, random forests, gradient boosting, model validation
xgboost / lightgbm	High-performance gradient boosting for forecasting
matplotlib / seaborn / plotly	Visualisation
geopandas / rasterio	Spatial analysis
windrose	Wind rose plotting

A worked Python example

import pandas as pd
import numpy as np
import pymannkendall as mk

# ---- 1. Load and structure -------------------------------------------
df = pd.read_csv("caaqms_station.csv", parse_dates=["timestamp"])
df["timestamp"] = df["timestamp"].dt.tz_localize("Asia/Kolkata")
df = df.set_index("timestamp").sort_index()

# ---- 2. Validation ----------------------------------------------------
df.loc[(df.PM25 < 0) | (df.PM25 > 1000), "PM25"] = np.nan
df.loc[(df.PM10 < 0) | (df.PM10 > 2000), "PM10"] = np.nan
bad = df.PM25 > df.PM10                       # physically impossible
df.loc[bad, ["PM25", "PM10"]] = np.nan
df.loc[df.wind_speed < 0.5, "wind_dir"] = np.nan   # stalled anemometer

# ---- 3. Data capture --------------------------------------------------
capture = (df.PM25.resample("MS").count() /
           df.PM25.resample("MS").size() * 100).round(1)
print(capture)

# ---- 4. Aggregation with a 75% capture rule ---------------------------
def daily_mean(series, thresh=0.75):
    counts = series.resample("D").count()
    means  = series.resample("D").mean()
    return means.where(counts >= 24 * thresh)

pm25_daily = daily_mean(df.PM25)

# ---- 5. Rolling 8-hour maximum, the correct way for O3 and CO ---------
o3_8h = df.O3.rolling("8h", min_periods=6).mean()
o3_8h_max_daily = o3_8h.resample("D").max()

# ---- 6. Compliance statistics -----------------------------------------
valid = pm25_daily.dropna()
print({
    "valid_days":  len(valid),
    "annual_mean": round(valid.mean(), 1),
    "exceed_days": int((valid > 60).sum()),
    "p98":         round(valid.quantile(0.98), 1),
})

# ---- 7. Trend test with seasonality and autocorrelation handled -------
monthly = pm25_daily.resample("MS").mean().dropna()
result = mk.seasonal_test(monthly, period=12)
print(f"trend={result.trend}, p={result.p:.4f}, slope={result.slope:.2f}")

AQI computation in Python

The AQI is worth scripting once and reusing, because the breakpoint tables and averaging rules are exactly the kind of thing that gets copied wrong in spreadsheets. The methodology is set out in how to calculate the AQI in India.

BREAKPOINTS = {                       # CPCB National AQI
    "PM25": [(0,30,0,50),(30,60,51,100),(60,90,101,200),
             (90,120,201,300),(120,250,301,400),(250,1000,401,500)],
    "PM10": [(0,50,0,50),(50,100,51,100),(100,250,101,200),
             (250,350,201,300),(350,430,301,400),(430,2000,401,500)],
    "NO2":  [(0,40,0,50),(40,80,51,100),(80,180,101,200),
             (180,280,201,300),(280,400,301,400),(400,1000,401,500)],
    "SO2":  [(0,40,0,50),(40,80,51,100),(80,380,101,200),
             (380,800,201,300),(800,1600,301,400),(1600,3000,401,500)],
    "CO":   [(0,1,0,50),(1,2,51,100),(2,10,101,200),
             (10,17,201,300),(17,34,301,400),(34,100,401,500)],
    "O3":   [(0,50,0,50),(50,100,51,100),(100,168,101,200),
             (168,208,201,300),(208,748,301,400),(748,1500,401,500)],
    "NH3":  [(0,200,0,50),(200,400,51,100),(400,800,101,200),
             (800,1200,201,300),(1200,1800,301,400),(1800,3000,401,500)],
}

def sub_index(pollutant, conc):
    if conc is None or np.isnan(conc):
        return None
    for bp_lo, bp_hi, i_lo, i_hi in BREAKPOINTS[pollutant]:
        if bp_lo <= conc <= bp_hi:
            return round(i_lo + (i_hi - i_lo) / (bp_hi - bp_lo) * (conc - bp_lo))
    return 500

def aqi(concentrations):
    """concentrations: dict of pollutant -> value at the correct averaging period."""
    subs = {p: sub_index(p, c) for p, c in concentrations.items()}
    subs = {p: v for p, v in subs.items() if v is not None}
    # CPCB rule: >= 3 pollutants, at least one of PM2.5 / PM10
    if len(subs) < 3 or not ({"PM25", "PM10"} & set(subs)):
        return None, None
    prominent = max(subs, key=subs.get)
    return subs[prominent], prominent

Note the minimum-data rule enforced in the function itself. That is the point of scripting it: the rule cannot be forgotten because it is not a step someone has to remember.

Reproducibility Practices That Matter

Correct analysis is not enough; it has to be re-runnable and reviewable.

1. Never edit the raw data. Read it, transform it in code, write outputs elsewhere. The raw file stays pristine.
2. Script everything, including the plots. A chart produced by clicking cannot be regenerated when the data is corrected.
3. Use version control. Git on the analysis code, always. It is the difference between "the numbers changed" and "here is the commit where the numbers changed and why."
4. Pin your environment. renv in R, a lock file or requirements.txt in Python. Package updates silently change results.
5. Use literate reporting. R Markdown, Quarto or Jupyter, so the narrative and the code producing the numbers live in the same document. Every figure in the report is generated from the data at build time, which eliminates the pasted-stale-chart failure mode entirely.
6. Write tests for critical functions. The AQI function above should have a test asserting that 111 µg/m³ PM2.5 returns 268. Once written, it can never silently break.
7. Log the data-capture rate in the output, not in a side note, so the reader always sees how much data supports each number.
8. Separate configuration from code. Thresholds, station lists and breakpoints belong in a config file, not scattered through scripts.

These practices are covered from a general statistics perspective in our guides to R programming for statistics and statistical data analysis services.

Common Coding Errors in Air Quality Analysis

• Timezone mishandling. Aggregating to "days" in UTC when the standard requires IST calendar days shifts every daily average by 5.5 hours.
• Fixed 8-hour blocks instead of a rolling maximum for ozone and CO, which systematically under-reports.
• Aggregating without a capture threshold, so a day with three valid hours becomes a daily mean.
• Silent NA propagation — using mean(x, na.rm = TRUE) everywhere without checking how many values survived.
• Averaging AQI values rather than concentrations, which is invalid because the AQI transform is non-linear.
• Unit confusion, particularly CO in mg/m³ against everything else in µg/m³.
• Resampling before validating, which bakes bad values into every downstream statistic.
• Random train-test splits on time series, which leaks future information into the training set and produces optimistic model scores.

How DSTechnoverse Helps

DSTechnoverse builds and runs scripted air quality analysis pipelines for consultancies, industries, municipal programmes and research groups.

• Analysis pipeline development — reproducible R or Python pipelines covering ingestion, validation, flagging, aggregation, statistics and figure generation, delivered as code you own and can run yourself.
• Custom tooling — AQI engines with the CPCB rules correctly enforced, automated NAAQS compliance reporting, exceedance trackers, and dashboard back-ends built on validated data.
• Analysis on demand — where you have data and need results rather than a pipeline, we run the full analysis and deliver the report plus the scripts behind it.
• Code review and audit — we review existing analysis code and spreadsheets, and routinely find timezone errors, missing capture thresholds and incorrect rolling-window logic in systems already producing published numbers.
• Training — practical sessions for environmental teams on openair, pandas and reproducible reporting, using your own data.

The principle throughout is that you receive both the answer and the machinery, so next year's dataset takes an afternoon rather than a month. The same approach underpins our statistical and research data analysis service and our air quality analysis service.

---

Need air quality data turned into a defensible answer? DSTechnoverse provides end-to-end air quality analysis services — monitoring data validation, statistical analysis, source apportionment, dispersion modelling, GIS mapping, and regulator-ready reporting. Talk to our air quality team about your dataset or project.

Frequently Asked Questions

Should I use R or Python for air quality data analysis?

Use R with openair for descriptive analysis, diagnostic plots and trend testing, where purpose-built functions do in one line what would take a page of general-purpose code. Use Python for machine-learning forecasting, production data pipelines and web integration. Many teams use both, and the file formats are interchangeable.

What is the openair package?

openair is an R package built specifically for air quality data analysis. It provides functions for polar plots, wind and pollution roses, diurnal and seasonal profiles, calendar heatmaps, Theil-Sen trend analysis with Mann-Kendall significance, correct time averaging with data-capture thresholds, and back-trajectory clustering.

How do I calculate a rolling 8-hour maximum correctly?

Compute the rolling mean over every 8-hour window in the day, requiring a minimum number of valid hours in each window, then take the maximum of those rolling means. Using three fixed 8-hour blocks instead is a common error that systematically under-reports the true maximum.

Why is timezone handling important in air quality analysis?

Because regulatory averaging periods are defined on local calendar days. Aggregating in UTC when the standard requires IST shifts every daily average by five and a half hours, changing daily means, exceedance counts and the 98th-percentile value. Always localise timestamps explicitly at import.

Can air quality analysis be done in Excel?

Small compliance checks can. Anything involving a full year of hourly data across multiple parameters should not be, because rolling windows, capture thresholds and repeated aggregation are error-prone to build by hand, impossible to review, and have to be rebuilt each time new data arrives. Scripted analysis is reproducible and auditable.

What does reproducible analysis mean in practice?

It means that anyone with the raw data and the code can regenerate every number and figure in the report exactly. In practice that requires never editing raw data, scripting all transformations and plots, using version control, pinning package versions, and writing the report in a literate format where figures are generated at build time.

How do I enforce the 75% data capture rule in code?

Count valid records per averaging period alongside computing the mean, then set the mean to missing wherever the count falls below the threshold. In R, openair's timeAverage has a data.thresh argument that does this directly; in Python it is a few lines combining resample().count() and resample().mean().

What is the most common mistake in air quality analysis code?

Aggregating before validating. Once invalid records — calibration periods, negative spikes, stuck values — are averaged into daily means, they are effectively invisible and contaminate every downstream statistic including trends, exceedance counts and AQI values. Validation must come first, always.