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)

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