Subtree-merge thermograph-backend era5-seed into backend/

backend/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

backend/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)

backend/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