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.
This commit is contained in:
Emi Griffith
parent
db86277103
commit
b4d8d67825
5 changed files with 300 additions and 1 deletions
|
|
@ -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",
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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:
|
||||
|
|
|
|||
193
data/meteostat.py
Normal file
193
data/meteostat.py
Normal file
|
|
@ -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))
|
||||
|
|
@ -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"
|
||||
|
||||
|
|
|
|||
96
tests/data/test_meteostat.py
Normal file
96
tests/data/test_meteostat.py
Normal file
|
|
@ -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]
|
||||
Loading…
Reference in a new issue