From 8d90dcb456b6c6ab2b61eefb8321bbf808f8d5d0 Mon Sep 17 00:00:00 2001 From: Emi Griffith Date: Fri, 10 Jul 2026 17:29:47 -0700 Subject: [PATCH] Initial commit: app + VPS deploy pipeline --- app.py | 444 +++++++++++++++++++++++++++++++++++++++ audit.py | 130 ++++++++++++ climate.py | 527 +++++++++++++++++++++++++++++++++++++++++++++++ grading.py | 369 +++++++++++++++++++++++++++++++++ grid.py | 51 +++++ requirements.txt | 6 + 6 files changed, 1527 insertions(+) create mode 100644 app.py create mode 100644 audit.py create mode 100644 climate.py create mode 100644 grading.py create mode 100644 grid.py create mode 100644 requirements.txt diff --git a/app.py b/app.py new file mode 100644 index 0000000..8ae8175 --- /dev/null +++ b/app.py @@ -0,0 +1,444 @@ +"""Thermograph API — grade recent local weather against ~45 years of climatology.""" +import datetime +import os +import time + +import pandas as pd +from fastapi import APIRouter, FastAPI, HTTPException, Query +from fastapi.responses import FileResponse, RedirectResponse +from fastapi.staticfiles import StaticFiles + +import audit +import climate +import grading +import grid + +FRONTEND_DIR = os.path.join(os.path.dirname(__file__), "..", "frontend") + +# Everything (pages, assets, API) is served under this base path so the app can +# live at https:///thermograph/ behind a shared domain. Override with the +# THERMOGRAPH_BASE env var. The frontend uses base-relative URLs, so it follows +# this automatically without hardcoding the prefix. +BASE = "/" + os.environ.get("THERMOGRAPH_BASE", "/thermograph").strip("/") + +# 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 +# be absent on an older cache, so only carry the ones present. +OBS_COLS = ("tmax", "tmin", "precip", "feels", "humid", "wind", "gust") + + +def _obs_from_row(row) -> dict: + return {k: row[k] for k in OBS_COLS if k in row} + + +def _weather_fetch_error(e) -> HTTPException: + """A clean, retryable message for a rate limit; the raw error otherwise.""" + low = str(e).lower() + reason = (climate._rate_limit_reason(e) or "").lower() + blob = low + " " + reason + if climate._is_rate_limit(e) or "429" in low or "rate-limited" in low or "limit" in reason: + if "daily" in blob or "tomorrow" in blob: # daily quota exhausted — resets tomorrow + return HTTPException( + status_code=503, + detail="Open-Meteo's daily request limit is exhausted — new locations will work again " + "tomorrow. Places you've already viewed still work.", + ) + return HTTPException( + status_code=503, + detail="The weather service is rate-limited right now — please try again in a minute.", + ) + return HTTPException(status_code=502, detail=f"weather data fetch failed: {e}") + + +def _grade_rows(history, rows_df) -> list[dict]: + """Grade each daily row against its own ±7-day climatology window.""" + return [grading.grade_day(history, row["date"], _obs_from_row(row)) + for _, row in rows_df.iterrows()] + + +def _attach_dry_streaks(graded: list[dict], *precip_frames) -> None: + """Attach `dsr` (days since last measurable rain) to each graded day, computed + over a combined, de-duplicated precip series so streaks stay continuous across + the history / recent / forecast sources.""" + frames = [f[["date", "precip"]] for f in precip_frames if f is not None and not f.empty] + if not frames: + return + combined = (pd.concat(frames) + .drop_duplicates(subset="date", keep="last") + .sort_values("date")) + dsr_map = grading.dry_streaks(combined["date"].values, combined["precip"].values) + for g in graded: + g["dsr"] = dsr_map.get(g["date"]) + + +app = FastAPI(title="Thermograph", version="0.2.0") + + +@app.middleware("http") +async def no_store_static(request, call_next): + """Serve the frontend with no-store so phones/browsers always pick up the + latest HTML/JS/CSS instead of a stale cached copy. This is a LAN dev tool, so + skipping the cache entirely is fine and avoids "my change isn't showing" bugs.""" + response = await call_next(request) + path = request.url.path + pages = (BASE, f"{BASE}/", f"{BASE}/calendar", f"{BASE}/day", f"{BASE}/legend") + if path.endswith((".js", ".css", ".html")) or path in pages: + response.headers["Cache-Control"] = "no-store, must-revalidate" + return response + + +def api_geocode(q: str = Query(..., min_length=1)): + try: + return {"results": climate.geocode(q)} + except Exception as e: # noqa: BLE001 - surface upstream failures to the client + raise HTTPException(status_code=502, detail=f"geocoding failed: {e}") + + +def api_grade( + lat: float = Query(..., ge=-90, le=90), + lon: float = Query(..., ge=-180, le=180), + date: str | None = Query(None, description="target date YYYY-MM-DD (default today)"), + days: int = Query(14, ge=1, le=60, description="how many recent days to grade"), +): + if not grid.in_north_america(lat, lon): + raise HTTPException( + status_code=400, + detail="Location is outside the supported US + Canada region.", + ) + + target = ( + pd.Timestamp(date).normalize() + if date + else pd.Timestamp(datetime.date.today()) + ) + cell = grid.snap(lat, lon) + + with audit.RunAudit( + endpoint="grade", + lat=round(lat, 4), + lon=round(lon, 4), + target_date=target.date().isoformat(), + recent_days=days, + cell_id=cell["id"], + ) as run: + try: + with run.phase("history"): + history, cache_meta = climate.get_history(cell) + with run.phase("recent"): + recent = climate.get_recent_forecast(cell) + except Exception as e: # noqa: BLE001 + raise _weather_fetch_error(e) + + if history.empty: + raise HTTPException(status_code=404, detail="No historical data for this cell.") + + with run.phase("grading"): + # Only grade observed days up to and including the target date (last `days`). + recent = recent[recent["date"] <= target].sort_values("date").tail(days) + graded = _grade_rows(history, recent) + _attach_dry_streaks(graded, history, recent) + graded.reverse() # newest first for display + climo = grading.climatology(history, int(target.dayofyear)) + + with run.phase("reverse_geocode"): + place = climate.reverse_geocode(cell["center_lat"], cell["center_lon"]) + + # Summarize what this run actually covered. A fresh history fetch pulls the + # full ~45-year archive (run_type "full"); a cache hit only fetched the + # recent window (run_type "partial"). + years = pd.to_datetime(history["date"]).dt.year + full = not cache_meta.get("cached", False) + run.set( + run_type="full" if full else "partial", + history_source="fetch" if full else "cache", + history_rows=int(len(history)), + history_years=int(years.max() - years.min() + 1), + history_span=[int(years.min()), int(years.max())], + cache_age_days=cache_meta.get("cache_age_days"), + graded_days=len(graded), + place_found=place is not None, + ) + + return { + "cell": cell, + "place": place, + "target_date": target.date().isoformat(), + "cache": cache_meta, + "climatology": climo, + "recent": graded, + } + + +CAL_MAX_SPAN_DAYS = 732 # ~2 years — the most a single calendar request will grade + # (the frontend splits longer spans into 2-year chunks) + +# Short-lived in-memory cache of full calendar responses, keyed by (cell, span). +# Lets the home page prefetch the calendar in the background so navigating to it is +# instant, and makes repeat calendar loads free. Bounded + TTL'd. +CAL_CACHE_TTL = 600 # seconds (10 min) +_CAL_CACHE: dict[tuple, tuple[float, dict]] = {} + + +def _cal_cache_get(key): + hit = _CAL_CACHE.get(key) + if hit and (time.time() - hit[0]) < CAL_CACHE_TTL: + return hit[1] + return None + + +def _cal_cache_put(key, resp): + _CAL_CACHE[key] = (time.time(), resp) + if len(_CAL_CACHE) > 64: # evict the oldest entry to stay bounded + _CAL_CACHE.pop(min(_CAL_CACHE, key=lambda k: _CAL_CACHE[k][0]), None) + + +def api_calendar( + lat: float = Query(..., ge=-90, le=90), + lon: float = Query(..., ge=-180, le=180), + start: str | None = Query(None, description="first date YYYY-MM-DD (overrides months)"), + end: str | None = Query(None, description="last date YYYY-MM-DD (default = latest available)"), + months: int = Query(24, ge=1, le=24, description="months back to grade when no start given"), +): + """Grade every day over a date range (default: last 2 years) for the calendar view. + + A custom [start, end] range is supported and capped at ~2 years per request; + the frontend splits longer spans into successive 2-year chunks. Without a + start it falls back to `months` back from `end`. Grades against the already- + cached history record (which reaches to ~6 days ago), so no extra fetch is + needed. Returns a compact per-day shape (see grading.grade_range). New in API v2. + """ + if not grid.in_north_america(lat, lon): + raise HTTPException( + status_code=400, + detail="Location is outside the supported US + Canada region.", + ) + cell = grid.snap(lat, lon) + + with audit.RunAudit( + endpoint="calendar", + lat=round(lat, 4), + lon=round(lon, 4), + recent_days=months * 31, # approximate span graded + cell_id=cell["id"], + ) as run: + try: + with run.phase("history"): + history, cache_meta = climate.get_history(cell) + except Exception as e: # noqa: BLE001 + raise _weather_fetch_error(e) + if history.empty: + raise HTTPException(status_code=404, detail="No historical data for this cell.") + + first = pd.Timestamp(history["date"].min()).normalize() + last = pd.Timestamp(history["date"].max()).normalize() + end_ts = min(pd.Timestamp(end).normalize(), last) if end else last + if start: + start_ts = pd.Timestamp(start).normalize() + else: + # Default: whole months back landing on the SAME day-of-month as the end, + # so the two range endpoints line up (e.g. 2024-06-29 → 2026-06-29). + start_ts = (end_ts - pd.DateOffset(months=months)).normalize() + # Cap the graded span at ~2 years, and keep it inside the available record. + min_start = end_ts - pd.Timedelta(days=CAL_MAX_SPAN_DAYS - 1) + start_ts = max(start_ts, min_start, first) + start_ts = min(start_ts, end_ts) + + cache_key = (cell["id"], start_ts.date().isoformat(), end_ts.date().isoformat()) + cached = _cal_cache_get(cache_key) + if cached is not None: + run.set(run_type="cache", history_source="cache", graded_days=len(cached["days"])) + return cached + + with run.phase("grading"): + days = grading.grade_range(history, start_ts, end_ts) + + with run.phase("reverse_geocode"): + place = climate.reverse_geocode(cell["center_lat"], cell["center_lon"]) + + full = not cache_meta.get("cached", False) + run.set( + run_type="full" if full else "partial", + history_source="fetch" if full else "cache", + graded_days=len(days), + ) + + resp = { + "api_version": "v2", + "cell": cell, + "place": place, + "range": {"start": start_ts.date().isoformat(), "end": end_ts.date().isoformat()}, + "months": months, + "days": days, + } + _cal_cache_put(cache_key, resp) + return resp + + +def api_day( + lat: float = Query(..., ge=-90, le=90), + lon: float = Query(..., ge=-180, le=180), + date: str | None = Query(None, description="target date YYYY-MM-DD (default = latest available)"), +): + """Full percentile breakdown for a single day: the value at every tier boundary + in that day-of-year's ±7-day window, plus where the observed values land. + + Powers the single-day detail page. Grades against the cached history record; + for a date newer than the cache it pulls the recent window for the observation. + New in API v2. + """ + if not grid.in_north_america(lat, lon): + raise HTTPException( + status_code=400, + detail="Location is outside the supported US + Canada region.", + ) + cell = grid.snap(lat, lon) + + with audit.RunAudit( + endpoint="day", + lat=round(lat, 4), + lon=round(lon, 4), + cell_id=cell["id"], + ) as run: + with run.phase("history"): + history, cache_meta = climate.get_history(cell) + if history.empty: + raise HTTPException(status_code=404, detail="No historical data for this cell.") + + last = pd.Timestamp(history["date"].max()).normalize() + target = pd.Timestamp(date).normalize() if date else last + run.set(target_date=target.date().isoformat()) + + # Observed values for the target date: prefer the cached history record; if + # the date is newer than the cache (last few days), fetch the recent window. + obs = None + hit = history[history["date"] == target] + if not hit.empty: + obs = _obs_from_row(hit.iloc[0]) + elif target > last: + try: + with run.phase("recent"): + recent = climate.get_recent_forecast(cell) + rr = recent[recent["date"] == target] + if not rr.empty: + obs = _obs_from_row(rr.iloc[0]) + except Exception: # noqa: BLE001 - detail page still works from climatology alone + obs = None + + with run.phase("detail"): + detail = grading.day_detail(history, target, obs) + + with run.phase("reverse_geocode"): + place = climate.reverse_geocode(cell["center_lat"], cell["center_lon"]) + + full = not cache_meta.get("cached", False) + run.set( + run_type="full" if full else "partial", + history_source="fetch" if full else "cache", + has_observation=obs is not None, + ) + + return { + "api_version": "v2", + "cell": cell, + "place": place, + "latest": last.date().isoformat(), + "detail": detail, + } + + +def api_forecast( + lat: float = Query(..., ge=-90, le=90), + lon: float = Query(..., ge=-180, le=180), + days: int = Query(7, ge=1, le=14, description="how many forecast days ahead to grade"), +): + """Grade the next `days` of forecast weather against local climatology. + + Same shape as /grade (so the frontend renders it identically), but `recent` + holds the forward forecast — furthest-out day first. The forecast is fetched + fresh and cached only ~1 hour (see climate.get_recent_forecast) to track updates. + New in API v2. + """ + if not grid.in_north_america(lat, lon): + raise HTTPException( + status_code=400, + detail="Location is outside the supported US + Canada region.", + ) + cell = grid.snap(lat, lon) + + with audit.RunAudit(endpoint="forecast", lat=round(lat, 4), lon=round(lon, 4), + recent_days=days, cell_id=cell["id"]) as run: + try: + with run.phase("history"): + history, _ = climate.get_history(cell) + with run.phase("forecast"): + fc = climate.get_recent_forecast(cell) + except Exception as e: # noqa: BLE001 + raise _weather_fetch_error(e) + + if history.empty: + raise HTTPException(status_code=404, detail="No historical data for this cell.") + + today = pd.Timestamp(datetime.date.today()) + # Strictly future days (tomorrow onward), earliest→latest, capped at `days`. + future = fc[fc["date"] > today].sort_values("date").head(days) + + with run.phase("grading"): + graded = _grade_rows(history, future) + _attach_dry_streaks(graded, history, fc) + graded.reverse() # furthest-out first, so the chart reverses to L→R chronological + climo = grading.climatology(history, int(today.dayofyear)) + + with run.phase("reverse_geocode"): + place = climate.reverse_geocode(cell["center_lat"], cell["center_lon"]) + + run.set(run_type="partial", graded_days=len(graded), place_found=place is not None) + + return { + "api_version": "v2", + "cell": cell, + "place": place, + "target_date": today.date().isoformat(), + "forecast": True, + "climatology": climo, + "recent": graded, + } + + +# --- API versioning -------------------------------------------------------- +# Non-backward-compatible additions land in a new version; every version stays +# mounted and served simultaneously. v1 is the original contract (grade/geocode); +# v2 adds the calendar. The unversioned /api/* paths are kept as aliases of v1 so +# existing clients keep working. +v1 = APIRouter(tags=["v1"]) +v1.add_api_route("/geocode", api_geocode, methods=["GET"]) +v1.add_api_route("/grade", api_grade, methods=["GET"]) + +v2 = APIRouter(tags=["v2"]) +v2.add_api_route("/geocode", api_geocode, 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("/forecast", api_forecast, methods=["GET"]) + +app.include_router(v1, prefix=f"{BASE}/api") # legacy unversioned == v1 (backward compatible) +app.include_router(v1, prefix=f"{BASE}/api/v1") +app.include_router(v2, prefix=f"{BASE}/api/v2") + + +# --- static frontend (served under BASE) ----------------------------------- +def _page(name): + return FileResponse(os.path.join(FRONTEND_DIR, name)) + + +# Bare domain and the un-slashed base both land on the app's index. The trailing +# slash matters: the frontend's relative asset URLs resolve against BASE/ . +app.add_api_route("/", lambda: RedirectResponse(url=f"{BASE}/"), methods=["GET"], include_in_schema=False) +app.add_api_route(BASE, lambda: RedirectResponse(url=f"{BASE}/"), methods=["GET"], include_in_schema=False) +app.add_api_route(f"{BASE}/", lambda: _page("index.html"), methods=["GET"], include_in_schema=False) +app.add_api_route(f"{BASE}/calendar", lambda: _page("calendar.html"), methods=["GET"], include_in_schema=False) +app.add_api_route(f"{BASE}/day", lambda: _page("day.html"), methods=["GET"], include_in_schema=False) +app.add_api_route(f"{BASE}/legend", lambda: _page("legend.html"), methods=["GET"], include_in_schema=False) + +# Everything else under BASE (app.js, style.css, nav.js, …) is a static asset. +# Registered last so the explicit page routes above win. +app.mount(BASE, StaticFiles(directory=FRONTEND_DIR), name="static") diff --git a/audit.py b/audit.py new file mode 100644 index 0000000..29cadda --- /dev/null +++ b/audit.py @@ -0,0 +1,130 @@ +"""Lightweight file-based audit + error logging for Thermograph. + +Two log streams, written as JSON Lines (one self-contained JSON object per line): + +* ``logs/audit/audit-.jsonl`` — one line per graded run, recording the + total wall-clock duration, a per-phase timing breakdown, how many days the run + covered (the recent-window size and the historical span), and a ``run_type`` + tag of **full** (a fresh ~45-year history fetch happened) or **partial** (the + history came from cache, so only the recent days were fetched). + +* ``logs/errors/errors-.jsonl`` — a *separate folder* for anything that + went wrong: each upstream failure is tagged ``retry`` (another attempt will be + made) or ``error`` (final failure), and whole-run failures are tagged ``error``. + +Everything is best-effort: logging never raises into the request path. +""" +import contextlib +import contextvars +import datetime +import json +import os +import threading +import time +import uuid + +LOG_ROOT = os.path.join(os.path.dirname(__file__), "..", "logs") +AUDIT_DIR = os.path.join(LOG_ROOT, "audit") +ERROR_DIR = os.path.join(LOG_ROOT, "errors") + +_lock = threading.Lock() + +# Set for the duration of a run so error/retry logs from deep in the fetch layer +# can be correlated back to the run that triggered them. +current_run_id: "contextvars.ContextVar[str | None]" = contextvars.ContextVar( + "thermograph_run_id", default=None +) + + +def _now_iso() -> str: + return datetime.datetime.now().astimezone().isoformat(timespec="milliseconds") + + +def _today() -> str: + return datetime.date.today().isoformat() + + +def _append(path: str, record: dict) -> None: + try: + os.makedirs(os.path.dirname(path), exist_ok=True) + line = json.dumps(record, default=str) + with _lock: + with open(path, "a", encoding="utf-8") as f: + f.write(line + "\n") + except Exception: # noqa: BLE001 - logging must never break the request + pass + + +def log_event(tag: str, record: dict) -> dict: + """Record an error or retry in the separate errors folder. + + ``tag`` is conventionally ``"retry"`` or ``"error"``. + """ + rec = {"ts": _now_iso(), "tag": tag, "run_id": current_run_id.get(), **record} + _append(os.path.join(ERROR_DIR, f"errors-{_today()}.jsonl"), rec) + return rec + + +class RunAudit: + """Times a graded run and writes one audit line when the ``with`` block exits. + + Usage:: + + with audit.RunAudit(endpoint="grade", lat=..., lon=...) as run: + with run.phase("history"): + ... + run.set(run_type="full", history_rows=n) + return response + """ + + def __init__(self, **params): + self.run_id = uuid.uuid4().hex[:12] + self.fields = {"run_id": self.run_id, **params} + self.phases: dict[str, float] = {} + self._t0 = time.perf_counter() + self._token = None + + def __enter__(self) -> "RunAudit": + self._token = current_run_id.set(self.run_id) + return self + + def set(self, **kw) -> "RunAudit": + """Attach extra fields (e.g. run_type, history span) to the audit record.""" + self.fields.update(kw) + return self + + @contextlib.contextmanager + def phase(self, name: str): + """Time a named part of the run; the duration (ms) lands in ``phases``.""" + start = time.perf_counter() + try: + yield + finally: + self.phases[name] = round((time.perf_counter() - start) * 1000, 1) + + def __exit__(self, exc_type, exc, tb) -> bool: + total = round((time.perf_counter() - self._t0) * 1000, 1) + status = "ok" + if exc is not None: + status = "error" + self.fields["error"] = repr(exc) + code = getattr(exc, "status_code", None) + if code is not None: + self.fields["http_status"] = code + # Mirror whole-run failures into the separate errors folder. + log_event("error", { + "phase": "run", + "error": repr(exc), + **{k: self.fields.get(k) for k in ("endpoint", "lat", "lon", "cell_id")}, + }) + record = { + "ts": _now_iso(), + "status": status, + "duration_ms": total, + "phases": self.phases, + **self.fields, + } + _append(os.path.join(AUDIT_DIR, f"audit-{_today()}.jsonl"), record) + if self._token is not None: + current_run_id.reset(self._token) + return False # never suppress the exception diff --git a/climate.py b/climate.py new file mode 100644 index 0000000..c1266ab --- /dev/null +++ b/climate.py @@ -0,0 +1,527 @@ +"""Fetch historical + recent daily weather from Open-Meteo and cache to parquet. + +The full multi-decade daily record for a grid cell is fetched once and stored as +a zstd-compressed parquet file keyed by cell id. Percentiles are computed on the +fly from this raw record (see grading.py), which keeps the cache small and lets +us handle ties (e.g. many zero-precip days) correctly. +""" +import datetime +import os +import threading +import time + +import httpx +import numpy as np +import pandas as pd + +import audit + +CACHE_DIR = os.path.join(os.path.dirname(__file__), "..", "data", "cache") +MAX_ATTEMPTS = 3 # per upstream call, before giving up +START_DATE = "1980-01-01" # ERA5 reaches back to 1940; 1980 = 45 yrs, fast + robust +ARCHIVE_LATENCY_DAYS = 6 # ERA5 archive lags real time by a few days +# History rarely changes, so the full multi-decade archive is fetched once and +# cached indefinitely (refetched only to add new metric columns). Only the recent +# tail is refreshed — a small incremental fetch, at most hourly. +HISTORY_TOPUP_INTERVAL = 3600 # seconds between recent-tail refresh attempts +FORECAST_TTL_HOURS = 1 # refetch the forward forecast hourly to track updates + +ARCHIVE_URL = "https://archive-api.open-meteo.com/v1/archive" +FORECAST_URL = "https://api.open-meteo.com/v1/forecast" +# Backup archive when Open-Meteo is unavailable (e.g. its daily rate limit). NASA +# POWER is free + keyless, global, daily from 1981. It lacks gusts + apparent temp, +# so gusts read as unavailable and "feels like" is computed from heat index/chill. +NASA_POWER_URL = "https://power.larc.nasa.gov/api/temporal/daily/point" +NASA_START = "19810101" +NASA_FILL = -900.0 # POWER's missing-value sentinel is ~-999 + +DAILY_VARS = ( + "temperature_2m_max,temperature_2m_min,precipitation_sum," + "wind_speed_10m_max,wind_gusts_10m_max," + "apparent_temperature_max,apparent_temperature_min," + "relative_humidity_2m_mean" +) + +# Columns added after the original tmax/tmin/precip schema. A cached parquet +# missing any of these predates the wind/feels/humidity features (or the split +# apparent high/low) and is refetched so the new metrics show up immediately +# instead of after the 30-day cache TTL. +NEW_COLS = ("wind", "gust", "feels", "humid", "fmax", "fmin") + +# Thermoneutral baseline (°F). The combined "feels like" metric reports whichever +# apparent-temperature extreme — the heat-index-driven daily max or the +# wind-chill-driven daily min — sits further from this comfort point, so a single +# daily value captures both heat index (hot side) and wind chill (cold side). +COMFORT_F = 65.0 + + +def _cache_path(cell_id: str) -> str: + return os.path.join(CACHE_DIR, f"{cell_id}.parquet") + + +# One lock per cell so concurrent requests/refreshes don't each pull the full +# multi-decade archive (which quickly trips the upstream rate limit). +_LOCKS_GUARD = threading.Lock() +_CELL_LOCKS: dict[str, threading.Lock] = {} + +# When the archive rate-limits us (429), back off globally for a bit rather than +# re-hitting it on every refresh — that only prolongs the limit. +ARCHIVE_COOLDOWN = 120 # seconds +_archive_cooldown_until = 0.0 + + +_archive_limit_daily = False # is the current cooldown a daily-quota exhaustion? + + +def _is_rate_limit(e) -> bool: + return getattr(getattr(e, "response", None), "status_code", None) == 429 + + +def _rate_limit_reason(e) -> str: + """Upstream's human-readable 429 reason, if present in the response body.""" + resp = getattr(e, "response", None) + if resp is not None: + try: + j = resp.json() + if isinstance(j, dict) and j.get("reason"): + return str(j["reason"]) + except Exception: # noqa: BLE001 + pass + return "" + + +def _seconds_to_utc_reset() -> float: + """Seconds until just after the next UTC midnight (when daily quotas reset).""" + now = datetime.datetime.now(datetime.timezone.utc) + reset = (now + datetime.timedelta(days=1)).replace(hour=0, minute=10, second=0, microsecond=0) + return max((reset - now).total_seconds(), 600) + + +def _note_rate_limit(e) -> None: + """Record a rate limit: back off ~2 min, or until tomorrow for a daily quota.""" + global _archive_cooldown_until, _archive_limit_daily + reason = _rate_limit_reason(e).lower() + _archive_limit_daily = "daily" in reason or "tomorrow" in reason + _archive_cooldown_until = time.time() + ( + _seconds_to_utc_reset() if _archive_limit_daily else ARCHIVE_COOLDOWN) + + +def _cooldown_message() -> str: + if _archive_limit_daily: + return ("Open-Meteo's daily request limit is exhausted — new locations will work again " + "tomorrow. Places you've already viewed still work.") + return "The weather archive is rate-limited right now — please try again in a minute." + + +def _cell_lock(cell_id: str) -> threading.Lock: + with _LOCKS_GUARD: + lk = _CELL_LOCKS.get(cell_id) + if lk is None: + lk = _CELL_LOCKS[cell_id] = threading.Lock() + return lk + + +def _request(url, params, timeout, *, phase, headers=None, attempts=MAX_ATTEMPTS): + """GET with bounded retries; every retry and the final failure are logged to + the errors folder (tagged ``retry`` / ``error``). A 429 (rate limit) fails fast + without retrying, so we don't hammer the limit and make it worse.""" + last = None + for attempt in range(1, attempts + 1): + try: + r = httpx.get(url, params=params, timeout=timeout, headers=headers) + r.raise_for_status() + return r + except Exception as e: # noqa: BLE001 - upstream/network failures are expected + last = e + status = getattr(getattr(e, "response", None), "status_code", None) + rate_limited = status == 429 + final = rate_limited or attempt == attempts + audit.log_event( + "error" if final else "retry", + {"phase": phase, "attempt": attempt, "max_attempts": attempts, + "url": url, "status": status, "error": repr(e)}, + ) + if final: + break + time.sleep(min(0.5 * 2 ** (attempt - 1), 4.0)) + raise last + + +def _derive_humidity(df: pd.DataFrame) -> pd.DataFrame: + """Replace the raw mean *relative* humidity column with *absolute* humidity + (grams of water vapor per m³), in place. + + Absolute humidity is derived from the day's mean RH and its mean temperature + ((tmax+tmin)/2) via the Magnus saturation-vapor-pressure formula. It's a far + more informative "how muggy was it" signal than relative humidity, which mostly + tracks the day/night temperature swing (cold nights read ~100% RH regardless of + actual moisture). Applied at the read boundary so the parquet cache keeps the + raw RH the archive returns — no refetch needed for existing cells.""" + if "humid" not in df.columns: + return df + tmean_c = ((df["tmax"] + df["tmin"]) / 2.0 - 32.0) * 5.0 / 9.0 + rh = pd.to_numeric(df["humid"], errors="coerce") + es = 6.112 * np.exp(17.67 * tmean_c / (tmean_c + 243.5)) # sat. vapor pressure, hPa + df["humid"] = (es * rh * 2.1674 / (273.15 + tmean_c)).round(1) # absolute humidity, g/m³ + return df + + +def _combined_feels(amax: pd.Series, amin: pd.Series) -> pd.Series: + """One daily "feels like" value: the apparent-temperature extreme furthest from + the comfort baseline — the heat-index high on warm days, the wind-chill low on + cold ones. Falls back to whichever side is present if one is missing.""" + hot = amax - COMFORT_F + cold = COMFORT_F - amin + feels = amax.where(hot >= cold, amin) # NaN comparisons pick the min side + return feels.fillna(amax).fillna(amin) + + +def _to_frame(daily: dict) -> pd.DataFrame: + n = len(daily["time"]) + # Older upstream responses (or a narrowed variable set) may omit a series; fall + # back to an all-null column of the right length so the frame shape is stable. + def col(key): + vals = daily.get(key) + return vals if vals else [None] * n + + df = pd.DataFrame( + { + "date": pd.to_datetime(daily["time"]), + "tmax": daily["temperature_2m_max"], + "tmin": daily["temperature_2m_min"], + "precip": daily["precipitation_sum"], + "wind": col("wind_speed_10m_max"), + "gust": col("wind_gusts_10m_max"), + "humid": col("relative_humidity_2m_mean"), + } + ) + # The apparent (felt) high and low are kept as their own columns — graded + # independently — alongside the combined `feels` (whichever side is further + # from the fixed comfort baseline), which the weekly/day views still use. + amax = pd.to_numeric(pd.Series(col("apparent_temperature_max")), errors="coerce") + amin = pd.to_numeric(pd.Series(col("apparent_temperature_min")), errors="coerce") + df["fmax"] = amax + df["fmin"] = amin + df["feels"] = _combined_feels(amax, amin) + # Need a valid high/low to be a usable climate day; the rest may be NaN and are + # simply ungraded for that day (the percentile helpers return None). + df = df.dropna(subset=["tmax", "tmin"]).reset_index(drop=True) + df["doy"] = df["date"].dt.dayofyear.astype("int16") + return df + + +def _fetch_history(cell: dict) -> pd.DataFrame: + end = (datetime.date.today() - datetime.timedelta(days=ARCHIVE_LATENCY_DAYS)).isoformat() + params = { + "latitude": cell["center_lat"], + "longitude": cell["center_lon"], + "start_date": START_DATE, + "end_date": end, + "daily": DAILY_VARS, + "timezone": "auto", + "temperature_unit": "fahrenheit", + "precipitation_unit": "inch", + "wind_speed_unit": "mph", + } + r = _request(ARCHIVE_URL, params, 180, phase="history_fetch") + return _to_frame(r.json()["daily"]) + + +def _heat_index(t_f, rh): + """NWS heat index (°F); below ~80°F it's just the air temperature.""" + T, R = t_f, rh + hi = (-42.379 + 2.04901523 * T + 10.14333127 * R - 0.22475541 * T * R + - 0.00683783 * T * T - 0.05481717 * R * R + 0.00122874 * T * T * R + + 0.00085282 * T * R * R - 0.00000199 * T * T * R * R) + return np.where(np.asarray(T, dtype="float64") >= 80, hi, T) + + +def _wind_chill(t_f, v_mph): + """NWS wind chill (°F); applies only when cold + breezy, else the air temp.""" + T = np.asarray(t_f, dtype="float64") + V = np.clip(np.asarray(v_mph, dtype="float64"), 0, None) + Vp = np.power(V, 0.16) + wc = 35.74 + 0.6215 * T - 35.75 * Vp + 0.4275 * T * Vp + return np.where((T <= 50) & (V >= 3), wc, T) + + +def _nasa_to_frame(param: dict) -> pd.DataFrame: + """Map a NASA POWER daily response to our (tmax/tmin/precip/wind/gust/humid/feels) + schema, converting units (°C→°F, mm→in, m/s→mph) and computing feels-like.""" + dates = sorted(param.get("T2M_MAX", {}).keys()) + + def col(key, transform): + d = param.get(key, {}) + return [np.nan if (v is None or v <= NASA_FILL) else transform(v) + for v in (d.get(k) for k in dates)] + + c2f = lambda c: c * 9.0 / 5.0 + 32.0 + df = pd.DataFrame({ + "date": pd.to_datetime(dates, format="%Y%m%d"), + "tmax": col("T2M_MAX", c2f), + "tmin": col("T2M_MIN", c2f), + "precip": col("PRECTOTCORR", lambda mm: mm / 25.4), + "wind": col("WS10M_MAX", lambda ms: ms * 2.2369362920544), + "gust": [np.nan] * len(dates), # POWER has no gusts + "humid": col("RH2M", lambda v: v), + }) + # POWER has no apparent temperature; approximate the felt high with the NWS + # heat index and the felt low with NWS wind chill (each falls back to the air + # temperature outside its regime), mirroring the Open-Meteo columns. + df["fmax"] = pd.Series(_heat_index(df["tmax"], df["humid"]), index=df.index) + df["fmin"] = pd.Series(_wind_chill(df["tmin"], df["wind"]), index=df.index) + df["feels"] = _combined_feels(df["fmax"], df["fmin"]) + df = df.dropna(subset=["tmax", "tmin"]).reset_index(drop=True) + df["doy"] = df["date"].dt.dayofyear.astype("int16") + return df + + +def _fetch_history_nasa(cell: dict) -> pd.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") + 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, + "format": "JSON", + } + r = _request(NASA_POWER_URL, params, 180, phase="history_nasa") + return _nasa_to_frame(r.json()["properties"]["parameter"]) + + +def _fetch_history_range(cell: dict, start_date: str, end_date: str) -> pd.DataFrame: + """Fetch just a date range of archive history (used to top up the recent tail).""" + params = { + "latitude": cell["center_lat"], + "longitude": cell["center_lon"], + "start_date": start_date, + "end_date": end_date, + "daily": DAILY_VARS, + "timezone": "auto", + "temperature_unit": "fahrenheit", + "precipitation_unit": "inch", + "wind_speed_unit": "mph", + } + r = _request(ARCHIVE_URL, params, 60, phase="history_topup") + return _to_frame(r.json()["daily"]) + + +def _read_history_cache(path): + """Schema-complete cached frame (no doy) + its file age in seconds, or None. + + No age expiry — history is cached indefinitely; a pre-wind/humidity schema is + the only reason to refetch (to add the new metric columns).""" + if not os.path.exists(path): + return None + df = pd.read_parquet(path) + if not all(c in df.columns for c in NEW_COLS): + return None + return df, time.time() - os.path.getmtime(path) + + +def _topup_tail(cell: dict, df: pd.DataFrame, path: str) -> pd.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 + with _cell_lock(cell["id"]): + fresh = _read_history_cache(path) # another thread may have just refreshed it + if fresh is not None: + 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 = pd.Timestamp(df["date"].max()).date() + if cached_max >= expected: + try: os.utime(path, None) # tail already current — reset the hourly timer + except OSError: pass + return df + try: + recent = _fetch_history_range( + 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): + _note_rate_limit(e) + return df + merged = (pd.concat([df, recent.drop(columns=["doy"], errors="ignore")]) + .drop_duplicates(subset="date", keep="last") + .sort_values("date").reset_index(drop=True)) + merged.to_parquet(path, compression="zstd", index=False) + return merged + + +def get_history(cell: dict) -> tuple[pd.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) + _derive_humidity(df) + return df, meta + + +def _load_history(cell: dict) -> tuple[pd.DataFrame, dict]: + """Return (daily history frame, cache metadata) for a cell. + + The full archive is cached indefinitely (fetched once); only the recent tail is + topped up, at most hourly. Concurrent callers are serialized so only one archive + fetch happens; on an upstream failure (e.g. 429) any existing cache is served.""" + path = _cache_path(cell["id"]) + + hit = _read_history_cache(path) + if hit is not None: + df, age_s = hit + if age_s > HISTORY_TOPUP_INTERVAL: + df = _topup_tail(cell, df, path) # refresh just the recent days + df = df.copy() + df["doy"] = df["date"].dt.dayofyear.astype("int16") + return df, {"cached": True, "cache_age_days": round(age_s / 86400.0, 1)} + + global _archive_cooldown_until + + with _cell_lock(cell["id"]): + hit = _read_history_cache(path) # another thread may have populated it while we waited + if hit is not None: + df, age_s = hit + df["doy"] = df["date"].dt.dayofyear.astype("int16") + return df, {"cached": True, "cache_age_days": round(age_s / 86400.0, 1)} + + stale = pd.read_parquet(path) if os.path.exists(path) else None + + def _serve_stale(): + stale["doy"] = stale["date"].dt.dayofyear.astype("int16") + return 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. + df = None + source = None + if time.time() >= _archive_cooldown_until: + try: + df = _fetch_history(cell) + 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() + raise RuntimeError(_cooldown_message()) + + os.makedirs(CACHE_DIR, exist_ok=True) + # Store the raw record; percentiles are derived at request time. + df.drop(columns=["doy"]).to_parquet(path, compression="zstd", index=False) + return df, {"cached": False, "cache_age_days": 0, "source": source} + + +RECENT_PAST_DAYS = 25 # recent observations window (covers the ~2-week graded view) +FORECAST_DAYS = 8 # today + 7 days ahead + + +def _rf_cache_path(cell_id: str) -> str: + return os.path.join(CACHE_DIR, f"{cell_id}_rf.parquet") + + +def get_recent_forecast(cell: dict) -> pd.DataFrame: + """Recent observations + forward forecast, with humidity as absolute humidity + (g/m³). Thin wrapper over the raw loader (see below).""" + df = _load_recent_forecast(cell) + _derive_humidity(df) + return df + + +def _load_recent_forecast(cell: dict) -> pd.DataFrame: + """Recent observations AND the forward forecast in ONE forecast-API call. + + Both the recent (past) view and the forecast (future) view slice from this, so + a cell needs just one upstream forecast request per hour (plus the ~monthly + archive fetch). Cached per cell for one hour so it still follows model updates. + """ + path = _rf_cache_path(cell["id"]) + if os.path.exists(path): + age_h = (time.time() - os.path.getmtime(path)) / 3600.0 + if age_h < FORECAST_TTL_HOURS: + df = pd.read_parquet(path) + df["doy"] = df["date"].dt.dayofyear.astype("int16") + return df + + params = { + "latitude": cell["center_lat"], + "longitude": cell["center_lon"], + "daily": DAILY_VARS, + "timezone": "auto", + "temperature_unit": "fahrenheit", + "precipitation_unit": "inch", + "wind_speed_unit": "mph", + "past_days": RECENT_PAST_DAYS, + "forecast_days": FORECAST_DAYS, + } + r = _request(FORECAST_URL, params, 60, phase="recent_forecast_fetch") + df = _to_frame(r.json()["daily"]) + os.makedirs(CACHE_DIR, exist_ok=True) + df.drop(columns=["doy"]).to_parquet(path, compression="zstd", index=False) + return df + + +_REVGEO_CACHE: dict[tuple[float, float], str | None] = {} + + +def reverse_geocode(lat: float, lon: float) -> str | None: + """Best-effort "City, State" label for a point (OpenStreetMap Nominatim). + + Cached per ~cell so panning around doesn't hammer the service, and failures + return None so the caller can fall back to bare coordinates. + """ + key = (round(lat, 3), round(lon, 3)) + if key in _REVGEO_CACHE: + return _REVGEO_CACHE[key] + label = None + try: + r = _request( + "https://nominatim.openstreetmap.org/reverse", + {"lat": lat, "lon": lon, "format": "jsonv2", "zoom": 10, + "addressdetails": 1}, + 15, + phase="reverse_geocode", + headers={"User-Agent": "Thermograph/0.1 (local weather grading app)"}, + ) + a = r.json().get("address", {}) or {} + city = (a.get("city") or a.get("town") or a.get("village") + or a.get("hamlet") or a.get("suburb") or a.get("county")) + region = a.get("state") or a.get("province") or a.get("region") + label = ", ".join(p for p in (city, region) if p) or None + except Exception: # noqa: BLE001 - reverse geocoding is a nicety, never fatal + label = None + _REVGEO_CACHE[key] = label + return label + + +def geocode(name: str, count: int = 5) -> list[dict]: + """Look up places by name (US/Canada biased) 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"), + } + for g in results + ] diff --git a/grading.py b/grading.py new file mode 100644 index 0000000..7d9f57f --- /dev/null +++ b/grading.py @@ -0,0 +1,369 @@ +"""Turn a raw daily-weather record into day-of-year climatology and letter grades. + +For a given day of the year we build the reference distribution from every +historical day whose day-of-year is within +/- `HALF_WINDOW` days of it (wrapping +around the year end). An observed value is then placed on that distribution as an +empirical percentile and mapped to a human-readable grade. +""" +import numpy as np +import pandas as pd + +HALF_WINDOW = 7 # +/- 7 days -> a 15-day seasonal window +RAIN_THRESHOLD = 0.01 # inches; a day with < this much precip counts as "dry" +DSR_LOOKBACK_MIN_DAYS = 14 # min history required before a window to seed dry streaks + +# Metrics graded on the diverging temperature scale (percentile -> TEMP_BANDS): +# the two air temps, the combined "feels like" (heat index / wind chill) plus its +# separate apparent high/low sides (fmax/fmin — the calendar recombines them +# client-side around a user-chosen comfort temperature), wind speed and wind +# gust. Each is a single daily scalar placed on its own ±7-day historical +# distribution, so calm/low reads as the cool (blue) end and windy/high as the +# hot (red) end — same coloring + categories as temperature. +TEMP_METRICS = ("tmax", "tmin", "feels", "fmax", "fmin", "humid", "wind", "gust") +# All metrics that get a percentile summary (temperature-like + precipitation). +CLIMO_METRICS = TEMP_METRICS + ("precip",) + +# Percentile -> (grade label, css class) for temperature. Higher percentile = warmer. +# 9 symmetric tiers around the middle 40-60% "Normal", escalating to "Near Record" +# at both edges. Each entry's number is the tier's LOWER bound; a tier spans +# [lower, next-lower) — i.e. lower-inclusive, upper-exclusive (a p90 day is "Very +# High", not "High"). Boundaries sit on multiples of 5 (plus the 1/99 record edges). +TEMP_BANDS = [ + (99, "Near Record", "rec-hot"), # >=99 extreme high (danger) + (90, "Very High", "very-hot"), # 90-99 + (75, "High", "hot"), # 75-90 + (60, "Above Normal", "warm"), # 60-75 + (40, "Normal", "normal"), # 40-60 (the middle) + (25, "Below Normal", "cool"), # 25-40 + (10, "Low", "cold"), # 10-25 + (1, "Very Low", "very-cold"), # 1-10 + (0, "Near Record", "rec-cold"), # <1 extreme low (danger) +] + +# Precipitation is graded ONLY among days that actually rained (>= RAIN_THRESHOLD) +# in the seasonal window — a "rain percentile". Rain is one-directional (heavier = +# more extreme), so these 9 tiers are sequential light->heavy, using the SAME cut +# points as temperature. Dry days are handled separately (the "dry" class, colored +# by dry streak in the UI). Same [lower, upper) convention as TEMP_BANDS. +RAIN_BANDS = [ + (99, "Extreme", "wet-9"), # heaviest rain for the season (darkest) + (90, "Very Heavy", "wet-8"), # 90-99 + (75, "Heavy", "wet-7"), # 75-90 + (60, "Mod–Heavy", "wet-6"), # 60-75 + (40, "Moderate", "wet-5"), # 40-60 + (25, "Light–Mod", "wet-4"), # 25-40 + (10, "Light", "wet-3"), # 10-25 + (1, "Very Light", "wet-2"), # 1-10 + (0, "Trace", "wet-1"), # <1 lightest measurable rain +] + + +# Tier ladders for the single-day detail view: each tier and the two percentile +# marks that bound it. Kept aligned with TEMP_BANDS / RAIN_BANDS (same labels and +# css classes) so the UI names and colors stay consistent across every view. +# (class, label, percentile-range label, lower-bound pct, upper-bound pct) +_TEMP_LADDER = [ + ("rec-hot", "Near Record", ">99%", 99, None), + ("very-hot", "Very High", "90–99%", 90, 99), + ("hot", "High", "75–90%", 75, 90), + ("warm", "Above Normal", "60–75%", 60, 75), + ("normal", "Normal", "40–60%", 40, 60), + ("cool", "Below Normal", "25–40%", 25, 40), + ("cold", "Low", "10–25%", 10, 25), + ("very-cold", "Very Low", "1–10%", 1, 10), + ("rec-cold", "Near Record", "<1%", None, 1), +] +# Rain tiers are on the rain-day-only percentile scale (see _grade_precip). The +# lightest tier's lower bound is the smallest measured rain day (marked with 0). +_RAIN_LADDER = [ + ("wet-9", "Extreme", ">99%", 99, None), + ("wet-8", "Very Heavy", "90–99%", 90, 99), + ("wet-7", "Heavy", "75–90%", 75, 90), + ("wet-6", "Mod–Heavy", "60–75%", 60, 75), + ("wet-5", "Moderate", "40–60%", 40, 60), + ("wet-4", "Light–Mod", "25–40%", 25, 40), + ("wet-3", "Light", "10–25%", 10, 25), + ("wet-2", "Very Light", "1–10%", 1, 10), + ("wet-1", "Trace", "<1%", 0, 1), +] + + +def _band(pct: float, bands) -> tuple[str, str]: + # The top tier is strict (pct > its threshold): "Near Record" high means + # strictly beyond the 99th percentile — the top <1% — mirroring the + # strictly-below-1st bottom tier (its lower neighbor already catches pct >= 1). + # Everything in between stays lower-inclusive, upper-exclusive. + top_thr, top_label, top_css = bands[0] + if pct > top_thr: + return top_label, top_css + for threshold, label, css in bands[1:]: + if pct >= threshold: + return label, css + return bands[-1][1], bands[-1][2] + + +def window_mask(doys: np.ndarray, target_doy: int, half: int = HALF_WINDOW) -> np.ndarray: + diff = np.abs(doys.astype(int) - int(target_doy)) + circular = np.minimum(diff, 366 - diff) + return circular <= half + + +def empirical_percentile(samples: np.ndarray, value) -> float | None: + """Mid-rank percentile of `value` within `samples` (handles ties correctly).""" + n = samples.size + if n == 0 or value is None or (isinstance(value, float) and np.isnan(value)): + return None + less = int(np.sum(samples < value)) + equal = int(np.sum(samples == value)) + return round(100.0 * (less + 0.5 * equal) / n, 1) + + +def climatology(df: pd.DataFrame, target_doy: int) -> dict: + """Summarize the +/-7 day historical distribution for one day of the year.""" + doys = df["doy"].values + sub = df[window_mask(doys, target_doy)] + years = pd.to_datetime(sub["date"]).dt.year + out = { + "target_doy": int(target_doy), + "n_samples": int(len(sub)), + "year_range": [int(years.min()), int(years.max())] if len(sub) else None, + } + for var in CLIMO_METRICS: + if var not in sub.columns: + out[var] = None + continue + v = sub[var].dropna().values + if v.size == 0: + out[var] = None + continue + p10, p40, p50, p60, p90 = np.percentile(v, [10, 40, 50, 60, 90]) + out[var] = { + "min": round(float(np.min(v)), 2), + "p10": round(float(p10), 1), + "p40": round(float(p40), 1), # Normal band is now 40-60 (see TEMP_BANDS) + "p50": round(float(p50), 1), + "p60": round(float(p60), 1), + "p90": round(float(p90), 1), + "max": round(float(np.max(v)), 2), + "mean": round(float(np.mean(v)), 1), + } + return out + + +def _band_stats(samples: np.ndarray) -> dict | None: + if samples.size == 0: + return None + # Percentiles for the chart's nested "normal" fan, matching the 9 tier bounds: + # p40-p60 is the Normal band; p25/p75, p10/p90 and p1/p99 mark the successive + # Below/Above Normal, Low/High, Very Low/High and Near-Record edges. + p1, p10, p25, p40, p50, p60, p75, p90, p99 = np.percentile( + samples, [1, 10, 25, 40, 50, 60, 75, 90, 99] + ) + return { + "p1": round(float(p1), 1), + "p10": round(float(p10), 1), + "p25": round(float(p25), 1), + "p40": round(float(p40), 1), + "p50": round(float(p50), 1), + "p60": round(float(p60), 1), + "p75": round(float(p75), 1), + "p90": round(float(p90), 1), + "p99": round(float(p99), 1), + } + + +def _grade_value(samples: np.ndarray, value, bands) -> dict | None: + """Grade a temperature value by its empirical percentile in the window.""" + pct = empirical_percentile(samples, value) + if pct is None: + return None + label, css = _band(pct, bands) + return {"value": round(float(value), 2), "percentile": pct, "grade": label, "class": css} + + +def _grade_precip(samples: np.ndarray, value) -> dict | None: + """Grade precipitation by its "rain percentile" — the rank of the day's rainfall + among *rain days only* (>= RAIN_THRESHOLD) in the window. Dry days get the "dry" + class with no percentile (the UI colors them by dry streak instead).""" + if value is None or (isinstance(value, float) and np.isnan(value)): + return None + value = float(value) + if value < RAIN_THRESHOLD: + return {"value": round(value, 2), "percentile": None, "grade": "Dry", "class": "dry"} + rain = samples[samples >= RAIN_THRESHOLD] + pct = empirical_percentile(rain, value) + if pct is None: # no historical rain days in window (extremely rare) + pct = 100.0 + label, css = _band(pct, RAIN_BANDS) + return {"value": round(value, 2), "percentile": pct, "grade": label, "class": css} + + +def dry_streaks(dates, precips) -> dict[str, int]: + """Map each date (ISO string) to days since the last measurable rain, walking a + chronological precip series. NaN precip counts as a dry day (compares False).""" + out: dict[str, int] = {} + streak = 0 + for d, p in zip(dates, precips): + wet = p is not None and not (isinstance(p, float) and np.isnan(p)) and p >= RAIN_THRESHOLD + streak = 0 if wet else streak + 1 + out[pd.Timestamp(d).date().isoformat()] = streak + return out + + +def grade_range(df: pd.DataFrame, start, end) -> list[dict]: + """Grade every historical day in [start, end] against its own ±7-day window. + + Powers the calendar view. Reuses one window per day-of-year across the whole + range, so a 2-year span computes at most ~366 windows (not one per day). Each + day is returned in a compact shape: value (v), percentile (pct), css class (c), + grade label (g). + """ + start, end = pd.Timestamp(start), pd.Timestamp(end) + full = df.sort_values("date") + doys_all = full["doy"].values + cols = {v: full[v].values for v in CLIMO_METRICS if v in full.columns} + cache: dict[tuple[int, str], np.ndarray] = {} + + def samples(doy: int, var: str) -> np.ndarray: + key = (doy, var) + if key not in cache: + arr = cols[var][window_mask(doys_all, doy)] + cache[key] = arr[~np.isnan(arr)] + return cache[key] + + # Days since last measurable rain. Computed over the ENTIRE record that precedes + # the window — not just the window, nor a fixed N-day buffer — so the streak on + # the first shown day is exact even when a dry spell straddles the window start. + # A fixed 14-day lookback would be the bare minimum but still undercounts longer + # droughts (the record has 25-day dry streaks); using the full history (already + # cached per cell) is both correct for any streak length and free. NaN precip + # compares False against the threshold, so it counts as a dry day. + if full["date"].min() > start - pd.Timedelta(days=DSR_LOOKBACK_MIN_DAYS): + # Should never happen (history starts in 1980); guards against a future + # change that trims history and would silently truncate streaks. + import warnings + warnings.warn("dry-streak lookback shorter than the 14-day minimum", stacklevel=2) + dsr_map: dict[str, int] = {} + streak = 0 + for d, p in zip(full["date"].values, cols["precip"]): + streak = 0 if p >= RAIN_THRESHOLD else streak + 1 + dsr_map[pd.Timestamp(d).date().isoformat()] = streak + + sub = full[(full["date"] >= start) & (full["date"] <= end)] + out = [] + for _, row in sub.iterrows(): + doy = int(row["doy"]) + date = pd.Timestamp(row["date"]).date().isoformat() + rec = {"date": date, "dsr": dsr_map.get(date)} + for var in TEMP_METRICS: + if var not in cols: + rec[var] = None + continue + g = _grade_value(samples(doy, var), row[var], TEMP_BANDS) + rec[var] = {"v": g["value"], "pct": g["percentile"], "c": g["class"], "g": g["grade"]} if g else None + gp = _grade_precip(samples(doy, "precip"), row["precip"]) + rec["precip"] = {"v": gp["value"], "pct": gp["percentile"], "c": gp["class"], "g": gp["grade"]} if gp else None + out.append(rec) + return out + + +def _temp_ladder(samples: np.ndarray) -> dict | None: + """Value at each temperature tier boundary within the ±7-day window.""" + if samples.size == 0: + return None + marks = {m: round(float(np.percentile(samples, m)), 1) for m in (1, 10, 25, 40, 50, 60, 75, 90, 99)} + tiers = [ + {"c": c, "label": label, "range": rng, + "lo": marks[lo] if lo is not None else None, + "hi": marks[hi] if hi is not None else None} + for c, label, rng, lo, hi in _TEMP_LADDER + ] + return {"tiers": tiers, "median": marks[50], + "min": round(float(samples.min()), 1), "max": round(float(samples.max()), 1)} + + +def _precip_ladder(samples: np.ndarray) -> dict | None: + """Value at each rain-day tier boundary, plus how often the window is dry.""" + if samples.size == 0: + return None + rain = samples[samples >= RAIN_THRESHOLD] + n = int(samples.size) + tiers = [] + if rain.size: + marks = {m: round(float(np.percentile(rain, m)), 2) for m in (1, 10, 25, 40, 60, 75, 90, 99)} + rmin = round(float(rain.min()), 2) + tiers = [ + {"c": c, "label": label, "range": rng, + "lo": rmin if lo == 0 else marks[lo], + "hi": marks[hi] if hi is not None else None} + for c, label, rng, lo, hi in _RAIN_LADDER + ] + tiers.append({"c": "dry", "label": "Dry", "range": f"< {RAIN_THRESHOLD}\"", "lo": 0.0, "hi": None}) + return {"tiers": tiers, "dry_pct": round(100.0 * (n - rain.size) / n, 1), + "rain_days": int(rain.size), + "min": round(float(samples.min()), 2), "max": round(float(samples.max()), 2)} + + +def day_detail(df: pd.DataFrame, date: pd.Timestamp, obs: dict | None) -> dict: + """Full percentile breakdown for one day: the value at every tier boundary in + its own ±7-day window, plus where the observed values (if any) land. + + Powers the single-day detail page. `obs` may be None when the date isn't yet + in the record — then only the climatological ladders are returned.""" + ts = pd.Timestamp(date) + doy = int(ts.dayofyear) + sub = df[window_mask(df["doy"].values, doy)] + years = pd.to_datetime(sub["date"]).dt.year + + metrics = {} + for var in TEMP_METRICS: + if var not in sub.columns: + continue + vals = sub[var].dropna().values + metrics[var] = { + "ladder": _temp_ladder(vals), + "obs": _grade_value(vals, obs.get(var) if obs else None, TEMP_BANDS), + } + metrics["precip"] = { + "ladder": _precip_ladder(sub["precip"].dropna().values), + "obs": _grade_precip(sub["precip"].dropna().values, obs.get("precip") if obs else None), + } + return { + "date": ts.date().isoformat(), + "doy": doy, + "n_samples": int(len(sub)), + "year_range": [int(years.min()), int(years.max())] if len(sub) else None, + "metrics": metrics, + } + + +def grade_day(df: pd.DataFrame, date: pd.Timestamp, obs: dict) -> dict: + """Grade one observed day against its own day-of-year +/-7 window.""" + doy = int(pd.Timestamp(date).dayofyear) + doys = df["doy"].values + sub = df[window_mask(doys, doy)] + + result = {"date": pd.Timestamp(date).date().isoformat(), "doy": doy} + for var in TEMP_METRICS: + result[var] = ( + _grade_value(sub[var].dropna().values, obs.get(var), TEMP_BANDS) + if var in sub.columns else None + ) + result["precip"] = _grade_precip(sub["precip"].dropna().values, obs.get("precip")) + + # Per-day normal band (this day-of-year's ±7 window) so the trend chart can + # draw the "normal" envelope each actual value is compared against. + result["normals"] = { + var: _band_stats(sub[var].dropna().values) + for var in CLIMO_METRICS if var in sub.columns + } + + # A single "departure" score: how far the day strayed from the median (50th pct), + # taking the most extreme of high/low. 0 = perfectly normal, 50 = record extreme. + departures = [ + abs(result[v]["percentile"] - 50) + for v in ("tmax", "tmin") + if result[v] is not None + ] + result["departure"] = round(max(departures), 1) if departures else None + return result diff --git a/grid.py b/grid.py new file mode 100644 index 0000000..4291a67 --- /dev/null +++ b/grid.py @@ -0,0 +1,51 @@ +"""Snap an arbitrary lat/lon to a stable ~4-square-mile grid cell. + +The grid is defined by fixed latitude rows (~2 miles tall). Within each row the +longitude step is scaled by cos(latitude) so cells stay roughly square (~4 sq mi) +at every latitude instead of getting skinny toward the poles. Cell ids are +deterministic, so the same physical location always maps to the same cache file. +""" +import math + +# 1 degree of latitude ~= 69 miles. ~2 miles -> ~0.029 deg gives a ~4 sq mi cell. +LAT_STEP = 1.0 / 34.5 # ~= 0.02899 deg (~2.0 miles) + + +def _lon_step(center_lat: float) -> float: + """Longitude degrees that span ~2 miles at the given latitude.""" + c = math.cos(math.radians(center_lat)) + c = max(c, 0.05) # clamp near the poles to avoid a blow-up + return LAT_STEP / c + + +def snap(lat: float, lon: float) -> dict: + """Return the grid cell (id + center + span) containing (lat, lon).""" + i = math.floor(lat / LAT_STEP) + center_lat = (i + 0.5) * LAT_STEP + lon_step = _lon_step(center_lat) + j = math.floor(lon / lon_step) + center_lon = (j + 0.5) * lon_step + + # Approximate cell dimensions in miles for display. + height_mi = LAT_STEP * 69.0 + width_mi = lon_step * 69.0 * math.cos(math.radians(center_lat)) + + return { + "id": f"{i}_{j}", + "center_lat": round(center_lat, 5), + "center_lon": round(center_lon, 5), + "lat_step": LAT_STEP, + "lon_step": lon_step, + "bounds": { + "south": round(i * LAT_STEP, 5), + "north": round((i + 1) * LAT_STEP, 5), + "west": round(j * lon_step, 5), + "east": round((j + 1) * lon_step, 5), + }, + "area_sq_mi": round(height_mi * width_mi, 2), + } + + +def in_north_america(lat: float, lon: float) -> bool: + """Rough bounding box for the US (incl. Alaska/Hawaii) and Canada.""" + return 14.0 <= lat <= 84.0 and -172.0 <= lon <= -52.0 diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..022d732 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,6 @@ +fastapi==0.115.6 +uvicorn[standard]==0.34.0 +httpx==0.28.1 +pandas==2.2.3 +pyarrow==18.1.0 +numpy==2.2.1