From db862771037840f31676e3eb37bc9e739674d780 Mon Sep 17 00:00:00 2001 From: Emi Griffith Date: Wed, 22 Jul 2026 21:29:26 -0700 Subject: [PATCH 1/6] Move forward geocoding off Open-Meteo to local index + Nominatim MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit /suggest is now served purely from the local GeoNames index — no external geocoder call per keystroke. /geocode answers from the local index first and falls back to OSM Nominatim /search only on a miss (or while the index is still loading), covering the neighbourhood/postcode/tiny-village/native-name long tail the cities dump lacks. Nominatim shares the reverse geocoder's lock and ~1/sec pacing since both hit the same host, and only the low-volume /geocode miss path reaches it. Removes the Open-Meteo geocoding dependency entirely; the "geocode" metrics phase now maps to the nominatim source. --- core/metrics.py | 4 ++- data/climate.py | 72 +++++++++++++++++++++++++++----------- data/places.py | 25 +++++++------ tests/core/test_metrics.py | 2 +- tests/web/test_api.py | 53 +++++++++++++++++++++++++--- web/app.py | 31 ++++++++++------ 6 files changed, 138 insertions(+), 49 deletions(-) diff --git a/core/metrics.py b/core/metrics.py index ff6e4ee..d1355d6 100644 --- a/core/metrics.py +++ b/core/metrics.py @@ -44,7 +44,9 @@ _PHASE_SOURCE = { "history_fetch": "open-meteo-archive", "history_topup": "open-meteo-archive", "recent_forecast_fetch": "open-meteo-forecast", - "geocode": "open-meteo-geocoding", + # Forward geocoding moved off Open-Meteo to Nominatim (same host as reverse) — + # /geocode now backstops the local GeoNames index on a miss. + "geocode": "nominatim", "history_nasa": "nasa-power", "forecast_metno": "met-norway", "reverse_geocode": "nominatim", diff --git a/data/climate.py b/data/climate.py index f88fc43..6a434f9 100644 --- a/data/climate.py +++ b/data/climate.py @@ -910,24 +910,54 @@ def reverse_geocode(lat: float, lon: float) -> str | None: return label -def geocode(name: str, count: int = 5) -> list[dict]: - """Look up places by name worldwide via Open-Meteo's geocoder.""" - r = _request( - "https://geocoding-api.open-meteo.com/v1/search", - {"name": name, "count": count, "language": "en", "format": "json"}, - 30, - phase="geocode", - ) - results = r.json().get("results", []) or [] - return [ - { - "name": g.get("name"), - "admin1": g.get("admin1"), - "country": g.get("country"), - "country_code": g.get("country_code"), - "lat": g.get("latitude"), - "lon": g.get("longitude"), - "population": g.get("population"), - } - for g in results - ] +def geocode_nominatim(name: str, count: int = 5) -> list[dict]: + """Forward place-name lookup via OpenStreetMap Nominatim (keyless). + + Replaces the former Open-Meteo geocoder. Nominatim covers the long tail the + local GeoNames index can't — neighbourhoods, postcodes, sub-1000-population + villages, and alternate/native-language spellings — so /geocode falls back to + it on a local miss. It carries no population, so results keep Nominatim's own + relevance order; name/admin/country map straight across. + + Shares the reverse geocoder's lock and ~1/sec pacing: both hit the same host, + so serializing them together keeps total Nominatim traffic under the usage + policy. Only the low-volume /geocode miss path reaches here — autocomplete + (/suggest) is served purely from the local index and never calls out. + """ + global _revgeo_last + with _REVGEO_LOCK: + wait = _REVGEO_MIN_INTERVAL - (time.monotonic() - _revgeo_last) + if wait > 0: + time.sleep(wait) + try: + r = _request( + "https://nominatim.openstreetmap.org/search", + {"q": name, "format": "jsonv2", "addressdetails": 1, + "limit": count, "accept-language": "en"}, + 30, + phase="geocode", + headers={"User-Agent": "Thermograph/0.1 (local weather grading app)"}, + ) + rows = r.json() or [] + finally: + _revgeo_last = time.monotonic() + out = [] + for g in rows: + a = g.get("address", {}) or {} + # jsonv2 gives a `name` for named features; for addresses/areas fall back to + # the most specific place component, then the head of display_name. + label = (g.get("name") + or a.get("city") or a.get("town") or a.get("village") + or a.get("hamlet") or a.get("suburb") or a.get("neighbourhood") + or (g.get("display_name") or "").split(",")[0].strip() or None) + cc = (a.get("country_code") or "").upper() or None + out.append({ + "name": label, + "admin1": a.get("state") or a.get("province") or a.get("region"), + "country": a.get("country"), + "country_code": cc, + "lat": float(g["lat"]) if g.get("lat") is not None else None, + "lon": float(g["lon"]) if g.get("lon") is not None else None, + "population": None, # Nominatim has no population; order is relevance-based + }) + return out diff --git a/data/places.py b/data/places.py index c2faadb..ffd8af9 100644 --- a/data/places.py +++ b/data/places.py @@ -1,22 +1,25 @@ """Typo-tolerant place-name index for search suggestions. -Open-Meteo's geocoder (climate.geocode) only matches exact spellings — one -mistyped letter and the query returns nothing. This module keeps a local index -of world places (a GeoNames cities dump, downloaded once into data/geonames/ -and reused across restarts) so /suggest can answer instantly and tolerate a -single-letter typo — a substituted, missing, or extra letter, or two adjacent -letters swapped — anywhere in the query, including the first character. +This module keeps a local index of world places (a GeoNames cities dump, +downloaded once into data/geonames/ and reused across restarts) so /suggest can +answer instantly, offline, and tolerate a single-letter typo — a substituted, +missing, or extra letter, or two adjacent letters swapped — anywhere in the +query, including the first character. It also builds a vocabulary of place-name tokens so a typo'd word in a multi-word query can be respelled against words the index knows ("pest -seattle" → "west seattle"); the endpoint verifies those candidates against the -upstream geocoder, which covers neighbourhood-level places the cities dump -lacks. +seattle" → "west seattle"). + +`suggest()` still takes an injected ``upstream`` callable so the module stays +free of the fetch layer, but /suggest now passes a no-op: autocomplete is served +purely from this index and makes no per-keystroke network call. (/geocode, the +explicit search, is the path that falls back to Nominatim on a local miss for +neighbourhoods/postcodes/tiny places the cities dump lacks — see web/app.py.) The index loads in a background thread (start_loading()). Until it's ready — or forever, if the download fails — search()/corrections() return None/empty -and /suggest degrades to the plain upstream geocoder, so the app never needs -this data to boot or serve. +and /suggest simply returns no results, so the app never needs this data to +boot or serve. """ import bisect import heapq diff --git a/tests/core/test_metrics.py b/tests/core/test_metrics.py index 5a4a743..3561f74 100644 --- a/tests/core/test_metrics.py +++ b/tests/core/test_metrics.py @@ -14,7 +14,7 @@ def test_phase_source_map_covers_every_source(): assert m.source_for_phase("history_fetch") == "open-meteo-archive" assert m.source_for_phase("history_topup") == "open-meteo-archive" assert m.source_for_phase("recent_forecast_fetch") == "open-meteo-forecast" - assert m.source_for_phase("geocode") == "open-meteo-geocoding" + 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("reverse_geocode") == "nominatim" diff --git a/tests/web/test_api.py b/tests/web/test_api.py index 9ac1e6b..f6893fa 100644 --- a/tests/web/test_api.py +++ b/tests/web/test_api.py @@ -6,6 +6,7 @@ from fastapi.testclient import TestClient from web import app as appmod from data import climate +from data import places @pytest.fixture @@ -182,11 +183,53 @@ def test_cell_bundle_matches_per_view_payloads(client): assert r304.status_code == 304 -def test_suggest_falls_back_to_upstream_geocoder(client, monkeypatch): - upstream = [{"name": "Seattle", "admin1": "Washington", "country": "United States", - "country_code": "US", "lat": 47.6, "lon": -122.33, "population": 737015}] - monkeypatch.setattr(climate, "geocode", lambda q, count=5: list(upstream)) - r = client.get("/thermograph/api/v2/suggest", params={"q": "seattle-fallback-probe"}) +_LOCAL_HIT = [{"name": "Seattle", "admin1": "Washington", "country": "United States", + "country_code": "US", "lat": 47.6, "lon": -122.33, "population": 737015, + "match": "prefix"}] + + +def test_geocode_local_first_skips_network(client, monkeypatch): + """A local GeoNames hit answers /geocode with no outbound Nominatim call, and + the internal prefix/fuzzy tag is stripped from the wire shape.""" + monkeypatch.setattr(places, "search", lambda q, limit=5: list(_LOCAL_HIT)) + + def _boom(*a, **k): + raise AssertionError("Nominatim must not be called when the local index hits") + monkeypatch.setattr(climate, "geocode_nominatim", _boom) + + r = client.get("/thermograph/api/v2/geocode", params={"q": "seattle"}) + assert r.status_code == 200 + res = r.json()["results"] + assert res[0]["name"] == "Seattle" + assert "match" not in res[0] + + +def test_geocode_falls_back_to_nominatim_on_local_miss(client, monkeypatch): + """A local miss (or a still-loading index) falls back to Nominatim /search — + the path that covers neighbourhoods/postcodes the cities dump lacks.""" + monkeypatch.setattr(places, "search", lambda q, limit=5: []) # genuine miss + nom = [{"name": "West Seattle", "admin1": "Washington", "country": "United States", + "country_code": "US", "lat": 47.57, "lon": -122.38, "population": None}] + called = [] + monkeypatch.setattr(climate, "geocode_nominatim", + lambda q, count=5: (called.append(q), nom)[1]) + + r = client.get("/thermograph/api/v2/geocode", params={"q": "west seattle"}) + assert r.status_code == 200 + assert r.json()["results"][0]["name"] == "West Seattle" + assert called == ["west seattle"] + + +def test_suggest_is_local_only_no_network(client, monkeypatch): + """Autocomplete is served purely from the local index — even a reachable + Nominatim must never be hit per keystroke.""" + monkeypatch.setattr(places, "search", lambda q, limit=5: list(_LOCAL_HIT)) + + def _boom(*a, **k): + raise AssertionError("/suggest must not make an outbound geocoder call") + monkeypatch.setattr(climate, "geocode_nominatim", _boom) + + r = client.get("/thermograph/api/v2/suggest", params={"q": "seattle"}) assert r.status_code == 200 body = r.json() assert body["results"][0]["name"] == "Seattle" diff --git a/web/app.py b/web/app.py index 990a122..5f3e19b 100644 --- a/web/app.py +++ b/web/app.py @@ -2,7 +2,6 @@ import asyncio import contextlib import datetime -import functools import hashlib import hmac import json @@ -399,9 +398,22 @@ def _cached_response(request, run, kind, cell, key, token, build, cache_placeles # --- endpoints ---------------------------------------------------------------- +_GEOCODE_LIMIT = 5 + + def api_geocode(q: str = Query(..., min_length=1)): + """Forward place search. The local GeoNames index answers first (no network); + only a genuine miss — or a query the index can't hold (neighbourhoods, + postcodes, tiny villages, native-language names) — falls back to Nominatim. + `search` returns None while the index is still loading, which also falls + through to the network so search keeps working during warmup.""" + local = places.search(q, _GEOCODE_LIMIT) + if local: + # Strip the internal prefix/fuzzy tag so the wire shape matches the + # Nominatim branch (and the former Open-Meteo contract) exactly. + return {"results": [{k: v for k, v in r.items() if k != "match"} for r in local]} try: - return {"results": climate.geocode(q)} + return {"results": climate.geocode_nominatim(q, _GEOCODE_LIMIT)} except Exception as e: # noqa: BLE001 - surface upstream failures to the client raise HTTPException(status_code=502, detail=f"geocoding failed: {e}") @@ -409,21 +421,20 @@ def api_geocode(q: str = Query(..., min_length=1)): _SUGGEST_LIMIT = 5 -@functools.lru_cache(maxsize=1024) -def _suggest_upstream(q: str) -> tuple: - """Open-Meteo lookup for /suggest, memoized per query string — type-ahead - re-asks the same prefixes constantly (backspacing, retyping). Failures - raise and are not cached, so a transient upstream error doesn't stick.""" - return tuple(climate.geocode(q, count=_SUGGEST_LIMIT)) +def _no_upstream(q: str) -> tuple: + """Autocomplete is served purely from the local GeoNames index — no per-keystroke + call to any external geocoder. Typo-correction still runs against the local + index, so respellings ("pest seattle" → "west seattle") keep working offline.""" + return () def api_suggest(q: str = Query(..., min_length=1)): """Type-ahead location suggestions: the top 5 places for a (possibly typo'd) query prefix; `corrected` reports the respelling that produced the results ("pest seattle" → "west seattle"), or null. The ranking/typo policy - lives in places.suggest.""" + lives in places.suggest, now driven entirely by the local index.""" try: - results, corrected = places.suggest(q, _suggest_upstream, _SUGGEST_LIMIT) + results, corrected = places.suggest(q, _no_upstream, _SUGGEST_LIMIT) except Exception as e: # noqa: BLE001 - nothing at all to serve raise HTTPException(status_code=502, detail=f"geocoding failed: {e}") return {"results": results, "corrected": corrected} -- 2.45.2 From b4d8d6782507bcb4987ca418d1cb8a13644dced2 Mon Sep 17 00:00:00 2001 From: Emi Griffith Date: Wed, 22 Jul 2026 21:39:09 -0700 Subject: [PATCH 2/6] 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] -- 2.45.2 From 30153dbc64fa31c16a07ecf1bfd5bc7e02b50543 Mon Sep 17 00:00:00 2001 From: Emi Griffith Date: Wed, 22 Jul 2026 22:02:30 -0700 Subject: [PATCH 3/6] Flip history primary to NASA POWER, Open-Meteo to fallback MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The live history path now tries NASA POWER first (keyless, independent of Open-Meteo) and falls back to the Open-Meteo archive only when NASA is unavailable — the reverse of before. Gusts NASA lacks are filled from Meteostat inside _fetch_history_nasa (added in the prior change). Both sources must still return a plausibly-full span (MIN_ARCHIVE_DAYS) before being cached as complete; a short/partial response is rejected and the other source is tried. _fetch_history_nasa now accepts a start/end range so it also serves the recent tail top-up: _topup_tail fetches NASA-first via the new _fetch_history_tail helper (Open-Meteo range as fallback, skipped during its rate-limit cooldown), removing the last per-active-cell Open-Meteo call from the history path. The Open-Meteo cooldown no longer gates the top-up, since NASA is not subject to it. Open-Meteo remains wired as the fallback (deletion is a later step). Tests updated to the new source order plus tail-fetch coverage. --- data/climate.py | 61 ++++++++++++++++++++------------ tests/data/test_climate.py | 71 ++++++++++++++++++++++++++++++-------- 2 files changed, 95 insertions(+), 37 deletions(-) diff --git a/data/climate.py b/data/climate.py index dda822f..3fa340b 100644 --- a/data/climate.py +++ b/data/climate.py @@ -512,16 +512,22 @@ def _nasa_to_frame(param: dict) -> pl.DataFrame: return _finalize_approximated(df) -def _fetch_history_nasa(cell: dict) -> pl.DataFrame: - """Backup history fetch from NASA POWER (used when Open-Meteo is unavailable).""" - end = (datetime.date.today() - datetime.timedelta(days=ARCHIVE_LATENCY_DAYS)).strftime("%Y%m%d") +def _fetch_history_nasa(cell: dict, start: str = NASA_START, + end: "str | None" = None) -> pl.DataFrame: + """Primary history fetch from NASA POWER (keyless, independent of Open-Meteo). + + `start`/`end` accept NASA's YYYYMMDD or ISO YYYY-MM-DD; `end` defaults to the + archive-latency cutoff. Serves both the full multi-decade pull and the tail + top-up range. NASA carries no gusts, so they're filled from Meteostat below.""" + if end is None: + end = (datetime.date.today() - datetime.timedelta(days=ARCHIVE_LATENCY_DAYS)).strftime("%Y%m%d") params = { "parameters": "T2M_MAX,T2M_MIN,PRECTOTCORR,WS10M_MAX,RH2M", "community": "RE", "latitude": cell["center_lat"], "longitude": cell["center_lon"], - "start": NASA_START, - "end": end, + "start": start.replace("-", ""), + "end": end.replace("-", ""), "format": "JSON", } r = _request(NASA_POWER_URL, params, 180, phase="history_nasa") @@ -615,10 +621,21 @@ def _read_history_cache(path): return df, time.time() - os.path.getmtime(path) +def _fetch_history_tail(cell: dict, start_date: str, end_date: str) -> pl.DataFrame: + """Fetch a recent history range to top up the cached tail — NASA POWER primary, + Open-Meteo archive fallback (mirroring the full-history source order). The + Open-Meteo fallback is skipped while it's in a rate-limit cooldown.""" + try: + return _fetch_history_nasa(cell, start=start_date, end=end_date) + except Exception: # noqa: BLE001 - NASA unavailable; try Open-Meteo unless it's cooling down + if time.time() < _archive_cooldown_until: + raise + return _fetch_history_range(cell, start_date, end_date) + + def _topup_tail(cell: dict, df: pl.DataFrame) -> pl.DataFrame: """Append newly-available archive days to a cached record. Best-effort and serialized per cell; a small incremental fetch, not the full multi-decade pull.""" - global _archive_cooldown_until cell_id = cell["id"] with _cell_lock(cell_id): fresh = _read_history_backed(cell_id) # another thread may have just refreshed it @@ -626,15 +643,13 @@ def _topup_tail(cell: dict, df: pl.DataFrame) -> pl.DataFrame: df, age_s = fresh if age_s < HISTORY_TOPUP_INTERVAL: return df - if time.time() < _archive_cooldown_until: - return df expected = datetime.date.today() - datetime.timedelta(days=ARCHIVE_LATENCY_DAYS) cached_max = df["date"].max() if cached_max >= expected: _touch_history_backed(cell_id) # tail already current — reset the hourly timer return df try: - recent = _fetch_history_range( + recent = _fetch_history_tail( cell, (cached_max + datetime.timedelta(days=1)).isoformat(), expected.isoformat()) except Exception as e: # noqa: BLE001 - keep the cached record on failure if is_rate_limit(e): @@ -699,30 +714,30 @@ def _load_history(cell: dict) -> tuple[pl.DataFrame, dict]: def _serve_stale(): return _with_doy(stale), {"cached": True, "stale": True, "cache_age_days": None} - # Fetch. Open-Meteo is primary (richer: gusts + apparent temp); NASA POWER is - # the backup when Open-Meteo is unavailable (network error or its daily limit). - # Skip Open-Meteo entirely while it's in a rate-limit cooldown. + # Fetch. NASA POWER is the primary source (keyless, independent of Open-Meteo); + # the gusts it lacks are filled from Meteostat inside _fetch_history_nasa. + # Open-Meteo is the fallback when NASA is unavailable — still guarded by its + # rate-limit cooldown. Both sources must return a plausibly-full span before + # being accepted and cached as a complete record (a short/partial response is + # rejected rather than held indefinitely). df = None source = None - if time.time() >= _archive_cooldown_until: + try: + fetched = _fetch_history_nasa(cell) + if fetched.height >= MIN_ARCHIVE_DAYS: + df = fetched + source = "nasa-power" + except Exception: # noqa: BLE001 - primary unavailable; fall through to Open-Meteo + pass + if df is None and time.time() >= _archive_cooldown_until: try: fetched = _fetch_history(cell) - # Guard a self-hosted archive that isn't backfilled yet: an all-null or - # recent-only response reduces to a near-empty frame, which would - # otherwise be cached indefinitely as a complete record and never - # fall to NASA. Require a plausibly-full span before accepting it. if fetched.height >= MIN_ARCHIVE_DAYS: df = fetched source = "open-meteo" except Exception as e: # noqa: BLE001 if is_rate_limit(e): _note_rate_limit(e) - if df is None: - try: - df = _fetch_history_nasa(cell) - source = "nasa-power" - except Exception: # noqa: BLE001 - backup unavailable too - df = None if df is None: if stale is not None: return _serve_stale() diff --git a/tests/data/test_climate.py b/tests/data/test_climate.py index e8e5782..cfc0475 100644 --- a/tests/data/test_climate.py +++ b/tests/data/test_climate.py @@ -292,35 +292,78 @@ def _full_history_frame(n=4000): })) -def test_short_archive_is_rejected_and_falls_back_to_nasa(monkeypatch, tmp_path): - """A self-hosted archive that isn't backfilled yet returns too few days; the load - must reject it (not cache a near-empty history as complete) and use NASA instead.""" +def test_nasa_is_primary_history_source(monkeypatch, tmp_path): + """NASA POWER is the primary history source: a full NASA frame is accepted and + Open-Meteo is never consulted.""" monkeypatch.setattr(climate, "CACHE_DIR", str(tmp_path)) monkeypatch.setattr(climate, "_archive_cooldown_until", 0.0) - cell = {"id": "short_cell", "center_lat": 47.6, "center_lon": -122.3} + cell = {"id": "nasa_primary", "center_lat": 47.6, "center_lon": -122.3} - short = climate._to_frame(_om_daily(3)) # 2 usable days « threshold - full = _full_history_frame() # >= MIN_ARCHIVE_DAYS - assert short.height < climate.MIN_ARCHIVE_DAYS <= full.height - monkeypatch.setattr(climate, "_fetch_history", lambda c: short) + full = _full_history_frame() monkeypatch.setattr(climate, "_fetch_history_nasa", lambda c: full) + def _om_should_not_run(c): raise AssertionError("Open-Meteo should not be called") + monkeypatch.setattr(climate, "_fetch_history", _om_should_not_run) df, meta = climate._load_history(cell) - assert meta["source"] == "nasa-power" # rejected the short archive + assert meta["source"] == "nasa-power" assert df.height == full.height -def test_full_archive_is_accepted(monkeypatch, tmp_path): - """A backfilled archive (enough days) is accepted and served as open-meteo.""" +def test_falls_back_to_open_meteo_when_nasa_unavailable(monkeypatch, tmp_path): + """When NASA POWER fails, the load falls back to the Open-Meteo archive.""" monkeypatch.setattr(climate, "CACHE_DIR", str(tmp_path)) monkeypatch.setattr(climate, "_archive_cooldown_until", 0.0) - cell = {"id": "full_cell", "center_lat": 47.6, "center_lon": -122.3} + cell = {"id": "om_fallback", "center_lat": 47.6, "center_lon": -122.3} full = _full_history_frame() + def _nasa_down(c): raise RuntimeError("NASA POWER outage") + monkeypatch.setattr(climate, "_fetch_history_nasa", _nasa_down) monkeypatch.setattr(climate, "_fetch_history", lambda c: full) - def _nasa_should_not_run(c): raise AssertionError("NASA should not be called") - monkeypatch.setattr(climate, "_fetch_history_nasa", _nasa_should_not_run) df, meta = climate._load_history(cell) assert meta["source"] == "open-meteo" assert df.height == full.height + + +def test_short_nasa_is_rejected_and_falls_back_to_open_meteo(monkeypatch, tmp_path): + """A too-short NASA response (a partial/unavailable point) is rejected rather than + cached as a complete record, and Open-Meteo serves instead.""" + monkeypatch.setattr(climate, "CACHE_DIR", str(tmp_path)) + monkeypatch.setattr(climate, "_archive_cooldown_until", 0.0) + cell = {"id": "short_nasa", "center_lat": 47.6, "center_lon": -122.3} + + short = climate._to_frame(_om_daily(3)) # 2 usable days « threshold + full = _full_history_frame() # >= MIN_ARCHIVE_DAYS + assert short.height < climate.MIN_ARCHIVE_DAYS <= full.height + monkeypatch.setattr(climate, "_fetch_history_nasa", lambda c: short) + monkeypatch.setattr(climate, "_fetch_history", lambda c: full) + + df, meta = climate._load_history(cell) + assert meta["source"] == "open-meteo" # rejected the short NASA response + assert df.height == full.height + + +def test_history_tail_prefers_nasa(monkeypatch): + """The tail top-up fetches NASA POWER first; Open-Meteo's range fetch is not + touched when NASA succeeds.""" + cell = {"center_lat": 1.0, "center_lon": 2.0} + sentinel = object() + monkeypatch.setattr(climate, "_fetch_history_nasa", + lambda c, start=None, end=None: sentinel) + def _om_should_not_run(c, s, e): raise AssertionError("Open-Meteo tail should not run") + monkeypatch.setattr(climate, "_fetch_history_range", _om_should_not_run) + + assert climate._fetch_history_tail(cell, "2026-06-01", "2026-06-10") is sentinel + + +def test_history_tail_falls_back_to_open_meteo(monkeypatch): + """When NASA fails and Open-Meteo isn't cooling down, the tail falls back to the + Open-Meteo archive range.""" + monkeypatch.setattr(climate, "_archive_cooldown_until", 0.0) + cell = {"center_lat": 1.0, "center_lon": 2.0} + def _nasa_down(c, start=None, end=None): raise RuntimeError("NASA down") + monkeypatch.setattr(climate, "_fetch_history_nasa", _nasa_down) + sentinel = object() + monkeypatch.setattr(climate, "_fetch_history_range", lambda c, s, e: sentinel) + + assert climate._fetch_history_tail(cell, "2026-06-01", "2026-06-10") is sentinel -- 2.45.2 From 24d8cd9ba1e11516b31fa2a8944c54e0bb881fe7 Mon Sep 17 00:00:00 2001 From: Emi Griffith Date: Wed, 22 Jul 2026 22:07:00 -0700 Subject: [PATCH 4/6] Add one-time ERA5 seed for curated-city cells (keyless Icechunk) seed_era5.py backfills the curated-city cells with true ERA5 data from the Earthmover Icechunk ERA5 archive on AWS Open Data (anonymous, keyless), so high-traffic pages keep ERA5 fidelity now that NASA POWER is the live primary. It aggregates the hourly ERA5 point series to the daily store schema in polars (no pandas), derives feels-like via the existing heat-index/wind-chill path, and writes straight to the history store via climate._write_history_backed, bypassing the live fetch. ERA5 carries real gusts (i10fg), so seeded cells keep measured gusts rather than the Meteostat fallback. The icechunk/xarray/zarr stack is isolated in requirements-seed.txt (not the app image or CI) and lazy-imported, so the module and its unit-tested hourly->daily transform load without those deps. The S3 access layer targets a young API and is verified via `--dry-run` on the seed box, not in tests. Idempotent by default (skips already-cached cells; --overwrite to reseed). --- requirements-seed.txt | 9 +++ seed_era5.py | 174 ++++++++++++++++++++++++++++++++++++++++ tests/test_seed_era5.py | 72 +++++++++++++++++ 3 files changed, 255 insertions(+) create mode 100644 requirements-seed.txt create mode 100644 seed_era5.py create mode 100644 tests/test_seed_era5.py diff --git a/requirements-seed.txt b/requirements-seed.txt new file mode 100644 index 0000000..5791fa4 --- /dev/null +++ b/requirements-seed.txt @@ -0,0 +1,9 @@ +# Seed-only dependencies for seed_era5.py (the one-time ERA5 backfill). NOT installed +# in the app image or CI — kept out of requirements.txt so the runtime stays lean. +# Install on the machine that runs the seed: pip install -r requirements-seed.txt +# Versions intentionally unpinned: the icechunk/xarray/zarr stack moves fast; install +# the current compatible set on the seed box and verify with `python seed_era5.py --dry-run`. +-r requirements.txt +icechunk +xarray +zarr diff --git a/seed_era5.py b/seed_era5.py new file mode 100644 index 0000000..03be936 --- /dev/null +++ b/seed_era5.py @@ -0,0 +1,174 @@ +"""One-time ERA5 seed for the curated-city cells (keeps ERA5 fidelity without Open-Meteo). + +Now that the live history primary is NASA POWER, this backfills the high-traffic +curated-city cells (backend/cities.json) with true ERA5 data pulled from a public, +KEYLESS source — the Earthmover Icechunk ERA5 archive on AWS Open Data — so the +pages people actually visit keep ERA5-quality history. It writes straight into the +history store (bypassing the live fetch), so a seeded cell is served from cache and +never falls to NASA on first request. + + # on a machine with network + the seed deps (see requirements-seed.txt): + pip install -r requirements-seed.txt + python seed_era5.py --dry-run # fetch+transform one cell, print, no write + python seed_era5.py [--limit N] [--overwrite] + +Idempotent by default: a cell that already has a cached history is skipped (pass +--overwrite to reseed). Not run in CI or the app image — the icechunk/xarray/zarr +stack lives only in requirements-seed.txt. + +IMPORTANT — the S3/Icechunk *access* layer (_open_store / fetch_point) targets a +young, fast-moving API and the store's exact variable/coordinate names; it is NOT +exercised by the test suite. VERIFY it with `--dry-run` on the seed box before a +full run. The hourly→daily transform (hourly_to_daily) IS unit-tested and is the +part that must be numerically correct. + +Source: https://registry.opendata.aws/earthmover-era5/ (group="temporal", anonymous) +""" +import os +import sys + +import polars as pl + +from data import cities +from data import climate +from data import grid + +# Match the live archive span (climate.START_DATE / ARCHIVE_LATENCY_DAYS). +START_DATE = climate.START_DATE # "1980-01-01" + +# Anonymous AWS Open Data Icechunk store; overridable in case the registry moves. +ERA5_BUCKET = os.environ.get("THERMOGRAPH_ERA5_BUCKET", "earthmover-icechunk-era5") +ERA5_PREFIX = os.environ.get("THERMOGRAPH_ERA5_PREFIX", "") # verify on the seed box +ERA5_REGION = os.environ.get("THERMOGRAPH_ERA5_REGION", "us-east-1") +ERA5_GROUP = os.environ.get("THERMOGRAPH_ERA5_GROUP", "temporal") # time-contiguous layout + +# ERA5 store variable name -> our short name. ERA5 gusts (i10fg) are real here, so the +# seed keeps genuine measured gusts rather than the Meteostat/estimate fallback. +ERA5_VARS = { + "2m_temperature": "t2m", # K + "2m_dewpoint_temperature": "d2m", # K (→ RH with t2m) + "total_precipitation": "tp", # m, hourly accumulation + "10m_u_component_of_wind": "u10", # m/s + "10m_v_component_of_wind": "v10", # m/s + "instantaneous_10m_wind_gust": "i10fg", # m/s +} + +_M_TO_IN = 39.37007874 +_MS_TO_MPH = 2.2369362920544 + + +def hourly_to_daily(hourly: pl.DataFrame) -> pl.DataFrame: + """Aggregate an hourly ERA5 point frame to the daily schema climate.py stores, + in the app's units (°F, inches, mph, %RH). Input columns: time (Datetime, UTC), + t2m/d2m (K), tp (m, per-hour accumulation), u10/v10/i10fg (m/s). + + Days are UTC calendar days — the reanalysis is UTC, so this is a small offset + from a local-timezone daily boundary, immaterial for percentile climatology. + Relative humidity is the Magnus-formula RH from temperature and dewpoint.""" + tc = pl.col("t2m") - 273.15 + tdc = pl.col("d2m") - 273.15 + es = (17.625 * tc / (243.04 + tc)).exp() + e = (17.625 * tdc / (243.04 + tdc)).exp() + df = hourly.with_columns( + pl.col("time").dt.date().alias("date"), + ((pl.col("t2m") - 273.15) * 9.0 / 5.0 + 32.0).alias("t_f"), + ((pl.col("u10") ** 2 + pl.col("v10") ** 2).sqrt() * _MS_TO_MPH).alias("wind_mph"), + (pl.col("i10fg") * _MS_TO_MPH).alias("gust_mph"), + (pl.col("tp") * _M_TO_IN).alias("precip_in"), + (100.0 * e / es).clip(0.0, 100.0).alias("rh"), + ) + return ( + df.group_by("date") + .agg( + pl.col("t_f").max().alias("tmax"), + pl.col("t_f").min().alias("tmin"), + pl.col("precip_in").sum().alias("precip"), + pl.col("wind_mph").max().alias("wind"), + pl.col("gust_mph").max().alias("gust"), + pl.col("rh").mean().alias("humid"), + ) + .sort("date") + ) + + +def _default_end() -> str: + import datetime + return (datetime.date.today() + - datetime.timedelta(days=climate.ARCHIVE_LATENCY_DAYS)).isoformat() + + +def _open_store(): + """Open the anonymous Icechunk ERA5 store (temporal group). Lazy-imports the + seed-only deps. VERIFY against the current icechunk/xarray API on the seed box — + this young API and the store's exact paths are not covered by tests.""" + import icechunk # noqa: PLC0415 - seed-only dep (requirements-seed.txt) + import xarray as xr # noqa: PLC0415 + + storage = icechunk.s3_storage( + bucket=ERA5_BUCKET, prefix=ERA5_PREFIX, region=ERA5_REGION, anonymous=True) + repo = icechunk.Repository.open(storage) + session = repo.readonly_session("main") + return xr.open_zarr(session.store, group=ERA5_GROUP, consolidated=False) + + +def fetch_point(ds, lat: float, lon: float, start: str, end: str) -> pl.DataFrame: + """Pull the hourly ERA5 series for the nearest grid point over [start, end] into a + polars frame (no pandas). ERA5 longitudes are 0..360. Coordinate/variable names + may differ per store revision — verify with --dry-run.""" + sub = (ds[list(ERA5_VARS)] + .sel(latitude=lat, longitude=lon % 360.0, method="nearest") + .sel(time=slice(start, end)) + .load()) + data = {"time": sub["time"].values} + for src, dst in ERA5_VARS.items(): + data[dst] = sub[src].values + return pl.DataFrame(data) + + +def seed_cell(ds, cell: dict, start: str, end: str) -> int: + """Fetch → daily → finalize → write one cell's ERA5 history. Returns row count.""" + daily = hourly_to_daily(fetch_point(ds, cell["center_lat"], cell["center_lon"], start, end)) + # ERA5 has no apparent-temperature field; approximate fmax/fmin from the NWS heat + # index / wind chill (same as the NASA/MET paths). _finalize_approximated adds + # feels + doy and folds NaN→null; the write path strips doy. + frame = climate._finalize_approximated(daily) + climate._write_history_backed(cell["id"], frame) + return frame.height + + +def main(limit: "int | None" = None, overwrite: bool = False, + dry_run: bool = False, start: str = START_DATE, end: "str | None" = None) -> None: + end = end or _default_end() + ds = _open_store() + todo = cities.all_cities() + if limit: + todo = todo[:limit] + seeded = skipped = failed = 0 + for i, c in enumerate(todo, 1): + cell = grid.snap(c["lat"], c["lon"]) + if not overwrite: + cached = climate.load_cached_history(cell) + if cached is not None and not cached.is_empty(): + skipped += 1 + continue + try: + if dry_run: + daily = hourly_to_daily( + fetch_point(ds, cell["center_lat"], cell["center_lon"], start, end)) + print(daily.head()) + print(f"[dry-run] {c['slug']} ({cell['id']}): {daily.height} days " + f"({start}..{end}); nothing written") + return + n = seed_cell(ds, cell, start, end) + seeded += 1 + print(f"[{i}/{len(todo)}] seeded {c['slug']} ({cell['id']}): {n} days") + except Exception as e: # noqa: BLE001 - keep going; a failed cell self-heals via NASA + failed += 1 + print(f"[{i}/{len(todo)}] FAILED {c['slug']}: {e}") + print(f"done: seeded={seeded} skipped(cached)={skipped} failed={failed}") + + +if __name__ == "__main__": + argv = sys.argv[1:] + lim = int(argv[argv.index("--limit") + 1]) if "--limit" in argv else None + main(limit=lim, overwrite="--overwrite" in argv, dry_run="--dry-run" in argv) diff --git a/tests/test_seed_era5.py b/tests/test_seed_era5.py new file mode 100644 index 0000000..b323e58 --- /dev/null +++ b/tests/test_seed_era5.py @@ -0,0 +1,72 @@ +"""Unit tests for the ERA5 seed's hourly→daily transform (hermetic — no network, no +icechunk/xarray). The S3 access layer (_open_store/fetch_point) is deliberately not +tested here; it must be verified with `seed_era5.py --dry-run` on the seed box.""" +import datetime +import math + +import polars as pl +import pytest + +import seed_era5 + +_MS_TO_MPH = 2.2369362920544 +_M_TO_IN = 39.37007874 + + +def _rh(t_c, td_c): + es = math.exp(17.625 * t_c / (243.04 + t_c)) + e = math.exp(17.625 * td_c / (243.04 + td_c)) + return min(100.0, 100.0 * e / es) + + +def _hourly(rows): + """rows: list of (time, t2m_K, d2m_K, tp_m, u10, v10, i10fg).""" + cols = ["time", "t2m", "d2m", "tp", "u10", "v10", "i10fg"] + return pl.DataFrame({c: [r[i] for r in rows] for i, c in enumerate(cols)}) + + +def test_hourly_to_daily_aggregates_and_converts(): + d1 = datetime.datetime(2020, 1, 1) + d2 = datetime.datetime(2020, 1, 2) + hourly = _hourly([ + # day 1: 10°C sat, then 20°C / 10°C dew; winds 5 m/s then calm; gusts 10/12 + (d1.replace(hour=0), 283.15, 283.15, 0.001, 3.0, 4.0, 10.0), + (d1.replace(hour=1), 293.15, 283.15, 0.002, 0.0, 0.0, 12.0), + # day 2: 0°C sat; wind 10 m/s; gust 20 + (d2.replace(hour=0), 273.15, 273.15, 0.0, 6.0, 8.0, 20.0), + ]) + daily = seed_era5.hourly_to_daily(hourly).sort("date") + assert daily["date"].to_list() == [datetime.date(2020, 1, 1), datetime.date(2020, 1, 2)] + + r0 = daily.row(0, named=True) + assert r0["tmax"] == pytest.approx(68.0) # 20°C + assert r0["tmin"] == pytest.approx(50.0) # 10°C + assert r0["precip"] == pytest.approx(0.003 * _M_TO_IN) + assert r0["wind"] == pytest.approx(5.0 * _MS_TO_MPH) # max(5, 0) + assert r0["gust"] == pytest.approx(12.0 * _MS_TO_MPH) # max(10, 12) + assert r0["humid"] == pytest.approx((_rh(10, 10) + _rh(20, 10)) / 2) + + r1 = daily.row(1, named=True) + assert r1["tmin"] == pytest.approx(32.0) # 0°C + assert r1["wind"] == pytest.approx(10.0 * _MS_TO_MPH) # sqrt(6²+8²) + + +def test_saturated_air_reads_full_humidity(): + hourly = _hourly([(datetime.datetime(2021, 6, 1), 293.15, 293.15, 0.0, 1.0, 1.0, 2.0)]) + daily = seed_era5.hourly_to_daily(hourly) + assert daily.row(0, named=True)["humid"] == pytest.approx(100.0) + + +def test_finalize_produces_the_store_schema(): + """The daily frame, run through climate._finalize_approximated (as seed_cell does), + yields the persisted columns incl. derived feels.""" + from data import climate + hourly = _hourly([ + (datetime.datetime(2019, 7, 1, 0), 300.0, 290.0, 0.0, 2.0, 2.0, 5.0), + (datetime.datetime(2019, 7, 1, 1), 305.0, 291.0, 0.001, 3.0, 1.0, 7.0), + ]) + frame = climate._finalize_approximated(seed_era5.hourly_to_daily(hourly)) + for col in ("date", "tmax", "tmin", "precip", "wind", "gust", "humid", + "fmax", "fmin", "feels"): + assert col in frame.columns + assert frame.height == 1 -- 2.45.2 From 5569bcbc0a974ab799a04865e9499beeccdba785 Mon Sep 17 00:00:00 2001 From: Emi Griffith Date: Thu, 23 Jul 2026 06:32:31 -0700 Subject: [PATCH 5/6] Flip recent/forecast leg off Open-Meteo (NASA past + MET Norway forward) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The recent-observations + forward-forecast bundle is now built without Open-Meteo: the recent observed window comes from a NASA POWER range (measured, via the range fetch added for the history flip) and the forward days from MET Norway, merged by date. Meteostat fills the gusts neither source carries — measured for the observed days, estimated from wind for the forecast days. Open-Meteo's forecast API is demoted to the fallback, still gated by its existing _forecast_cooldown_until rate-limit cooldown; then a stale cache. MET Norway switched to /compact (same fields, smaller payload) and the bundle TTL relaxed 1h -> 4h. NASA's near-real-time lag leaves a 1-2 day recent-edge gap that grades as missing, bracketed by history behind and forecast ahead. Tests updated to the new source order (NASA+MET merge, Open-Meteo fallback, and the forecast cooldown now gating that fallback); deletion of Open-Meteo is not done — it stays as the dormant fallback. --- backend/data/climate.py | 133 +++++++++++++++++------------ backend/tests/data/test_climate.py | 88 +++++++++++-------- 2 files changed, 129 insertions(+), 92 deletions(-) diff --git a/backend/data/climate.py b/backend/data/climate.py index 909ad24..fc5e411 100644 --- a/backend/data/climate.py +++ b/backend/data/climate.py @@ -60,7 +60,8 @@ MIN_ARCHIVE_DAYS = 3650 # ~10 yrs: far above any sync window, far bel # cached indefinitely (refetched only to add new metric columns). Only the recent # tail is refreshed — a small incremental fetch, at most hourly. HISTORY_TOPUP_INTERVAL = 3600 # seconds between recent-tail refresh attempts -FORECAST_TTL_HOURS = 1 # refetch the forward forecast hourly to track updates +FORECAST_TTL_HOURS = 4 # refetch the recent+forecast bundle every ~4h (lower IO; + # daily-resolution forecast normals move slowly) # The archive (historical) endpoint. Overridable so it can point at a self-hosted # Open-Meteo instance (ERA5 in object storage) instead of the rate-limited public @@ -80,12 +81,13 @@ NASA_POWER_URL = "https://power.larc.nasa.gov/api/temporal/daily/point" NASA_START = "19810101" NASA_FILL = -900.0 # POWER's missing-value sentinel is ~-999 -# Backup forward forecast when Open-Meteo's forecast API is unavailable. MET Norway -# (yr.no) is free + keyless + global, mirroring NASA POWER's role for history. It -# returns a sub-daily timeseries in metric units with no gusts or apparent temp, so -# we aggregate to daily, convert units, and derive feels-like like the NASA path. -# It is forecast-only — no recent past days — so it's a degraded-but-working fallback. -METNO_URL = "https://api.met.no/weatherapi/locationforecast/2.0/complete" +# Forward-forecast source (primary for the forward days). MET Norway (yr.no) is free +# + keyless + global. It returns a sub-daily timeseries in metric units with no gusts +# or apparent temp, so we aggregate to daily, convert units, derive feels-like like +# the NASA path, and fill gusts from Meteostat. It is forecast-only (no recent past), +# so the recent observed days come from a NASA POWER range instead. The /compact +# endpoint carries every field we use at a smaller payload than /complete. +METNO_URL = "https://api.met.no/weatherapi/locationforecast/2.0/compact" # MET Norway's ToS requires an identifying User-Agent (a missing/generic one is # 403'd); include the app and a contact URL so they can reach us if usage misbehaves. METNO_UA = "Thermograph/0.2 (+https://thermograph.org)" @@ -824,6 +826,10 @@ def _load_history(cell: dict) -> tuple[pl.DataFrame, dict]: RECENT_PAST_DAYS = 25 # recent observations window (covers the ~2-week graded view) FORECAST_DAYS = 8 # today + 7 days ahead +# NASA POWER near-real-time lags a couple of days, so the recent observed window ends +# here; the small today-1/today-2 gap before the MET Norway forecast begins grades as +# missing (drop_nulls), bracketed by history behind and forecast ahead. +RECENT_END_LAG_DAYS = 2 def _rf_cache_path(cell_id: str) -> str: @@ -865,20 +871,27 @@ def load_cached_recent_forecast(cell: dict) -> "pl.DataFrame | None": return None -def _load_recent_forecast(cell: dict) -> pl.DataFrame: - """Recent observations AND the forward forecast in ONE forecast-API call. +def _fetch_recent_forecast(cell: dict) -> pl.DataFrame: + """Recent observations + forward forecast WITHOUT Open-Meteo: the recent observed + window from a NASA POWER range (measured), the forward days from MET Norway, + merged by date. Gusts — which neither source carries — are filled from Meteostat: + measured for the observed days, estimated from wind for the forecast days.""" + today = datetime.date.today() + start = (today - datetime.timedelta(days=RECENT_PAST_DAYS)).isoformat() + end = (today - datetime.timedelta(days=RECENT_END_LAG_DAYS)).isoformat() + past = _fetch_history_nasa(cell, start=start, end=end) # measured + Meteostat gusts + fwd = meteostat.fill_gusts( + cell["center_lat"], cell["center_lon"], _fetch_forecast_metno(cell)) + # Forecast wins on any overlapping day (freshest model value for today); the NASA + # observed days fill the recent past MET Norway lacks. + return (pl.concat([past, fwd], how="diagonal_relaxed") + .unique(subset="date", keep="last", maintain_order=True) + .sort("date")) - Both the recent (past) view and the forecast (future) view slice from this, so - a cell needs just one upstream forecast request per hour (plus the ~monthly - archive fetch). Cached per cell for one hour so it still follows model updates. - """ - cell_id = cell["id"] - hit = _read_recent_backed(cell_id) - if hit is not None: - cached, age_s = hit - if age_s / 3600.0 < FORECAST_TTL_HOURS: - return _with_doy(cached) +def _fetch_recent_forecast_om(cell: dict) -> pl.DataFrame: + """Fallback recent+forecast bundle from the Open-Meteo forecast API (the former + primary): recent past + forward days in one call.""" params = { "latitude": cell["center_lat"], "longitude": cell["center_lon"], @@ -890,49 +903,59 @@ def _load_recent_forecast(cell: dict) -> pl.DataFrame: "past_days": RECENT_PAST_DAYS, "forecast_days": FORECAST_DAYS, } - # Skip Open-Meteo's forecast endpoint entirely while it's in its own - # rate-limit cooldown (see _note_forecast_rate_limit) -- mirrors - # _load_history's archive-cooldown check, so a forecast brownout doesn't - # retry-storm the endpoint from every subscribed cell and every live request - # independently; it falls straight through to the MET Norway backup instead. + r = _request(FORECAST_URL, params, 60, phase="recent_forecast_fetch") + return _to_frame(r.json()["daily"]) + + +def _load_recent_forecast(cell: dict) -> pl.DataFrame: + """Recent observations AND the forward forecast for a cell. + + The primary source is keyless and Open-Meteo-free: a NASA POWER recent-past range + plus the MET Norway forward forecast (see _fetch_recent_forecast). The Open-Meteo + forecast API is the fallback (skipped while it's in its own rate-limit cooldown, + see _note_forecast_rate_limit), then a stale cache. Cached per cell for + FORECAST_TTL_HOURS so it follows model updates without per-hour upstream IO. + """ + cell_id = cell["id"] + hit = _read_recent_backed(cell_id) + if hit is not None: + cached, age_s = hit + if age_s / 3600.0 < FORECAST_TTL_HOURS: + return _with_doy(cached) + df = None - primary_error = None - if time.time() >= _forecast_cooldown_until: + try: + df = _fetch_recent_forecast(cell) # NASA recent + MET forward (primary) + except Exception: # noqa: BLE001 - keyless primary unavailable; try Open-Meteo + pass + + forecast_error = None + if df is None and time.time() >= _forecast_cooldown_until: + # Fall back to the Open-Meteo forecast API, skipped while it's in its own + # rate-limit cooldown so a brownout doesn't retry-storm the endpoint. try: - r = _request(FORECAST_URL, params, 60, phase="recent_forecast_fetch") - df = _to_frame(r.json()["daily"]) + df = _fetch_recent_forecast_om(cell) except Exception as e: # noqa: BLE001 if is_rate_limit(e): _note_forecast_rate_limit(e) - primary_error = e + forecast_error = e if df is None: - # Open-Meteo forecast unavailable (rate limit, cooldown, or outage). Fall - # back to MET Norway (yr.no) — global + keyless — mirroring the archive's - # NASA POWER backup. MET Norway is forecast-only (no recent past days), - # so this keeps the forecast / day-ahead views working in a degraded form. - try: - df = _fetch_forecast_metno(cell) - except Exception: # noqa: BLE001 - backup unavailable too - # Last resort: serve the stale cache if we have one. An hours-old bundle - # (which still carries the recent observed days MET Norway lacks) beats a - # hard failure, mirroring _load_history's stale-serve. Return WITHOUT - # rewriting it, so recent_stamp stays old and the next request still - # retries upstream first rather than serving this as if it were fresh. - if hit is not None: - return _with_doy(hit[0]) - # No cache either: surface the original error, classifying an Open-Meteo - # rate limit as the typed, daily-aware WeatherUnavailable (the archive - # path classifies its own inside _load_history). primary_error is None - # when the primary fetch was skipped outright (cooldown active) -- - # report the cooldown itself in that case. - if primary_error is not None and is_rate_limit(primary_error): - daily = "daily" in _rate_limit_reason(primary_error).lower() - raise WeatherUnavailable(limit_message(daily), daily=daily) from primary_error - if primary_error is not None: - raise primary_error - raise WeatherUnavailable(limit_message(_forecast_limit_daily), - daily=_forecast_limit_daily) + # Last resort: serve the stale cache if we have one, WITHOUT rewriting it, so + # recent_stamp stays old and the next request retries upstream first (mirrors + # _load_history's stale-serve). + if hit is not None: + return _with_doy(hit[0]) + # No cache either: surface the error, classifying an Open-Meteo rate limit as + # the typed, daily-aware WeatherUnavailable. forecast_error is None when the + # fallback was skipped outright (cooldown active) — report the cooldown then. + if forecast_error is not None and is_rate_limit(forecast_error): + daily = "daily" in _rate_limit_reason(forecast_error).lower() + raise WeatherUnavailable(limit_message(daily), daily=daily) from forecast_error + if forecast_error is not None: + raise forecast_error + raise WeatherUnavailable(limit_message(_forecast_limit_daily), + daily=_forecast_limit_daily) _write_recent_backed(cell_id, df) return df diff --git a/backend/tests/data/test_climate.py b/backend/tests/data/test_climate.py index 6fab033..79c3d71 100644 --- a/backend/tests/data/test_climate.py +++ b/backend/tests/data/test_climate.py @@ -137,34 +137,53 @@ def test_metno_to_frame_aggregates_daily_and_converts_units(): assert round(d2["precip"], 4) == round(6.0 / 25.4, 4) # only the 6-hour block -def test_recent_forecast_falls_back_to_metno(monkeypatch, tmp_path): - """When Open-Meteo's forecast API fails, the MET Norway backup serves the frame.""" +def test_recent_forecast_merges_nasa_past_and_metno_forward(monkeypatch, tmp_path): + """The recent/forecast bundle is built from a NASA POWER recent-past range plus the + MET Norway forward forecast, merged by date (Open-Meteo not consulted).""" monkeypatch.setattr(climate, "CACHE_DIR", str(tmp_path)) - cell = {"id": "fallback_cell", "center_lat": 47.6062, "center_lon": -122.3321} - met = {"timeseries": [ - _metno_step("2026-07-16T00:00:00Z", 18.0, 55.0, 3.0, p1=0.0), - _metno_step("2026-07-17T00:00:00Z", 22.0, 45.0, 5.0, p6=1.0), - ]} + cell = {"id": "rf_primary", "center_lat": 47.6, "center_lon": -122.3} - class Resp: - def json(self): return {"properties": met} + past = climate._finalize_approximated(pl.DataFrame({ + "date": [datetime.date(2026, 7, 10), datetime.date(2026, 7, 11)], + "tmax": [70.0, 72.0], "tmin": [50.0, 51.0], "precip": [0.0, 0.1], + "wind": [5.0, 6.0], "humid": [60.0, 55.0], + })) + fwd = climate._finalize_approximated(pl.DataFrame({ + "date": [datetime.date(2026, 7, 16), datetime.date(2026, 7, 17)], + "tmax": [80.0, 82.0], "tmin": [60.0, 61.0], "precip": [0.0, 0.0], + "wind": [3.0, 4.0], "humid": [40.0, 45.0], + })) + monkeypatch.setattr(climate, "_fetch_history_nasa", lambda c, start, end: past) + monkeypatch.setattr(climate, "_fetch_forecast_metno", lambda c: fwd) + monkeypatch.setattr(climate.meteostat, "fill_gusts", lambda lat, lon, df: df) + def _om_should_not_run(c): raise AssertionError("Open-Meteo forecast should not run") + monkeypatch.setattr(climate, "_fetch_recent_forecast_om", _om_should_not_run) - def fake_request(url, params, timeout, *, phase, headers=None, attempts=climate.MAX_ATTEMPTS): - if url == climate.FORECAST_URL: - raise RuntimeError("open-meteo forecast outage") - assert url == climate.METNO_URL - assert headers and headers.get("User-Agent"), "MET Norway needs a User-Agent" - return Resp() - - monkeypatch.setattr(climate, "_request", fake_request) df = climate._load_recent_forecast(cell) - assert len(df) == 2 and df["date"].max() == datetime.date(2026, 7, 17) + assert df["date"].to_list() == [ + datetime.date(2026, 7, 10), datetime.date(2026, 7, 11), + datetime.date(2026, 7, 16), datetime.date(2026, 7, 17)] + + +def test_recent_forecast_falls_back_to_open_meteo(monkeypatch, tmp_path): + """When the NASA + MET Norway primary fails, the Open-Meteo forecast API serves.""" + monkeypatch.setattr(climate, "CACHE_DIR", str(tmp_path)) + monkeypatch.setattr(climate, "_forecast_cooldown_until", 0.0) + cell = {"id": "rf_fallback", "center_lat": 47.6, "center_lon": -122.3} + + def _primary_down(c): raise RuntimeError("NASA + MET Norway down") + monkeypatch.setattr(climate, "_fetch_recent_forecast", _primary_down) + om = climate._to_frame(_om_daily(3)) + monkeypatch.setattr(climate, "_fetch_recent_forecast_om", lambda c: om) + + df = climate._load_recent_forecast(cell) + assert df.height == om.height def test_recent_forecast_serves_stale_cache_when_all_sources_fail(monkeypatch, tmp_path): - """With Open-Meteo AND MET Norway both down, an existing (stale) cache is served - rather than failing — and it is NOT rewritten, so its mtime stays old and the - next request still retries upstream first.""" + """With every source down (the NASA + MET Norway primary and the Open-Meteo + fallback), an existing (stale) cache is served rather than failing — and it is NOT + rewritten, so its mtime stays old and the next request still retries upstream first.""" import os import time monkeypatch.setattr(climate, "CACHE_DIR", str(tmp_path)) @@ -370,28 +389,23 @@ def test_history_tail_falls_back_to_open_meteo(monkeypatch): # --- forecast cooldown (mirrors the archive path's, tracked separately) ----- -def test_forecast_cooldown_skips_open_meteo_and_falls_to_metno(monkeypatch, tmp_path): - """While the forecast cooldown is active, the forecast fetch must not call - Open-Meteo at all -- straight to the MET Norway backup, mirroring - _load_history's archive-cooldown skip.""" +def test_forecast_cooldown_skips_the_open_meteo_fallback(monkeypatch, tmp_path): + """While the forecast cooldown is active, the Open-Meteo FALLBACK is skipped: with + the keyless primary (NASA + MET Norway) also down and no cache, the load surfaces + WeatherUnavailable without ever calling Open-Meteo.""" import time as time_mod monkeypatch.setattr(climate, "CACHE_DIR", str(tmp_path)) monkeypatch.setattr(climate, "_forecast_cooldown_until", time_mod.time() + 60) cell = {"id": "fc_cooldown_cell", "center_lat": 47.6, "center_lon": -122.3} - met = {"timeseries": [_metno_step("2026-07-16T00:00:00Z", 18.0, 55.0, 3.0, p1=0.0)]} + def _primary_down(c): raise RuntimeError("NASA + MET Norway down") + monkeypatch.setattr(climate, "_fetch_recent_forecast", _primary_down) + def _om_should_not_run(c): + raise AssertionError("Open-Meteo fallback must not run during cooldown") + monkeypatch.setattr(climate, "_fetch_recent_forecast_om", _om_should_not_run) - class Resp: - def json(self): return {"properties": met} - - def fake_request(url, params, timeout, *, phase, headers=None, attempts=climate.MAX_ATTEMPTS): - assert url != climate.FORECAST_URL, "must not call Open-Meteo during cooldown" - assert url == climate.METNO_URL - return Resp() - - monkeypatch.setattr(climate, "_request", fake_request) - df = climate._load_recent_forecast(cell) - assert len(df) == 1 + with pytest.raises(climate.WeatherUnavailable): + climate._load_recent_forecast(cell) def test_forecast_429_sets_its_own_cooldown(monkeypatch, tmp_path): -- 2.45.2 From abf37e6376c75e2bdf321319395dc946d0ffa9f5 Mon Sep 17 00:00:00 2001 From: Emi Griffith Date: Thu, 23 Jul 2026 06:33:54 -0700 Subject: [PATCH 6/6] Add value-drift check comparing NASA vs Open-Meteo (ERA5) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit backend/drift_check.py quantifies the migration's data impact by fetching both live history sources (NASA POWER primary + the Open-Meteo/ERA5 fallback) for a sample of cells and reporting, per variable: mean-absolute-difference / bias / max (raw drift) and grade-band divergence (the share of days whose graded band actually changes). Since Open-Meteo is ERA5, this doubles as a fidelity check for the ERA5 seed. Because grades are percentiles within each source's own distribution, a uniform bias moves MAD but not grades — so the two numbers are read together. Open-Meteo is kept as a dormant fallback (not deleted), so both sources stay fetchable for this comparison. The comparison logic is unit-tested; the per-cell fetch runs on a networked box (python drift_check.py [--limit N] [--days N]). --- backend/drift_check.py | 136 ++++++++++++++++++++++++++++++ backend/tests/test_drift_check.py | 55 ++++++++++++ 2 files changed, 191 insertions(+) create mode 100644 backend/drift_check.py create mode 100644 backend/tests/test_drift_check.py diff --git a/backend/drift_check.py b/backend/drift_check.py new file mode 100644 index 0000000..b91df04 --- /dev/null +++ b/backend/drift_check.py @@ -0,0 +1,136 @@ +"""Value-drift check: quantify how much the migration moved the data, per variable +and by grade band, by comparing the two live history sources head-to-head. + +Open-Meteo is kept as a dormant fallback, so both sources are still fetchable: this +pulls NASA POWER (the new primary) and the Open-Meteo archive (ERA5) for a sample of +cells and reports, per variable: + - mean-absolute-difference / bias / max — the RAW value drift, and + - grade-band divergence — the share of (day, metric) whose graded band actually + changes between sources. +Since Open-Meteo IS ERA5, this doubles as a fidelity check for the ERA5 seed. + + python drift_check.py [--limit N] [--days N] # on a networked box + +Note a real property of the grading model worth reading the two numbers together: +grades are percentiles within each source's OWN distribution, so a uniform bias +between sources shifts MAD but leaves grades unchanged. High MAD with low grade +divergence means "different absolute values, same story" — usually fine. High grade +divergence is the signal that actually matters. + +The comparison logic (compare_values / grade_band_divergence) is unit-tested; the +fetch is not. +""" +import datetime +import sys + +import polars as pl + +from data import cities +from data import climate +from data import grading +from data import grid + +# Variables compared for raw drift (feels is derived; precip/wind/gust/humid optional). +_METRICS = ("tmax", "tmin", "precip", "wind", "gust", "humid", "feels") + + +def compare_values(a: pl.DataFrame, b: pl.DataFrame, metrics=_METRICS) -> dict: + """Per-variable raw-drift stats over the dates the two frames share: n aligned + days, mean-absolute-difference, mean bias (a - b), and max absolute difference. + Days where either source is null for a metric are dropped from that metric.""" + j = a.join(b, on="date", how="inner", suffix="_b") + out = {} + for m in metrics: + if m not in a.columns or m not in b.columns: + continue + d = j.select((pl.col(m) - pl.col(f"{m}_b")).alias("d")).drop_nulls() + if d.height == 0: + out[m] = {"n": 0, "mad": None, "bias": None, "max_abs": None} + continue + diff = d["d"] + out[m] = { + "n": d.height, + "mad": round(float(diff.abs().mean()), 3), + "bias": round(float(diff.mean()), 3), + "max_abs": round(float(diff.abs().max()), 3), + } + return out + + +def grade_band_divergence(a: pl.DataFrame, b: pl.DataFrame, start, end) -> dict: + """Per-metric share of days whose graded band differs between the two sources over + [start, end]. Each source's days are graded against its OWN ±7-day climatology + (grading.grade_range), so this measures whether the migration changes the story a + day tells, not just its raw value.""" + ga = {d["date"]: d for d in grading.grade_range(a, start, end)} + gb = {d["date"]: d for d in grading.grade_range(b, start, end)} + per: dict[str, list[int]] = {} + for date in ga.keys() & gb.keys(): + da, db = ga[date], gb[date] + for m in grading.CLIMO_METRICS: + va, vb = da.get(m), db.get(m) + if va and vb: + slot = per.setdefault(m, [0, 0]) # [compared, differ] + slot[0] += 1 + slot[1] += int(va["g"] != vb["g"]) + return {m: {"compared": c, "differ": d, "pct": round(100.0 * d / c, 1) if c else None} + for m, (c, d) in sorted(per.items())} + + +def check_cell(cell: dict, days: int) -> dict: + """Fetch both sources for a cell and compare them over the most recent `days`.""" + nasa = climate._fetch_history_nasa(cell) # new primary (+ Meteostat gusts) + om = climate._fetch_history(cell) # Open-Meteo archive (ERA5) + end = min(nasa["date"].max(), om["date"].max()) + start = end - datetime.timedelta(days=days) + return {"values": compare_values(nasa, om), + "grades": grade_band_divergence(nasa, om, start, end)} + + +def _mean(xs: list) -> "float | None": + xs = [x for x in xs if x is not None] + return round(sum(xs) / len(xs), 3) if xs else None + + +def main(limit: "int | None" = None, days: int = 365) -> None: + todo = cities.all_cities() + if limit: + todo = todo[:limit] + mad_by_metric: dict[str, list] = {} + grade_by_metric: dict[str, list] = {} + checked = failed = 0 + for i, c in enumerate(todo, 1): + cell = grid.snap(c["lat"], c["lon"]) + try: + rep = check_cell(cell, days) + except Exception as e: # noqa: BLE001 - one cell's outage shouldn't abort the sweep + failed += 1 + print(f"[{i}/{len(todo)}] FAILED {c['slug']}: {e}") + continue + checked += 1 + vals = ", ".join(f"{m} mad={s['mad']}" for m, s in rep["values"].items() + if s.get("mad") is not None) + grds = ", ".join(f"{m} {s['pct']}%" for m, s in rep["grades"].items()) + print(f"[{i}/{len(todo)}] {c['slug']} ({cell['id']})") + print(f" values: {vals}") + print(f" grade drift: {grds}") + for m, s in rep["values"].items(): + mad_by_metric.setdefault(m, []).append(s.get("mad")) + for m, s in rep["grades"].items(): + grade_by_metric.setdefault(m, []).append(s.get("pct")) + + print(f"\n=== summary over {checked} cells (failed={failed}) ===") + print("mean MAD (raw drift):") + for m in _METRICS: + if m in mad_by_metric: + print(f" {m:7s} {_mean(mad_by_metric[m])}") + print("mean grade-band divergence (share of days a grade changes):") + for m in sorted(grade_by_metric): + print(f" {m:7s} {_mean(grade_by_metric[m])}%") + + +if __name__ == "__main__": + argv = sys.argv[1:] + lim = int(argv[argv.index("--limit") + 1]) if "--limit" in argv else None + dys = int(argv[argv.index("--days") + 1]) if "--days" in argv else 365 + main(limit=lim, days=dys) diff --git a/backend/tests/test_drift_check.py b/backend/tests/test_drift_check.py new file mode 100644 index 0000000..698dc95 --- /dev/null +++ b/backend/tests/test_drift_check.py @@ -0,0 +1,55 @@ +"""Unit tests for the drift-check comparison logic (hermetic — no network). The +source fetch (check_cell) is not tested here; it runs on a networked box.""" +import datetime +import math + +import polars as pl +import pytest + +import drift_check +from data import climate + +_D = [datetime.date(2024, 1, d) for d in (1, 2, 3)] + + +def test_compare_values_stats(): + a = pl.DataFrame({"date": _D, "tmax": [70.0, 80.0, 90.0], "wind": [5.0, 6.0, 7.0]}) + b = pl.DataFrame({"date": _D, "tmax": [72.0, 78.0, 90.0], "wind": [5.0, 6.0, 7.0]}) + out = drift_check.compare_values(a, b, metrics=("tmax", "wind")) + assert out["tmax"] == {"n": 3, "mad": pytest.approx(4 / 3, abs=1e-3), + "bias": pytest.approx(0.0), "max_abs": pytest.approx(2.0)} + assert out["wind"] == {"n": 3, "mad": 0.0, "bias": 0.0, "max_abs": 0.0} + + +def test_compare_values_skips_nulls_and_missing_metric(): + a = pl.DataFrame({"date": _D, "tmax": [70.0, 80.0, 90.0]}) + b = pl.DataFrame({"date": _D, "tmax": [71.0, None, 90.0]}) + out = drift_check.compare_values(a, b, metrics=("tmax", "gust")) + assert out["tmax"]["n"] == 2 # the null day is dropped + assert out["tmax"]["mad"] == pytest.approx(0.5) # |70-71| and |90-90| + assert "gust" not in out # metric absent from both frames + + +def _gradeable(n=800): + """A finalized, gradeable multi-year frame (seasonal temp swing so windows have a + real distribution).""" + start = datetime.date(2019, 1, 1) + dates = [start + datetime.timedelta(days=i) for i in range(n)] + tmax = [70.0 + 15.0 * math.sin(i / 58.0) for i in range(n)] + tmin = [t - 15.0 for t in tmax] + return climate._finalize_frame(pl.DataFrame({ + "date": dates, "tmax": tmax, "tmin": tmin, "precip": [0.0] * n, + "wind": [5.0] * n, "gust": [8.0] * n, "humid": [60.0] * n, + "fmax": tmax, "fmin": tmin, + })) + + +def test_grade_band_divergence_identical_is_zero(): + f = _gradeable() + end = f["date"].max() + start = end - datetime.timedelta(days=60) + div = drift_check.grade_band_divergence(f, f, start, end) + assert div # some metrics were compared + assert any(s["compared"] > 0 for s in div.values()) + assert all(s["differ"] == 0 for s in div.values()) # identical frames never diverge + assert all(s["pct"] == 0.0 for s in div.values()) -- 2.45.2