thermograph/grading.py

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