From b4d8d6782507bcb4987ca418d1cb8a13644dced2 Mon Sep 17 00:00:00 2001 From: Emi Griffith Date: Wed, 22 Jul 2026 21:39:09 -0700 Subject: [PATCH] Add Meteostat gust supplier for the gust-less backup sources NASA POWER (history) and MET Norway (forecast) carry no wind gusts, but gust is a graded metric. This adds data/meteostat.py, which finds the nearest Meteostat station to a cell and reads its daily peak gust (wpgt) from the keyless gzipped bulk endpoints, filling the gust column on the NASA POWER history frame. Where no station is within ~100 km it estimates from sustained wind (wind * GUST_FACTOR). Bulk files (station list + per-station daily) are cached on disk so steady-state network IO is near zero, and any lookup/fetch failure degrades to pure estimation so a history fetch never fails. A constant estimate factor makes an estimated gust redundant with wind, so real signal comes only where a station backs it. Activates once NASA POWER becomes the primary history source; Open-Meteo already carries its own gusts. New history_gust / meteostat_stations metrics phases map to the meteostat source. --- core/metrics.py | 4 + data/climate.py | 6 +- data/meteostat.py | 193 +++++++++++++++++++++++++++++++++++ tests/core/test_metrics.py | 2 + tests/data/test_meteostat.py | 96 +++++++++++++++++ 5 files changed, 300 insertions(+), 1 deletion(-) create mode 100644 data/meteostat.py create mode 100644 tests/data/test_meteostat.py diff --git a/core/metrics.py b/core/metrics.py index d1355d6..36b085b 100644 --- a/core/metrics.py +++ b/core/metrics.py @@ -49,6 +49,10 @@ _PHASE_SOURCE = { "geocode": "nominatim", "history_nasa": "nasa-power", "forecast_metno": "met-norway", + # Wind gusts for the gust-less backup sources (NASA POWER / MET Norway), read + # from Meteostat's keyless bulk endpoints (station list + per-station daily). + "history_gust": "meteostat", + "meteostat_stations": "meteostat", "reverse_geocode": "nominatim", } diff --git a/data/climate.py b/data/climate.py index 6a434f9..dda822f 100644 --- a/data/climate.py +++ b/data/climate.py @@ -18,6 +18,7 @@ from core import audit from core import metrics import paths from data import climate_store +from data import meteostat from data import store CACHE_DIR = os.path.join(paths.DATA_DIR, "cache") @@ -524,7 +525,10 @@ def _fetch_history_nasa(cell: dict) -> pl.DataFrame: "format": "JSON", } r = _request(NASA_POWER_URL, params, 180, phase="history_nasa") - return _nasa_to_frame(r.json()["properties"]["parameter"]) + df = _nasa_to_frame(r.json()["properties"]["parameter"]) + # NASA POWER carries no gusts; fill from the nearest Meteostat station, falling + # back to an estimate from sustained wind where no station is in range. + return meteostat.fill_gusts(cell["center_lat"], cell["center_lon"], df) def _metno_to_frame(props: dict) -> pl.DataFrame: diff --git a/data/meteostat.py b/data/meteostat.py new file mode 100644 index 0000000..10df786 --- /dev/null +++ b/data/meteostat.py @@ -0,0 +1,193 @@ +"""Wind-gust supplier via Meteostat's free, keyless bulk data. + +The archive sources that back Thermograph when Open-Meteo is unavailable — NASA +POWER (history) and MET Norway (forecast) — carry no wind gusts, yet gust is a +graded metric. This module fills that gap: it finds the nearest Meteostat weather +station to a grid cell and reads its daily peak-gust record (``wpgt``) straight +from Meteostat's gzipped bulk endpoints (no API key, no per-query quota). + +Where no station is within range — oceans, remote cells — it falls back to an +*estimate* from that day's sustained wind (``wind * GUST_FACTOR``). Note that a +constant factor makes an estimated gust move in lockstep with wind, so the gust +grade carries independent signal only where a real station backs it; estimation +just keeps the metric populated rather than leaving a hole. + +Everything degrades safely: any station-lookup or fetch failure falls through to +estimation, so a Meteostat outage can never break a history fetch. Bulk files are +downloaded once and cached on disk (like the GeoNames dump), so steady-state +network IO is near zero. + +Wired in at the NASA-POWER history path (data/climate.py); it activates for real +once NASA POWER becomes the primary history source. Open-Meteo already carries its +own gusts, so nothing here runs while Open-Meteo is primary. +""" +import datetime +import gzip +import json +import math +import os +import threading + +import httpx +import numpy as np +import polars as pl + +from core import audit +from core import metrics +import paths + +META_DIR = os.path.join(paths.DATA_DIR, "meteostat") +STATIONS_URL = "https://bulk.meteostat.net/v2/stations/lite.json.gz" +DAILY_URL = "https://bulk.meteostat.net/v2/daily/{}.csv.gz" + +# How far a station may sit from a cell's center before we treat it as "no +# station" and estimate. ~100 km keeps gusts representative of the same synoptic +# conditions without demanding a station in every remote cell. +MAX_STATION_KM = 100.0 +# Sustained-wind → gust ratio for the estimate used where no station is in range. +# ~1.4 is a typical over-land daily-max gust factor; see the module note on why an +# estimated gust is redundant with wind (only measured gusts add real signal). +GUST_FACTOR = 1.4 +KMH_TO_MPH = 0.6213711922 + +# Meteostat daily bulk CSV is headerless; columns are positional. wpgt (peak gust, +# km/h) is column 8; the date (YYYY-MM-DD) is column 0. +_CSV_DATE = 0 +_CSV_WPGT = 8 + +_STATIONS: tuple[list[str], np.ndarray, np.ndarray] | None = None +_STATIONS_LOCK = threading.Lock() + + +def _fetch_bytes(url: str, phase: str, timeout: float = 60.0) -> bytes: + """GET raw bytes, gunzipping a gzip payload (bulk files are ``.gz``; httpx only + auto-decodes Content-Encoding, not a gzip *file*). Records an outbound metric.""" + r = httpx.get(url, timeout=timeout, follow_redirects=True) + r.raise_for_status() + metrics.record_outbound(phase, "ok") + data = r.content + if data[:2] == b"\x1f\x8b": # gzip magic + return gzip.decompress(data) + return data + + +def _parse_stations(raw: bytes) -> tuple[list[str], np.ndarray, np.ndarray]: + """Parse the bulk stations JSON into (ids, lat[], lon[]). Entries without an id + or coordinates are skipped.""" + ids: list[str] = [] + lats: list[float] = [] + lons: list[float] = [] + for e in json.loads(raw): + loc = e.get("location") or {} + sid = e.get("id") + lat = loc.get("latitude") + lon = loc.get("longitude") + if sid is None or lat is None or lon is None: + continue + ids.append(str(sid)) + lats.append(float(lat)) + lons.append(float(lon)) + return ids, np.asarray(lats), np.asarray(lons) + + +def _load_stations() -> tuple[list[str], np.ndarray, np.ndarray]: + """Load the station index (downloaded once into META_DIR, cached forever).""" + global _STATIONS + if _STATIONS is not None: + return _STATIONS + with _STATIONS_LOCK: + if _STATIONS is not None: + return _STATIONS + os.makedirs(META_DIR, exist_ok=True) + cache = os.path.join(META_DIR, "stations.json") + if os.path.exists(cache): + with open(cache, "rb") as fh: + raw = fh.read() + else: + raw = _fetch_bytes(STATIONS_URL, phase="meteostat_stations") + with open(cache, "wb") as fh: + fh.write(raw) + _STATIONS = _parse_stations(raw) + return _STATIONS + + +def nearest_station(lat: float, lon: float) -> str | None: + """Meteostat station id nearest to (lat, lon) within MAX_STATION_KM, or None.""" + ids, lats, lons = _load_stations() + if not ids: + return None + # Great-circle distance, vectorized over all stations. + rlat = math.radians(lat) + dlat = np.radians(lats - lat) + dlon = np.radians(lons - lon) + a = (np.sin(dlat / 2) ** 2 + + math.cos(rlat) * np.cos(np.radians(lats)) * np.sin(dlon / 2) ** 2) + km = 6371.0 * 2 * np.arcsin(np.sqrt(a)) + i = int(np.argmin(km)) + return ids[i] if km[i] <= MAX_STATION_KM else None + + +def _parse_daily_csv(text: str) -> dict[datetime.date, float]: + """Map a station's daily bulk CSV to {date: peak_gust_mph}, skipping days with + no recorded gust.""" + out: dict[datetime.date, float] = {} + for line in text.splitlines(): + cols = line.split(",") + if len(cols) <= _CSV_WPGT: + continue + raw = cols[_CSV_WPGT].strip() + if not raw: + continue + try: + out[datetime.date.fromisoformat(cols[_CSV_DATE])] = float(raw) * KMH_TO_MPH + except ValueError: + continue + return out + + +def daily_gusts(station_id: str) -> dict[datetime.date, float]: + """Measured daily peak gusts (mph) for a station, keyed by date. Cached on disk + per station (bulk history barely changes; the recent tail staleness is immaterial + against a cache that is itself refetched only on a cold history pull).""" + os.makedirs(META_DIR, exist_ok=True) + cache = os.path.join(META_DIR, f"daily_{station_id}.csv") + if os.path.exists(cache): + with open(cache, "r", encoding="utf-8") as fh: + return _parse_daily_csv(fh.read()) + text = _fetch_bytes(DAILY_URL.format(station_id), phase="history_gust").decode("utf-8") + with open(cache, "w", encoding="utf-8") as fh: + fh.write(text) + return _parse_daily_csv(text) + + +def fill_gusts(lat: float, lon: float, df: pl.DataFrame) -> pl.DataFrame: + """Populate an all-null ``gust`` column from the nearest Meteostat station, + estimating from ``wind`` where the station has no value (or no station is in + range). A no-op if the frame already carries gusts or lacks the needed columns. + + Never raises: any station-lookup / fetch failure degrades to pure estimation, + so the caller's history fetch always succeeds.""" + if df.is_empty() or "gust" not in df.columns or "wind" not in df.columns: + return df + if df["gust"].drop_nulls().len() > 0: # source already provided gusts + return df + measured: dict[datetime.date, float] = {} + try: + sid = nearest_station(lat, lon) + if sid: + measured = daily_gusts(sid) + except Exception as e: # noqa: BLE001 - gusts are a nicety; never fail history + audit.log_event("error", {"phase": "history_gust", "error": repr(e)}) + measured = {} + dates = df["date"].to_list() + wind = df["wind"].to_list() + gust = [] + for d, w in zip(dates, wind): + m = measured.get(d) + if m is not None: + gust.append(m) + elif w is not None: + gust.append(w * GUST_FACTOR) + else: + gust.append(None) + return df.with_columns(pl.Series("gust", gust, dtype=pl.Float64)) diff --git a/tests/core/test_metrics.py b/tests/core/test_metrics.py index 3561f74..010f246 100644 --- a/tests/core/test_metrics.py +++ b/tests/core/test_metrics.py @@ -17,6 +17,8 @@ def test_phase_source_map_covers_every_source(): assert m.source_for_phase("geocode") == "nominatim" # forward geocoding moved off Open-Meteo assert m.source_for_phase("history_nasa") == "nasa-power" assert m.source_for_phase("forecast_metno") == "met-norway" + assert m.source_for_phase("history_gust") == "meteostat" + assert m.source_for_phase("meteostat_stations") == "meteostat" assert m.source_for_phase("reverse_geocode") == "nominatim" assert m.source_for_phase("unknown-phase") == "other" diff --git a/tests/data/test_meteostat.py b/tests/data/test_meteostat.py new file mode 100644 index 0000000..1182de4 --- /dev/null +++ b/tests/data/test_meteostat.py @@ -0,0 +1,96 @@ +"""Unit tests for the Meteostat gust supplier — all hermetic (no network): the +station index and per-station daily fetch are monkeypatched, and the pure parse +helpers are fed literal payloads.""" +import datetime +import json + +import numpy as np +import polars as pl +import pytest + +from data import meteostat + + +def _frame(gust=(None, None, None), wind=(10.0, 20.0, None)): + return pl.DataFrame( + { + "date": [datetime.date(2020, 1, 1), datetime.date(2020, 1, 2), + datetime.date(2020, 1, 3)], + "wind": list(wind), + "gust": pl.Series(list(gust), dtype=pl.Float64), + } + ) + + +def test_parse_daily_csv_converts_kmh_to_mph_and_skips_blanks(): + csv = "\n".join([ + "2020-01-01,,,,,,,,50.0,,", # wpgt (col 8) = 50 km/h + "2020-01-02,,,,,,,,,,", # no gust -> skipped + "2020-01-03,,,,,,,,80.0,,", # wpgt = 80 km/h + "bad-date,,,,,,,,50.0,,", # unparseable date -> skipped + "x,y", # too few columns -> skipped + ]) + out = meteostat._parse_daily_csv(csv) + assert set(out) == {datetime.date(2020, 1, 1), datetime.date(2020, 1, 3)} + assert out[datetime.date(2020, 1, 1)] == pytest.approx(50.0 * meteostat.KMH_TO_MPH) + assert out[datetime.date(2020, 1, 3)] == pytest.approx(80.0 * meteostat.KMH_TO_MPH) + + +def test_parse_stations_skips_entries_without_id_or_coords(): + raw = json.dumps([ + {"id": "A", "location": {"latitude": 47.6, "longitude": -122.3}}, + {"id": "B", "location": {"latitude": 51.5, "longitude": -0.1}}, + {"id": "C", "location": {}}, # no coords -> skip + {"location": {"latitude": 1, "longitude": 2}}, # no id -> skip + ]).encode() + ids, lats, lons = meteostat._parse_stations(raw) + assert ids == ["A", "B"] + assert lats.tolist() == [47.6, 51.5] + assert lons.tolist() == [-122.3, -0.1] + + +def test_nearest_station_picks_closest_within_range(monkeypatch): + monkeypatch.setattr(meteostat, "_STATIONS", + (["seattle", "london"], + np.array([47.6, 51.5]), np.array([-122.3, -0.1]))) + assert meteostat.nearest_station(47.61, -122.31) == "seattle" + # Middle of the Pacific — both stations are far beyond MAX_STATION_KM. + assert meteostat.nearest_station(0.0, -160.0) is None + + +def test_fill_gusts_prefers_measured_then_estimates(monkeypatch): + monkeypatch.setattr(meteostat, "nearest_station", lambda lat, lon: "s1") + monkeypatch.setattr(meteostat, "daily_gusts", + lambda sid: {datetime.date(2020, 1, 1): 33.0}) + out = meteostat.fill_gusts(47.6, -122.3, _frame()) + gust = out["gust"].to_list() + assert gust[0] == pytest.approx(33.0) # measured wins + assert gust[1] == pytest.approx(20.0 * meteostat.GUST_FACTOR) # estimated from wind + assert gust[2] is None # no wind -> no estimate + + +def test_fill_gusts_estimates_everywhere_with_no_station(monkeypatch): + monkeypatch.setattr(meteostat, "nearest_station", lambda lat, lon: None) + out = meteostat.fill_gusts(0.0, 0.0, _frame()) + gust = out["gust"].to_list() + assert gust[0] == pytest.approx(10.0 * meteostat.GUST_FACTOR) + assert gust[1] == pytest.approx(20.0 * meteostat.GUST_FACTOR) + assert gust[2] is None + + +def test_fill_gusts_degrades_on_fetch_error(monkeypatch): + def _boom(lat, lon): + raise RuntimeError("meteostat down") + monkeypatch.setattr(meteostat, "nearest_station", _boom) + # Must not raise — falls back to pure estimation. + out = meteostat.fill_gusts(47.6, -122.3, _frame()) + assert out["gust"].to_list()[1] == pytest.approx(20.0 * meteostat.GUST_FACTOR) + + +def test_fill_gusts_noop_when_source_has_gusts(monkeypatch): + def _fail(*a, **k): + raise AssertionError("must not consult Meteostat when gusts already present") + monkeypatch.setattr(meteostat, "nearest_station", _fail) + df = _frame(gust=(40.0, 41.0, 42.0)) + out = meteostat.fill_gusts(47.6, -122.3, df) + assert out["gust"].to_list() == [40.0, 41.0, 42.0]