climatology.py

"""
climatology.py
==============
Per-zone day-of-year climatology and anomaly (z-score) computation.

 This module turns absolute real-world readings into the z-score
 anomalies the rest of the pipeline was designed around.

DATA SOURCES (two tiers, matching the codebase's degrade-safely philosophy)
----------------------------------------------------------------------------
  1. Real: Open-Meteo historical archive API (free, no API key -- the same
     endpoint era5_data_pipeline._fetch_openmeteo already uses). Daily
     precipitation_sum + temperature_2m_mean are well-established archive
     variables. soil_moisture_0_to_7cm_mean is requested as documented in
     the Open-Meteo archive docs at write time; if the API rejects it or
     returns nothing, soil climatology degrades to the precip-tracked
     model below, logged -- never silently zeroed.
  2. Synthetic: a deterministic, latitude-aware maritime-continent monsoon
     model (see _synthetic_climatology). It is a HEURISTIC, calibrated to
     the broad shape of the Indonesian wet/dry season (SH monsoon: wet
     Dec-Mar, dry Jun-Sep; weaker/bimodal near the equator; shifted peak
     for northern Sumatra). It exists so the pipeline keeps producing
     sensible anomalies offline and in tests. It is NOT a retrieval --
     treat its absolute values as plausible shapes, not measurements.
"""

from __future__ import annotations

import json
import logging
import math
import os
from dataclasses import dataclass, field
from datetime import datetime, timedelta, timezone
from pathlib import Path
from typing import Any, Dict, List, Optional, Tuple

import zone_observation as _zo

assert _zo.SCHEMA_VERSION == 3, (
    f"climatology: zone_observation schema mismatch "
    f"(expected 3, got {_zo.SCHEMA_VERSION})"
)

from zone_observation import DataSource, ZoneObs, _clip, _stable_seed

logger = logging.getLogger(__name__)

try:
    import requests
    _REQUESTS_AVAILABLE = True
except ImportError:
    _REQUESTS_AVAILABLE = False
    logger.info("requests not installed -- climatology will use the synthetic model")


# ---------------------------------------------------------------------------
# Constants
# ---------------------------------------------------------------------------

_OPENMETEO_ARCHIVE_URL = "https://archive-api.open-meteo.com/v1/archive"
_TIMEOUT_S = int(os.environ.get("WEATHER_HTTP_TIMEOUT", "60"))

_CACHE_DIR = Path(os.environ.get("WEATHER_CACHE_DIR", ".cache/era5")) / "climatology"
_CACHE_DIR.mkdir(parents=True, exist_ok=True)
_CACHE_TTL_DAYS = 90  # climatology drifts slowly; quarterly refresh is ample

_DAYS_PER_YEAR = 365.25
_TABLE_LEN = 366  # DOY table indexed doy-1; DOY 60 = Feb 29 (leap mapping below)

# Std floors: prevent division blow-ups in convectively uniform seasons.
_PRECIP_STD_FLOOR = 1.5   # mm/day
_TEMP_STD_FLOOR   = 0.4   # deg C
_SOIL_STD_FLOOR   = 2.0   # percent
_RH_STD_FLOOR     = 3.0   # percent -- RH is bounded [0,100] and often near-
                          # saturated in the tropics, so day-to-day variance
                          # is naturally small; floor prevents z-score blowup
_PRECIP_WINDOW_DAYS = 30


# ---------------------------------------------------------------------------
# Day-of-year helpers (fixed 366-entry table regardless of leap years)
# ---------------------------------------------------------------------------

def _is_leap(year: int) -> bool:
    return year % 4 == 0 and (year % 100 != 0 or year % 400 == 0)


def doy_index(dt: datetime) -> int:
    doy = dt.timetuple().tm_yday
    if not _is_leap(dt.year) and doy >= 60:
        doy += 1
    return min(doy, _TABLE_LEN)


def _circular_smooth(values: List[float], half_window: int = 7) -> List[float]:
    n = len(values)
    out = []
    for i in range(n):
        acc = 0.0
        cnt = 0
        for j in range(-half_window, half_window + 1):
            acc += values[(i + j) % n]
            cnt += 1
        out.append(acc / cnt)
    return out


# ---------------------------------------------------------------------------
# ZoneClimatology
# ---------------------------------------------------------------------------

@dataclass
class ZoneClimatology:
    zone_id: str
    source: str                       # 'openmeteo_archive' | 'synthetic_model' | 'mixed'
    n_years: int
    period_start_year: int
    period_end_year: int
    precip_mean_mm: List[float] = field(default_factory=list)   # mm/day
    precip_std_mm:  List[float] = field(default_factory=list)
    temp_mean_c:    List[float] = field(default_factory=list)
    temp_std_c:     List[float] = field(default_factory=list)
    soil_mean_pct:  List[float] = field(default_factory=list)
    soil_std_pct:   List[float] = field(default_factory=list)
    rh_mean_pct:    List[float] = field(default_factory=list)
    rh_std_pct:     List[float] = field(default_factory=list)

    def __post_init__(self) -> None:
        for name in ("precip_mean_mm", "precip_std_mm", "temp_mean_c",
                     "temp_std_c", "soil_mean_pct", "soil_std_pct",
                     "rh_mean_pct", "rh_std_pct"):
            v = getattr(self, name)
            if len(v) != _TABLE_LEN:
                raise ValueError(
                    f"ZoneClimatology('{self.zone_id}'): {name} has length "
                    f"{len(v)}, expected {_TABLE_LEN}"
                )

    # --- Window statistics ------------------------------------------------
    def window_precip_stats(self, dt: datetime, window_days: int) -> Tuple[float, float]:
        idx0 = doy_index(dt) - 1
        mean = 0.0
        var = 0.0
        for k in range(window_days):
            i = (idx0 - k) % _TABLE_LEN
            mean += self.precip_mean_mm[i]
            var += self.precip_std_mm[i] ** 2
        return mean, max(math.sqrt(var), _PRECIP_STD_FLOOR)

    def daily_temp_stats(self, dt: datetime) -> Tuple[float, float]:
        i = doy_index(dt) - 1
        return self.temp_mean_c[i], max(self.temp_std_c[i], _TEMP_STD_FLOOR)

    def daily_soil_stats(self, dt: datetime) -> Tuple[float, float]:
        i = doy_index(dt) - 1
        return self.soil_mean_pct[i], max(self.soil_std_pct[i], _SOIL_STD_FLOOR)

    def daily_rh_stats(self, dt: datetime) -> Tuple[float, float]:
        i = doy_index(dt) - 1
        return self.rh_mean_pct[i], max(self.rh_std_pct[i], _RH_STD_FLOOR)

    # --- Serialisation (JSON cache) ---------------------------------------
    def to_dict(self) -> Dict[str, Any]:
        return {
            "zone_id": self.zone_id,
            "source": self.source,
            "n_years": self.n_years,
            "period_start_year": self.period_start_year,
            "period_end_year": self.period_end_year,
            "precip_mean_mm": self.precip_mean_mm,
            "precip_std_mm": self.precip_std_mm,
            "temp_mean_c": self.temp_mean_c,
            "temp_std_c": self.temp_std_c,
            "soil_mean_pct": self.soil_mean_pct,
            "soil_std_pct": self.soil_std_pct,
            "rh_mean_pct": self.rh_mean_pct,
            "rh_std_pct": self.rh_std_pct,
        }

    @classmethod
    def from_dict(cls, d: Dict[str, Any]) -> "ZoneClimatology":
        rh_mean = d.get("rh_mean_pct")
        rh_std = d.get("rh_std_pct")
        if rh_mean is None or len(rh_mean) != _TABLE_LEN:
            rh_mean = [85.0] * _TABLE_LEN
        if rh_std is None or len(rh_std) != _TABLE_LEN:
            rh_std = [_RH_STD_FLOOR] * _TABLE_LEN
        return cls(
            zone_id=d["zone_id"],
            source=d.get("source", "unknown"),
            n_years=int(d.get("n_years", 0)),
            period_start_year=int(d.get("period_start_year", 0)),
            period_end_year=int(d.get("period_end_year", 0)),
            precip_mean_mm=[float(v) for v in d["precip_mean_mm"]],
            precip_std_mm=[float(v) for v in d["precip_std_mm"]],
            temp_mean_c=[float(v) for v in d["temp_mean_c"]],
            temp_std_c=[float(v) for v in d["temp_std_c"]],
            soil_mean_pct=[float(v) for v in d["soil_mean_pct"]],
            soil_std_pct=[float(v) for v in d["soil_std_pct"]],
            rh_mean_pct=[float(v) for v in rh_mean],
            rh_std_pct=[float(v) for v in rh_std],
        )


# ---------------------------------------------------------------------------
# Tier 2 -- deterministic synthetic maritime-continent climatology
# ---------------------------------------------------------------------------

def _synthetic_climatology(zone_id: str, lat: float, n_years: int = 0) -> ZoneClimatology:
    seed = _stable_seed(f"clim_{zone_id}")
    # Deterministic jitter in [0.9, 1.1] from the seed's low bits.
    jitter = 0.9 + 0.2 * ((seed % 1000) / 1000.0)

    abs_lat = abs(lat)
    # Wet-season peak: late Jan in the south, shifting earlier north of ~1N.
    peak_doy = 30.0 if lat <= 1.0 else max(300.0, 30.0 - 12.0 * lat)
    amp_scale = _clip(abs_lat / 8.0, 0.35, 1.0)

    precip_base = max(3.0, (7.0 - 0.25 * abs_lat) * jitter)   # mm/day
    precip_amp  = 4.5 * amp_scale * jitter                    # mm/day
    temp_base   = 27.0 - 0.30 * abs_lat
    # SH zones: coolest around DOY ~200 (mid-Jul). NH: weaker, reversed.
    temp_amp    = 1.3 if lat < 0.0 else -0.5

    precip_mean, precip_std = [], []
    temp_mean, temp_std = [], []
    soil_mean, soil_std = [], []
    rh_mean, rh_std = [], []

    daily_precip_for_soil: List[float] = []
    for doy in range(1, _TABLE_LEN + 1):
        phase = 2.0 * math.pi * (doy - peak_doy) / _DAYS_PER_YEAR
        p = precip_base + precip_amp * math.cos(phase)
        p = max(0.8, p)
        daily_precip_for_soil.append(p)
        precip_mean.append(p)
        precip_std.append(max(_PRECIP_STD_FLOOR, 0.9 * p))

        t_phase = 2.0 * math.pi * (doy - 200.0) / _DAYS_PER_YEAR
        t = temp_base - temp_amp * math.cos(t_phase)
        temp_mean.append(t)
        temp_std.append(0.7)

        r = 86.0 + amp_scale * 6.0 * math.cos(phase)
        rh_mean.append(_clip(r, 78.0, 94.0))
        rh_std.append(max(_RH_STD_FLOOR, 3.5))

    # Soil tracks precip with a 20-day lag.
    for doy in range(1, _TABLE_LEN + 1):
        lagged = daily_precip_for_soil[(doy - 1 - 20) % _TABLE_LEN]
        s = _clip(16.0 + 2.4 * lagged, 8.0, 52.0)
        soil_mean.append(s)
        soil_std.append(max(_SOIL_STD_FLOOR, 4.0))

    return ZoneClimatology(
        zone_id=zone_id,
        source="synthetic_model",
        n_years=n_years,
        period_start_year=0,
        period_end_year=0,
        precip_mean_mm=_circular_smooth(precip_mean),
        precip_std_mm=precip_std,
        temp_mean_c=_circular_smooth(temp_mean),
        temp_std_c=temp_std,
        soil_mean_pct=_circular_smooth(soil_mean),
        soil_std_pct=soil_std,
        rh_mean_pct=_circular_smooth(rh_mean),
        rh_std_pct=rh_std,
    )


# ---------------------------------------------------------------------------
# Tier 1 -- real climatology from the Open-Meteo archive
# ---------------------------------------------------------------------------

def _fetch_openmeteo_climatology(
    zone_id: str,
    lat: float,
    lon: float,
    years: int,
    end_year: Optional[int] = None,
) -> ZoneClimatology:
    if not _REQUESTS_AVAILABLE:
        raise RuntimeError("requests not installed")

    last_full_year = (end_year if end_year is not None
                      else datetime.now(timezone.utc).year - 1)
    start_year = last_full_year - years + 1

    params = {
        "latitude": lat,
        "longitude": lon,
        "start_date": f"{start_year}-01-01",
        "end_date": f"{last_full_year}-12-31",
        "daily": ",".join([
            "precipitation_sum",
            "temperature_2m_mean",
            "soil_moisture_0_to_7cm_mean",
            "relative_humidity_2m_mean",
        ]),
        "timezone": "UTC",
    }

    resp = requests.get(_OPENMETEO_ARCHIVE_URL, params=params, timeout=_TIMEOUT_S)
    resp.raise_for_status()
    data = resp.json()
    daily = data.get("daily", {})
    dates = daily.get("time", [])
    if not dates:
        raise RuntimeError(f"Open-Meteo archive returned no daily rows for {zone_id}")

    precip_series = daily.get("precipitation_sum", [])
    temp_series = daily.get("temperature_2m_mean", [])
    soil_series = daily.get("soil_moisture_0_to_7cm_mean", [])
    rh_series = daily.get("relative_humidity_2m_mean", [])

    p_sum = [0.0] * _TABLE_LEN
    p_sq = [0.0] * _TABLE_LEN
    p_n = [0] * _TABLE_LEN
    t_sum = [0.0] * _TABLE_LEN
    t_sq = [0.0] * _TABLE_LEN
    t_n = [0] * _TABLE_LEN
    s_sum = [0.0] * _TABLE_LEN
    s_sq = [0.0] * _TABLE_LEN
    s_n = [0] * _TABLE_LEN
    r_sum = [0.0] * _TABLE_LEN
    r_sq = [0.0] * _TABLE_LEN
    r_n = [0] * _TABLE_LEN

    def _val(series: List[Any], i: int) -> Optional[float]:
        if i >= len(series):
            return None
        v = series[i]
        if v is None:
            return None
        try:
            return float(v)
        except (TypeError, ValueError):
            return None

    for i, date_str in enumerate(dates):
        try:
            dt = datetime.fromisoformat(date_str).replace(tzinfo=timezone.utc)
        except ValueError:
            continue
        k = doy_index(dt) - 1

        p = _val(precip_series, i)
        if p is not None:
            p_sum[k] += p
            p_sq[k] += p * p
            p_n[k] += 1
        t = _val(temp_series, i)
        if t is not None:
            t_sum[k] += t
            t_sq[k] += t * t
            t_n[k] += 1
        s = _val(soil_series, i)
        if s is not None:
            s_pct = s * 100.0  # m3/m3 -> % (matches ZoneObs.soil_moisture_pct)
            s_sum[k] += s_pct
            s_sq[k] += s_pct * s_pct
            s_n[k] += 1
        r = _val(rh_series, i)
        if r is not None:
            r_sum[k] += r
            r_sq[k] += r * r
            r_n[k] += 1

    if sum(p_n) < 300 * years or sum(t_n) < 300 * years:
        raise RuntimeError(
            f"Open-Meteo archive coverage too thin for {zone_id}: "
            f"precip_days={sum(p_n)} temp_days={sum(t_n)} over {years}y"
        )

    def _mean_std(sums, sqs, ns, floor):
        means, stds = [], []
        for k in range(_TABLE_LEN):
            n = ns[k]
            if n == 0:
                # Should not happen with full-year coverage; guard anyway.
                means.append(0.0)
                stds.append(floor)
                continue
            m = sums[k] / n
            var = max(0.0, sqs[k] / n - m * m)
            means.append(m)
            stds.append(max(floor, math.sqrt(var)))
        return means, stds

    precip_mean, precip_std = _mean_std(p_sum, p_sq, p_n, _PRECIP_STD_FLOOR)
    temp_mean, temp_std = _mean_std(t_sum, t_sq, t_n, _TEMP_STD_FLOOR)

    soil_days = sum(s_n)
    soil_thin = soil_days < 300 * years
    if not soil_thin:
        soil_mean, soil_std = _mean_std(s_sum, s_sq, s_n, _SOIL_STD_FLOOR)
    else:
        logger.warning(
            "climatology: soil_moisture_0_to_7cm_mean coverage thin for %s "
            "(%d days over %dy) -- deriving soil tables from the real precip "
            "series (lagged mapping). Provenance marked 'mixed'.",
            zone_id, soil_days, years,
        )
        soil_mean, soil_std = [], []
        for doy in range(1, _TABLE_LEN + 1):
            lagged = precip_mean[(doy - 1 - 20) % _TABLE_LEN]
            soil_mean.append(_clip(16.0 + 2.4 * lagged, 8.0, 52.0))
            soil_std.append(max(_SOIL_STD_FLOOR, 4.0))

    rh_days = sum(r_n)
    rh_thin = rh_days < 300 * years
    if not rh_thin:
        rh_mean, rh_std = _mean_std(r_sum, r_sq, r_n, _RH_STD_FLOOR)
    else:
        logger.warning(
            "climatology: relative_humidity_2m_mean coverage thin for %s "
            "(%d days over %dy) -- deriving RH tables from the real precip "
            "series (wet-season correlation, same phase). Provenance marked "
            "'mixed'.",
            zone_id, rh_days, years,
        )
        rh_mean, rh_std = [], []
        p_min, p_max = min(precip_mean), max(precip_mean)
        p_span = max(p_max - p_min, 1e-6)
        for doy in range(1, _TABLE_LEN + 1):
            p_frac = (precip_mean[doy - 1] - p_min) / p_span  # 0..1
            rh_mean.append(_clip(78.0 + 12.0 * p_frac, 78.0, 94.0))
            rh_std.append(max(_RH_STD_FLOOR, 3.5))

    source = "openmeteo_archive" if not (soil_thin or rh_thin) else "mixed"

    return ZoneClimatology(
        zone_id=zone_id,
        source=source,
        n_years=years,
        period_start_year=start_year,
        period_end_year=last_full_year,
        precip_mean_mm=_circular_smooth(precip_mean),
        precip_std_mm=precip_std,
        temp_mean_c=_circular_smooth(temp_mean),
        temp_std_c=temp_std,
        soil_mean_pct=_circular_smooth(soil_mean),
        soil_std_pct=soil_std,
        rh_mean_pct=_circular_smooth(rh_mean),
        rh_std_pct=rh_std,
    )


# ---------------------------------------------------------------------------
# Public API: cached climatology + anomaly application
# ---------------------------------------------------------------------------

def _cache_path(zone_id: str, years: int, end_year: Optional[int]) -> Path:
    key = _stable_seed(f"{zone_id}|{years}|{end_year}")
    return _CACHE_DIR / f"{zone_id}_{key}.json"


def get_zone_climatology(
    zone_id: str,
    lat: float,
    lon: float,
    years: int = 10,
    end_year: Optional[int] = None,
    prefer_real: bool = True,
    use_cache: bool = True,
) -> ZoneClimatology:
    path = _cache_path(zone_id, years, end_year)
    if use_cache and path.exists():
        age_days = (
            datetime.now(timezone.utc)
            - datetime.fromtimestamp(path.stat().st_mtime, tz=timezone.utc)
        ).days
        if age_days < _CACHE_TTL_DAYS:
            try:
                with open(path) as f:
                    return ZoneClimatology.from_dict(json.load(f))
            except Exception as e:
                logger.warning("climatology: cache read failed (%s) -- rebuilding", e)

    clim: Optional[ZoneClimatology] = None
    if prefer_real and _REQUESTS_AVAILABLE:
        try:
            clim = _fetch_openmeteo_climatology(zone_id, lat, lon, years, end_year)
            logger.info(
                "climatology: built real %dy climatology for %s (%d-%d)",
                years, zone_id, clim.period_start_year, clim.period_end_year,
            )
        except Exception as e:
            logger.warning(
                "climatology: real fetch failed for %s (%s) -- synthetic model",
                zone_id, e,
            )
            clim = None

    if clim is None:
        clim = _synthetic_climatology(zone_id, lat, n_years=years)

    if use_cache:
        try:
            with open(path, "w") as f:
                json.dump(clim.to_dict(), f)
        except Exception as e:
            logger.warning("climatology: cache write failed (%s) -- continuing", e)

    return clim


def apply_climatology_anomalies(obs: ZoneObs, clim: ZoneClimatology) -> ZoneObs:
    if obs.source == DataSource.SYNTHETIC:
        logger.debug(
            "climatology: %s source is SYNTHETIC -- anomalies left as injected",
            obs.zone_id,
        )
        return obs

    d = obs.to_dict()
    d.pop("_schema_version", None)
    
    precip_values = (obs.precip_30d_mm, obs.precip_14d_mm,
                     obs.precip_7d_mm, obs.precip_24h_mm)
    has_precip_data = any(v != 0.0 for v in precip_values) or obs.precip_anomaly_idx != 0.0
    if has_precip_data:
        mean_w, std_w = clim.window_precip_stats(obs.valid_time, _PRECIP_WINDOW_DAYS)
        d["precip_anomaly_idx"] = _clip(
            (obs.precip_30d_mm - mean_w) / std_w, -5.0, 5.0
        )

    if obs.temp_mean_c != 0.0:
        t_mean, t_std = clim.daily_temp_stats(obs.valid_time)
        d["temp_anomaly_idx"] = _clip((obs.temp_mean_c - t_mean) / t_std, -5.0, 5.0)

    if obs.soil_moisture_pct > 0.0:
        s_mean, s_std = clim.daily_soil_stats(obs.valid_time)
        d["soil_moisture_anom"] = _clip(
            (obs.soil_moisture_pct - s_mean) / s_std, -5.0, 5.0
        )

    if obs.rh_mean_pct > 0.0:
        r_mean, r_std = clim.daily_rh_stats(obs.valid_time)
        d["rh_anomaly_idx"] = _clip(
            (obs.rh_mean_pct - r_mean) / r_std, -5.0, 5.0
        )

    return ZoneObs.from_dict(d)


def apply_anomalies_by_zone_id(
    obs: ZoneObs,
    lat: float,
    lon: float,
    years: int = 10,
    prefer_real: bool = True,
) -> ZoneObs:
    clim = get_zone_climatology(obs.zone_id, lat, lon, years=years,
                                prefer_real=prefer_real)
    return apply_climatology_anomalies(obs, clim)


# ---------------------------------------------------------------------------
# Self-test (python climatology.py) -- fully offline
# ---------------------------------------------------------------------------

if __name__ == "__main__":
    import sys
    from datetime import timezone as _tz
    from zone_observation import make_synthetic_zone_obs

    logging.basicConfig(level=logging.WARNING)
    print("climatology.py self-test (offline: prefer_real=False)\n")
    failures: List[str] = []

    def _assert(cond: bool, msg: str) -> None:
        if not cond:
            failures.append(msg)
            print(f"  FAIL: {msg}")

    LAT, LON = -6.3, 107.3  # Karawang, West Java

    # 1. Synthetic climatology: shape, length, round-trip
    clim = _synthetic_climatology("test_zone", LAT)
    _assert(len(clim.precip_mean_mm) == _TABLE_LEN, "precip table length")
    d = clim.to_dict()
    clim2 = ZoneClimatology.from_dict(d)
    _assert(clim2.zone_id == clim.zone_id, "ZoneClimatology round-trip zone_id")
    _assert(abs(clim2.precip_mean_mm[100] - clim.precip_mean_mm[100]) < 1e-12,
            "ZoneClimatology round-trip values")

    # 2. Seasonality: Java should be much wetter in Jan than in Aug
    jan_mean = sum(clim.precip_mean_mm[0:31]) / 31.0
    aug_mean = sum(clim.precip_mean_mm[212:243]) / 31.0
    _assert(jan_mean > aug_mean * 1.3,
            f"monsoon shape wrong: Jan={jan_mean:.1f} vs Aug={aug_mean:.1f} mm/day")
    print(f"  Seasonality OK: Jan {jan_mean:.1f} mm/day vs Aug {aug_mean:.1f} mm/day")

    # 3. doy_index leap mapping: Mar 1 maps to the same entry in every year
    d1 = doy_index(datetime(2023, 3, 1, tzinfo=_tz.utc))
    d2 = doy_index(datetime(2024, 3, 1, tzinfo=_tz.utc))
    _assert(d1 == d2 == 61, f"Mar 1 mapping inconsistent: {d1} vs {d2}")
    _assert(doy_index(datetime(2024, 2, 29, tzinfo=_tz.utc)) == 60, "Feb 29 mapping")
    print(f"  doy_index OK (Mar 1 -> {d1}, Feb 29 -> 60)")

    # 4. Window stats: 30-day wet-season aggregate exceeds dry-season
    wet_dt = datetime(2024, 1, 31, tzinfo=_tz.utc)
    dry_dt = datetime(2024, 8, 31, tzinfo=_tz.utc)
    wet_mean, wet_std = clim.window_precip_stats(wet_dt, 30)
    dry_mean, _ = clim.window_precip_stats(dry_dt, 30)
    _assert(wet_mean > dry_mean, "window aggregate seasonality wrong")
    _assert(wet_std >= _PRECIP_STD_FLOOR, "window std floor violated")
    print(f"  Window stats OK: wet30={wet_mean:.0f}mm dry30={dry_mean:.0f}mm")

    # 5. apply_climatology_anomalies: real-source obs gets anomalies
    obs = make_synthetic_zone_obs("realish_zone", seed=1)
    od = obs.to_dict()
    od.pop("_schema_version", None)
    od["source"] = DataSource.OPENMETEO_LIVE.value   # pretend real
    # A real fetcher leaves anomaly fields at 0.0 ("unset") -- mirror that so
    # this test measures exactly what this module adds.
    od["precip_anomaly_idx"] = od["temp_anomaly_idx"] = od["soil_moisture_anom"] = 0.0
    real_obs = ZoneObs.from_dict(od)
    pre_precip_z = real_obs.precip_anomaly_idx
    out = apply_climatology_anomalies(real_obs, clim)
    _assert(out is not real_obs, "apply should return a NEW object")
    _assert(real_obs.precip_anomaly_idx == pre_precip_z, "input obs was mutated!")
    # The obs built by make_synthetic_zone_obs has neutral-ish aggregates;
    # anomaly must be finite and within clip range.
    _assert(-5.0 <= out.precip_anomaly_idx <= 5.0, "anomaly outside clip")
    _assert(-5.0 <= out.temp_anomaly_idx <= 5.0, "temp anomaly outside clip")
    _assert(-5.0 <= out.soil_moisture_anom <= 5.0, "soil anomaly outside clip")
    print(f"  Anomaly application OK: precip_z={out.precip_anomaly_idx:+.2f} "
          f"temp_z={out.temp_anomaly_idx:+.2f} soil_z={out.soil_moisture_anom:+.2f}")

    # 6. Drought/wet extremes produce correctly-signed anomalies
    # (keep aggregates monotonic: 24h <= 7d <= 14d <= 30d)
    dry_obs_d = dict(od)
    dry_obs_d["precip_24h_mm"] = 0.0
    dry_obs_d["precip_7d_mm"] = 0.01 * wet_mean / 4.0
    dry_obs_d["precip_14d_mm"] = 0.02 * wet_mean / 2.0
    dry_obs_d["precip_30d_mm"] = 0.05 * wet_mean   # 5% of wet climatology
    dry_obs_d["valid_time"] = wet_dt.isoformat()
    dry_out = apply_climatology_anomalies(ZoneObs.from_dict(dry_obs_d), clim)
    _assert(dry_out.precip_anomaly_idx < -1.0,
            f"dry obs should get negative anomaly, got {dry_out.precip_anomaly_idx}")
    wet_obs_d = dict(od)
    wet_obs_d["precip_24h_mm"] = 2.5 * dry_mean / 30.0
    wet_obs_d["precip_7d_mm"] = 2.5 * dry_mean / 4.0
    wet_obs_d["precip_14d_mm"] = 2.5 * dry_mean / 2.0
    wet_obs_d["precip_30d_mm"] = 2.5 * dry_mean
    wet_obs_d["valid_time"] = dry_dt.isoformat()
    wet_out = apply_climatology_anomalies(ZoneObs.from_dict(wet_obs_d), clim)
    _assert(wet_out.precip_anomaly_idx > 1.0,
            f"wet obs should get positive anomaly, got {wet_out.precip_anomaly_idx}")
    print(f"  Sign check OK: dry_z={dry_out.precip_anomaly_idx:+.2f} "
          f"wet_z={wet_out.precip_anomaly_idx:+.2f}")

    # 7. SYNTHETIC-source obs is skipped unchanged
    syn = make_synthetic_zone_obs("syn_zone", drought=True, seed=2)
    syn_out = apply_climatology_anomalies(syn, clim)
    _assert(syn_out.precip_anomaly_idx == syn.precip_anomaly_idx,
            "SYNTHETIC obs anomaly was modified (should be skipped)")
    print("  SYNTHETIC skip OK")

    # 8. Zero-precip obs (smap-style real fetch: anomaly fields unset at 0.0)
    # must NOT read as a catastrophic false drought.
    zero_d = dict(od)
    zero_d["precip_24h_mm"] = zero_d["precip_7d_mm"] = 0.0
    zero_d["precip_14d_mm"] = zero_d["precip_30d_mm"] = 0.0
    zero_d["precip_anomaly_idx"] = 0.0   # real fetchers leave this unset
    zero_out = apply_climatology_anomalies(ZoneObs.from_dict(zero_d), clim)
    _assert(zero_out.precip_anomaly_idx == 0.0,
            "precip-less fetch should keep anomaly 0.0 (no false drought)")
    print("  Zero-precip guard OK")

    # 9. Cache round-trip through get_zone_climatology (offline path)
    clim_c = get_zone_climatology("cache_zone", LAT, LON, years=5,
                                  prefer_real=False, use_cache=True)
    clim_c2 = get_zone_climatology("cache_zone", LAT, LON, years=5,
                                   prefer_real=False, use_cache=True)
    _assert(clim_c.precip_mean_mm == clim_c2.precip_mean_mm,
            "cached climatology not identical")
    print("  Cache round-trip OK")

    print()
    if failures:
        print(f"FAILED  {len(failures)} test(s):")
        for f in failures:
            print(f"  - {f}")
        sys.exit(1)
    else:
        print("All 9 test groups passed.")