"""Mapping NO2 air pollution from space: Sentinel-5P TROPOMI over the UK and north-west Europe. Searches the Microsoft Planetary Computer STAC API for every daytime Sentinel-5P TROPOMI Level-2 NO2 file of 2025 (filling any gaps from the AWS Open Data copy of the same archive), streams only the swath rows over north-west Europe from each NetCDF file, keeps pixels with qa_value > 0.75, and grids them onto a regular 0.05 degree raster. It then builds annual, seasonal and monthly composites, measures how much cloud removes, ranks cities, compares weekdays with Sundays and saves the annual mean as a GeoTIFF (used by the kriging tutorial). Run with: python s5p_no2_air_pollution.py """ import platform import re import time from concurrent.futures import ThreadPoolExecutor, as_completed from importlib.metadata import version from pathlib import Path import fsspec import geopandas as gpd import h5py import matplotlib matplotlib.use("Agg") import matplotlib.dates as mdates import matplotlib.pyplot as plt import numpy as np import pandas as pd import planetary_computer import pystac_client import rasterio import requests from rasterio.transform import from_origin HERE = Path(__file__).parent OUT, DATA = HERE / "output", HERE / "data" GRIDS = DATA / "grids" for d in (OUT, DATA, GRIDS): d.mkdir(exist_ok=True) plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9}) pd.set_option("display.width", 140) print("=== versions ===") print(f"python {platform.python_version()} | pystac-client {pystac_client.__version__} | " f"planetary-computer {version('planetary-computer')} | h5py {h5py.__version__} | " f"fsspec {fsspec.__version__} | numpy {np.__version__} | pandas {pd.__version__}") def section(title: str) -> None: print(f"\n=== {title} ===") # --- Step 1: search the Sentinel-5P archive with STAC ------------------------------- STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1" YEAR = 2025 W, S, E, N = -11.0, 48.0, 10.0, 60.5 # north-west Europe: Ireland to the Ruhr, Paris to Shetland RES = 0.05 # output grid in degrees (about 3.3 km E-W x 5.6 km N-S here) NX, NY = round((E - W) / RES), round((N - S) / RES) catalog = pystac_client.Client.open(STAC_URL, modifier=planetary_computer.sign_inplace) items = [] for month in range(1, 13): # one search per month keeps each response small start = pd.Timestamp(YEAR, month, 1) end = start + pd.offsets.MonthEnd(1) search = catalog.search(collections=["sentinel-5p-l2-netcdf"], bbox=[W, S, E, N], datetime=f"{start:%Y-%m-%d}/{end:%Y-%m-%d}", query={"s5p:product_name": {"eq": "no2"}, "s5p:processing_mode": {"eq": "OFFL"}}) items += list(search.items()) def describe(href, source): """One row per file: the start and end times and orbit number are in the file name.""" parts = href.split("?")[0].split("/")[-1].split("_") start, end, orbit, proc = parts[-6], parts[-5], parts[-4], parts[-2] t0, t1 = pd.Timestamp(start), pd.Timestamp(end) return {"id": f"S5P_L2_NO2____{start}_{end}_{orbit}", "time": t0 + (t1 - t0) / 2, "orbit": int(orbit), "processor": proc, "href": href.split("?")[0], "source": source} # The box also catches night-time orbits that clip the Arctic: keep the daytime overpasses. files = pd.DataFrame([describe(it.assets["no2"].href, "Planetary Computer") for it in items]) files = files[files["time"].dt.hour.between(9, 15)].drop_duplicates("id") section("STAC search") print(f"{len(items)} OFFL NO2 files intersect the box in {YEAR}; {len(files)} are daytime orbits") calendar = pd.date_range(f"{YEAR}-01-01", f"{YEAR}-12-31") missing = calendar.difference(pd.DatetimeIndex(files["time"].dt.normalize().unique())) print("days with no file on the Planetary Computer, by month:", missing.month.value_counts().sort_index().to_dict()) # Fill the gaps from the AWS Open Data copy of the same OFFL archive (bucket meeo-s5p, no login). AWS = "https://meeo-s5p.s3.amazonaws.com" filled = [] for d in missing: xml = requests.get(AWS, params={"list-type": "2", "prefix": f"OFFL/L2__NO2___/{d:%Y/%m/%d}/"}, timeout=60).text for key in re.findall(r"([^<]+\.nc)", xml): row = describe(f"{AWS}/{key}", "AWS Open Data") if 9 <= row["time"].hour <= 15: filled.append(row) files = pd.concat([files, pd.DataFrame(filled)]).sort_values("time").reset_index(drop=True) print(f"added {len(filled)} daytime files for {len(missing)} missing days from the AWS mirror") per_day = files["time"].dt.normalize().value_counts() print(f"{len(files)} files on {per_day.size} days; daytime orbits per day {per_day.min()} to {per_day.max()}") print("files by source:", files["source"].value_counts().to_dict()) print("processor versions:", files["processor"].value_counts().sort_index().to_dict()) # --- Step 2: look inside one file --------------------------------------------------- def open_remote(href): """Open the NetCDF4/HDF5 file over HTTPS; h5py fetches only the chunks we read.""" if "blob.core.windows.net" in href: href = planetary_computer.sign(href) # short-lived SAS token, no account needed f = fsspec.filesystem("https").open(href, block_size=2**20, cache_type="blockcache") return h5py.File(f, "r") sample = files[(files["time"].dt.normalize() == f"{YEAR}-04-08") & files["time"].dt.hour.between(12, 14)].iloc[0] section("one file") print(sample["id"]) with open_remote(sample["href"]) as h: p = h["PRODUCT"] for name in ["nitrogendioxide_tropospheric_column", "qa_value", "latitude"]: d, a = p[name], p[name].attrs attrs = {k: (a[k].decode() if isinstance(a[k], bytes) else float(np.ravel(a[k])[0])) for k in ["units", "scale_factor", "_FillValue", "multiplication_factor_to_convert_to_molecules_percm2"] if k in a} print(f"{name:36s} {d.dtype} shape {d.shape} chunks {d.chunks} {attrs}") # --- Step 3: read the rows over the box, filter on qa_value and grid ------------------ QA_MIN = 75 # qa_value is stored as 0-100 (scale factor 0.01): keep qa_value > 0.75 TO_UMOL = 1e6 # mol m-2 -> umol m-2 def read_swath(href): """Return lon, lat, NO2 (umol m-2) and raw qa (0-100) for the swath rows that cross the box.""" with open_remote(href) as h: p = h["PRODUCT"] lat_mid = p["latitude"][0, :, 225] # centre of the swath, to find the rows we need rows = np.where((lat_mid > S - 1) & (lat_mid < N + 1))[0] if rows.size == 0: return None r0, r1 = rows.min(), rows.max() + 1 lat = p["latitude"][0, r0:r1, :] lon = p["longitude"][0, r0:r1, :] no2 = p["nitrogendioxide_tropospheric_column"][0, r0:r1, :].astype("float64") qa = p["qa_value"][0, r0:r1, :] no2[no2 > 1e30] = np.nan # _FillValue 9.97e36 return lon, lat, no2 * TO_UMOL, qa def grid_swath(lon, lat, no2, qa): """Bin the centres of pixels that pass the QA filter into the 0.05 degree grid: sum and count.""" ok = (qa > QA_MIN) & np.isfinite(no2) & (lon >= W) & (lon < E) & (lat >= S) & (lat < N) ix = np.minimum(((lon[ok] - W) / RES).astype(int), NX - 1) # guard against float rounding at the edges iy = np.minimum(((N - lat[ok]) / RES).astype(int), NY - 1) cell = iy * NX + ix total = np.bincount(cell, weights=no2[ok], minlength=NX * NY).reshape(NY, NX) count = np.bincount(cell, minlength=NX * NY).reshape(NY, NX) return total.astype("float32"), count.astype("uint16"), int(ok.sum()) lon, lat, no2, qa = read_swath(sample["href"]) inbox = (lon >= W) & (lon < E) & (lat >= S) & (lat < N) & np.isfinite(no2) good = inbox & (qa > QA_MIN) section("QA filter on the sample orbit") print(f"{inbox.sum():,} pixels in the box; {good.sum():,} ({good.sum() / inbox.sum():.0%}) pass qa_value > 0.75") print(f"good pixels: median {np.median(no2[good]):.1f}, 99th percentile {np.percentile(no2[good], 99):.1f} umol/m2; " f"{(no2[good] < 0).mean():.1%} are negative") print(f"rejected pixels: median {np.median(no2[inbox & ~good]):.1f} umol/m2") # --- Step 4: process the whole year (cached per orbit) ------------------------------ def process(row): dest = GRIDS / f"{row.id}.npz" if dest.exists(): return dest for attempt in range(4): try: swath = read_swath(row.href) break except Exception as exc: # flaky HTTPS reads: retry with backoff if attempt == 3: raise RuntimeError(f"{row.id}: {exc}") from exc time.sleep(5 * (attempt + 1)) if swath is None: total, count, n = np.zeros((NY, NX), "float32"), np.zeros((NY, NX), "uint16"), 0 else: total, count, n = grid_swath(*swath) np.savez_compressed(dest, total=total, count=count, n=n) return dest t0 = time.time() with ThreadPoolExecutor(max_workers=8) as pool: futures = [pool.submit(process, row) for row in files.itertuples()] for i, fut in enumerate(as_completed(futures), 1): fut.result() if i % 250 == 0: print(f" {i}/{len(files)} orbits gridded", flush=True) dates = pd.DatetimeIndex(sorted(files["time"].dt.normalize().unique())) day_total = np.zeros((len(dates), NY, NX), "float32") day_count = np.zeros((len(dates), NY, NX), "uint16") n_valid = 0 for row in files.itertuples(): g = np.load(GRIDS / f"{row.id}.npz") k = dates.get_loc(row.time.normalize()) day_total[k] += g["total"] day_count[k] += g["count"] n_valid += int(g["n"]) section("year gridded") print(f"{len(files)} orbits on {len(dates)} days: {n_valid:,} good-quality pixels on a {NX} x {NY} grid of {RES} deg") def composite(day_mask, min_obs=5): """Mean of all good pixels on the selected days; NaN where there are fewer than min_obs.""" tot, cnt = day_total[day_mask].sum(0), day_count[day_mask].sum(0).astype("float64") with np.errstate(invalid="ignore", divide="ignore"): mean = tot / cnt mean[cnt < min_obs] = np.nan return mean, cnt annual, annual_n = composite(np.ones(len(dates), bool)) lons = W + RES * (np.arange(NX) + 0.5) lats = N - RES * (np.arange(NY) + 0.5) # --- Step 5: save the annual mean as a GeoTIFF -------------------------------------- tif = OUT / f"s5p_no2_{YEAR}_mean.tif" with rasterio.open(tif, "w", driver="GTiff", width=NX, height=NY, count=2, dtype="float32", crs="EPSG:4326", transform=from_origin(W, N, RES, RES), nodata=np.nan, compress="deflate") as dst: dst.write(annual.astype("float32"), 1) dst.write(annual_n.astype("float32"), 2) dst.set_band_description(1, f"NO2 tropospheric column, umol m-2, qa_value > 0.75, mean of {YEAR}") dst.set_band_description(2, "number of good-quality TROPOMI pixels") section("annual mean") print(f"saved {tif.name}: band 1 = mean column (umol/m2), band 2 = pixel count") print(f"good pixels per 0.05 deg cell in {YEAR}: median {np.median(annual_n):.0f}, " f"5th percentile {np.percentile(annual_n, 5):.0f}, minimum {annual_n.min():.0f}") print(f"annual mean over the box: median {np.nanmedian(annual):.1f}, max {np.nanmax(annual):.1f} umol/m2 " f"({np.nanmax(annual) * 6.02214e-2:.1f} x 10^15 molecules/cm2)") iy, ix = np.unravel_index(np.nanargmax(annual), annual.shape) print(f"highest cell: {lats[iy]:.3f} N, {lons[ix]:.3f} E") # --- Step 6: rank cities and compare weekdays with weekends -------------------------- CITIES = { "London": (51.507, -0.128), "Paris": (48.857, 2.352), "Birmingham": (52.486, -1.890), "Manchester": (53.481, -2.243), "Leeds": (53.800, -1.549), "Liverpool": (53.408, -2.991), "Sheffield": (53.381, -1.470), "Newcastle": (54.978, -1.618), "Bristol": (51.455, -2.588), "Glasgow": (55.864, -4.252), "Edinburgh": (55.953, -3.188), "Cardiff": (51.481, -3.179), "Belfast": (54.597, -5.930), "Dublin": (53.350, -6.260), "Brussels": (50.850, 4.352), "Antwerp": (51.219, 4.402), "Amsterdam": (52.370, 4.895), "Rotterdam": (51.924, 4.478), "Lille": (50.629, 3.057), "Cologne": (50.938, 6.960), "Essen (Ruhr)": (51.456, 7.012), } RADIUS_KM = 10 LON2, LAT2 = np.meshgrid(lons, lats) def city_mask(lat0, lon0): dx = (LON2 - lon0) * 111.32 * np.cos(np.radians(lat0)) dy = (LAT2 - lat0) * 110.57 return np.hypot(dx, dy) <= RADIUS_KM masks = {c: city_mask(*ll) for c, ll in CITIES.items()} def city_means(day_mask): tot, cnt = day_total[day_mask].sum(0), day_count[day_mask].sum(0) return {c: tot[m].sum() / cnt[m].sum() for c, m in masks.items()} weekday = dates.dayofweek < 5 saturday = dates.dayofweek == 5 sunday = dates.dayofweek == 6 rank = pd.DataFrame({"annual": city_means(np.ones(len(dates), bool)), "weekday": city_means(weekday), "saturday": city_means(saturday), "sunday": city_means(sunday)}) rank["sat_vs_wkday_%"] = 100 * (rank["saturday"] / rank["weekday"] - 1) rank["sun_vs_wkday_%"] = 100 * (rank["sunday"] / rank["weekday"] - 1) rank = rank.sort_values("annual", ascending=False) section(f"city ranking: mean column within {RADIUS_KM} km of the centre, {YEAR} (umol/m2)") print(rank.round(1).to_string()) UK = ["London", "Birmingham", "Manchester", "Leeds", "Liverpool", "Sheffield", "Newcastle", "Bristol", "Glasgow", "Edinburgh", "Cardiff", "Belfast"] print(f"UK cities (median): Saturday {rank.loc[UK, 'sat_vs_wkday_%'].median():+.0f}%, " f"Sunday {rank.loc[UK, 'sun_vs_wkday_%'].median():+.0f}% vs weekdays") bg = (LAT2 > 56.8) & (LAT2 < 57.6) & (LON2 > -5.5) & (LON2 < -4.0) # rural Scottish Highlands print(f"rural background (Scottish Highlands): {np.nanmean(annual[bg]):.1f} umol/m2") # --- Step 7: seasonal and monthly composites ----------------------------------------- SEASONS = {"Winter (DJF)": [12, 1, 2], "Spring (MAM)": [3, 4, 5], "Summer (JJA)": [6, 7, 8], "Autumn (SON)": [9, 10, 11]} seasonal = {s: composite(dates.month.isin(m)) for s, m in SEASONS.items()} section("seasonal composites (umol/m2)") for s, (mean, cnt) in seasonal.items(): print(f"{s:13s} London {np.nanmean(mean[masks['London']]):6.1f} | Manchester " f"{np.nanmean(mean[masks['Manchester']]):6.1f} | box median {np.nanmedian(mean):5.1f} | " f"cells with >= 5 good pixels {np.isfinite(mean).mean():.0%}") cover = (day_count > 0).mean(axis=(1, 2)) # share of the box with a good pixel each day monthly = pd.DataFrame({"cover_%": 100 * cover}, index=dates).resample("MS").mean() for c in ["London", "Manchester", "Glasgow"]: m = masks[c] tot = pd.Series(day_total[:, m].sum(1), index=dates).resample("MS").sum() cnt = pd.Series(day_count[:, m].sum(1).astype(float), index=dates).resample("MS").sum() monthly[c] = tot / cnt monthly["London_days"] = pd.Series((day_count[:, masks["London"]].sum(1) > 0).astype(int), index=dates).resample("MS").sum() monthly.index = monthly.index.strftime("%b") section("monthly composites (city columns in umol/m2)") print(monthly.round(1).to_string()) # --- Figures ------------------------------------------------------------------------ ne = DATA / "countries.gpkg" if not ne.exists(): gpd.read_file("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_0_countries.zip").to_file(ne) land = gpd.read_file(ne).clip((W - 1, S - 1, E + 1, N + 1)) CMAP = "Spectral_r" EXTENT = (W, E, S, N) LABEL = "NO$_2$ tropospheric column (µmol/m²)" def basemap(ax, extent=EXTENT): land.boundary.plot(ax=ax, color="#333333", lw=0.4) ax.set_xlim(extent[0], extent[1]) ax.set_ylim(extent[2], extent[3]) ax.set_aspect(1 / np.cos(np.radians(54))) ax.set_xlabel("Longitude (°)", fontsize=7) ax.set_ylabel("Latitude (°)", fontsize=7) ax.tick_params(labelsize=7) # Figure 1: one orbit before and after the QA filter. fig, axes = plt.subplots(1, 2, figsize=(11, 5), constrained_layout=True) for ax, title, keep in [(axes[0], "All retrievals", inbox), (axes[1], "qa_value > 0.75", good)]: sc = ax.scatter(lon[keep], lat[keep], c=no2[keep], s=0.5, cmap=CMAP, vmin=0, vmax=150, rasterized=True) basemap(ax) ax.set_title(title, fontsize=10) fig.colorbar(sc, ax=axes, shrink=0.75, label=LABEL, extend="both") fig.suptitle(f"One Sentinel-5P overpass, {sample['time']:%d %B %Y} {sample['time']:%H:%M} UTC " f"(orbit {sample['orbit']})", fontsize=11) fig.savefig(OUT / "01-single-orbit-qa-filter.png") plt.close(fig) # Figure 2: annual mean with labelled cities. fig, ax = plt.subplots(figsize=(10, 7.6)) im = ax.imshow(annual, extent=EXTENT, cmap=CMAP, vmin=0, vmax=100, interpolation="nearest") basemap(ax) for c in ["London", "Paris", "Birmingham", "Manchester", "Leeds", "Glasgow", "Brussels", "Rotterdam", "Essen (Ruhr)", "Dublin", "Newcastle", "Belfast"]: la, lo = CITIES[c] ax.plot(lo, la, "o", ms=4, mfc="none", mec="black", mew=0.8) ax.annotate(c, (lo, la), xytext=(4, 3), textcoords="offset points", fontsize=7.5, bbox=dict(boxstyle="round,pad=0.15", fc="white", ec="none", alpha=0.8)) fig.colorbar(im, ax=ax, shrink=0.75, label=LABEL, extend="max") ax.set_title(f"Sentinel-5P TROPOMI NO$_2$, mean of {YEAR} (qa_value > 0.75, {RES}° grid)", fontsize=11) fig.savefig(OUT / "02-annual-mean-no2.png") plt.close(fig) # Figure 3: seasonal composites over Great Britain and Ireland. UKX = (-10.8, 2.2, 49.8, 59.2) fig, axes = plt.subplots(1, 4, figsize=(13, 5.4), constrained_layout=True) for ax, (s, (mean, _)) in zip(axes, seasonal.items()): im = ax.imshow(mean, extent=EXTENT, cmap=CMAP, vmin=0, vmax=120, interpolation="nearest") basemap(ax, UKX) ax.set_title(s, fontsize=10) fig.colorbar(im, ax=axes, shrink=0.6, label=LABEL, extend="max", orientation="horizontal", aspect=40) fig.suptitle(f"Seasonal composites, {YEAR} (winter = January, February and December)", fontsize=11) fig.savefig(OUT / "03-seasonal-composites.png") plt.close(fig) # Figure 4: monthly coverage and city series. mdx = pd.date_range(f"{YEAR}-01-01", periods=12, freq="MS") fig, (a1, a2) = plt.subplots(2, 1, figsize=(9, 6.4), sharex=True, constrained_layout=True) a1.bar(mdx, monthly["cover_%"], width=20, color="#3288bd") a1.set_ylabel("Grid cells with a good\npixel on an average day (%)") a1.set_title("How much of the map survives the cloud and quality filter each day", fontsize=10) for c, col in [("London", "#d53e4f"), ("Manchester", "#f46d43"), ("Glasgow", "#5e4fa2")]: a2.plot(mdx, monthly[c].to_numpy(), "o-", color=col, label=c, lw=1.6, ms=4) a2.set_ylabel("Monthly mean column (µmol/m²)") a2.set_title(f"Monthly composites within {RADIUS_KM} km of each city centre", fontsize=10) a2.legend(frameon=False) a2.xaxis.set_major_formatter(mdates.DateFormatter("%b")) for a in (a1, a2): a.grid(axis="y", alpha=0.3) a.spines[["top", "right"]].set_visible(False) fig.savefig(OUT / "04-monthly-coverage-and-cities.png") plt.close(fig) # Figure 5: weekday vs Sunday over England and the Low Countries. wk, _ = composite(weekday) su, _ = composite(sunday) ZX = (-4.5, 7.5, 50.3, 54.6) fig, axes = plt.subplots(3, 1, figsize=(8, 11), constrained_layout=True) for ax, arr, title in [(axes[0], wk, "Monday to Friday"), (axes[1], su, "Sunday")]: im = ax.imshow(arr, extent=EXTENT, cmap=CMAP, vmin=0, vmax=120, interpolation="nearest") basemap(ax, ZX) ax.set_title(f"{title} mean, {YEAR}", fontsize=10) fig.colorbar(im, ax=axes[:2], shrink=0.8, label="µmol/m²", extend="max") im = axes[2].imshow(su - wk, extent=EXTENT, cmap="RdBu_r", vmin=-40, vmax=40, interpolation="nearest") basemap(axes[2], ZX) axes[2].set_title("Sunday minus weekday", fontsize=10) fig.colorbar(im, ax=axes[2], shrink=0.8, label="µmol/m²", extend="both") fig.savefig(OUT / "05-weekday-vs-sunday.png") plt.close(fig) # Figure 6: city ranking with weekday and Sunday bars. r = rank.sort_values("annual") fig, ax = plt.subplots(figsize=(8, 7.5)) y = np.arange(len(r)) ax.barh(y + 0.2, r["weekday"], height=0.4, color="#d53e4f", label="Monday to Friday") ax.barh(y - 0.2, r["sunday"], height=0.4, color="#3288bd", label="Sunday") for yi, (w_, s_, pct) in enumerate(zip(r["weekday"], r["sunday"], r["sun_vs_wkday_%"])): ax.text(max(w_, s_) + 1.5, yi, f"{pct:+.0f}%", va="center", fontsize=7.5) ax.set_yticks(y, r.index) ax.set_xlabel(f"Mean NO$_2$ column within {RADIUS_KM} km of the centre, {YEAR} (µmol/m²)") ax.set_title("Cities ranked by annual mean; label = Sunday change vs weekdays", fontsize=10) ax.legend(frameon=False, loc="lower right") ax.spines[["top", "right"]].set_visible(False) ax.grid(axis="x", alpha=0.3) fig.savefig(OUT / "06-city-ranking-weekday-sunday.png") plt.close(fig) rank.round(2).to_csv(OUT / "city_ranking.csv") monthly.round(3).to_csv(OUT / "monthly_composites.csv") print("\nfigures and CSVs written to", OUT.name)