Showcase — wildfire smoke chokes New York (OpenAQ, June 2023)¶
In early June 2023 smoke from the record Canadian wildfires drifted south and gave New York City the worst air quality of any major city in the world — the AQI hit ~484 ('hazardous') on 7 June, a level unseen since the 1960s. This notebook uses the earthlens OpenAQ backend (a global network of ground air-quality monitors) to pull the PM2.5 record across the event and compare it against health guidelines.
Needs a free OpenAQ API key in
OPENAQ_API_KEY; the query skips gracefully without it.
What this notebook does¶
- Query OpenAQ for PM2.5 at New York monitors across early June 2023.
- Plot the city time series of fine-particulate pollution.
- Compare the peak against the WHO guideline, the US standard, and the EPA AQI hazard bands.
PM2.5 (particles < 2.5 µm) penetrates deep into the lungs; the WHO 24-hour guideline is 15 µg/m³ and the US 24-hour standard is 35 µg/m³.
Setup¶
First the imports and a bit of noise-suppression. earthlens provides the unified EarthLens entry point; pandas, numpy, and matplotlib handle the tidy table and the plot.
import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from loguru import logger
from earthlens import EarthLens
logger.disable("earthlens")
logger.disable("pyramids")
1 · Query OpenAQ for PM2.5¶
OpenAQ serves ground-monitor readings as a tidy table, so earthlens treats it as a tabular backend: load() returns a DataFrame with station_id, datetime_utc, value (µg/m³), lat, and lon. We build the request — source, pollutant (variables), date window, and a bounding box over New York City — then load() the readings into memory. The whole fetch is guarded so it skips gracefully without an OPENAQ_API_KEY.
pm = None
if not os.environ.get("OPENAQ_API_KEY"):
print("set OPENAQ_API_KEY to run this query")
else:
try:
nyc = EarthLens(
data_source="openaq",
variables=["pm25"],
start="2023-06-01",
end="2023-06-12",
aoi=[-74.1, 40.5, -73.7, 41.0], # New York City (W, S, E, N)
)
pm = nyc.load(progress_bar=False)
pm["datetime_utc"] = pd.to_datetime(pm["datetime_utc"], utc=True)
print(
f"{len(pm)} PM2.5 records from {pm.station_id.nunique()} stations; "
f"peak {pm.value.max():.0f} µg/m³"
)
except Exception as exc: # noqa: BLE001
print("OpenAQ query skipped:", exc)
156 PM2.5 records from 13 stations; peak 200 µg/m³
2 · The PM2.5 spike against health guidelines¶
The city-wide mean is plotted with each monitor faintly behind it, and the WHO / US thresholds and EPA AQI hazard bands overlaid. PM2.5 leaps from a clean ~10 µg/m³ to many times the standards as the smoke arrives on 6–7 June.
Clean and aggregate¶
First drop obviously-bad readings (negative or implausibly large) and collapse the individual monitors into a city-wide daily mean.
if pm is not None and len(pm):
pm = pm[(pm.value >= 0) & (pm.value < 1000)]
city = pm.groupby(pm.datetime_utc.dt.floor("D")).value.mean()
Plot the time series¶
Each station is drawn faintly in grey behind the heavy black city-mean line, with the EPA AQI hazard bands shaded across the background and the WHO / US guideline lines overlaid.
if pm is not None and len(pm):
fig, ax = plt.subplots(figsize=(11, 5))
for sid, g in pm.groupby("station_id"):
ax.plot(g.datetime_utc, g.value, color="0.8", lw=0.8, zorder=1)
ax.plot(
city.index,
city.values,
color="black",
lw=2,
marker="o",
label="NYC mean PM2.5",
zorder=3,
)
bands = [
(0, 12, "#00e400", "Good"),
(12, 35.4, "#ffff00", "Moderate"),
(35.4, 55.4, "#ff7e00", "Unhealthy (sensitive)"),
(55.4, 150.4, "#ff0000", "Unhealthy"),
(150.4, 250.4, "#8f3f97", "Very unhealthy"),
(250.4, 500, "#7e0023", "Hazardous"),
]
for lo, hi, col, lab in bands:
ax.axhspan(lo, hi, color=col, alpha=0.12)
ax.axhline(15, color="tab:green", ls="--", lw=1, label="WHO 24h guideline (15)")
ax.axhline(35, color="tab:blue", ls="--", lw=1, label="US 24h standard (35)")
ax.set(
ylabel="PM2.5 (µg/m³)",
title="New York City PM2.5 — Canadian wildfire smoke, June 2023",
)
ax.legend(loc="upper left", fontsize=8)
fig.autofmt_xdate()
plt.show()
else:
print("no data — set OPENAQ_API_KEY and rerun")
How far over the limits¶
Express the peak reading as a multiple of the WHO guideline and the US standard to put the smoke event in context.
if pm is not None and len(pm):
peak = pm.value.max()
print(
f"peak PM2.5 {peak:.0f} µg/m³ = {peak/15:.0f}x the WHO 24h guideline, "
f"{peak/35:.1f}x the US standard"
)
peak PM2.5 200 µg/m³ = 13x the WHO 24h guideline, 5.7x the US standard
Recap¶
The earthlens OpenAQ backend turns a global ground-monitor network into a tidy DataFrame, so a public-health air-quality event — smoke from fires 1,000 km away — becomes a few lines of pandas and a time-series plot against the WHO/US thresholds.
Try it yourself¶
- Re-point the bbox at Chicago/Detroit (also smoke-hit), or query
no2/o3. - Compute exposure (person-hours above 'unhealthy') or compare years.
- See the OpenAQ backend reference.