Jinemi/thermograph
4
0
Fork 0
Code Issues 2 Pull requests 9 Projects Releases Packages Wiki Activity Actions 1

Initial commit: app +

Browse source 
VPS deploy pipeline
This commit is contained in:
Emi Griffith 2026-07-10 17:29:47 -07:00
commit 8d90dcb456
6 changed files with 1527 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

444
app.py Normal file
View file

@ -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://<host>/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")

130
audit.py Normal file
View file

@ -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-<date>.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-<date>.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

527
climate.py Normal file
View file

@ -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
]

369
grading.py Normal file
View file

@ -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

51
grid.py Normal file
View file

@ -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

6
requirements.txt Normal file
View file

@ -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