From db862771037840f31676e3eb37bc9e739674d780 Mon Sep 17 00:00:00 2001 From: Emi Griffith Date: Wed, 22 Jul 2026 21:29:26 -0700 Subject: [PATCH 1/4] 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/4] 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/4] 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/4] 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