"""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))