MediterraneanWillForge Pipeline Data Report¶
Project: MediterraneanWillForge
Stack: Open-Meteo air quality and ERA5 weather + OpenAQ v3 + WAQI → Delta Lake (Backblaze B2) → Gold marts + Anomaly Detection
Observability: Grafana Cloud (mohamedwillforge.grafana.net)
This notebook reads the Gold tables on Backblaze B2 and builds the daily report: geographic and temporal coverage, pollutant levels, WHO guideline exceedance, wildfire risk, anomalies, seasonal weather and source coverage.
Published artifacts are generated from the retained real-source lakehouse. Test fixtures are never used by report generation.
In [1]:
import sys
from pathlib import Path
REPO_ROOT = Path.cwd().parent
if str(REPO_ROOT) not in sys.path:
sys.path.insert(0, str(REPO_ROOT))
import os
import warnings
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
import pandas as pd
from deltalake import DeltaTable
from dotenv import load_dotenv
from data.reporting.analytics import (
anomaly_daily_rates,
filter_anomaly_model_sources,
filter_report_countries,
latest_reporting_dates,
reporting_dates,
write_readiness_diagnostics,
)
warnings.filterwarnings("ignore")
# Load B2 credentials from .env (MINIO_ENDPOINT, MINIO_ACCESS_KEY, MINIO_SECRET_KEY)
load_dotenv(dotenv_path=os.path.join(os.path.dirname(os.getcwd()), ".env"))
load_dotenv() # fallback: current dir .env
plt.rcParams.update({
"figure.facecolor": "white",
"axes.facecolor": "#f8f8f8",
"axes.grid": True,
"grid.color": "white",
"grid.linewidth": 1.2,
"font.family": "DejaVu Sans",
"axes.spines.top": False,
"axes.spines.right": False,
})
PALETTE = ["#2196F3", "#FF5722", "#4CAF50", "#FFC107", "#9C27B0", "#00BCD4"]
MAX_DATE_TICKS = 12
def date_tick_positions(count: int, max_ticks: int = MAX_DATE_TICKS) -> list[int]:
if count <= 0:
return []
stride = max(1, int(np.ceil(count / max_ticks)))
ticks = list(range(0, count, stride))
if ticks[-1] != count - 1:
ticks.append(count - 1)
return ticks
print("Environment ready.")
Environment ready.
In [2]:
def storage_options() -> dict:
endpoint = os.environ["MINIO_ENDPOINT"]
return {
"AWS_ENDPOINT_URL": endpoint,
"AWS_ACCESS_KEY_ID": os.environ["MINIO_ACCESS_KEY"],
"AWS_SECRET_ACCESS_KEY": os.environ["MINIO_SECRET_KEY"],
"AWS_REGION": os.environ.get("AWS_REGION", "eu-central-003"),
"AWS_ALLOW_HTTP": "true" if endpoint.startswith("http://") else "false",
"AWS_S3_ALLOW_UNSAFE_RENAME": "true",
}
GOLD = os.environ.get("MINIO_BUCKET_GOLD", "med-ops-mohamed-gold")
opts = storage_options()
def read_gold(table: str) -> pd.DataFrame:
path = f"s3://{GOLD}/{table}"
df = DeltaTable(path, storage_options=opts).to_pandas()
if df.empty:
raise RuntimeError(f"Gold table {table} is empty; refusing to publish report artifacts")
print(f" {table}: {len(df):,} rows")
return df
print("Loading Gold tables from B2...")
summary = filter_report_countries(read_gold("daily_country_summary"))
risk = filter_report_countries(read_gold("wildfire_risk_index"))
anomalies = filter_anomaly_model_sources(
filter_report_countries(read_gold("anomaly_alerts"))
)
weather = filter_report_countries(read_gold("daily_country_weather"))
READY_DATES = reporting_dates(summary)
DISPLAY_DATES = latest_reporting_dates(summary)
readiness = write_readiness_diagnostics(summary, anomalies, "reporting_readiness.csv")
print(f"Reporting-ready dates: {len(READY_DATES):,} of {summary['partition_date'].nunique() if not summary.empty else 0:,}")
print(f"Dense chart window: latest {len(DISPLAY_DATES):,} reporting-ready dates")
print("Latest date diagnostics:")
print(readiness.tail(5).to_string(index=False))
print("Done.")
Loading Gold tables from B2...
daily_country_summary: 3,263 rows
wildfire_risk_index: 14,472 rows
anomaly_alerts: 11,696 rows
daily_country_weather: 2,108 rows
Reporting-ready dates: 184 of 191
Dense chart window: latest 45 reporting-ready dates
Latest date diagnostics:
partition_date station_days countries source_count sources coverage_baseline_station_days coverage_ratio is_reporting_ready anomaly_flags rows anomaly_rate_pct
2026-10-07 83 11 3 openaq, openmeteo, waqi 82.0 1.012 True 4 69 5.80
2026-10-08 82 11 3 openaq, openmeteo, waqi 82.0 1.000 True 7 68 10.29
2026-10-09 82 11 3 openaq, openmeteo, waqi 82.0 1.000 True 7 68 10.29
2026-10-10 77 11 3 openaq, openmeteo, waqi 82.0 0.939 True 6 63 9.52
2026-10-11 59 11 3 openaq, openmeteo, waqi 82.0 0.720 True 4 45 8.89
Done.
In [3]:
# ── Summary stats ──────────────────────────────────────────────────────────────
dates = sorted(summary["partition_date"].unique())
countries = sorted(summary["country_code"].dropna().unique())
sources = sorted(summary["source"].dropna().unique())
print(f"Date range : {dates[0]} → {dates[-1]}")
print(f"Dates total : {len(dates)}")
print(f"Countries : {len(countries)} - {', '.join(countries)}")
print(f"Sources : {', '.join(sources)}")
print(f"\nSummary rows : {len(summary):,}")
print(f"Risk index rows : {len(risk):,}")
print(f"Anomaly rows : {len(anomalies):,}")
print(f"Weather rows : {len(weather):,}")
total_stations = risk["station_id"].nunique() if not risk.empty else "N/A"
print(f"Unique stations : {total_stations}")
Date range : 2026-03-31 → 2026-10-11 Dates total : 191 Countries : 11 - DZ, EG, ES, FR, GR, IT, LB, LY, MA, TN, TR Sources : openaq, openmeteo, waqi Summary rows : 3,163 Risk index rows : 13,171 Anomaly rows : 11,076 Weather rows : 2,108 Unique stations : 128
In [4]:
# Station-day heatmap per country, limited to the latest reporting-ready dates.
public_summary = summary[summary["partition_date"].isin(DISPLAY_DATES)].copy()
pivot = (
public_summary
.groupby(["partition_date", "country_code"])["station_count"]
.sum()
.unstack("country_code", fill_value=0)
)
pivot.index = pd.to_datetime(pivot.index)
fig, ax = plt.subplots(figsize=(14, max(4, len(pivot.columns) * 0.45)))
im = ax.imshow(
pivot.values.T,
aspect="auto",
cmap="YlOrRd",
interpolation="nearest",
)
ax.set_yticks(range(len(pivot.columns)))
ax.set_yticklabels(pivot.columns, fontsize=9)
ticks = date_tick_positions(len(pivot.index))
ax.set_xticks(ticks)
ax.set_xticklabels(
[pivot.index[t].strftime("%b %d") for t in ticks],
rotation=60, ha="right", fontsize=7
)
plt.colorbar(im, ax=ax, label="Station-days")
ax.set_title(
f"Station Coverage per Country - Latest {len(DISPLAY_DATES)} Stable Dates",
fontsize=13,
fontweight="bold",
pad=12,
)
ax.set_facecolor("white")
plt.tight_layout()
plt.savefig("coverage_heatmap.png", dpi=150, bbox_inches="tight")
plt.show()
2. Pollutant Concentrations by Country¶
Mean PM2.5 and O3 across all collected dates, split by data source.
In [5]:
# ── Mean PM2.5 by country × source ────────────────────────────────────────────
pm25_by_country = (
summary
.groupby(["country_code", "source"])["mean_pm2_5"]
.mean()
.round(2)
.unstack("source", fill_value=np.nan)
.sort_values(by=summary["source"].mode()[0], ascending=False)
)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# PM2.5
ax = axes[0]
x = np.arange(len(pm25_by_country))
w = 0.35
for i, (src, col) in enumerate(zip(pm25_by_country.columns, PALETTE)):
vals = pm25_by_country[src].values
bars = ax.barh(x + i * w - w/2, vals, height=w, label=src, color=col, alpha=0.85)
ax.axvline(15, color="red", linestyle="--", linewidth=1, label="WHO PM2.5 guideline (15 µg/m³)")
ax.set_yticks(x)
ax.set_yticklabels(pm25_by_country.index)
ax.set_xlabel("Mean PM2.5 (µg/m³)")
ax.set_title("Mean PM2.5 by Country", fontweight="bold")
ax.legend(fontsize=8)
# O3
o3_by_country = (
summary
.groupby(["country_code", "source"])["mean_o3"]
.mean()
.round(2)
.unstack("source", fill_value=np.nan)
.reindex(pm25_by_country.index)
)
ax = axes[1]
for i, (src, col) in enumerate(zip(o3_by_country.columns, PALETTE)):
vals = o3_by_country[src].values
ax.barh(x + i * w - w/2, vals, height=w, label=src, color=col, alpha=0.85)
ax.axvline(100, color="red", linestyle="--", linewidth=1, label="WHO O3 guideline (100 µg/m³)")
ax.set_yticks(x)
ax.set_yticklabels(pm25_by_country.index)
ax.set_xlabel("Mean O3 (µg/m³)")
ax.set_title("Mean O3 by Country", fontweight="bold")
ax.legend(fontsize=8)
plt.suptitle("Air Pollutant Concentrations - All Dates Combined", fontsize=13, fontweight="bold", y=1.02)
plt.tight_layout()
plt.savefig("pollutants_by_country.png", dpi=150, bbox_inches="tight")
plt.show()
3. WHO Guideline Exceedance¶
Percentage of station-days where concentrations exceeded WHO 2021 guidelines.
In [6]:
who_cols = {
"who_pm25_exceed_pct": "PM2.5",
"who_pm10_exceed_pct": "PM10",
"who_no2_exceed_pct": "NO2",
"who_o3_exceed_pct": "O3",
}
who_cols = {k: v for k, v in who_cols.items() if k in summary.columns}
exceed = (
summary
.groupby("country_code")[list(who_cols.keys())]
.mean()
.round(1)
.rename(columns=who_cols)
)
exceed = exceed.sort_values(by=list(who_cols.values())[0], ascending=True)
fig, ax = plt.subplots(figsize=(11, max(4, len(exceed) * 0.55)))
x = np.arange(len(exceed))
n = len(who_cols)
w = 0.18
for i, (col, color) in enumerate(zip(exceed.columns, PALETTE)):
ax.barh(
x + (i - n / 2 + 0.5) * w,
exceed[col],
height=w,
label=col,
color=color,
alpha=0.85,
)
ax.set_yticks(x)
ax.set_yticklabels(exceed.index)
ax.set_xlabel("% of station-days exceeding WHO guideline")
ax.set_title("WHO Guideline Exceedance Rate by Country", fontsize=13, fontweight="bold")
ax.xaxis.set_major_formatter(mticker.PercentFormatter())
ax.legend(title="Pollutant", loc="lower right", fontsize=8)
plt.tight_layout()
plt.savefig("who_exceedance.png", dpi=150, bbox_inches="tight")
plt.show()
print("\nCountries with highest PM2.5 exceedance rate:")
print(exceed[["PM2.5"]].sort_values("PM2.5", ascending=False).head(5).to_string())
Countries with highest PM2.5 exceedance rate:
PM2.5
country_code
LB 71.2
TR 59.4
EG 59.3
ES 38.0
GR 36.6
4. Wildfire Risk Index¶
Composite score (0–100) per station × day: 60% normalised O3 + 40% normalised PM2.5.
Levels: low (<25) · moderate (25–50) · high (50–75) · extreme (>75)
In [7]:
if risk.empty:
print("No wildfire risk data available.")
else:
level_order = ["low", "moderate", "high", "extreme"]
level_colors = {"low": "#4CAF50", "moderate": "#FFC107", "high": "#FF5722", "extreme": "#B71C1C"}
# ── Distribution of risk levels by country ─────────────────────────────────
risk_counts = (
risk
.groupby(["country_code", "risk_level"])
.size()
.unstack("risk_level", fill_value=0)
.reindex(columns=[l for l in level_order if l in risk["risk_level"].unique()], fill_value=0)
)
risk_pct = risk_counts.div(risk_counts.sum(axis=1), axis=0) * 100
risk_pct = risk_pct.sort_values(by=[c for c in ["extreme", "high"] if c in risk_pct.columns], ascending=False)
fig, axes = plt.subplots(1, 2, figsize=(14, max(4, len(risk_pct) * 0.55)))
# Stacked bar: risk level distribution
ax = axes[0]
left = np.zeros(len(risk_pct))
for level in risk_pct.columns:
vals = risk_pct[level].values
ax.barh(risk_pct.index, vals, left=left, label=level, color=level_colors.get(level, "grey"), alpha=0.9)
left += vals
ax.set_xlabel("% of station-days")
ax.set_title("Risk Level Distribution by Country", fontweight="bold")
ax.xaxis.set_major_formatter(mticker.PercentFormatter())
ax.legend(title="Risk level", bbox_to_anchor=(1.01, 1), loc="upper left", fontsize=8)
# Scatter: O3 vs PM2.5 coloured by risk level
ax = axes[1]
for level in level_order:
mask = risk["risk_level"] == level
if mask.any():
ax.scatter(
risk.loc[mask, "pm2_5"],
risk.loc[mask, "ozone"],
c=level_colors[level],
label=level,
alpha=0.5,
s=12,
)
ax.set_xlabel("PM2.5 (µg/m³)")
ax.set_ylabel("O3 (µg/m³)")
ax.set_title("PM2.5 vs O3 - Risk Level", fontweight="bold")
ax.legend(title="Risk level", fontsize=8)
plt.suptitle("Wildfire Risk Index", fontsize=13, fontweight="bold", y=1.01)
plt.tight_layout()
plt.savefig("wildfire_risk.png", dpi=150, bbox_inches="tight")
plt.show()
print("\nTop 10 highest-risk station-days:")
cols = [c for c in ["partition_date", "country_code", "station_name", "risk_index", "risk_level", "pm2_5", "ozone"] if c in risk.columns]
print(risk.nlargest(10, "risk_index")[cols].to_string(index=False))
Top 10 highest-risk station-days:
partition_date country_code station_name risk_index risk_level pm2_5 ozone
2026-08-29 TR Selimiye, Turkey 65.3 high 121.0000 99.2000
2026-08-27 TR Selimiye, Turkey 62.7 high 111.0000 99.2000
2026-08-28 TR Selimiye, Turkey 60.3 high 102.0000 99.2000
2026-07-22 IT LT TASSO 60.0 high NaN 235.0000
2026-08-11 IT LT TASSO 53.0 high NaN 159.0000
2026-08-22 ES EL MAILLO 52.9 high 58.3500 112.0000
2026-08-13 ES Casa De Campo, Madrid, Spain 51.3 high 61.0000 105.1000
2026-04-17 EG Alexandria 50.6 high 83.1375 85.3333
2026-08-23 ES EL MAILLO 50.6 high 44.8042 116.0000
2026-06-23 ES Barcelona, Catalunya, Spain 50.3 high 158.0000 30.9000
5. Anomaly Detection¶
WAQI is excluded from anomaly detection because its iaqi pollutant fields are index values rather than raw concentrations. The model uses Open-Meteo and OpenAQ concentration-compatible rows only.
Isolation Forest trained on the latest 60 days of Silver (PM2.5, O3, NO2).
Contamination rate: 5%.
In [8]:
public_anomalies = anomalies[anomalies["partition_date"].isin(DISPLAY_DATES)].copy()
if public_anomalies.empty:
print("No anomaly data available for the latest reporting-ready dates.")
else:
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# Anomaly rate per date. Raw counts are misleading when latest coverage changes.
ax = axes[0]
anom_by_date = anomaly_daily_rates(public_anomalies)
x = range(len(anom_by_date))
ax.bar(x, anom_by_date["anomaly_rate_pct"], color="#FF5722", alpha=0.8)
ticks = date_tick_positions(len(anom_by_date))
ax.set_xticks(ticks)
ax.set_xticklabels(
[str(anom_by_date.iloc[t]["partition_date"])[-5:] for t in ticks],
rotation=70, ha="right", fontsize=7
)
ax.set_ylabel("Anomaly rate")
ax.set_title("Anomaly Rate per Date", fontweight="bold")
ax.yaxis.set_major_formatter(mticker.PercentFormatter())
# Anomaly rate by country
ax = axes[1]
anom_by_country = (
public_anomalies
.groupby("country_code")["is_anomaly"]
.agg(["sum", "count"])
.rename(columns={"sum": "anomalies", "count": "total"})
)
anom_by_country["rate"] = (anom_by_country["anomalies"] / anom_by_country["total"] * 100).round(1)
anom_by_country = anom_by_country.sort_values("rate", ascending=True)
ax.barh(anom_by_country.index, anom_by_country["rate"], color="#9C27B0", alpha=0.8)
ax.set_xlabel("Anomaly rate")
ax.set_title("Anomaly Rate by Country", fontweight="bold")
ax.xaxis.set_major_formatter(mticker.PercentFormatter())
# Anomaly score distribution
ax = axes[2]
normal = public_anomalies.loc[public_anomalies["is_anomaly"] == 0, "anomaly_score"]
flagged = public_anomalies.loc[public_anomalies["is_anomaly"] == 1, "anomaly_score"]
ax.hist(normal, bins=30, color="#4CAF50", alpha=0.7, label=f"Normal ({len(normal):,})")
ax.hist(flagged, bins=15, color="#FF5722", alpha=0.85, label=f"Anomaly ({len(flagged):,})")
ax.set_xlabel("Isolation Forest score (lower = more anomalous)")
ax.set_ylabel("Count")
ax.set_title("Score Distribution", fontweight="bold")
ax.legend(fontsize=8)
plt.suptitle(
f"Anomaly Detection - Latest {len(DISPLAY_DATES)} Stable Dates",
fontsize=13,
fontweight="bold",
y=1.01,
)
plt.tight_layout()
plt.savefig("anomaly_detection.png", dpi=150, bbox_inches="tight")
plt.show()
print("\nTop 10 most anomalous readings in chart window (lowest score = most deviant):")
cols = [c for c in ["partition_date", "country_code", "station_name", "pm2_5", "ozone", "nitrogen_dioxide", "anomaly_score"] if c in public_anomalies.columns]
print(public_anomalies[public_anomalies["is_anomaly"] == 1].nsmallest(10, "anomaly_score")[cols].to_string(index=False))
Top 10 most anomalous readings in chart window (lowest score = most deviant):
partition_date country_code station_name pm2_5 ozone nitrogen_dioxide anomaly_score
2026-10-11 TR Istanbul 43.0458 31.6250 60.7833 -0.758853
2026-10-08 TR Istanbul 31.9333 28.3750 50.9333 -0.732370
2026-10-09 TR Istanbul 32.9625 38.7083 52.3792 -0.725616
2026-10-10 TR Istanbul 30.6542 37.2500 53.9417 -0.724391
2026-10-10 ES ERANDIO 57.0000 NaN 33.0000 -0.705097
2026-09-08 ES Sant Adrià de Besòs 24.0000 7.0000 34.7000 -0.701841
2026-08-31 GR PIREAUS-1 9.7500 15.0000 97.6000 -0.699781
2026-10-08 GR KORDELIO NaN 10.0000 107.0000 -0.697782
2026-08-27 GR PIREAUS-1 13.1000 22.0000 81.1000 -0.697734
2026-08-31 DZ Algiers 36.5417 100.2083 18.8542 -0.696011
6. Top Anomaly of the Latest Date¶
The single most anomalous reading on the most recent reporting-ready date, shown against that day's spread across every station. The marked line is the flagged station; the histogram is every other station reporting that pollutant, so the reason the Isolation Forest singled it out is visible rather than asserted.
A station with no sensor for a pollutant is labelled as having no reading and is not plotted as zero.
In [9]:
# Spotlight the single most anomalous reading on the latest reporting-ready date.
# The three panels show that day's spread across every station, with the flagged
# station marked, so the reason it was flagged shows on the chart.
TOP_DATE = max(DISPLAY_DATES) if len(DISPLAY_DATES) else None
day = anomalies[anomalies["partition_date"] == TOP_DATE].copy() if TOP_DATE else pd.DataFrame()
if day.empty:
print("No anomaly data for the latest reporting-ready date.")
else:
top = day.nsmallest(1, "anomaly_score").iloc[0]
POLLUTANTS = [
("pm2_5", "PM2.5", "#2196F3"),
("ozone", "Ozone", "#4CAF50"),
("nitrogen_dioxide", "NO2", "#9C27B0"),
]
fig, axes = plt.subplots(1, 3, figsize=(16, 4.6))
for ax, (col, label, colour) in zip(axes, POLLUTANTS):
values = pd.to_numeric(day[col], errors="coerce").dropna()
station_value = pd.to_numeric(pd.Series([top.get(col)]), errors="coerce").iloc[0]
if values.empty:
ax.text(0.5, 0.5, f"No {label} readings\non this date",
ha="center", va="center", fontsize=10, color="#777")
ax.set_xticks([]); ax.set_yticks([])
ax.set_title(label, fontweight="bold")
continue
ax.hist(values, bins=min(20, max(5, len(values) // 2)),
color=colour, alpha=0.55, edgecolor="white")
if pd.isna(station_value):
# A station with no sensor for this pollutant reports nothing at all.
ax.set_title(f"{label} - station has no reading", fontweight="bold", fontsize=10)
else:
pct = (values < station_value).mean() * 100
ax.axvline(station_value, color="#FF5722", linewidth=2.5, zorder=5)
# Label on whichever side of the line has room, so it never clips.
lo, hi = ax.get_xlim()
on_right = station_value > lo + 0.6 * (hi - lo)
ax.annotate(
f"{station_value:,.1f}\n{pct:.0f}th pct",
xy=(station_value, ax.get_ylim()[1] * 0.92),
xytext=(-7 if on_right else 7, 0), textcoords="offset points",
ha="right" if on_right else "left",
color="#FF5722", fontweight="bold", fontsize=9, va="top",
)
ax.set_title(f"{label} - median {values.median():,.1f}", fontweight="bold", fontsize=10)
ax.set_xlabel(f"{label} (ug/m3)")
ax.set_ylabel("Stations")
flagged = "flagged anomaly" if int(top.get("is_anomaly", 0)) == 1 else "most deviant reading"
plt.suptitle(
f"Top Anomaly {TOP_DATE} - {top['station_name']} ({top['country_code']}, {top['source']}) "
f"- score {top['anomaly_score']:.3f} - {flagged}",
fontsize=13, fontweight="bold", y=1.02,
)
plt.tight_layout()
plt.savefig("top_anomaly.png", dpi=150, bbox_inches="tight")
plt.show()
cols = [c for c in ["country_code", "station_name", "source", "pm2_5", "ozone",
"nitrogen_dioxide", "anomaly_score", "is_anomaly"] if c in day.columns]
print(f"\nMost anomalous readings on {TOP_DATE} (lowest score = most deviant):")
print(day.nsmallest(5, "anomaly_score")[cols].to_string(index=False))
Most anomalous readings on 2026-10-11 (lowest score = most deviant):
country_code station_name source pm2_5 ozone nitrogen_dioxide anomaly_score is_anomaly
TR Istanbul openmeteo 43.0458 31.6250 60.7833 -0.758853 1
LB Beirut openmeteo 29.7167 68.8333 15.5208 -0.639730 1
GR PATISION openaq NaN 5.4200 41.0000 -0.599652 1
EG Cairo openmeteo 19.1750 60.0000 41.1542 -0.582444 1
LY Tripoli openmeteo 2.7792 87.2917 1.1958 -0.541580 0
7. Seasonal Weather¶
Daily highs, lows, wind, rain and Saharan dust for the same 12 grid points the air-quality model uses, from Open-Meteo's ERA5 reanalysis.
The section follows the season of the latest date. Winter charts nightly lows
and cold waves. Spring charts Saharan dust, which peaks over the central and
eastern basin in spring. Summer charts daily highs and heatwaves. Autumn charts
rain and the days with 10 mm or more, the flash-flood season. An active
heatwave or cold wave takes over in any season. The mapping lives in
data/reporting/season.py.
A day counts as unusually hot when it clears two bars at once: 30 C in absolute terms, and the 90th percentile of that station's own previous 30 days. The relative bar stops an ordinary Mediterranean August from reading as one long heatwave. The absolute bar stops a mild April day from qualifying because the fortnight before it was milder still. Three such days in a row make a heatwave, and a heatwave whose high reaches 40 C is extreme. Cold alerts mirror the same shape against the 10th percentile and 5 C.
Below the chart sits a generated paragraph on the country leading the season's measure: what weather at that level does to people or to the ground, and whether it arrived gradually or as a swing. It is labelled as generated text. The pipeline computes every figure in it, including the size of each swing.
In [10]:
# The latest date sets the season, and the season sets what this section charts.
from IPython.display import Markdown, display
from data.reporting.season import focus_for, rank_countries, season_of
latest_weather_date = weather["partition_date"].max()
latest_weather = weather[weather["partition_date"] == latest_weather_date]
focus = focus_for(weather, latest_weather_date)
season = season_of(latest_weather_date)
display(Markdown(f"**{season.title()} focus: {focus.title}.** Latest date {latest_weather_date}."))
public_weather = weather[weather["partition_date"].isin(DISPLAY_DATES)].copy()
public_weather[focus.column] = pd.to_numeric(public_weather[focus.column], errors="coerce")
focus_pivot = (
public_weather
.pivot_table(
index="partition_date",
columns="country_code",
values=focus.column,
aggfunc="max",
)
.sort_index()
)
focus_pivot.index = pd.to_datetime(focus_pivot.index)
event_days = (
weather[focus.is_event(weather)]
.groupby("country_code")
.size()
.sort_values()
)
fig, axes = plt.subplots(
1, 2,
figsize=(15, max(4.5, len(focus_pivot.columns) * 0.5)),
gridspec_kw={"width_ratios": [2.2, 1]},
)
ax = axes[0]
im = ax.imshow(focus_pivot.values.T, aspect="auto", cmap=focus.cmap, interpolation="nearest")
ax.set_yticks(range(len(focus_pivot.columns)))
ax.set_yticklabels(focus_pivot.columns, fontsize=9)
ticks = date_tick_positions(len(focus_pivot.index))
ax.set_xticks(ticks)
ax.set_xticklabels(
[focus_pivot.index[t].strftime("%b %d") for t in ticks],
rotation=60, ha="right", fontsize=7,
)
plt.colorbar(im, ax=ax, label=focus.label)
ax.set_title(
f"{focus.label.split(' (')[0]} per Country - Latest {len(DISPLAY_DATES)} Stable Dates",
fontsize=12, fontweight="bold", pad=12,
)
ax.set_facecolor("white")
ax = axes[1]
if event_days.empty:
ax.text(0.5, 0.5, f"{focus.event_label}: none on record", ha="center",
va="center", transform=ax.transAxes, fontsize=10)
ax.set_axis_off()
else:
ax.barh(event_days.index, event_days.values, color=PALETTE[1], alpha=0.85)
ax.set_xlabel(focus.event_label)
ax.set_title(f"{focus.title} Days per Country", fontsize=12, fontweight="bold", pad=12)
ax.xaxis.set_major_locator(mticker.MaxNLocator(integer=True))
plt.tight_layout()
plt.savefig("seasonal_weather.png", dpi=150, bbox_inches="tight")
plt.show()
conditions = latest_weather.sort_values(
focus.column, ascending=focus.lowest_first, na_position="last"
)[[
"country_code", "temp_max_c", "temp_min_c", "precipitation_mm", "condition",
"wind_level", "dust_level", "heat_alert", "heat_streak_days",
"cold_alert", "cold_streak_days",
]]
conditions.columns = [
"Country", "High C", "Low C", "Rain mm", "Sky", "Wind", "Dust",
"Heat alert", "Heat streak", "Cold alert", "Cold streak",
]
print(f"Conditions on {latest_weather_date}:\n")
print(conditions.to_string(index=False))
heat_alerted = (latest_weather["heat_alert"] != "none").sum()
cold_alerted = (latest_weather["cold_alert"] != "none").sum()
print(f"\nUnder a heat alert: {heat_alerted} of {len(latest_weather)} countries")
print(f"Under a cold alert: {cold_alerted} of {len(latest_weather)} countries")
ranking = rank_countries(weather, focus, str(latest_weather_date))
if not ranking.empty:
print(f"Leader on the {focus.measure}: {ranking.index[0]} at {ranking.iloc[0]:.1f}")
Autumn focus: Rain. Latest date 2026-10-11.
Conditions on 2026-10-11:
Country High C Low C Rain mm Sky Wind Dust Heat alert Heat streak Cold alert Cold streak
GR 23.9 17.4 10.2 rain windy none none 0 none 0
MA 24.5 19.3 5.0 rain breezy none none 0 none 0
IT 23.8 14.5 0.2 drizzle breezy none none 0 none 0
DZ 26.5 14.7 0.0 cloudy breezy none none 0 none 0
EG 29.2 17.8 0.0 cloudy breezy none none 0 none 0
ES 22.0 11.4 0.0 cloudy breezy none none 0 none 0
FR 18.4 14.2 0.0 cloudy windy none none 0 none 0
LY 27.1 22.6 0.0 cloudy windy none none 0 none 0
TN 25.9 17.5 0.0 cloudy breezy none none 0 none 0
TR 25.2 15.8 0.0 cloudy calm none none 0 none 0
Under a heat alert: 0 of 10 countries
Under a cold alert: 0 of 10 countries
Leader on the highest average daily rain over the last 10 days: ES at 10.2
In [11]:
# The season's spotlight note, written by Gemini from the figures above. No key
# or a spent quota leaves it out of the run.
import json as _json
from pathlib import Path as _Path
_brief_path = _Path("ai_brief.json")
_spotlight = None
if _brief_path.exists():
_spotlight = _json.loads(_brief_path.read_text(encoding="utf-8")).get(
"weather_spotlight"
)
if not _spotlight:
print("No spotlight note for this run (no GEMINI_API_KEY, or the call failed).")
else:
_figs = _spotlight.get("figures") or {}
_bits = []
if _figs.get("high_c") is not None:
_bits.append(f"high {_figs['high_c']} C")
if _figs.get("low_c") is not None:
_bits.append(f"low {_figs['low_c']} C")
if _figs.get("window_rain_total_mm"):
_bits.append(
f"{_figs['window_rain_total_mm']} mm of rain over "
f"{_figs.get('days_of_history', '?')} days"
)
if _figs.get("dust_today") is not None:
_bits.append(f"dust {_figs['dust_today']} ug/m3")
for _alert in ("heat_alert", "cold_alert"):
if _figs.get(_alert) and _figs[_alert] != "none":
_bits.append(f"alert `{_figs[_alert]}`")
if _figs.get("biggest_day_to_day_change_c") is not None:
_bits.append(
f"sharpest day-to-day change {_figs['biggest_day_to_day_change_c']:+g} C"
)
_model = _spotlight.get("model") or "Gemini"
_date = _spotlight.get("date", "the latest")
_caption = ", ".join(_bits) if _bits else "figures above"
_season = (_spotlight.get("season") or "").title()
display(
Markdown(
f"""#### {_season} spotlight: {_spotlight['country_code']}
Led on the {_spotlight.get('measure', 'measure above')}: {_spotlight.get('value')} {_spotlight.get('unit', '')}.
{_spotlight['paragraph']}
*Generated by `{_model}` from {_date} Gold figures: {_caption}.*"""
)
)
Autumn spotlight: ES¶
Led on the highest average daily rain over the last 10 days: 10.2 mm.
Spain recorded the highest average daily rainfall over the last 10 days at 10.2 mm. Within this period, one day experienced 31.9 mm of rain, an amount that, in Mediterranean climates, can lead to runoff and contribute to flash flood conditions if the ground is dry. Temperatures across the window showed a spread of 15.1 C, ranging from a high of 25.4 C to a low of 10.3 C, reflecting varying conditions over the period. No heat or cold alerts are currently in effect.
Generated by gemini-2.5-flash from 2026-10-11 Gold figures: high 22.0 C, low 11.4 C, 102.5 mm of rain over 10 days, dust 0.0 ug/m3, sharpest day-to-day change -5.5 C.
8. AI Daily Brief¶
Written by Gemini from this run's Gold layer, covering both air quality and the weather above it. The anomaly fact-check is given the flagged reading, that day's distribution and that country's weather, and searches the web for a real-world explanation such as a wildfire, dust intrusion, heatwave or traffic event. It reports an absence of evidence as an absence of evidence. Most anomalies go uncorroborated, and the model does not invent a cause.
This section is generated text and is labelled as such. Every figure it cites comes from the pipeline; anything it claims about the wider world carries its source link.
In [12]:
import json
from pathlib import Path
from IPython.display import Markdown, display
brief_path = Path("ai_brief.json")
if not brief_path.exists():
print("No AI brief for this run (GEMINI_API_KEY unset, or the call failed).")
else:
brief = json.loads(brief_path.read_text(encoding="utf-8"))
lines = [
f"*Generated by `{brief.get('model') or 'Gemini'}` from the "
f"{brief.get('date', 'latest')} Gold layer. Figures come from the pipeline; "
"the fact-check additionally searched the web and cites its sources.*",
"",
]
fc = brief.get("fact_check")
if fc:
lines += [
"### Anomaly fact-check",
"",
f"**{fc['station']} ({fc['country_code']})**: {fc['verdict']}",
"",
]
if fc.get("sources"):
lines.append("Sources consulted:")
lines += [f"- [{s['title']}]({s['url']})" for s in fc["sources"]]
lines.append("")
elif fc.get("searched") is False:
lines += ["*No web sources were returned for this reading.*", ""]
briefings = brief.get("briefings") or []
if briefings:
lines += ["### Country briefings", ""]
for entry in sorted(briefings, key=lambda b: b.get("country_code", "")):
lines.append(f"- **{entry['country_code']}**: {entry['briefing']}")
lines.append("")
display(Markdown("\n".join(lines)))
print(f"AI brief: {len(briefings)} country briefing(s), "
f"fact-check={'yes' if fc else 'no'}, generated {brief.get('generated_at', '?')}")
Generated by gemini-3.7-flash from the 2026-10-11 Gold layer. Figures come from the pipeline; the fact-check additionally searched the web and cites its sources.
Anomaly fact-check¶
Istanbul (TR): The reading is plausible but highly anomalous, with PM2.5 at 43.0 ug/m3 (median 8.9 ug/m3) and NO2 at 60.8 ug/m3 (median 15.1 ug/m3) being significantly elevated, while Ozone at 31.6 ug/m3 (median 56.8 ug/m3) is considerably lower. The calm winds and cloudy conditions are consistent with local pollutant accumulation. However, as the date (2026-10-11) is in the future, no corroborating real-world events like industrial incidents or wildfires can be found to explain these specific anomalies.
No web sources were returned for this reading.
Country briefings¶
- DZ: Measurements from a single station showed NO2 averaging 29.6 ug/m3, above the WHO guideline of 25 ug/m3, while PM2.5 remained clean at 8.4 ug/m3. The day was cloudy and breezy with a high of 26.5 C.
- EG: Based on two stations, NO2 averaged 28.2 ug/m3, exceeding the WHO guideline of 25 ug/m3, and half of station-days were above the 15 ug/m3 threshold for PM2.5 with a mean of 13.6 ug/m3 and peak of 19.2 ug/m3. Weather was cloudy and breezy, reaching a high of 29.2 C with gusts up to 36.4 km/h.
- ES: All stations exceeded the WHO PM2.5 limit, averaging 39.3 ug/m3 and reaching a maximum of 59.0 ug/m3, while PM10 averaged 18.3 ug/m3. Skies were cloudy and breezy with a daily high of 22.0 C.
- FR: Reporting from a single station showed fine particulate matter at 4.2 ug/m3 and NO2 at 5.4 ug/m3, both well within WHO limits. Strong winds brought gusts up to 54.7 km/h under cloudy skies, with temperatures reaching 18.4 C.
- GR: Rain delivered 10.2 mm of precipitation amid windy conditions with gusts of 42.8 km/h and a high of 23.9 C. Daily PM2.5 averaged 15.8 ug/m3, slightly over the WHO 15 ug/m3 guideline, with 35.3% of station readings exceeding the threshold and a peak of 30.0 ug/m3.
- IT: Air quality remained within WHO guidelines across the network, with PM2.5 averaging 8.1 ug/m3 and NO2 at 15.9 ug/m3. Conditions were mild and breezy with light drizzle of 0.2 mm and an afternoon high of 23.8 C.
- LB: A single reporting station recorded PM2.5 at 29.7 ug/m3 and PM10 at 47.7 ug/m3, both exceeding their respective WHO guidelines of 15 ug/m3 and 45 ug/m3. NO2 averaged 15.5 ug/m3 and ozone reached 68.8 ug/m3; no weather data was reported.
- LY: Data from a single station showed particulate levels well below guidelines, with PM2.5 at 2.8 ug/m3 and PM10 at 6.6 ug/m3, while ozone averaged 87.3 ug/m3. Strong winds reached gusts of 40.7 km/h under cloudy skies, with a high of 27.1 C.
- MA: Based on a single station, PM2.5 averaged 11.6 ug/m3 and PM10 stood at 20.0 ug/m3, keeping within WHO guidelines. The day saw 5.0 mm of rain and breezy conditions, reaching a daytime high of 24.5 C.
- TN: Readings from a single station indicated clean air, with PM2.5 at 4.5 ug/m3 and NO2 at 3.3 ug/m3. The weather was cloudy and breezy with peak gusts of 37.8 km/h and a maximum temperature of 25.9 C.
- TR: Pollution levels were elevated across all stations, with PM2.5 averaging 79.9 ug/m3 (reaching 139.0 ug/m3), PM10 at 65.4 ug/m3, and NO2 at 45.7 ug/m3, all surpassing WHO guidelines. The day was calm and cloudy with an afternoon high of 25.2 C.
AI brief: 11 country briefing(s), fact-check=yes, generated 2026-10-12T14:45:27.930401+00:00
9. Source Coverage by Data Source¶
Station-days contributed by each source (Open-Meteo, OpenAQ, WAQI) per country. WAQI fills the gaps in LB and MA, where OpenAQ has few stations.
In [13]:
public_summary = summary[summary["partition_date"].isin(READY_DATES)].copy()
if public_summary.empty:
print("No summary data available for reporting-ready dates.")
else:
coverage = (
public_summary
.groupby(["country_code", "source"])["station_count"]
.sum()
.unstack("source", fill_value=0)
)
coverage = coverage.loc[coverage.sum(axis=1).sort_values().index]
fig, axes = plt.subplots(1, 2, figsize=(14, max(4, len(coverage) * 0.55)))
# Absolute station-days per source
ax = axes[0]
left = np.zeros(len(coverage))
for src, color in zip(coverage.columns, PALETTE):
vals = coverage[src].values
ax.barh(coverage.index, vals, left=left, label=src, color=color, alpha=0.85)
left += vals
ax.set_xlabel("Station-days")
ax.set_title("Station-Days per Country by Source", fontweight="bold")
ax.legend(title="Source", fontsize=8)
# Percentage share per source
coverage_pct = coverage.div(coverage.sum(axis=1), axis=0) * 100
ax = axes[1]
left = np.zeros(len(coverage_pct))
for src, color in zip(coverage_pct.columns, PALETTE):
vals = coverage_pct[src].values
ax.barh(coverage_pct.index, vals, left=left, label=src, color=color, alpha=0.85)
left += vals
ax.set_xlabel("% of station-days")
ax.set_title("Source Share per Country (%)", fontweight="bold")
ax.xaxis.set_major_formatter(mticker.PercentFormatter())
ax.legend(title="Source", fontsize=8)
plt.suptitle(
"Data Source Coverage - Station-Days per Country",
fontsize=13, fontweight="bold", y=1.01
)
plt.tight_layout()
plt.savefig("source_coverage.png", dpi=150, bbox_inches="tight")
plt.show()
print("\nStation-days per source per country:")
print(coverage.to_string())
Station-days per source per country: source openaq openmeteo waqi country_code FR 0 178 0 TN 0 181 0 LY 0 182 0 DZ 0 184 0 LB 0 184 0 MA 0 183 227 EG 0 366 294 TR 0 184 858 IT 2763 184 0 GR 3116 184 395 ES 2753 153 878
In [14]:
n_anomalies = int(anomalies["is_anomaly"].sum()) if not anomalies.empty else 0
n_anom_rows = len(anomalies)
anom_rate = f"{n_anomalies / n_anom_rows:.1%}" if n_anom_rows else "N/A"
top_risk_country = (
risk.groupby("country_code")["risk_index"].mean().idxmax()
if not risk.empty else "N/A"
)
hottest_country = (
weather.groupby("country_code")["temp_max_c"].max().idxmax()
if not weather.empty else "N/A"
)
coldest_country = (
weather.groupby("country_code")["temp_min_c"].min().idxmin()
if not weather.empty else "N/A"
)
cold_wave_country_days = int(
weather["cold_alert"].isin(["cold_wave", "severe_cold_wave"]).sum()
) if not weather.empty else 0
heatwave_country_days = int(
weather["heat_alert"].isin(["heatwave", "extreme_heatwave"]).sum()
) if not weather.empty else 0
highest_exceed_country = (
summary.groupby("country_code")["who_pm25_exceed_pct"].mean().idxmax()
if "who_pm25_exceed_pct" in summary.columns else "N/A"
)
report = pd.DataFrame([
["Date range", f"{dates[0]} → {dates[-1]}"],
["Total dates collected", len(dates)],
["Countries covered", f"{len(countries)} ({', '.join(countries)})"],
["Data sources", ", ".join(sources)],
["Unique stations", risk["station_id"].nunique() if not risk.empty else "N/A"],
["Gold summary rows", f"{len(summary):,}"],
["Gold risk index rows", f"{len(risk):,}"],
["Gold anomaly rows", f"{n_anom_rows:,}"],
["Gold weather rows", f"{len(weather):,}"],
["Anomalies flagged", f"{n_anomalies} ({anom_rate})"],
["Highest avg risk country", top_risk_country],
["Highest PM2.5 exceedance", highest_exceed_country],
["Hottest country", hottest_country],
["Heatwave country-days", f"{heatwave_country_days:,}"],
["Coldest country", coldest_country],
["Cold-wave country-days", f"{cold_wave_country_days:,}"],
["Season focus", f"{season} ({focus.title.lower()})"],
], columns=["Metric", "Value"])
print(report.to_string(index=False))
Metric Value
Date range 2026-03-31 → 2026-10-11
Total dates collected 191
Countries covered 11 (DZ, EG, ES, FR, GR, IT, LB, LY, MA, TN, TR)
Data sources openaq, openmeteo, waqi
Unique stations 128
Gold summary rows 3,163
Gold risk index rows 13,171
Gold anomaly rows 11,076
Gold weather rows 2,108
Anomalies flagged 561 (5.1%)
Highest avg risk country LB
Highest PM2.5 exceedance LB
Hottest country TN
Heatwave country-days 90
Coldest country TR
Cold-wave country-days 0
Season focus autumn (rain)
Regenerate the Report¶
Install requirements-report.txt, then run from the repository root:
cd docs
jupyter nbconvert --to notebook --execute pipeline_report.ipynb --output pipeline_report
jupyter nbconvert --to html pipeline_report.ipynb --output pipeline_report