From aec64b40582f4a3ef4b6ee0f970e71614f7dd853 Mon Sep 17 00:00:00 2001 From: Emi Griffith Date: Sun, 19 Jul 2026 16:02:33 -0700 Subject: [PATCH] Add a climate-score page from recent-vs-baseline percentile divergence (#196) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Score how far a location's last 6 years have drifted from its full 45-year record. For each metric and percentile category (p10/p25/p50/p75/p90), the recent-years value is placed on the baseline distribution and the gap from the expected percentile is the divergence — unit-free, so metrics compare directly. Scored per meteorological season plus annual, weighted into per-metric and overall scores (temps, humidity and feels-like weighted heaviest). - backend/scoring.py: divergence math, seasonal slicing, precip zero-inflation split (wet-day frequency + amount), tier mapping onto the existing temp scale. - climate.py: derive a wet-bulb column (Stull 2011) at the read boundary, before the humidity column is converted to absolute — via a shared _derive_metrics wrapper at all four read sites. - api/v2/score endpoint + build_score payload, cached on the history token with a scoring-version key. - frontend score page: overall hero, per-metric cards, by-season chips, and a button-revealed summary (sentences + metrics×season table + per-percentile detail). Score nav link across all headers. - Tests for the scoring math, wet-bulb formula, payload shape and route. --- app.py | 36 +++++- climate.py | 42 ++++++- scoring.py | 244 +++++++++++++++++++++++++++++++++++++++++ templates/base.html.j2 | 1 + tests/test_api.py | 48 ++++++++ tests/test_climate.py | 44 ++++++++ tests/test_scoring.py | 181 ++++++++++++++++++++++++++++++ tests/test_views.py | 43 ++++++++ views.py | 28 +++++ 9 files changed, 661 insertions(+), 6 deletions(-) create mode 100644 scoring.py create mode 100644 tests/test_scoring.py diff --git a/app.py b/app.py index 8775b03..95f6605 100644 --- a/app.py +++ b/app.py @@ -177,8 +177,8 @@ async def revalidate_static(request, call_next): "path": path, "status": response.status_code, "cat": cat}) except Exception: # noqa: BLE001 - never let instrumentation break a response pass - pages = (BASE, f"{BASE}/", f"{BASE}/calendar", f"{BASE}/day", f"{BASE}/compare", - f"{BASE}/legend", f"{BASE}/alerts", f"{BASE}/privacy") + pages = (BASE, f"{BASE}/", f"{BASE}/calendar", f"{BASE}/day", f"{BASE}/score", + f"{BASE}/compare", f"{BASE}/legend", f"{BASE}/alerts", f"{BASE}/privacy") if path.endswith((".js", ".css", ".html")) or path in pages: response.headers["Cache-Control"] = "no-cache" return response @@ -430,6 +430,36 @@ def api_day( build) +def api_score( + request: Request, + lat: float = Query(..., ge=-90, le=90), + lon: float = Query(..., ge=-180, le=180), +): + """Climate-drift score for a location: how far the last 6 years have moved + from the full 45-year baseline, per metric, percentile category and season, + rolled into per-metric and overall scores. + + Derived purely from the archive record, so it is cached until the record's + tail advances (the hourly top-up); the scoring-math version is baked into the + cache key so a math change invalidates only score rows. New in API v2. + """ + cell = grid.snap(lat, lon) + + with audit.RunAudit(endpoint="score", lat=round(lat, 4), lon=round(lon, 4), + cell_id=cell["id"]) as run: + history, cache_meta, _ = _fetch_history(run, cell) + + def build(place): + payload = views.build_score(cell, history, place, run) + full = not cache_meta.get("cached", False) + run.set(run_type="full" if full else "partial", + history_source="fetch" if full else "cache") + return payload + + return _cached_response(request, run, "score", cell, + views.score_key(), views.history_token(history), build) + + def api_forecast( request: Request, lat: float = Query(..., ge=-90, le=90), @@ -639,6 +669,7 @@ v2.add_api_route("/place", api_place, methods=["GET"]) v2.add_api_route("/grade", api_grade, methods=["GET"]) v2.add_api_route("/calendar", api_calendar, methods=["GET"]) v2.add_api_route("/day", api_day, methods=["GET"]) +v2.add_api_route("/score", api_score, methods=["GET"]) v2.add_api_route("/forecast", api_forecast, methods=["GET"]) v2.add_api_route("/cell", api_cell, methods=["GET"]) v2.add_api_route("/metrics", api_metrics, methods=["GET"], include_in_schema=False) @@ -709,6 +740,7 @@ if BASE: # readers as real HTML. app.add_api_route(f"{BASE}/calendar", _page("calendar.html"), methods=["GET", "HEAD"], include_in_schema=False) app.add_api_route(f"{BASE}/day", _page("day.html"), methods=["GET", "HEAD"], include_in_schema=False) +app.add_api_route(f"{BASE}/score", _page("score.html"), methods=["GET", "HEAD"], include_in_schema=False) app.add_api_route(f"{BASE}/compare", _page("compare.html"), methods=["GET", "HEAD"], include_in_schema=False) app.add_api_route(f"{BASE}/legend", _page("legend.html"), methods=["GET", "HEAD"], include_in_schema=False) app.add_api_route(f"{BASE}/alerts", _page("subscriptions.html"), methods=["GET", "HEAD"], include_in_schema=False) diff --git a/climate.py b/climate.py index 088a1a8..884391a 100644 --- a/climate.py +++ b/climate.py @@ -215,6 +215,40 @@ def _derive_humidity(df: pl.DataFrame) -> pl.DataFrame: (es * rh * 2.1674 / (273.15 + tmean_c)).round(1).alias("humid")) # abs. humidity, g/m³ +def _derive_wetbulb(df: pl.DataFrame) -> pl.DataFrame: + """Add a `wetbulb` column (°F): the daytime-peak wet-bulb temperature via the + Stull (2011) single-value approximation from the day's high (`tmax`) and mean + relative humidity. + + Wet bulb is the temperature a parcel reaches by evaporative cooling to + saturation — the ceiling on how much the body can shed heat by sweating, so it + is the sharper heat-stress signal than dry-bulb temperature or humidity alone. + Must run BEFORE ``_derive_humidity`` (which replaces the raw RH column with + absolute humidity, the input this needs). Read-time only, like the humidity + derivation, so the parquet cache is untouched. Stull's fit is valid for + RH 5-99% and −20..50 °C; days outside that range grade as null.""" + if "humid" not in df.columns: + return df + t = (pl.col("tmax") - 32.0) * 5.0 / 9.0 # daily high, °C + rh = pl.col("humid").cast(pl.Float64, strict=False) + tw = (t * (0.151977 * (rh + 8.313659).sqrt()).arctan() + + (t + rh).arctan() + - (rh - 1.676331).arctan() + + 0.00391838 * rh.pow(1.5) * (0.023101 * rh).arctan() + - 4.686035) # wet bulb, °C + tw_f = (tw * 9.0 / 5.0 + 32.0).round(1) + valid = (rh >= 5) & (rh <= 99) & (t >= -20) & (t <= 50) + return df.with_columns(pl.when(valid).then(tw_f).otherwise(None).alias("wetbulb")) + + +def _derive_metrics(df: pl.DataFrame) -> pl.DataFrame: + """Read-boundary derivations that depend on the raw relative-humidity column. + Wet bulb must be computed before ``_derive_humidity`` replaces raw RH with + absolute humidity, so both live behind this single wrapper to keep the order + right at every read site.""" + return _derive_humidity(_derive_wetbulb(df)) + + def _with_doy(df: pl.DataFrame) -> pl.DataFrame: """(Re)attach the int16 day-of-year column the grading windows key on.""" return df.with_columns(pl.col("date").dt.ordinal_day().cast(pl.Int16).alias("doy")) @@ -492,7 +526,7 @@ def get_history(cell: dict) -> tuple[pl.DataFrame, dict]: """Return (daily history frame, cache metadata) for a cell, with humidity as absolute humidity (g/m³). Thin wrapper over the raw loader (see below).""" df, meta = _load_history(cell) - return _derive_humidity(df), meta + return _derive_metrics(df), meta def load_cached_history(cell: dict) -> pl.DataFrame | None: @@ -503,7 +537,7 @@ def load_cached_history(cell: dict) -> pl.DataFrame | None: hit = _read_history_cache(_cache_path(cell["id"])) if hit is None: return None - return _derive_humidity(_with_doy(hit[0])) + return _derive_metrics(_with_doy(hit[0])) def _load_history(cell: dict) -> tuple[pl.DataFrame, dict]: @@ -584,7 +618,7 @@ def recent_stamp(cell_id: str) -> int: def get_recent_forecast(cell: dict) -> pl.DataFrame: """Recent observations + forward forecast, with humidity as absolute humidity (g/m³). Thin wrapper over the raw loader (see below).""" - return _derive_humidity(_load_recent_forecast(cell)) + return _derive_metrics(_load_recent_forecast(cell)) def load_cached_recent_forecast(cell: dict) -> "pl.DataFrame | None": @@ -597,7 +631,7 @@ def load_cached_recent_forecast(cell: dict) -> "pl.DataFrame | None": if not os.path.exists(path): return None try: - return _derive_humidity(_with_doy(_normalize_read(pl.read_parquet(path)))) + return _derive_metrics(_with_doy(_normalize_read(pl.read_parquet(path)))) except Exception: # noqa: BLE001 - a truncated/corrupt cache file is a miss, not a crash return None diff --git a/scoring.py b/scoring.py new file mode 100644 index 0000000..141e1ce --- /dev/null +++ b/scoring.py @@ -0,0 +1,244 @@ +"""Climate-shift scoring: how far a cell's recent record (the last ``RECENT_YEARS`` +years) has drifted from its full multi-decade baseline. + +For each metric and each percentile category ``q`` we take the recent-years value +at that percentile and ask where it lands in the baseline distribution. The gap +between where it lands and ``q`` — in percentile points, so it is unit-free and +comparable across metrics — is the divergence. A metric's score is the mean +absolute divergence over the categories, mapped to 0-100; its sign (bias) says +which way it drifted. Metrics are computed per meteorological season and pooled +into an annual view, then weighted (temps / humidity / feels-like heaviest) into +one overall score. + +Baseline is the ENTIRE record, including the recent years — the same all-years +climatology every other view grades against. The recent window's ~13% overlap +mildly attenuates the divergence, uniformly for every cell; it is flagged in the +payload (``baseline_overlaps_recent``). +""" +import numpy as np +import polars as pl + +import grading + +RECENT_YEARS = 6 +QS = (10, 25, 50, 75, 90) # the scored percentile categories +SEASONS = {"djf": (12, 1, 2), "mam": (3, 4, 5), + "jja": (6, 7, 8), "son": (9, 10, 11)} +SLICES = ("annual", "djf", "mam", "jja", "son") + +# All scored metrics. wetbulb is derived at the read boundary (climate.py); the +# rest are raw daily columns. Order here is the canonical compute order; the +# frontend renders in its own display order (Precip first). +SCORE_METRICS = ("tmax", "tmin", "feels", "humid", "wetbulb", "wind", "gust", "precip") + +# Temps / humidity / feels-like dominate; wet bulb medium; wind + precip lightest. +WEIGHTS = {"feels": 2.0, "humid": 2.0, "tmax": 1.5, "tmin": 1.5, + "wetbulb": 1.25, "precip": 0.75, "wind": 0.5, "gust": 0.5} + +# Metrics whose signed drift feeds the headline "warmer / cooler" (so wind and +# precip signs never cancel the temperature story in the overall bias). +TEMP_DIR_METRICS = ("tmax", "tmin", "feels", "wetbulb") + +MIN_BASELINE_YEARS = 25 # below this, no scoring at all (payload carries a reason) +MIN_SLICE_SAMPLES = 300 # recent-slice floor (~540 expected per 6-yr season) +MIN_WET_DAYS = 30 # recent wet-day floor for the precip amount component +SATURATION = 15.0 # mean |divergence| (pct points) that maps to score 100 + +METRIC_LABELS = { + "tmax": "High temp", "tmin": "Low temp", "feels": "Feels like", + "humid": "Humidity", "wetbulb": "Wet bulb", "wind": "Wind", + "gust": "Gusts", "precip": "Precip", +} + +# (word for a positive drift, word for a negative drift). +DIRECTION_WORDS = { + "tmax": ("warmer", "cooler"), "tmin": ("warmer", "cooler"), + "feels": ("warmer", "cooler"), "wetbulb": ("warmer", "cooler"), + "humid": ("muggier", "drier"), "wind": ("windier", "calmer"), + "gust": ("gustier", "calmer"), "precip": ("wetter", "drier"), +} + +# score threshold -> (label, css for positive bias, css for negative bias). The +# css classes are the existing 9-step temperature scale, so no new color tokens +# are needed — the frontend's TIER_COLORS resolves them for free. Lower-bound +# inclusive, checked high-to-low. +TIERS = [ + (80, "Extreme shift", "rec-hot", "rec-cold"), + (60, "Strong shift", "very-hot", "very-cold"), + (35, "Notable shift", "hot", "cold"), + (15, "Mild shift", "warm", "cool"), + (0, "Steady", "normal", "normal"), +] + + +def _minus_years(d, years: int): + """`d` shifted back `years` calendar years, clamping Feb 29 to Feb 28.""" + try: + return d.replace(year=d.year - years) + except ValueError: + return d.replace(year=d.year - years, day=28) + + +def recent_baseline_split(df: pl.DataFrame): + """(recent, baseline): recent = the last ``RECENT_YEARS`` years, baseline = the + full record (the recent years included — see module docstring).""" + latest = df["date"].max() + cutoff = _minus_years(latest, RECENT_YEARS) + return df.filter(pl.col("date") > cutoff), df + + +def _slice(df: pl.DataFrame, key: str) -> pl.DataFrame: + """A slice of the record: the whole thing for ``"annual"``, else the pooled + days of one meteorological season (DJF wraps the year end, but the samples are + pooled across all years so the boundary is irrelevant).""" + if key == "annual": + return df + return df.filter(pl.col("date").dt.month().is_in(list(SEASONS[key]))) + + +def divergence(base: np.ndarray, rec: np.ndarray) -> dict | None: + """Core math for one metric within one slice. For each category ``q`` in + ``QS``: the recent value at that percentile, where it lands in the baseline, + and the signed gap ``d = landed − q``. Returns per-category detail plus the + mean-absolute-divergence (``mad``, drives the score) and mean signed + divergence (``bias``, drives the direction). ``None`` when the recent slice is + too thin to trust.""" + if rec.size < MIN_SLICE_SAMPLES or base.size == 0: + return None + per_q, deltas = [], [] + for q in QS: + v6 = float(np.percentile(rec, q)) + pct = grading.empirical_percentile(base, v6) + if pct is None: + continue + per_q.append({"q": q, "v6": round(v6, 2), "pct": pct, "d": round(pct - q, 1)}) + deltas.append(pct - q) + if not per_q: + return None + a = np.asarray(deltas) + return {"per_q": per_q, + "mad": round(float(np.mean(np.abs(a))), 1), + "bias": round(float(np.mean(a)), 1)} + + +def precip_divergence(base: np.ndarray, rec: np.ndarray) -> dict | None: + """Precip drift, handling its zero-inflation as two parts: how the wet-day + *amounts* shifted (percentile divergence over rain days only) and how the + wet-day *frequency* shifted (points of the wet-day share). They average into + one mad/bias; the frequency alone carries it when rain days are too few for a + stable amount distribution.""" + if rec.size < MIN_SLICE_SAMPLES or base.size == 0: + return None + thr = grading.RAIN_THRESHOLD + f6 = float(np.mean(rec >= thr)) * 100.0 + f45 = float(np.mean(base >= thr)) * 100.0 + freq_d = f6 - f45 + wet_rec = rec[rec >= thr] + amount = divergence(base[base >= thr], wet_rec) if wet_rec.size >= MIN_WET_DAYS else None + if amount is not None: + mad = 0.5 * amount["mad"] + 0.5 * abs(freq_d) + bias = 0.5 * amount["bias"] + 0.5 * freq_d + else: + mad, bias = abs(freq_d), freq_d + return {"per_q": amount["per_q"] if amount else [], + "mad": round(mad, 1), "bias": round(bias, 1), + "freq": {"f6": round(f6, 1), "f45": round(f45, 1), "d": round(freq_d, 1)}} + + +def score_of(mad: float) -> int: + """Mean-absolute-divergence (percentile points) -> 0-100 score. 0 = matches + the baseline; 100 = a drift of ``SATURATION`` points or more.""" + return int(round(100.0 * min(mad / SATURATION, 1.0))) + + +def tier_of(score: int, bias: float): + """(label, css class) for a score, the class picked from the warm or cool + ladder by the sign of ``bias``.""" + for thr, label, hot, cold in TIERS: + if score >= thr: + return label, (hot if bias >= 0 else cold) + return TIERS[-1][1], TIERS[-1][2] + + +def _direction(metric: str, bias: float) -> str: + up, down = DIRECTION_WORDS[metric] + return up if bias >= 0 else down + + +def _null_entry(metric: str, reason: str) -> dict: + return {"key": metric, "label": METRIC_LABELS[metric], "score": None, + "weight": WEIGHTS[metric], "reason": reason} + + +def _metric_entry(metric: str, base: np.ndarray, rec: np.ndarray) -> dict: + div = precip_divergence(base, rec) if metric == "precip" else divergence(base, rec) + if div is None: + return _null_entry(metric, "not enough recent data") + score = score_of(div["mad"]) + tier, css = tier_of(score, div["bias"]) + direction = _direction(metric, div["bias"]) + entry = { + "key": metric, "label": METRIC_LABELS[metric], + "score": score, "mad": div["mad"], "bias": div["bias"], + "tier": tier, "class": css, "direction": direction, + "grade": tier if score < 15 else f"{tier} — {direction}", + "weight": WEIGHTS[metric], "per_q": div["per_q"], + "n_recent": int(rec.size), "n_base": int(base.size), + } + if "freq" in div: + entry["freq"] = div["freq"] + return entry + + +def _overall(entries: dict) -> dict | None: + """Weighted roll-up of the present metrics in one slice. Weights renormalize + over whatever is present (a metric missing for the source drops out cleanly). + The headline bias comes only from the temperature-direction metrics so the + overall reads 'warmer / cooler' rather than being muddied by wind/precip.""" + present = [e for e in entries.values() if e.get("score") is not None] + if not present: + return None + wsum = sum(e["weight"] for e in present) + mad = sum(e["weight"] * e["mad"] for e in present) / wsum + dirs = [e for e in present if e["key"] in TEMP_DIR_METRICS] + bias = (sum(e["weight"] * e["bias"] for e in dirs) / sum(e["weight"] for e in dirs) + if dirs else 0.0) + score = score_of(mad) + tier, css = tier_of(score, bias) + direction = "warmer" if bias >= 0 else "cooler" + return {"score": score, "mad": round(mad, 1), "bias": round(bias, 1), + "tier": tier, "class": css, "direction": direction, + "grade": tier if score < 15 else f"{tier} — {direction}"} + + +def build_scores(history: pl.DataFrame) -> dict: + """Full score payload for a cell's history frame. Returns an ``{"unavailable": + reason}`` shape when the record is too short to score.""" + years = history["date"].dt.year() + span = int(years.max() - years.min() + 1) + if span < MIN_BASELINE_YEARS: + return {"unavailable": f"Only {span} years of record here — climate drift " + f"needs at least {MIN_BASELINE_YEARS}."} + recent, baseline = recent_baseline_split(history) + latest = history["date"].max() + + slices = {} + for key in SLICES: + base_s, rec_s = _slice(baseline, key), _slice(recent, key) + entries = {} + for m in SCORE_METRICS: + if m not in history.columns: + entries[m] = _null_entry(m, "not available for this location") + continue + entries[m] = _metric_entry(m, grading._finite(base_s[m]), grading._finite(rec_s[m])) + slices[key] = {"overall": _overall(entries), "metrics": entries} + + return { + "recent_years": RECENT_YEARS, + "recent_range": [recent["date"].min().isoformat(), latest.isoformat()], + "baseline_range": [history["date"].min().isoformat(), latest.isoformat()], + "n_recent": int(recent.height), + "n_baseline": int(history.height), + "baseline_overlaps_recent": True, + "slices": slices, + } diff --git a/templates/base.html.j2 b/templates/base.html.j2 index c310b4a..841150b 100644 --- a/templates/base.html.j2 +++ b/templates/base.html.j2 @@ -51,6 +51,7 @@ Calendar Day Detail Compare + Score Climate Alerts diff --git a/tests/test_api.py b/tests/test_api.py index b1ab604..a5228cd 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -231,3 +231,51 @@ def test_warm_cell_never_fetches_upstream(client, monkeypatch): monkeypatch.setattr(climate, "get_recent_forecast", boom) monkeypatch.setattr(climate, "load_cached_history", lambda cell: None) appmod._warm_cell({"id": "1_1", "center_lat": 0.03, "center_lon": 0.03}) # cold: no-op + + +def _score_history(years=45, seed=5): + import datetime + import numpy as np + import polars as pl + end = datetime.date(2026, 7, 11) + start = datetime.date(end.year - years, end.month, end.day) + dates = [start + datetime.timedelta(days=i) for i in range((end - start).days + 1)] + n = len(dates) + rng = np.random.default_rng(seed) + doy = np.array([d.timetuple().tm_yday for d in dates]) + tmax = 55 + 30 * np.sin((doy - 100) / 366.0 * 2 * np.pi) + rng.normal(0, 8, n) + return pl.DataFrame({ + "date": dates, "tmax": np.round(tmax, 1), "tmin": np.round(tmax - 15, 1), + "feels": np.round(tmax + 1, 1), "humid": np.round(np.clip(12 + rng.normal(0, 3, n), 1, None), 1), + "wetbulb": np.round(tmax - 12, 1), "wind": np.round(np.clip(8 + rng.normal(0, 3, n), 0, None), 1), + "gust": np.round(np.clip(16 + rng.normal(0, 5, n), 0, None), 1), + "precip": np.where(rng.random(n) < 0.3, 0.2, 0.0), + }).with_columns(pl.col("date").dt.ordinal_day().cast(pl.Int16).alias("doy")) + + +@pytest.fixture +def score_client(monkeypatch): + hist = _score_history() + monkeypatch.setattr(climate, "get_history", + lambda cell: (hist.clone(), {"cached": True, "cache_age_days": 3})) + monkeypatch.setattr(climate, "reverse_geocode", lambda lat, lon: "Testville, Washington") + return TestClient(appmod.app) + + +def test_score_shape_and_conditional_revalidation(score_client): + r = score_client.get("/thermograph/api/v2/score", params=Q) + assert r.status_code == 200 + body = r.json() + assert body["place"] == "Testville, Washington" + s = body["scores"] + assert set(s["slices"]) == {"annual", "djf", "mam", "jja", "son"} + assert s["slices"]["annual"]["overall"]["score"] is not None + + etag = r.headers["etag"] + r304 = score_client.get("/thermograph/api/v2/score", params=Q, + headers={"If-None-Match": etag}) + assert r304.status_code == 304 and r304.headers["etag"] == etag + + # Second GET replays the exact same bytes from the derived store. + r2 = score_client.get("/thermograph/api/v2/score", params=Q) + assert r2.status_code == 200 and r2.content == r.content diff --git a/tests/test_climate.py b/tests/test_climate.py index 283abe3..8091594 100644 --- a/tests/test_climate.py +++ b/tests/test_climate.py @@ -3,6 +3,7 @@ no network).""" import datetime import polars as pl +import pytest import climate @@ -172,3 +173,46 @@ def test_write_cache_strips_the_derived_doy(tmp_path): stored = pl.read_parquet(path) assert "doy" not in stored.columns assert len(stored) == len(df) + + +def _wb_frame(tmax, humid, tmin=None): + """Frame with raw RH (`humid`) for the wet-bulb / humidity derivations.""" + n = len(tmax) + return pl.DataFrame({ + "date": [datetime.date(2026, 7, 1) + datetime.timedelta(days=i) for i in range(n)], + "tmax": [float(x) for x in tmax], + "tmin": [float(t) for t in (tmin or [x - 15 for x in tmax])], + "humid": [float(h) for h in humid], + }) + + +def test_wetbulb_stull_spot_values(): + # Stull (2011): 20 °C / 50% -> ~13.7 °C; 30 °C / 80% -> ~27.2 °C. In °F here. + df = climate._derive_wetbulb(_wb_frame([68.0, 86.0], [50.0, 80.0])) + wb = df["wetbulb"].to_list() + assert wb[0] == pytest.approx((13.7 * 9 / 5) + 32, abs=0.6) + assert wb[1] == pytest.approx((27.2 * 9 / 5) + 32, abs=0.6) + + +def test_wetbulb_never_exceeds_dry_bulb(): + df = climate._derive_wetbulb(_wb_frame([40.0, 60.0, 85.0, 100.0], [20.0, 55.0, 70.0, 95.0])) + for tmax, wb in zip(df["tmax"], df["wetbulb"]): + assert wb is not None and wb <= tmax + 0.05 + + +def test_wetbulb_null_outside_validity_range(): + # RH 3% (< 5) and a scorching 130 °F (> 50 °C) both fall outside Stull's fit. + df = climate._derive_wetbulb(_wb_frame([70.0, 130.0], [3.0, 40.0])) + assert df["wetbulb"].to_list() == [None, None] + + +def test_wetbulb_noop_without_humidity(): + df = pl.DataFrame({"tmax": [70.0], "tmin": [55.0]}) + assert "wetbulb" not in climate._derive_wetbulb(df).columns + + +def test_derive_metrics_adds_wetbulb_and_absolute_humidity(): + out = climate._derive_metrics(_wb_frame([68.0], [50.0])) + assert "wetbulb" in out.columns + # _derive_humidity replaces raw RH (50) with absolute humidity (g/m³, ~8-9). + assert out["humid"][0] < 30 and out["humid"][0] != 50.0 diff --git a/tests/test_scoring.py b/tests/test_scoring.py new file mode 100644 index 0000000..1cd5453 --- /dev/null +++ b/tests/test_scoring.py @@ -0,0 +1,181 @@ +"""Divergence-scoring tests: synthetic records with known shifts drive known +scores. The recent window is always the last 6 years OF THE SAME frame, so a +"shift" here means perturbing only those trailing years.""" +import datetime + +import numpy as np +import polars as pl +import pytest + +import scoring + +ALL_METRICS = ("tmax", "tmin", "feels", "humid", "wetbulb", "wind", "gust", "precip") + + +def _years_before(d: datetime.date, years: int) -> datetime.date: + try: + return d.replace(year=d.year - years) + except ValueError: + return d.replace(year=d.year - years, day=28) + + +def make_frame(years: int = 45, seed: int = 3, drop: tuple = (), end: str | None = None): + """A full multi-metric daily record, every metric i.i.d. across years around a + seasonal cycle — so the last 6 years are statistically identical to the rest + until a test perturbs them. Columns in ``drop`` are omitted (to exercise the + missing-metric path).""" + end_ts = datetime.date.fromisoformat(end) if end else datetime.date(2026, 7, 11) + start = _years_before(end_ts, years) + dates = [start + datetime.timedelta(days=i) for i in range((end_ts - start).days + 1)] + n = len(dates) + rng = np.random.default_rng(seed) + doy = np.array([d.timetuple().tm_yday for d in dates]) + seasonal = 55 + 30 * np.sin((doy - 100) / 366.0 * 2 * np.pi) + tmax = seasonal + rng.normal(0, 8, n) + tmin = tmax - 15 + rng.normal(0, 3, n) + cols = { + "date": dates, + "tmax": np.round(tmax, 1), + "tmin": np.round(tmin, 1), + "feels": np.round(tmax + rng.normal(0, 2, n), 1), + "humid": np.round(np.clip(12 + rng.normal(0, 3, n), 1, None), 1), + "wetbulb": np.round(tmax - 12 + rng.normal(0, 2, n), 1), + "wind": np.round(np.clip(8 + rng.normal(0, 3, n), 0, None), 1), + "gust": np.round(np.clip(16 + rng.normal(0, 5, n), 0, None), 1), + "precip": np.where(rng.random(n) < 0.3, np.round(rng.gamma(1.5, 0.2, n), 2), 0.0), + } + for d in drop: + cols.pop(d) + df = pl.DataFrame(cols) + return df.with_columns(pl.col("date").dt.ordinal_day().cast(pl.Int16).alias("doy")) + + +def _recent_mask(df: pl.DataFrame) -> np.ndarray: + cutoff = scoring._minus_years(df["date"].max(), scoring.RECENT_YEARS) + return (df["date"] > cutoff).to_numpy() + + +def shift_metric(df, metric, amount, months=None): + """Add ``amount`` to ``metric`` for the recent-6-years rows (optionally only + within ``months``) — the perturbation the scorer should detect.""" + mask = _recent_mask(df) + if months is not None: + mask = mask & df["date"].dt.month().is_in(list(months)).to_numpy() + vals = df[metric].to_numpy().copy() + vals[mask] = vals[mask] + amount + return df.with_columns(pl.Series(metric, np.round(vals, 2))) + + +# --- baseline / no-shift ------------------------------------------------------ + +def test_no_shift_is_steady(): + out = scoring.build_scores(make_frame()) + ann = out["slices"]["annual"] + assert ann["overall"]["score"] < 15 # Steady + for k in ("tmax", "tmin", "feels", "humid"): + m = ann["metrics"][k] + assert m["score"] is not None + assert abs(m["bias"]) <= 5 # within sampling noise + assert m["mad"] <= 5 + + +def test_payload_metadata(): + out = scoring.build_scores(make_frame(years=45)) + assert out["recent_years"] == 6 + assert out["baseline_overlaps_recent"] is True + assert out["n_recent"] < out["n_baseline"] + assert set(out["slices"]) == set(scoring.SLICES) + assert out["recent_range"][1] == out["baseline_range"][1] + + +# --- known directional shift -------------------------------------------------- + +def test_warming_shift_scores_up_and_warmer(): + base = scoring.build_scores(make_frame())["slices"]["annual"]["metrics"]["tmax"] + hot = scoring.build_scores(shift_metric(make_frame(), "tmax", 9.0)) + m = hot["slices"]["annual"]["metrics"]["tmax"] + assert m["score"] >= 30 and m["score"] > base["score"] + assert m["bias"] > 0 + assert m["direction"] == "warmer" + assert "warmer" in m["grade"] + assert all(q["d"] > 0 for q in m["per_q"]) # every category shifted up + # The overall headline follows the temperature drift. + assert hot["slices"]["annual"]["overall"]["direction"] == "warmer" + + +def test_cooling_shift_reads_cooler(): + out = scoring.build_scores(shift_metric(make_frame(), "tmin", -9.0)) + m = out["slices"]["annual"]["metrics"]["tmin"] + assert m["bias"] < 0 + assert m["direction"] == "cooler" + assert m["class"] in ("cool", "cold", "very-cold", "rec-cold") + + +# --- seasonal isolation ------------------------------------------------------- + +def test_summer_only_shift_isolated_to_jja(): + out = scoring.build_scores(shift_metric(make_frame(), "tmax", 12.0, months=(6, 7, 8))) + jja = out["slices"]["jja"]["metrics"]["tmax"] + djf = out["slices"]["djf"]["metrics"]["tmax"] + assert jja["mad"] > djf["mad"] + assert jja["score"] > djf["score"] + assert jja["bias"] > 0 + + +# --- precip zero-inflation ---------------------------------------------------- + +def test_precip_frequency_shift(): + df = make_frame() + mask = _recent_mask(df) + rng = np.random.default_rng(99) + p = df["precip"].to_numpy().copy() + # Double the wet-day frequency in the recent window at similar amounts. + extra = mask & (rng.random(len(p)) < 0.35) & (p < scoring.grading.RAIN_THRESHOLD) + p[extra] = 0.12 + df = df.with_columns(pl.Series("precip", np.round(p, 2))) + m = scoring.build_scores(df)["slices"]["annual"]["metrics"]["precip"] + assert m["freq"]["d"] > 8 # wetter days more often + assert m["direction"] == "wetter" + assert m["bias"] > 0 + + +# --- missing metric / weight renormalization ---------------------------------- + +def test_missing_gust_column_nulls_out_but_overall_survives(): + out = scoring.build_scores(make_frame(drop=("gust",))) + ann = out["slices"]["annual"] + g = ann["metrics"]["gust"] + assert g["score"] is None and "not available" in g["reason"] + assert ann["overall"]["score"] is not None # renormalized over present metrics + + +def test_all_null_metric_reports_not_enough_data(): + df = make_frame().with_columns(pl.lit(None, dtype=pl.Float64).alias("humid")) + m = scoring.build_scores(df)["slices"]["annual"]["metrics"]["humid"] + assert m["score"] is None + assert "not enough" in m["reason"] + + +# --- guards ------------------------------------------------------------------- + +def test_short_record_is_unavailable(): + out = scoring.build_scores(make_frame(years=10)) + assert "unavailable" in out + assert "slices" not in out + + +# --- unit helpers ------------------------------------------------------------- + +def test_score_and_tier_mapping(): + assert scoring.score_of(0) == 0 + assert scoring.score_of(scoring.SATURATION) == 100 + assert scoring.score_of(scoring.SATURATION * 2) == 100 # clamped + assert scoring.tier_of(90, 5)[1] == "rec-hot" + assert scoring.tier_of(90, -5)[1] == "rec-cold" + assert scoring.tier_of(5, 1)[0] == "Steady" + + +def test_divergence_none_below_min_samples(): + base = np.linspace(0, 100, 5000) + rec = np.linspace(0, 100, scoring.MIN_SLICE_SAMPLES - 1) + assert scoring.divergence(base, rec) is None diff --git a/tests/test_views.py b/tests/test_views.py index 9cb0c8f..7893f7c 100644 --- a/tests/test_views.py +++ b/tests/test_views.py @@ -141,3 +141,46 @@ def test_build_forecast_only_future_days(history, recent): assert days == sorted(days, reverse=True) # furthest-out first assert min(days) > today.isoformat() assert payload["forecast"] is True + + +# ---- build_score --------------------------------------------------------------- + +def _full_history(years=45, seed=5): + """A 45-year all-metric record — enough span for the climate score (the shared + 20-year `history` fixture is intentionally below MIN_BASELINE_YEARS).""" + import numpy as np + end = datetime.date(2026, 7, 11) + start = datetime.date(end.year - years, end.month, end.day) + dates = [start + datetime.timedelta(days=i) for i in range((end - start).days + 1)] + n = len(dates) + rng = np.random.default_rng(seed) + doy = np.array([d.timetuple().tm_yday for d in dates]) + tmax = 55 + 30 * np.sin((doy - 100) / 366.0 * 2 * np.pi) + rng.normal(0, 8, n) + return pl.DataFrame({ + "date": dates, "tmax": np.round(tmax, 1), "tmin": np.round(tmax - 15, 1), + "feels": np.round(tmax + 1, 1), "humid": np.round(np.clip(12 + rng.normal(0, 3, n), 1, None), 1), + "wetbulb": np.round(tmax - 12, 1), "wind": np.round(np.clip(8 + rng.normal(0, 3, n), 0, None), 1), + "gust": np.round(np.clip(16 + rng.normal(0, 5, n), 0, None), 1), + "precip": np.where(rng.random(n) < 0.3, 0.2, 0.0), + }).with_columns(pl.col("date").dt.ordinal_day().cast(pl.Int16).alias("doy")) + + +def test_build_score_shape(): + hist = _full_history() + payload = views.build_score(CELL, hist, "Testville") + assert payload["api_version"] == "v2" + assert payload["cell"] == CELL and payload["place"] == "Testville" + assert payload["latest"] == views.hist_end(hist) + s = payload["scores"] + assert set(s["slices"]) == {"annual", "djf", "mam", "jja", "son"} + ann = s["slices"]["annual"] + assert ann["overall"]["score"] is not None + assert ann["metrics"]["tmax"]["score"] is not None + assert ann["metrics"]["wetbulb"]["score"] is not None # derived metric flows through + + +def test_score_key_is_the_version(): + assert views.score_key() == views.SCORE_VER + # Score payloads are history-only, so they ride the plain history token. + hist = _full_history() + assert views.history_token(hist) == f"{views.PAYLOAD_VER}:{views.hist_end(hist)}" diff --git a/views.py b/views.py index ae246b3..bba612d 100644 --- a/views.py +++ b/views.py @@ -15,6 +15,7 @@ import polars as pl import climate import grading +import scoring # Observed values pulled from a daily record row for grading. Includes the # temperature-scale metrics (tmax/tmin/feels/wind/gust) plus precip; a column may @@ -103,6 +104,17 @@ def forecast_key(today, days: int) -> str: return f"{today.isoformat()}:{days}" +# The score is derived purely from the full archive record, so its validity token +# is the plain `history_token` (expires when the archive tail advances). The key +# just carries the scoring-math version, so a math change invalidates only score +# rows without disturbing any other kind. +SCORE_VER = "s1" + + +def score_key() -> str: + return SCORE_VER + + def _obs_from_row(row: dict) -> dict: return {k: row[k] for k in OBS_COLS if k in row} @@ -276,3 +288,19 @@ def build_forecast(cell, days, history, fc, today, place, run=None) -> dict: "climatology": climo, "recent": graded, } + + +def build_score(cell, history, place, run=None) -> dict: + """/score payload: how far the cell's last 6 years have drifted from its full + 45-year baseline, per metric and season (see scoring.build_scores).""" + run = run or NullRun() + with run.phase("scoring"): + scores = scoring.build_scores(history) + run.set(place_found=place is not None) + return { + "api_version": "v2", + "cell": cell, + "place": place, + "latest": hist_end(history), + "scores": scores, + }