How to Compute an Air Quality Index in Python

On this page

Problem statement

An AQI looks like a lookup table, and most implementations are exactly that: take a concentration, find its band, interpolate. The table is the easy part. The work is everything around it — each pollutant on its own averaging period, gases converted from µg/m³ to the ppb or ppm the US tables use, concentrations truncated to the table's precision before the lookup, days with too few hours thrown out rather than averaged, and the hourly "current AQI" that American sites publish computed with a NowCast rather than a 24-hour mean.

Get any of those wrong and the index is off by a band on a fraction of days — which is the only part of an index anybody looks at.

Quick answer

A vectorised US AQI for four pollutants, with the averaging and unit handling each one needs:

import io, requests, numpy as np, pandas as pd

BP = {   # (C_lo, C_hi, I_lo, I_hi); PM in µg/m³, NO2 in ppb, O3 in ppm
    "PM2.5": [(0.0, 9.0, 0, 50), (9.1, 35.4, 51, 100), (35.5, 55.4, 101, 150),
              (55.5, 125.4, 151, 200), (125.5, 225.4, 201, 300), (225.5, 325.4, 301, 500)],
    "PM10":  [(0, 54, 0, 50), (55, 154, 51, 100), (155, 254, 101, 150),
              (255, 354, 151, 200), (355, 424, 201, 300), (425, 604, 301, 500)],
    "NO2":   [(0, 53, 0, 50), (54, 100, 51, 100), (101, 360, 101, 150),
              (361, 649, 151, 200), (650, 1249, 201, 300), (1250, 2049, 301, 500)],
    "O3":    [(0.000, 0.054, 0, 50), (0.055, 0.070, 51, 100), (0.071, 0.085, 101, 150),
              (0.086, 0.105, 151, 200), (0.106, 0.200, 201, 300)],
}
DECIMALS = {"PM2.5": 1, "PM10": 0, "NO2": 0, "O3": 3}

def subindex(c, pollutant):
    """US AQI sub-index for an array of concentrations already in the table's unit."""
    c = np.floor(np.asarray(c, float) * 10 ** DECIMALS[pollutant]) / 10 ** DECIMALS[pollutant]
    out = np.full(c.shape, np.nan)
    for lo, hi, ilo, ihi in BP[pollutant]:
        m = (c >= lo) & (c <= hi)
        out[m] = np.round((ihi - ilo) / (hi - lo) * (c[m] - lo) + ilo)
    out[c > BP[pollutant][-1][1]] = 500
    return out

def aurn(site, year=2023):
    text = requests.get(f"https://uk-air.defra.gov.uk/datastore/data_files/site_data/"
                        f"{site}_{year}.csv", timeout=120).text
    d = pd.read_csv(io.StringIO(text), skiprows=4, low_memory=False)
    end = pd.to_datetime(d.Date, dayfirst=True) + pd.to_timedelta(d.time.str[:2].astype(int), "h")
    idx = pd.DatetimeIndex(end - pd.Timedelta("1h"), name="start")
    col = lambda name: pd.Series(pd.to_numeric(d[name], errors="coerce").to_numpy(), index=idx)
    return pd.DataFrame({
        "PM2.5": col("PM<sub>2.5</sub> particulate matter (Hourly measured)"),
        "PM10": col("PM<sub>10</sub> particulate matter (Hourly measured)"),
        "NO2": col("Nitrogen dioxide") / 1.88,          # µg/m³ → ppb at 25 °C
        "O3": col("Ozone") / 1.96 / 1000,               # µg/m³ → ppm at 25 °C
    })

def daily_aqi(h):
    enough = lambda s: s.resample("D").count() >= 18                 # 75% of the day
    pm25 = h["PM2.5"].resample("D").mean().where(enough(h["PM2.5"]))
    pm10 = h["PM10"].resample("D").mean().where(enough(h["PM10"]))
    no2 = h["NO2"].resample("D").max()                               # daily max 1-h
    o3 = h["O3"].rolling(8, min_periods=6).mean().resample("D").max() # daily max 8-h
    sub = pd.DataFrame({"PM2.5": subindex(pm25, "PM2.5"), "PM10": subindex(pm10, "PM10"),
                        "NO2": subindex(no2, "NO2"), "O3": subindex(o3, "O3")}, index=pm25.index)
    sub = sub.dropna(how="all")                       # a day with no valid pollutant has no AQI
    return sub.assign(AQI=sub.max(axis=1), dominant=sub.idxmax(axis=1))

h = aurn("MY1")
daily = daily_aqi(h)
cats = pd.cut(daily.AQI, [-1, 50, 100, 150, 200, 300, 500],
              labels=["Good", "Moderate", "USG", "Unhealthy", "Very unhealthy", "Hazardous"])
print(f"Marylebone Road 2023: {daily.AQI.notna().sum()} days with an AQI")
print(cats.value_counts().reindex(cats.cat.categories).to_string())
print("\nhighest days:")
print(daily.nlargest(3, "AQI")[["PM2.5", "PM10", "NO2", "O3", "AQI", "dominant"]].to_string())

Marylebone Road 2023: 364 days with an AQI
Good              129
Moderate          233
USG                 2
Unhealthy           0
Very unhealthy      0
Hazardous           0

highest days:
            PM2.5  PM10   NO2    O3    AQI dominant
start                                              
2023-02-14  111.0  49.0  46.0   4.0  111.0    PM2.5
2023-09-10  101.0  41.0  48.0  29.0  101.0    PM2.5
2023-02-15   99.0  44.0  42.0  11.0   99.0    PM2.5

The daily AQI is the maximum of four sub-indices, each computed from a different aggregation of the hourly data. On the three worst days it is PM2.5 — a February episode and a September one — with NO2 and ozone far behind.

Table units, own averaging periods, sufficiency, truncation

1. Put each pollutant in the table's unit

US breakpoints are in µg/m³ for particles, ppb for NO2 and ppm for ozone, at 25 °C. European data is in µg/m³ throughout. Convert once, on the hourly data, before any averaging.

2. Aggregate each pollutant on its own period

PM2.5 and PM10: 24-hour mean. NO2: the hourly value (daily reports use the day's maximum hour). Ozone: 8-hour running mean (daily reports use the day's maximum). Build four daily series, not one.

3. Refuse days without enough data

A 24-hour mean needs at least 18 hours; an 8-hour mean at least 6. A day with four hours of PM2.5 has no PM2.5 index, not an index from four hours.

4. Truncate before the lookup

The tables are defined at a fixed precision: 9.04 µg/m³ is 9.0 and "Good"; 9.1 is "Moderate". Rounding instead of truncating moves values across boundaries.

5. Interpolate within the band and take the maximum

Linear interpolation inside the band gives the sub-index; the highest sub-index is the AQI; its pollutant is the dominant pollutant. Report both.

6. Use NowCast for hourly PM reporting

An hourly "current" AQI for particles is computed from the NowCast — a weighted mean of the last 12 hours that responds faster than a 24-hour mean.

Four pollutants, four aggregations, one maximum.

NowCast and boundary precision

The data used on this page

Hourly measurements from the UK's Automatic Urban and Rural Network, downloaded by the page's own code from DEFRA's site-data archive: London Bloomsbury (CLL2), London North Kensington (KC1), London Marylebone Road (MY1). The files are ratified data; the hour in each row is the end of the measurement hour, in GMT.

import requests

print("CLL2_2023.csv:", requests.head("https://uk-air.defra.gov.uk/datastore/data_files/site_data/CLL2_2023.csv", timeout=60, allow_redirects=True).status_code)
print("KC1_2023.csv:", requests.head("https://uk-air.defra.gov.uk/datastore/data_files/site_data/KC1_2023.csv", timeout=60, allow_redirects=True).status_code)
print("MY1_2023.csv:", requests.head("https://uk-air.defra.gov.uk/datastore/data_files/site_data/MY1_2023.csv", timeout=60, allow_redirects=True).status_code)

CLL2_2023.csv: 200
KC1_2023.csv: 200
MY1_2023.csv: 200

Example 1 — NowCast for hourly PM2.5

def nowcast(series):
    """EPA NowCast over the last 12 hours; needs 2 of the 3 most recent hours."""
    values = series.to_numpy()
    out = np.full(len(values), np.nan)
    for i in range(11, len(values)):
        c = values[i - 11:i + 1][::-1]                         # most recent first
        if np.isnan(c[:3]).sum() > 1:
            continue
        valid = c[~np.isnan(c)]
        w = max(valid.min() / valid.max(), 0.5) if valid.max() > 0 else 1.0
        weights = w ** np.arange(12)
        ok = ~np.isnan(c)
        out[i] = (weights[ok] * c[ok]).sum() / weights[ok].sum()
    return pd.Series(out, index=series.index)

pm = h["PM2.5"]
nc = nowcast(pm)
jump = pm.diff().idxmax()
window = pd.DataFrame({"hourly": pm, "NowCast": nc,
                       "24-h mean": pm.rolling(24, min_periods=18).mean()}).loc[
    jump - pd.Timedelta("3h"): jump + pd.Timedelta("4h")].round(1)
print(f"sharpest hourly rise of the year, at {jump}:")
print(window.to_string())

sharpest hourly rise of the year, at 2023-01-25 14:00:00:
                     hourly  NowCast  24-h mean
start                                          
2023-01-25 11:00:00    31.0     33.1       28.0
2023-01-25 12:00:00    33.0     33.1       28.2
2023-01-25 13:00:00    35.0     33.7       29.0
2023-01-25 14:00:00    68.0     50.8       31.2
2023-01-25 15:00:00    34.0     42.4       32.0
2023-01-25 16:00:00    39.0     40.7       32.9
2023-01-25 17:00:00    36.0     38.4       33.7
2023-01-25 18:00:00    34.0     36.2       34.3

A one-hour spike from 35 to 68 µg/m³ lifts the NowCast by 17 at once, and it decays over the following hours; the 24-hour mean moves by 2. That responsiveness is the point of it — and the reason an hourly AQI and a daily AQI for the same day can disagree.

Example 2 — the precision rule matters at boundaries

for c in (8.95, 9.04, 9.06, 9.1, 35.44, 35.46):
    trunc = subindex([c], "PM2.5")[0]
    rounded = subindex([round(c, 1) + 1e-9], "PM2.5")[0]
    print(f"PM2.5 {c:>6}: truncated → {trunc:>4.0f}, rounded first → {rounded:>4.0f}")

PM2.5   8.95: truncated →   49, rounded first →   49
PM2.5   9.04: truncated →   50, rounded first →   50
PM2.5   9.06: truncated →   50, rounded first →   51
PM2.5    9.1: truncated →   51, rounded first →   51
PM2.5  35.44: truncated →  100, rounded first →  100
PM2.5  35.46: truncated →  100, rounded first →  101

Values between 9.05 and 9.1 are Good when truncated and Moderate when rounded. Over a year of daily means only a handful of days sit that close to a boundary, which is exactly why the discrepancy survives testing.

Example 3 — several sites, one table

rows = {}
for site in ("MY1", "KC1", "CLL2"):
    d = daily_aqi(aurn(site))
    rows[site] = {"days": int(d.AQI.notna().sum()),
                  "median AQI": d.AQI.median(),
                  "days > 100": int((d.AQI > 100).sum()),
                  "top pollutant": d.dominant.value_counts().idxmax()}
print(pd.DataFrame(rows).T.to_string())

     days median AQI days > 100 top pollutant
MY1   364       53.0          2         PM2.5
KC1   365       38.0         10         PM2.5
CLL2  365       38.0          4         PM2.5

The kerbside site has the highest median but the fewest days above 100: North Kensington, the background site, has ten. The table also hides which sub-indices each site could compute — a site without PM10 or ozone monitors computes its AQI from fewer pollutants and can look cleaner for that reason alone. Say which were available.

The NowCast jumps with a spike and decays over hours; the 24-hour mean barely notices.

Convert first, truncate, require enough data

Why the unit conversion happens first

Converting a daily mean in µg/m³ to ppb gives the same number as averaging converted hourly values, because the conversion is linear. The reason to do it first is bookkeeping: once every hourly column is in its table's unit, the aggregation code cannot mix them.

Why truncation, not rounding

The breakpoint tables have gaps — 9.0 then 9.1 — because they are defined at a fixed number of decimals. Truncating to that precision guarantees every value falls into exactly one band. Rounding would put 9.05 into a gap on some implementations and into the upper band on others.

Why an index needs a data-sufficiency rule

A 24-hour mean computed from a few morning hours can be far from the day's true mean, especially for traffic pollutants with strong diurnal cycles. Index values published on insufficient data are the most common source of spurious "unhealthy" days in automated reports.

The precision rule lives in the first line; forget it and boundary days change category.

One-hour ozone, off-table values and time zones

Ozone 1-hour: the US table has a separate 1-hour ozone scale for high values; the 8-hour one applies below 0.200 ppm.
Above the table: concentrations beyond the last breakpoint are reported as 500 (or "beyond the AQI").
Negative concentrations from instruments should be kept for means and clipped to zero only for the lookup.
Time zones: US daily AQI uses local days; UK and EU data in GMT or local standard time shifts which hours fall into which day.
Other indices (UK DAQI, EU EAQI) use their own tables, units and periods; the structure of the code is the same.

Internal links

Air quality indices explained — what an index is and why they differ
AQI averaged the wrong way — the fix guide for averaging index values
ppb and µg/m³ mixed up — the conversion in detail
How to download air quality station data in Python — the hourly inputs
How to analyse an air quality time series — exceedances without an index
Air quality data gaps and negative values — data sufficiency in depth

FAQ

How do I calculate the AQI from PM2.5?

Take the 24-hour mean in µg/m³, truncate to one decimal, find its band in the EPA table and interpolate linearly within the band.

What is the NowCast?

A weighted average of the last 12 hours used for hourly PM reporting in the US. It reacts faster than a 24-hour mean.

Which units do US AQI tables use?

µg/m³ for particles, ppb for NO2 and SO2, ppm for ozone and CO — at 25 °C.

How many hours do I need for a daily mean?

At least 18 of 24 (75%) is the usual rule; 6 of 8 for an 8-hour mean.

Should I round or truncate before the lookup?

Truncate to the table's precision; rounding moves values across band boundaries.

Can I compute an AQI if a site measures only some pollutants?

Yes, from the pollutants available — but say which, because missing pollutants can only lower the result.