era5_data_pipeline.py

"""
era5_data_pipeline.py
=====================
Causal observation fetchers for WeatherForecastEnv / historical caches.

Contract
--------
- Callers pass an explicit valid_time (the "as of" date). Fetchers build
  their own BACKWARD 30-day window ending on that date. There is no
  date_range tuple for observations.
- precip_24h / 7d / 14d / 30d are last-N-days totals through valid_time,
  never a forward sum from valid_time.
- Historical Open-Meteo dates use archive-api.open-meteo.com and are
  stamped DataSource.OPENMETEO_ARCHIVE when that enum member exists.
- Insufficient history raises; it does not zero-pad (zero-pad looks like
  drought).
- A force_data_source other than SYNTHETIC does not fall back to
  make_synthetic_zone_obs unless allow_synthetic_obs_fallback=True.
- This module does not issue real NWP forecasts. forecast_backend="baseline"
  is persistence of causal precip_30d and must not be reported as skill.
"""

from __future__ import annotations

import hashlib
import json
import logging
import math
import os
import random
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"Schema mismatch: expected 3, got {_zo.SCHEMA_VERSION}"
)

from zone_observation import (
    BasinContext,
    DataSource,
    EpisodeContext,
    ForecastConfig,
    GeoPolygon,
    ZoneObs,
    derive_helio_regime,
    make_synthetic_basin_context,
    make_synthetic_episode_context,
    make_synthetic_zone_obs,
    _stable_seed,
)

logger = logging.getLogger(__name__)


def _openmeteo_archive_source():
    """Prefer DataSource.OPENMETEO_ARCHIVE when zone_observation has it.

    Historical dates must not be stamped OPENMETEO_LIVE. If the enum member
    has not been added yet, fall back to LIVE and warn -- the fetch URL is
    still the archive endpoint; only the provenance tag is degraded.
    """
    src = getattr(DataSource, "OPENMETEO_ARCHIVE", None)
    if src is not None:
        return src
    logger.warning(
        "DataSource.OPENMETEO_ARCHIVE is missing from zone_observation.py. "
        "Add OPENMETEO_ARCHIVE = 'openmeteo_archive' next to OPENMETEO_LIVE. "
        "Stamping OPENMETEO_LIVE for archive fetches until that lands."
    )
    return DataSource.OPENMETEO_LIVE


# ---------------------------------------------------------------------------
# Optional dependencies
# ---------------------------------------------------------------------------
try:
    import requests
    REQUESTS_AVAILABLE = True
except ImportError:
    REQUESTS_AVAILABLE = False
    logger.warning("requests not installed — Open-Meteo unavailable")

try:
    import cdsapi
    CDSAPI_AVAILABLE = True
except ImportError:
    CDSAPI_AVAILABLE = False
    logger.info("cdsapi not installed — ERA5 will fall back to synthetic")

try:
    import numpy as np
    NUMPY_AVAILABLE = True
except ImportError:
    NUMPY_AVAILABLE = False

try:
    import ee  # Google Earth Engine — IMERG / CHIRPS / SMAP access
    EE_AVAILABLE = True
except ImportError:
    EE_AVAILABLE = False
    logger.info("earthengine-api not installed — satellite sources unavailable")

_EE_INITIALIZED = False


def _ensure_ee_initialized() -> None:
    """Lazily call ee.Initialize() once per process."""
    global _EE_INITIALIZED
    if _EE_INITIALIZED:
        return
    if not EE_AVAILABLE:
        raise RuntimeError("earthengine-api not installed")
    project = os.environ.get("EARTHENGINE_PROJECT")
    if project:
        ee.Initialize(project=project)
    else:
        ee.Initialize()
    _EE_INITIALIZED = True

# ---------------------------------------------------------------------------
# Zone registry
# ---------------------------------------------------------------------------
_ZONE_REGISTRY: Dict[str, GeoPolygon] = {}


def register_zone(polygon: GeoPolygon) -> None:
    """Register a sourcing zone polygon for lat/lon resolution."""
    _ZONE_REGISTRY[polygon.zone_id] = polygon
    logger.info(
        f"Registered zone {polygon.zone_id} centroid={polygon.centroid}"
    )


def _resolve_latlon(zone_id: str) -> Tuple[float, float]:
    if zone_id not in _ZONE_REGISTRY:
        raise KeyError(
            f"Zone '{zone_id}' not registered. "
            f"Call register_zone() before fetching data."
        )
    return _ZONE_REGISTRY[zone_id].centroid


# ---------------------------------------------------------------------------
# Config / cache
# ---------------------------------------------------------------------------
_CACHE_DIR = Path(os.environ.get("WEATHER_CACHE_DIR", ".cache/era5"))
_CACHE_DIR.mkdir(parents=True, exist_ok=True)

_ERA5_CACHE_DIR  = _CACHE_DIR / "era5_nc"
_ERA5_CACHE_DIR.mkdir(parents=True, exist_ok=True)

_OPENMETEO_URL   = "https://api.open-meteo.com/v1/forecast"
_OPENMETEO_ARCHIVE_URL = "https://archive-api.open-meteo.com/v1/archive"

_TIMEOUT_S       = int(os.environ.get("WEATHER_HTTP_TIMEOUT",    "30"))
_CACHE_TTL_DAYS  = int(os.environ.get("WEATHER_CACHE_TTL_DAYS",  "7"))
_ERA5_TTL_DAYS   = int(os.environ.get("WEATHER_ERA5_TTL_DAYS",   "30"))

# ERA5 bounding box padding in degrees around zone centroid
_ERA5_BOX_PAD    = float(os.environ.get("WEATHER_ERA5_BOX_PAD",  "0.5"))


# ---------------------------------------------------------------------------
# Cache utilities
# ---------------------------------------------------------------------------
def _cache_key(url: str, params: Dict[str, Any]) -> str:
    payload = f"{url}{json.dumps(params, sort_keys=True)}"
    return hashlib.sha256(payload.encode()).hexdigest()


def _cached_get(url: str, params: Dict[str, Any]) -> Dict[str, Any]:
    """HTTP GET with file-based JSON cache."""
    if not REQUESTS_AVAILABLE:
        raise RuntimeError("requests not installed — cannot fetch HTTP data")

    path = _CACHE_DIR / f"{_cache_key(url, params)}.json"

    if path.exists():
        age = (
            datetime.now(timezone.utc)
            - datetime.fromtimestamp(path.stat().st_mtime, tz=timezone.utc)
        ).days
        if age < _CACHE_TTL_DAYS:
            with open(path) as f:
                return json.load(f)

    resp = requests.get(url, params=params, timeout=_TIMEOUT_S)
    resp.raise_for_status()
    data = resp.json()

    with open(path, "w") as f:
        json.dump(data, f)

    return data


def _era5_cache_path(zone_id: str, date_range: Tuple[datetime, datetime]) -> Path:
    """Deterministic cache file path for an ERA5 download."""
    key = _stable_seed(
        zone_id
        + date_range[0].date().isoformat()
        + date_range[1].date().isoformat()
    )
    return _ERA5_CACHE_DIR / f"{zone_id}_{key}.nc"


# ---------------------------------------------------------------------------
# Timezone helper
# ---------------------------------------------------------------------------
def _ensure_utc(dt: datetime) -> datetime:
    """Normalise a datetime to UTC at the pipeline boundary."""
    if dt.tzinfo is None:
        return dt.replace(tzinfo=timezone.utc)
    return dt.astimezone(timezone.utc)


_ERA5_VARIABLES: List[str] = [
    "total_precipitation",
    "2m_temperature",
    "2m_dewpoint_temperature",
    "10m_u_component_of_wind",
    "10m_v_component_of_wind",
    "surface_pressure",
    "volumetric_soil_water_layer_1",
    "potential_evaporation",
]


def _dewpoint_to_rh(temp_c: float, dewpoint_c: float) -> float:
    """Magnus formula: relative humidity from temperature and dewpoint (%)."""
    a, b = 17.625, 243.04  # Magnus coefficients
    rh = 100.0 * math.exp(
        (a * dewpoint_c / (b + dewpoint_c))
        - (a * temp_c / (b + temp_c))
    )
    return max(0.0, min(100.0, rh))


def _build_era5_obs(
    zone_id:    str,
    start:      datetime,
    nc_path:    Path,
) -> ZoneObs:
    """Map a downloaded ERA5 NetCDF file to a ZoneObs.

    The NetCDF file covers a BACKWARD-anchored window ending at `start`
    (see _fetch_era5's request construction) -- so the most recent data is
    at the END of each variable's time axis, not the beginning. All
    indexing here is from the end ([-N:]), not the start ([:N]).
    """
    try:
        import netCDF4 as nc  # type: ignore
    except ImportError:
        try:
            import xarray as xr  # type: ignore
            return _build_era5_obs_xarray(zone_id, start, nc_path)
        except ImportError:
            raise ImportError(
                "netCDF4 or xarray required for ERA5 ingestion. "
                "Install with: pip install netCDF4  or  pip install xarray"
            )

    import numpy as np

    ds = nc.Dataset(str(nc_path), "r")
    try:
        def _daily_mean(var: str) -> float:
            if var not in ds.variables:
                return 0.0
            data = ds.variables[var][-24:].flatten()
            data = np.ma.filled(data, np.nan)
            valid = data[np.isfinite(data)]
            return float(np.mean(valid)) if len(valid) > 0 else 0.0

        def _daily_max(var: str) -> float:
            if var not in ds.variables:
                return 0.0
            data = ds.variables[var][-24:].flatten()
            data = np.ma.filled(data, np.nan)
            valid = data[np.isfinite(data)]
            return float(np.max(valid)) if len(valid) > 0 else 0.0

        def _daily_sum(var: str, scale: float = 1.0) -> float:
            if var not in ds.variables:
                return 0.0
            data = ds.variables[var][-24:].flatten()
            data = np.ma.filled(data, np.nan)
            valid = data[np.isfinite(data)]
            return float(np.sum(valid) * scale) if len(valid) > 0 else 0.0

        temp_mean_c = _daily_mean("t2m") - 273.15
        temp_max_c  = _daily_max("t2m")  - 273.15
        if "t2m" in ds.variables:
            t_data = np.ma.filled(ds.variables["t2m"][-24:].flatten(), np.nan)
            valid  = t_data[np.isfinite(t_data)]
            temp_min_c = float(np.min(valid)) - 273.15 if len(valid) > 0 else temp_mean_c
        else:
            temp_min_c = temp_mean_c

        dewpoint_mean_c = _daily_mean("d2m") - 273.15
        rh_mean = _dewpoint_to_rh(temp_mean_c, dewpoint_mean_c)
        dewpoint_max_c  = _daily_max("d2m") - 273.15
        rh_max  = _dewpoint_to_rh(temp_min_c, dewpoint_max_c)

        precip_24h = _daily_sum("tp", scale=1000.0)

        total_hours = len(ds.variables.get("tp", [])) if "tp" in ds.variables else 24
        MIN_HOURS_30D = 24 * 30
        if total_hours < MIN_HOURS_30D:
            raise RuntimeError(
                f"_build_era5_obs: only {total_hours} hours of 'tp' available "
                f"for {zone_id} at {start.date().isoformat()}, need "
                f"{MIN_HOURS_30D} (30 days) for backward precip accumulators. "
                f"Refusing to zero-pad or extrapolate -- drop this point or "
                f"widen the ERA5 request window."
            )

        def _window_sum(var: str, hours: int, scale: float = 1.0) -> float:
            if var not in ds.variables:
                return 0.0
            data = np.ma.filled(
                ds.variables[var][-hours:].flatten(), np.nan
            )
            valid = data[np.isfinite(data)]
            return float(np.sum(valid) * scale) if len(valid) > 0 else 0.0
        precip_7d  = _window_sum("tp", 24*7,  scale=1000.0)
        precip_14d = _window_sum("tp", 24*14, scale=1000.0)
        precip_30d = _window_sum("tp", 24*30, scale=1000.0)

        u = _daily_mean("u10")
        v = _daily_mean("v10")
        wind_mean = math.sqrt(u**2 + v**2)
        u_max = _daily_max("u10")
        v_max = _daily_max("v10")
        wind_max = math.sqrt(u_max**2 + v_max**2)

        soil_pct = _daily_mean("swvl1") * 100.0

        pe_m = _daily_sum("pev", scale=1.0)
        et0_mm = abs(pe_m) * 1000.0

        return ZoneObs(
            zone_id=zone_id,
            valid_time=start,
            source=DataSource.ERA5_REANALYSIS,
            precip_24h_mm=max(0.0, precip_24h),
            precip_7d_mm=max(0.0, precip_7d),
            precip_14d_mm=max(0.0, precip_14d),
            precip_30d_mm=max(0.0, precip_30d),
            temp_mean_c=temp_mean_c,
            temp_max_c=max(temp_mean_c, temp_max_c),
            temp_min_c=min(temp_mean_c, temp_min_c),
            temp_anomaly_idx=0.0,
            precip_anomaly_idx=0.0,
            evapotranspiration_mm=max(0.0, et0_mm),
            wind_speed_mean_ms=max(0.0, wind_mean),
            wind_speed_max_ms=max(0.0, wind_max),
            rh_mean_pct=rh_mean,
            rh_max_pct=max(rh_mean, rh_max),
            soil_moisture_pct=max(0.0, min(100.0, soil_pct)),
            soil_moisture_anom=0.0,
            quality_flag=0,
        )
    finally:
        ds.close()


def _build_era5_obs_xarray(
    zone_id: str,
    start:   datetime,
    nc_path: Path,
) -> ZoneObs:
    """xarray fallback for _build_era5_obs.

    This function's docstring previously claimed to mirror
    _build_era5_obs's backward-window, end-indexed convention, but only
    the precip window-sums (_window_sum_xr) actually did -- _mean/_max/
    _sum had no windowing at all and aggregated the ENTIRE downloaded
    cube (all 30 days) for every "daily" stat: temp_mean/max/min, RH,
    wind, soil moisture, ET0, and precip_24h. Fixed by giving every
    helper the same [-hours:] slicing _build_era5_obs uses, so this
    fallback path (taken whenever netCDF4 isn't installed but xarray is)
    can no longer silently diverge from the primary path's correctness.
    """
    import xarray as xr
    import numpy as np

    ds = xr.open_dataset(str(nc_path))
    try:
        def _windowed(var: str, hours: int) -> "np.ndarray":
            if var not in ds:
                return np.array([])
            vals = ds[var].values.flatten()[-hours:]
            return vals[np.isfinite(vals)]

        def _daily_mean(var: str) -> float:
            valid = _windowed(var, 24)
            return float(valid.mean()) if len(valid) > 0 else 0.0

        def _daily_max(var: str) -> float:
            valid = _windowed(var, 24)
            return float(valid.max()) if len(valid) > 0 else 0.0

        def _daily_min(var: str, fallback: float) -> float:
            valid = _windowed(var, 24)
            return float(valid.min()) if len(valid) > 0 else fallback

        def _daily_sum(var: str, scale: float = 1.0) -> float:
            valid = _windowed(var, 24)
            return float(valid.sum() * scale) if len(valid) > 0 else 0.0

        def _window_sum_xr(var: str, hours: int, scale: float = 1.0) -> float:
            valid = _windowed(var, hours)
            return float(valid.sum() * scale) if len(valid) > 0 else 0.0

        temp_mean_c = _daily_mean("t2m") - 273.15
        temp_max_c  = _daily_max("t2m")  - 273.15
        temp_min_c  = _daily_min("t2m", fallback=temp_mean_c + 273.15) - 273.15

        dewpoint_mean_c = _daily_mean("d2m") - 273.15
        dewpoint_max_c  = _daily_max("d2m")  - 273.15
        rh_mean = _dewpoint_to_rh(temp_mean_c, dewpoint_mean_c)
        rh_max  = _dewpoint_to_rh(temp_min_c, dewpoint_max_c)

        precip_24h = max(0.0, _daily_sum("tp", scale=1000.0))

        total_hours = len(ds["tp"].values.flatten()) if "tp" in ds else 0
        MIN_HOURS_30D = 24 * 30
        if total_hours < MIN_HOURS_30D:
            raise RuntimeError(
                f"_build_era5_obs_xarray: only {total_hours} hours of 'tp' "
                f"available for {zone_id} at {start.date().isoformat()}, need "
                f"{MIN_HOURS_30D} (30 days) for backward precip accumulators. "
                f"Refusing to zero-pad or extrapolate -- drop this point or "
                f"widen the ERA5 request window."
            )

        precip_7d  = max(0.0, _window_sum_xr("tp", 24 * 7,  scale=1000.0))
        precip_14d = max(0.0, _window_sum_xr("tp", 24 * 14, scale=1000.0))
        precip_30d = max(0.0, _window_sum_xr("tp", 24 * 30, scale=1000.0))

        u = _daily_mean("u10")
        v = _daily_mean("v10")
        wind_mean = math.sqrt(u**2 + v**2)
        wind_max  = math.sqrt(_daily_max("u10")**2 + _daily_max("v10")**2)
        soil_pct  = _daily_mean("swvl1") * 100.0
        et0_mm    = abs(_daily_sum("pev")) * 1000.0

        return ZoneObs(
            zone_id=zone_id,
            valid_time=start,
            source=DataSource.ERA5_REANALYSIS,
            precip_24h_mm=precip_24h,
            precip_7d_mm=precip_7d,
            precip_14d_mm=precip_14d,
            precip_30d_mm=precip_30d,
            temp_mean_c=temp_mean_c,
            temp_max_c=max(temp_mean_c, temp_max_c),
            temp_min_c=min(temp_mean_c, temp_min_c),
            evapotranspiration_mm=max(0.0, et0_mm),
            wind_speed_mean_ms=max(0.0, wind_mean),
            wind_speed_max_ms=max(0.0, wind_max),
            rh_mean_pct=rh_mean,
            rh_max_pct=max(rh_mean, rh_max),
            soil_moisture_pct=max(0.0, min(100.0, soil_pct)),
            quality_flag=0,
        )
    finally:
        ds.close()


def _fetch_openmeteo(zone_id: str, valid_time: datetime) -> ZoneObs:
    """Fetch from Open-Meteo API with full variable coverage.

    Takes an explicit valid_time (the anchor / "as of" date), not a
    (start, end) tuple. This function always requests and computes a
    BACKWARD 30-day window ending at valid_time.

    Why not a date_range tuple: this codebase had TWO callers of the
    fetch pipeline with OPPOSITE assumptions about which end of a
    (start, end) tuple held the anchor. build_continuous_historical.py
    called fetch_episode_context(zone_id, (day, day+30), cfg) treating
    index 0 as the anchor; backtest_indonesia.py called
    fetch_zone_obs(zone_id, (vt-35, vt), cfg) treating index 1 as the
    anchor. Under the old code (which used date_range[0] as both
    valid_time and the start of a forward sum), the first caller leaked
    up to 29 days of future precipitation into "last N days" fields, and
    the second caller silently set obs.valid_time to vt-35 instead of vt
    (stale by 35 days, in the opposite direction) -- two different bugs
    from one ambiguous parameter shape. An explicit valid_time parameter
    makes both classes of bug structurally impossible for any caller.
    """
    lat, lon = _resolve_latlon(zone_id)
    anchor   = _ensure_utc(valid_time)

    MIN_DAYS_30D = 30
    api_start = anchor - timedelta(days=MIN_DAYS_30D - 1)

    today = datetime.now(timezone.utc).date()
    use_archive = anchor.date() < today - timedelta(days=5)
    url = _OPENMETEO_ARCHIVE_URL if use_archive else _OPENMETEO_URL

    daily_vars = ",".join([
        "temperature_2m_mean",
        "temperature_2m_max",
        "temperature_2m_min",
        "precipitation_sum",
        "precipitation_hours",
        "wind_speed_10m_mean",
        "wind_speed_10m_max",
        "relative_humidity_2m_mean",
        "relative_humidity_2m_max",
        "et0_fao_evapotranspiration",
        "shortwave_radiation_sum",
        "soil_moisture_0_to_7cm_mean",
    ])

    params = {
        "latitude":   lat,
        "longitude":  lon,
        "daily":      daily_vars,
        "start_date": api_start.date().isoformat(),
        "end_date":   anchor.date().isoformat(),
        "timezone":   "UTC",
    }

    data  = _cached_get(url, params)
    daily = data.get("daily", {})

    # Pin the series to the anchor calendar date. Open-Meteo usually
    # honors end_date, but the live forecast endpoint can still return
    # days after `anchor`. Indexing "from the end" of that payload would
    # reintroduce a short forward leak. Slice through the last index
    # whose `time` label equals the anchor (or the last element if no
    # time axis is present).
    times = daily.get("time") or []
    anchor_iso = anchor.date().isoformat()
    if times:
        end_idx = None
        for i in range(len(times) - 1, -1, -1):
            if str(times[i])[:10] == anchor_iso:
                end_idx = i
                break
        if end_idx is None:
            raise RuntimeError(
                f"_fetch_openmeteo: anchor {anchor_iso} not present on "
                f"Open-Meteo time axis for {zone_id} "
                f"({times[0]} .. {times[-1]})."
            )
        daily = {
            k: (v[: end_idx + 1] if isinstance(v, list) else v)
            for k, v in daily.items()
        }

    # Chronological, ending at `anchor` -- index from the END of the
    # sliced series, so "last N days" means the N most recent days
    # relative to anchor, never days after it.
    precip_series = daily.get("precipitation_sum", [])
    n_valid_precip = sum(1 for v in precip_series if v is not None)
    if n_valid_precip < MIN_DAYS_30D:
        raise RuntimeError(
            f"_fetch_openmeteo: only {n_valid_precip} valid precipitation "
            f"days available for {zone_id} ending {anchor.date().isoformat()} "
            f"(need {MIN_DAYS_30D}). Refusing to zero-pad -- that would "
            f"fabricate a drought signal. Caller should drop this point or "
            f"fall back to another source."
        )

    def _safe_last(key: str) -> float:
        vals = daily.get(key, [])
        return float(vals[-1]) if vals and vals[-1] is not None else 0.0

    def _safe_sum_last(key: str, n: int) -> float:
        vals = daily.get(key, [])
        window = vals[-n:] if len(vals) >= n else vals
        return float(sum(v for v in window if v is not None))

    temp_mean = _safe_last("temperature_2m_mean")
    temp_max  = _safe_last("temperature_2m_max")
    temp_min  = _safe_last("temperature_2m_min")

    temp_max  = max(temp_mean, temp_max)
    temp_min  = min(temp_mean, temp_min)

    rh_mean = _safe_last("relative_humidity_2m_mean")
    rh_max  = _safe_last("relative_humidity_2m_max")
    rh_max  = max(rh_mean, rh_max)

    wind_mean = _safe_last("wind_speed_10m_mean")
    wind_max  = _safe_last("wind_speed_10m_max")
    wind_max  = max(wind_mean, wind_max)

    return ZoneObs(
        zone_id=zone_id,
        valid_time=anchor,
        source=_openmeteo_archive_source() if use_archive else DataSource.OPENMETEO_LIVE,
        precip_24h_mm=_safe_last("precipitation_sum"),
        precip_7d_mm=_safe_sum_last("precipitation_sum", 7),
        precip_14d_mm=_safe_sum_last("precipitation_sum", 14),
        precip_30d_mm=_safe_sum_last("precipitation_sum", 30),
        temp_mean_c=temp_mean,
        temp_max_c=temp_max,
        temp_min_c=temp_min,
        evapotranspiration_mm=_safe_last("et0_fao_evapotranspiration"),
        wind_speed_mean_ms=wind_mean / 3.6,
        wind_speed_max_ms=wind_max  / 3.6,
        rh_mean_pct=rh_mean,
        rh_max_pct=rh_max,
        soil_moisture_pct=max(0.0, min(100.0,
                               _safe_last("soil_moisture_0_to_7cm_mean") * 100.0)),
        quality_flag=0,
    )


def _fetch_era5(zone_id: str, valid_time: datetime) -> ZoneObs:
    """Explicit valid_time (anchor); a BACKWARD 30-day CDS request is
    constructed from it internally. See _fetch_openmeteo's docstring for
    why this takes valid_time rather than a date_range tuple."""
    if not CDSAPI_AVAILABLE:
        raise RuntimeError(
            f"_fetch_era5: cdsapi is not installed — cannot fetch ERA5 for "
            f"{zone_id} at {valid_time.date().isoformat()}. fetch_zone_obs "
            f"may fall back to another source; this function will not invent "
            f"synthetic weather itself."
        )

    lat, lon = _resolve_latlon(zone_id)
    anchor   = _ensure_utc(valid_time)
    start    = anchor - timedelta(days=29)
    end      = anchor

    bbox = [
        round(lat + _ERA5_BOX_PAD, 2),
        round(lon - _ERA5_BOX_PAD, 2),
        round(lat - _ERA5_BOX_PAD, 2),
        round(lon + _ERA5_BOX_PAD, 2),
    ]

    nc_path = _era5_cache_path(zone_id, (start, end))

    if nc_path.exists():
        age_days = (
            datetime.now(timezone.utc)
            - datetime.fromtimestamp(nc_path.stat().st_mtime, tz=timezone.utc)
        ).days
        if age_days < _ERA5_TTL_DAYS:
            logger.debug("ERA5 cache hit for %s", zone_id)
        else:
            nc_path.unlink()

    if not nc_path.exists():
        years:  set = set()
        months: set = set()
        days:   set = set()
        cursor = start.date()
        end_date = end.date()
        while cursor <= end_date:
            years.add(cursor.year)
            months.add(cursor.month)
            days.add(cursor.day)
            cursor += timedelta(days=1)

        logger.info(
            "Requesting ERA5 data for %s  bbox=%s  dates=%s to %s (backward, "
            "ending at anchor)",
            zone_id, bbox,
            start.date().isoformat(),
            end.date().isoformat(),
        )
        try:
            client = cdsapi.Client(quiet=True)
            client.retrieve(
                "reanalysis-era5-single-levels",
                {
                    "product_type": "reanalysis",
                    "variable":     _ERA5_VARIABLES,
                    "year":         sorted(str(y) for y in years),
                    "month":        sorted(f"{m:02d}" for m in months),
                    "day":          sorted(f"{d:02d}" for d in days),
                    "time":         [f"{h:02d}:00" for h in range(24)],
                    "area":         bbox,
                    "format":       "netcdf",
                },
                str(nc_path),
            )
            logger.info("ERA5 download complete: %s", nc_path)
        except Exception as e:
            if nc_path.exists():
                nc_path.unlink()
            raise RuntimeError(
                f"_fetch_era5: CDS request failed for {zone_id}: {e}"
            ) from e

    try:
        return _build_era5_obs(zone_id, anchor, nc_path)
    except Exception as e:
        raise RuntimeError(
            f"_fetch_era5: NetCDF parse failed for {zone_id}: {e}"
        ) from e


def _fetch_imerg(zone_id: str, valid_time: datetime) -> ZoneObs:
    """Explicit valid_time. Also fixes the same forward-window bug found
    in _fetch_openmeteo: the previous version queried
    [start, start+n_days) for "last n_days" precip -- forward, not
    backward. Now queries [start-n_days, start+1day).

    The "+1day" on the end bound matters: Earth Engine's
    ImageCollection.filterDate(start, end) is half-open -- [start, end).
    Passing end=start (valid_time) would silently exclude valid_time's
    own day's image(s), making every "last N days" window actually mean
    "last N days ending yesterday." Using end=start+1day makes it
    genuinely inclusive of valid_time, matching precip_24h/7d/14d/30d's
    documented meaning and _fetch_openmeteo's equivalent window.
    """
    _ensure_ee_initialized()

    lat, lon = _resolve_latlon(zone_id)
    start = _ensure_utc(valid_time)
    inclusive_end = start + timedelta(days=1)

    region = ee.Geometry.Point([lon, lat]).buffer(_ERA5_BOX_PAD * 111_000)

    def _window_sum_mm(n_days: int) -> float:
        window_start = start - timedelta(days=n_days - 1)
        coll = (
            ee.ImageCollection("NASA/GPM_L3/IMERG_V07")
            .filterDate(window_start.isoformat(), inclusive_end.isoformat())
            .filterBounds(region)
            .select("precipitation")
        )
        total_mm_image = coll.sum().multiply(0.5)
        stats = total_mm_image.reduceRegion(
            reducer=ee.Reducer.mean(), geometry=region, scale=11_000, bestEffort=True,
        ).getInfo()
        return float(stats.get("precipitation", 0.0) or 0.0)

    precip_24h = max(0.0, _window_sum_mm(1))
    precip_7d  = max(precip_24h, _window_sum_mm(7))
    precip_14d = max(precip_7d,  _window_sum_mm(14))
    precip_30d = max(precip_14d, _window_sum_mm(30))

    return ZoneObs(
        zone_id=zone_id,
        valid_time=start,
        source=DataSource.SATELLITE_PRECIP,
        precip_24h_mm=precip_24h,
        precip_7d_mm=precip_7d,
        precip_14d_mm=precip_14d,
        precip_30d_mm=precip_30d,
        precip_satellite_mm=precip_24h,
        quality_flag=0,
    )


def _fetch_smap(zone_id: str, valid_time: datetime) -> ZoneObs:
    """Single point-in-time satellite soil moisture read (not a
    multi-day accumulator, so the forward/backward window issue that
    affects the precip fetchers doesn't apply the same way here)."""
    _ensure_ee_initialized()

    lat, lon = _resolve_latlon(zone_id)
    start = _ensure_utc(valid_time)
    end = start + timedelta(days=1)

    region = ee.Geometry.Point([lon, lat]).buffer(_ERA5_BOX_PAD * 111_000)

    coll = (
        ee.ImageCollection("NASA/SMAP/SPL4SMGP/007")
        .filterDate(start.isoformat(), end.isoformat())
        .filterBounds(region)
        .select("sm_surface")
    )
    stats = coll.mean().reduceRegion(
        reducer=ee.Reducer.mean(), geometry=region, scale=9_000, bestEffort=True,
    ).getInfo()

    vwc = stats.get("sm_surface")
    if vwc is None:
        raise RuntimeError(f"SMAP: no data returned for zone={zone_id} date={start.date()}")
    soil_pct = max(0.0, min(100.0, float(vwc) * 100.0))

    return ZoneObs(
        zone_id=zone_id,
        valid_time=start,
        source=DataSource.SATELLITE_SOIL,
        soil_moisture_pct=soil_pct,
        soil_moisture_satellite_pct=soil_pct,
        quality_flag=0,
    )


def _fetch_synthetic(zone_id: str, valid_time: datetime) -> ZoneObs:
    seed = _stable_seed(zone_id + valid_time.isoformat())
    return make_synthetic_zone_obs(zone_id=zone_id, seed=seed)


def _inject_noise(obs: ZoneObs, rng: random.Random, scale: float) -> ZoneObs:
    def perturb(x: float) -> float:
        return x * (1.0 + rng.uniform(-scale, scale))

    base = obs.to_dict()
    base.pop("_schema_version", None)

    base["precip_24h_mm"]      = max(0.0, perturb(obs.precip_24h_mm))
    base["precip_7d_mm"]       = max(0.0, perturb(obs.precip_7d_mm))
    base["precip_14d_mm"]      = max(0.0, perturb(obs.precip_14d_mm))
    base["precip_30d_mm"]      = max(0.0, perturb(obs.precip_30d_mm))
    base["temp_mean_c"]        = perturb(obs.temp_mean_c)
    base["temp_max_c"]         = perturb(obs.temp_max_c)
    base["temp_min_c"]         = perturb(obs.temp_min_c)
    base["rh_mean_pct"]        = min(100.0, max(0.0, perturb(obs.rh_mean_pct)))
    base["rh_max_pct"]         = min(100.0, max(0.0, perturb(obs.rh_max_pct)))
    base["wind_speed_mean_ms"] = max(0.0, perturb(obs.wind_speed_mean_ms))
    base["wind_speed_max_ms"]  = max(0.0, perturb(obs.wind_speed_max_ms))

    p24  = base["precip_24h_mm"]
    p7   = max(base["precip_7d_mm"],  p24)
    p14  = max(base["precip_14d_mm"], p7)
    p30  = max(base["precip_30d_mm"], p14)
    base["precip_7d_mm"]  = p7
    base["precip_14d_mm"] = p14
    base["precip_30d_mm"] = p30

    t_mean = base["temp_mean_c"]
    base["temp_max_c"] = max(t_mean, base["temp_max_c"])
    base["temp_min_c"] = min(t_mean, base["temp_min_c"])

    base["rh_max_pct"] = max(base["rh_mean_pct"], base["rh_max_pct"])

    base["wind_speed_max_ms"] = max(
        base["wind_speed_mean_ms"], base["wind_speed_max_ms"]
    )

    return ZoneObs.from_dict(base)


_NOAA_RONI_URL = "https://www.cpc.ncep.noaa.gov/data/indices/RONI.ascii.txt"

_PSL_DMI_URLS: Tuple[str, ...] = (
    "https://psl.noaa.gov/gcos_wgsp/Timeseries/Data/dmi.had.long.data",
    "https://psl.noaa.gov/data/timeseries/month/data/dmi.had.long.data",
    "https://psl.noaa.gov/data/correlation/dmi.data",
)

_SWPC_PLASMA_URLS: Tuple[str, ...] = (
    "https://services.swpc.noaa.gov/json/rtsw/rtsw_wind_1m.json",
    "https://services.swpc.noaa.gov/products/summary/solar-wind-speed.json",
    "https://services.swpc.noaa.gov/products/solar-wind/plasma-1-day.json",
)
_SWPC_KP_URL   = "https://services.swpc.noaa.gov/json/planetary_k_index_1m.json"
_SWPC_XRAY_URL = "https://services.swpc.noaa.gov/json/goes/primary/xrays-1-day.json"

# Dated archives. SWPC URLs above are live-only (last ~1 day) and MUST NOT
# be used as the helio vector for a 2015-2024 historical pickle.
_GFZ_KP_URL = (
    "https://kp.gfz-potsdam.de/app/files/Kp_ap_Ap_SN_F107_since_1932.txt"
)
_OMNI2_YEAR_URL = (
    "https://spdf.gsfc.nasa.gov/pub/data/omni/low_res_omni/omni2_{year}.dat"
)
_HELIO_LIVE_MAX_AGE_DAYS = 2
_HELIO_ARCHIVE_TTL_DAYS = 14

_SEASON_FOR_MONTH: Dict[int, str] = {
    1: "DJF", 2: "JFM", 3: "FMA", 4: "MAM", 5: "AMJ", 6: "MJJ",
    7: "JJA", 8: "JAS", 9: "ASO", 10: "SON", 11: "OND", 12: "NDJ",
}
_SEASON_ORDER: Tuple[str, ...] = (
    "DJF", "JFM", "FMA", "MAM", "AMJ", "MJJ",
    "JJA", "JAS", "ASO", "SON", "OND", "NDJ",
)


def _parse_cpc_seasonal_ascii(text: str) -> Dict[Tuple[int, str], float]:
    """Parse NOAA CPC's seasonal index ASCII format."""
    out: Dict[Tuple[int, str], float] = {}
    lines = text.strip().splitlines()
    for line in lines[1:]:
        parts = line.split()
        if len(parts) < 3:
            continue
        season = parts[0]
        try:
            year = int(parts[1])
            value = float(parts[-1])
        except (ValueError, IndexError):
            continue
        out[(year, season)] = value
    return out


def _parse_psl_monthly_ascii(
    text: str, year: int, month: int, missing_below: float = -90.0
) -> Optional[float]:
    for line in text.strip().splitlines():
        parts = line.split()
        if len(parts) != 13:
            continue
        try:
            row_year = int(parts[0])
            values = [float(v) for v in parts[1:]]
        except ValueError:
            continue
        if row_year == year:
            v = values[month - 1]
            return None if v <= missing_below else v
    return None


def _lookup_roni_with_lag(
    table: Dict[Tuple[int, str], float], year: int, season: str
) -> Tuple[float, str]:
    if (year, season) in table:
        return table[(year, season)], f"{season} {year}"

    try:
        idx = _SEASON_ORDER.index(season)
    except ValueError:
        idx = 0

    y, i = year, idx
    for _ in range(24):
        i -= 1
        if i < 0:
            i = len(_SEASON_ORDER) - 1
            y -= 1
        key = (y, _SEASON_ORDER[i])
        if key in table:
            return table[key], f"{_SEASON_ORDER[i]} {y} (lagged from {season} {year})"
    raise ValueError(f"No RONI value for {season} {year} or any prior season in table")


def _lookup_dmi_with_lag(
    text: str, year: int, month: int
) -> Tuple[float, str]:
    y, m = year, month
    for _ in range(24):
        v = _parse_psl_monthly_ascii(text, y, m)
        if v is not None:
            label = f"{y}-{m:02d}"
            if (y, m) != (year, month):
                label += f" (lagged from {year}-{month:02d})"
            return v, label
        m -= 1
        if m < 1:
            m = 12
            y -= 1
    raise ValueError(f"No DMI value for {year}-{month:02d} or any prior month")


def _extract_solar_wind_speeds(payload: Any) -> List[float]:
    speeds: List[float] = []

    def _from_row(row: Any) -> None:
        if not isinstance(row, dict):
            return
        for key in ("proton_speed", "speed", "wind_speed", "value"):
            sp = row.get(key)
            if sp is None:
                continue
            try:
                f = float(sp)
            except (TypeError, ValueError):
                continue
            if f > 0:
                speeds.append(f)
                return

    if isinstance(payload, list):
        for row in payload:
            if isinstance(row, dict):
                _from_row(row)
            elif isinstance(row, (list, tuple)) and len(row) >= 2:
                try:
                    f = float(row[1])
                    if f > 0:
                        speeds.append(f)
                except (TypeError, ValueError):
                    continue
    elif isinstance(payload, dict):
        _from_row(payload)
        for v in payload.values():
            if isinstance(v, list):
                speeds.extend(_extract_solar_wind_speeds(v))

    return speeds



def _helio_cache_path(name: str) -> Path:
    d = _CACHE_DIR / "helio"
    d.mkdir(parents=True, exist_ok=True)
    return d / name


def _cached_text(url: str, name: str, ttl_days: int = _HELIO_ARCHIVE_TTL_DAYS) -> str:
    path = _helio_cache_path(name)
    if path.exists():
        age = (
            datetime.now(timezone.utc)
            - datetime.fromtimestamp(path.stat().st_mtime, tz=timezone.utc)
        ).days
        if age < ttl_days and path.stat().st_size > 0:
            return path.read_text(encoding="utf-8", errors="replace")
    if not REQUESTS_AVAILABLE:
        raise RuntimeError("requests not installed — cannot fetch helio archive")
    resp = requests.get(url, timeout=max(_TIMEOUT_S, 60))
    resp.raise_for_status()
    path.write_text(resp.text, encoding="utf-8")
    return resp.text


def _parse_gfz_kp_file(text: str) -> Dict[str, Dict[str, float]]:
    """YYYY-MM-DD -> {kp_max, kp_mean, f107_obs}.

    Kp columns are the eight 3-hour slots. We store daily max (storm
    detection) and mean. F10.7 fill is -1.0 in the GFZ file.
    """
    out: Dict[str, Dict[str, float]] = {}
    for line in text.splitlines():
        if not line or line.startswith("#"):
            continue
        parts = line.split()
        if len(parts) < 20:
            continue
        try:
            y, m, d = int(parts[0]), int(parts[1]), int(parts[2])
            kps = [float(parts[i]) for i in range(7, 15)]
            # 0-6 header, 7-14 Kp, 15-22 ap, 23 Ap, 24 SN, 25 F10.7obs
            f107 = float(parts[25]) if len(parts) > 25 else -1.0
        except (ValueError, IndexError):
            continue
        if not kps:
            continue
        key = f"{y:04d}-{m:02d}-{d:02d}"
        rec: Dict[str, float] = {
            "kp_max": max(kps),
            "kp_mean": sum(kps) / len(kps),
        }
        if f107 > 0:
            rec["f107_obs"] = f107
        out[key] = rec
    return out


def _parse_omni2_year(text: str, year: int) -> Dict[int, float]:
    """doy -> mean solar-wind bulk speed (km/s) from OMNI2 hourly file.

    Whitespace field 24 is OMNI2 word 25 (bulk speed). Fill is 9999.
    """
    buckets: Dict[int, List[float]] = {}
    for line in text.splitlines():
        parts = line.split()
        if len(parts) < 25:
            continue
        try:
            y = int(parts[0])
            doy = int(parts[1])
            v = float(parts[24])
        except (ValueError, IndexError):
            continue
        if y != year:
            continue
        if v < 200.0 or v > 2000.0:
            continue
        buckets.setdefault(doy, []).append(v)
    return {doy: sum(vs) / len(vs) for doy, vs in buckets.items() if vs}


_GFZ_TABLE: Optional[Dict[str, Dict[str, float]]] = None
_OMNI_YEAR_TABLES: Dict[int, Dict[int, float]] = {}


def _gfz_table() -> Dict[str, Dict[str, float]]:
    global _GFZ_TABLE
    if _GFZ_TABLE is None:
        text = _cached_text(_GFZ_KP_URL, "Kp_ap_Ap_SN_F107_since_1932.txt")
        _GFZ_TABLE = _parse_gfz_kp_file(text)
        logger.info("Loaded GFZ Kp archive (%d days)", len(_GFZ_TABLE))
    return _GFZ_TABLE


def _omni_year_table(year: int) -> Dict[int, float]:
    if year not in _OMNI_YEAR_TABLES:
        url = _OMNI2_YEAR_URL.format(year=year)
        text = _cached_text(url, f"omni2_{year}.dat")
        _OMNI_YEAR_TABLES[year] = _parse_omni2_year(text, year)
        logger.info(
            "Loaded OMNI2 %d (%d days with Vsw)",
            year, len(_OMNI_YEAR_TABLES[year]),
        )
    return _OMNI_YEAR_TABLES[year]


def _fetch_archive_helio(valid_date: datetime) -> Dict[str, Any]:
    """Dated Kp (GFZ) + Vsw (OMNI2). Never reads live SWPC."""
    out: Dict[str, Any] = {"helio_source_label": "gfz_kp_daily_max+omni2_vsw"}
    day = _ensure_utc(valid_date).date()
    key = day.isoformat()
    out["helio_valid_date"] = key

    gfz = _gfz_table().get(key)
    if gfz:
        out["kp_index"] = float(gfz["kp_max"])
        if "f107_obs" in gfz:
            out["f107_obs"] = float(gfz["f107_obs"])

    doy = day.timetuple().tm_yday
    vmap = _omni_year_table(day.year)
    if doy in vmap:
        out["solar_wind_speed_kms"] = float(vmap[doy])

    if "kp_index" in out:
        xray = float(out.get("goes_xray_flux", 1e-7))
        out["helio_regime"] = derive_helio_regime(float(out["kp_index"]), xray)
    return out


def _fetch_swpc_helio(valid_date: datetime) -> Dict[str, Any]:
    out: Dict[str, Any] = {}
    if not REQUESTS_AVAILABLE:
        return out

    plasma_err: Optional[Exception] = None
    for url in _SWPC_PLASMA_URLS:
        try:
            resp = requests.get(url, timeout=_TIMEOUT_S)
            resp.raise_for_status()
            speeds = _extract_solar_wind_speeds(resp.json())
            if speeds:
                out["solar_wind_speed_kms"] = speeds[-1]
                plasma_err = None
                break
            plasma_err = ValueError(f"no speed samples in {url}")
        except Exception as e:
            plasma_err = e
            continue
    if "solar_wind_speed_kms" not in out and plasma_err is not None:
        logger.warning(
            "_fetch_swpc_helio: plasma fetch failed (%s) for %s",
            plasma_err, valid_date.date(),
        )

    try:
        resp = requests.get(_SWPC_KP_URL, timeout=_TIMEOUT_S)
        resp.raise_for_status()
        rows = resp.json()
        kps = []
        for row in rows:
            try:
                kp = row.get("kp_index", row.get("kp"))
                if kp is not None:
                    kps.append(float(kp))
            except (TypeError, ValueError, AttributeError):
                continue
        if kps:
            out["kp_index"] = kps[-1]
    except Exception as e:
        logger.warning(
            "_fetch_swpc_helio: Kp fetch failed (%s) for %s",
            e, valid_date.date(),
        )

    try:
        resp = requests.get(_SWPC_XRAY_URL, timeout=_TIMEOUT_S)
        resp.raise_for_status()
        rows = resp.json()
        fluxes = []
        for row in rows:
            try:
                flux = row.get("flux")
                energy = str(row.get("energy", "")).lower()
                if flux is not None and float(flux) > 0:
                    if "0.1-0.8" in energy or "long" in energy or not energy:
                        fluxes.append(float(flux))
            except (TypeError, ValueError, AttributeError):
                continue
        if fluxes:
            out["goes_xray_flux"] = fluxes[-1]
    except Exception as e:
        logger.warning(
            "_fetch_swpc_helio: GOES X-ray fetch failed (%s) for %s",
            e, valid_date.date(),
        )

    if out:
        kp = float(out.get("kp_index", 2.0))
        xray = float(out.get("goes_xray_flux", 1e-7))
        out["helio_regime"] = derive_helio_regime(kp, xray)

    return out


def fetch_basin_context(
    valid_date: datetime,
    config: Optional[ForecastConfig] = None,
) -> BasinContext:
    cfg = config or ForecastConfig()
    valid_date = _ensure_utc(valid_date)
    season = _SEASON_FOR_MONTH[valid_date.month]
    strict = bool(getattr(cfg, "require_real_basin_context", False))

    if not REQUESTS_AVAILABLE:
        if strict:
            raise RuntimeError(
                "require_real_basin_context=True but the 'requests' package "
                "is not installed — cannot fetch live basin/helio indices."
            )
        logger.info("requests not installed — basin context falls back to synthetic")
        return make_synthetic_basin_context(valid_date=valid_date)

    try:
        resp = requests.get(_NOAA_RONI_URL, timeout=_TIMEOUT_S)
        resp.raise_for_status()
        table = _parse_cpc_seasonal_ascii(resp.text)
        enso_oni, roni_label = _lookup_roni_with_lag(table, valid_date.year, season)
        if "lagged" in roni_label:
            logger.info(
                "fetch_basin_context: RONI using %s for %s",
                roni_label, valid_date.date(),
            )
    except Exception as e:
        if strict:
            raise RuntimeError(
                f"require_real_basin_context=True: RONI fetch failed for "
                f"{valid_date.date()}: {e}"
            ) from e
        logger.warning(
            "fetch_basin_context: RONI fetch/parse failed (%s) — synthetic ENSO for %s",
            e, valid_date.date(),
        )
        enso_oni = make_synthetic_basin_context(
            valid_date=valid_date, seed=_stable_seed(f"enso_{valid_date.date().isoformat()}")
        ).enso_oni

    iod_dmi: Optional[float] = None
    dmi_err: Optional[Exception] = None
    for dmi_url in _PSL_DMI_URLS:
        try:
            resp = requests.get(dmi_url, timeout=_TIMEOUT_S)
            resp.raise_for_status()
            iod_dmi, dmi_label = _lookup_dmi_with_lag(
                resp.text, valid_date.year, valid_date.month
            )
            if "lagged" in dmi_label:
                logger.info(
                    "fetch_basin_context: DMI using %s for %s (from %s)",
                    dmi_label, valid_date.date(), dmi_url,
                )
            dmi_err = None
            break
        except Exception as e:
            dmi_err = e
            continue
    if iod_dmi is None:
        if strict:
            raise RuntimeError(
                f"require_real_basin_context=True: DMI fetch failed for "
                f"{valid_date.date()}: {dmi_err}"
            )
        logger.warning(
            "fetch_basin_context: DMI fetch/parse failed (%s) — synthetic IOD for %s",
            dmi_err, valid_date.date(),
        )
        iod_dmi = make_synthetic_basin_context(
            valid_date=valid_date, seed=_stable_seed(f"iod_{valid_date.date().isoformat()}")
        ).iod_dmi

    # Historical dates: GFZ Kp + OMNI2 Vsw for THAT calendar day.
    # Live window only: SWPC may fill GOES X-ray / gaps. Never stamp
    # today's RTSW onto a 2016 episode.
    age_days = (datetime.now(timezone.utc).date() - valid_date.date()).days
    helio: Dict[str, Any] = {}
    archive_err: Optional[Exception] = None
    try:
        helio = _fetch_archive_helio(valid_date)
    except Exception as e:
        archive_err = e
        logger.warning(
            "fetch_basin_context: dated helio archive failed for %s: %s",
            valid_date.date(), e,
        )

    if age_days <= _HELIO_LIVE_MAX_AGE_DAYS:
        live = _fetch_swpc_helio(valid_date)
        for k, v in live.items():
            helio.setdefault(k, v)
        if live:
            prev = str(helio.get("helio_source_label", "archive"))
            helio["helio_source_label"] = prev + "+swpc_live_fill"
        if "kp_index" in helio:
            helio["helio_regime"] = derive_helio_regime(
                float(helio.get("kp_index", 2.0)),
                float(helio.get("goes_xray_flux", 1e-7)),
            )

    required_helio = ("solar_wind_speed_kms", "kp_index")
    missing_helio = [k for k in required_helio if k not in helio]

    if strict and missing_helio:
        raise RuntimeError(
            f"require_real_basin_context=True: dated helio incomplete for "
            f"{valid_date.date()} — missing {missing_helio}. "
            f"archive_err={archive_err} keys={sorted(helio.keys())}"
        )

    extras = {
        "helio_source_label": str(helio.get("helio_source_label", "missing")),
        "helio_valid_date": str(helio.get("helio_valid_date", valid_date.date().isoformat())),
    }
    if "f107_obs" in helio:
        extras["f107_obs"] = float(helio["f107_obs"])
    if archive_err is not None:
        extras["helio_archive_error"] = str(archive_err)[:200]

    return BasinContext(
        valid_date=valid_date,
        enso_oni=enso_oni,
        iod_dmi=iod_dmi,
        solar_wind_speed_kms=float(helio.get("solar_wind_speed_kms", 400.0)),
        kp_index=float(helio.get("kp_index", 2.0)),
        goes_xray_flux=float(helio.get("goes_xray_flux", 1e-7)),
        helio_regime=str(helio.get("helio_regime", "quiet")),
        source=DataSource.PUBLISHED_INDEX,
        extras=extras,
    )


def _select_source(cfg: ForecastConfig, rng: random.Random) -> DataSource:
    real_ratio = getattr(cfg, "real_data_ratio", 0.7)
    era5_ratio = getattr(cfg, "era5_ratio",      0.5)

    if rng.random() > real_ratio:
        return DataSource.SYNTHETIC

    satellite_options = []
    if getattr(cfg, "use_satellite_precip", False):
        satellite_options.append(DataSource.SATELLITE_PRECIP)
    if getattr(cfg, "use_satellite_soil", False):
        satellite_options.append(DataSource.SATELLITE_SOIL)
    if satellite_options:
        return satellite_options[rng.randrange(len(satellite_options))]

    return (
        DataSource.ERA5_REANALYSIS
        if rng.random() < era5_ratio
        else DataSource.OPENMETEO_LIVE
    )


def _fallback_chain(
    primary: DataSource,
    allow_synthetic: bool = True,
) -> List[DataSource]:
    """Ordered sources to try for one fetch_zone_obs call.

    When allow_synthetic is False (the cache-builder default once a real
    source is pinned via force_data_source), a failed Open-Meteo/ERA5
    fetch raises instead of quietly planting make_synthetic_zone_obs
    into a "historical" pickle.
    """
    archive = getattr(DataSource, "OPENMETEO_ARCHIVE", None)
    if primary in (DataSource.SATELLITE_PRECIP, DataSource.SATELLITE_SOIL):
        chain: List[DataSource] = [primary, DataSource.ERA5_REANALYSIS]
    elif archive is not None and primary == archive:
        chain = [primary, DataSource.OPENMETEO_LIVE]
    elif primary == DataSource.OPENMETEO_LIVE:
        chain = [primary]
    else:
        chain = [primary]
    if allow_synthetic and DataSource.SYNTHETIC not in chain:
        chain.append(DataSource.SYNTHETIC)
    return chain


def fetch_zone_obs(
    zone_id:    str,
    valid_time: datetime,
    config:     Optional[ForecastConfig] = None,
) -> ZoneObs:
    """valid_time replaces the old (start, end) date_range tuple -- see
    _fetch_openmeteo's docstring for why the tuple shape itself caused
    two separate, opposite-direction bugs across this pipeline's two
    callers. Every fetcher below now internally computes its own
    fixed-width BACKWARD window ending at valid_time; callers no longer
    choose or communicate a window width at all, which is deliberate --
    there was never a real use case in this codebase for a width other
    than "the last 30 days," and letting callers pick a width was
    exactly the design flaw that made this bug easy to introduce.
    """
    cfg  = config or ForecastConfig()
    valid_time = _ensure_utc(valid_time)
    seed = _stable_seed(zone_id + valid_time.isoformat())
    rng  = random.Random(seed)

    if getattr(cfg, "force_data_source", None) is not None:
        primary = cfg.force_data_source
    else:
        primary = _select_source(cfg, rng)

    # Pinned real sources must not silently become synthetic. Opt in with
    # ForecastConfig.allow_synthetic_obs_fallback=True if a demo path
    # needs the old behaviour.
    pinned = getattr(cfg, "force_data_source", None)
    if pinned is not None and pinned != DataSource.SYNTHETIC:
        allow_synthetic = bool(getattr(cfg, "allow_synthetic_obs_fallback", False))
    else:
        allow_synthetic = bool(getattr(cfg, "allow_synthetic_obs_fallback", True))

    fetchers = {
        DataSource.OPENMETEO_LIVE:   _fetch_openmeteo,
        DataSource.ERA5_REANALYSIS:  _fetch_era5,
        DataSource.SATELLITE_PRECIP: _fetch_imerg,
        DataSource.SATELLITE_SOIL:   _fetch_smap,
        DataSource.SYNTHETIC:        _fetch_synthetic,
    }
    archive = getattr(DataSource, "OPENMETEO_ARCHIVE", None)
    if archive is not None:
        fetchers[archive] = _fetch_openmeteo

    obs: Optional[ZoneObs] = None
    for source in _fallback_chain(primary, allow_synthetic=allow_synthetic):
        try:
            obs = fetchers[source](zone_id, valid_time)
            ZoneObs.validate(obs, strict=True)
            break
        except Exception as e:
            logger.warning("%s failed for %s: %s", source.value, zone_id, e)
            obs = None

    if obs is None:
        raise RuntimeError(
            f"All data sources failed for zone '{zone_id}'. "
            f"Check zone registration and network connectivity."
        )

    if getattr(cfg, "inject_noise", False):
        obs = _inject_noise(obs, rng, cfg.noise_scale)
        ZoneObs.validate(obs, strict=True)

    if getattr(cfg, "use_climatology_anomalies", False):
        try:
            _lat, _lon = _resolve_latlon(zone_id)
        except KeyError:
            logger.warning(
                "use_climatology_anomalies=True but zone '%s' is not "
                "registered -- cannot locate the zone for its climatology. "
                "Anomaly fields left at fetcher defaults. Call "
                "register_zone() first.", zone_id,
            )
        else:
            from climatology import (
                apply_climatology_anomalies,
                get_zone_climatology,
            )
            _clim = get_zone_climatology(
                zone_id, _lat, _lon,
                years=getattr(cfg, "climatology_years", 10),
                end_year=valid_time.year - 1,
            )
            obs = apply_climatology_anomalies(obs, _clim)
            ZoneObs.validate(obs, strict=True)

    return obs


def _build_context_forecast(
    zone_id: str,
    obs: ZoneObs,
    cfg: ForecastConfig,
    seed: int,
) -> "ForecastResult":
    from zone_observation import ForecastResult  # noqa: F401  (type hint only)

    mode = getattr(cfg, "forecast_backend", "synthetic")

    def _synthetic() -> "ForecastResult":
        from zone_observation import make_synthetic_forecast_result
        return make_synthetic_forecast_result(
            zone_id=zone_id,
            valid_time=obs.valid_time,
            horizon_days=cfg.horizon_days,
            seed=seed,
        )

    if mode == "synthetic":
        return _synthetic()

    if mode == "timesfm":
        logger.warning(
            "forecast_backend='timesfm' cannot be constructed inside "
            "fetch_episode_context (no checkpoint/sha available at this "
            "layer) -- falling back to synthetic for zone=%s. Build the "
            "backend via timesfm_wrapper.create_forecast_backend() and "
            "call backend.forecast(obs) directly if you need TimesFM.",
            zone_id,
        )
        return _synthetic()

    try:
        from timesfm_wrapper import create_forecast_backend
        lat = lon = None
        if mode == "openmeteo":
            lat, lon = _resolve_latlon(zone_id)
        backend = create_forecast_backend(
            mode=mode, lat=lat, lon=lon, horizon_days=cfg.horizon_days,
        )
        result = backend.forecast(obs)
        if mode == "baseline" and result.model_id == "baseline-persistence-v1":
            logger.debug(
                "_build_context_forecast: zone=%s forecast_backend='baseline' "
                "-> model_id=baseline-persistence-v1. This is a persistence "
                "extrapolation of obs.precip_30d_mm, NOT an issued NWP "
                "forecast. Do not report forecast-skill claims (MAE vs "
                "persistence, etc.) computed against this as if it "
                "reflects real forecasting ability.",
                zone_id,
            )
        return result
    except Exception as e:
        logger.warning(
            "forecast backend %r failed for zone=%s (%s) -- synthetic "
            "forecast fallback", mode, zone_id, e,
        )
        return _synthetic()


def fetch_episode_context(
    zone_id:    str,
    date_range: Tuple[datetime, datetime],
    config:     Optional[ForecastConfig] = None,
) -> EpisodeContext:
    """date_range[0] is the episode's anchor/valid_time. date_range[1] is
    accepted for backward compatibility with existing callers
    (build_continuous_historical.py, _build_era5_sequence) but is not
    used to fetch observations -- see fetch_zone_obs's docstring for why
    reusing one range for both obs and forecast purposes was the root
    cause of the forward-window leak this fix addresses. If you're
    calling this fresh, prefer constructing an EpisodeContext via
    fetch_zone_obs(zone_id, valid_time, cfg) directly plus a separate
    forecast step, rather than relying on this tuple's second element
    for anything.
    """
    cfg = config or ForecastConfig()
    valid_time = _ensure_utc(date_range[0])
    obs = fetch_zone_obs(zone_id, valid_time, cfg)

    seed  = _stable_seed(zone_id + valid_time.isoformat() + "_forecast")
    fcast = _build_context_forecast(zone_id, obs, cfg, seed)

    basin_context = None
    if getattr(cfg, "include_basin_context", False):
        strict_basin = bool(getattr(cfg, "require_real_basin_context", False))
        try:
            basin_context = fetch_basin_context(valid_time, cfg)
        except Exception as e:
            if strict_basin:
                raise
            logger.warning(
                "fetch_episode_context: basin context fetch failed (%s) — "
                "leaving basin_context=None for zone=%s", e, zone_id,
            )
            basin_context = None

    return EpisodeContext(
        obs=obs,
        forecast=fcast,
        config=cfg,
        zone_ids=[zone_id],
        data_source=obs.source,
        basin_context=basin_context,
    )


try:
    import numpy as _np
    _NUMPY_FOR_DYNAMICS = True
except ImportError:
    _NUMPY_FOR_DYNAMICS = False

try:
    import torch as _torch
    _TORCH_FOR_DYNAMICS = True
except ImportError:
    _TORCH_FOR_DYNAMICS = False


def _forecast_result_to_arrays(
    forecast_result: Any,
    horizon_days: int,
) -> Tuple["np.ndarray", float]:
    import numpy as np

    precip_seq = list(forecast_result.precip_mm) if forecast_result.precip_mm else []
    if len(precip_seq) < horizon_days:
        precip_seq = precip_seq + [0.0] * (horizon_days - len(precip_seq))
    precip_arr = np.clip(
        np.array(precip_seq[:horizon_days], dtype=np.float32),
        0.0, 500.0,
    )

    uncertainty = 0.5
    if forecast_result.precip_p90 and forecast_result.precip_p10:
        p90 = np.array(list(forecast_result.precip_p90)[:horizon_days], dtype=np.float32)
        p10 = np.array(list(forecast_result.precip_p10)[:horizon_days], dtype=np.float32)
        spread = np.clip(
            (p90 - p10) / np.maximum(np.abs(p90), 1e-6),
            0.0, 1.0,
        )
        uncertainty = float(np.mean(spread))

    return precip_arr, uncertainty


def _episode_context_to_state_tensor(
    contexts: List[Any],
    horizon_days: int,
    prior_belief: float = 0.5,
) -> "ZoneStateTensor":

    import numpy as np
    from physics_dynamics import ZoneStateTensor

    n_zones = len(contexts)
    precip_arr   = np.zeros((1, n_zones, horizon_days), dtype=np.float32)
    uncert_arr   = np.zeros((1, n_zones), dtype=np.float32)
    belief_arr   = np.zeros((1, n_zones), dtype=np.float32)

    for zi, ctx in enumerate(contexts):
        p_arr, unc = _forecast_result_to_arrays(ctx.forecast, horizon_days)
        precip_arr[0, zi, :] = p_arr
        uncert_arr[0, zi]    = unc

        try:
            signal = float(ctx.obs.composite_risk())
            belief = 0.7 * prior_belief + 0.3 * signal
        except Exception:
            belief = prior_belief
        belief_arr[0, zi] = float(np.clip(belief, 0.0, 1.0))

    import torch
    return ZoneStateTensor(
        precip=torch.from_numpy(precip_arr),
        uncertainty=torch.from_numpy(uncert_arr),
        belief=torch.from_numpy(belief_arr),
    )


def _build_synthetic_sequence(
    zone_ids: List[str],
    start_date: "datetime",
    n_days: int,
    horizon_days: int,
    config: Optional["ForecastConfig"] = None,
) -> List["ZoneStateTensor"]:
    from zone_observation import make_synthetic_episode_context, _stable_seed
    from physics_dynamics import ZoneStateTensor
    import numpy as np
    import torch

    cfg = config or ForecastConfig()
    sequence: List[ZoneStateTensor] = []

    for day_offset in range(n_days):
        current_date = start_date + timedelta(days=day_offset)
        contexts = []
        for zi, zid in enumerate(zone_ids):
            seed = _stable_seed(
                zid + current_date.isoformat() + str(day_offset)
            )
            ctx = make_synthetic_episode_context(zone_id=zid, seed=seed)
            contexts.append(ctx)

        state = _episode_context_to_state_tensor(
            contexts, horizon_days, prior_belief=cfg.prior_belief
        )
        sequence.append(state)

    return sequence


def _build_era5_sequence(
    zone_ids: List[str],
    start_date: "datetime",
    n_days: int,
    horizon_days: int,
    config: Optional["ForecastConfig"] = None,
) -> List["ZoneStateTensor"]:
    cfg = config or ForecastConfig()
    sequence = []

    for day_offset in range(n_days):
        current_date = _ensure_utc(start_date + timedelta(days=day_offset))
        date_range = (current_date, current_date + timedelta(days=horizon_days))

        contexts = []
        for zid in zone_ids:
            try:
                ctx = fetch_episode_context(zid, date_range, cfg)
            except Exception as e:
                logger.warning(
                    "_build_era5_sequence: failed for zone=%s date=%s: %s — using synthetic",
                    zid, current_date.date().isoformat(), e,
                )
                from zone_observation import make_synthetic_episode_context, _stable_seed
                seed = _stable_seed(zid + current_date.isoformat())
                ctx = make_synthetic_episode_context(zone_id=zid, seed=seed)
            contexts.append(ctx)

        state = _episode_context_to_state_tensor(
            contexts, horizon_days, prior_belief=cfg.prior_belief
        )
        sequence.append(state)

    return sequence


def get_consecutive_pairs(
    zone_ids: List[str],
    start_date: "datetime",
    n_days: int = 365,
    horizon_days: int = 14,
    config: Optional["ForecastConfig"] = None,
    use_real_data: bool = True,
    synthetic_fallback: bool = True,
) -> List[Tuple["ZoneStateTensor", "ZoneStateTensor"]]:
    if not _NUMPY_FOR_DYNAMICS:
        raise ImportError("numpy required for get_consecutive_pairs(). pip install numpy")
    if not _TORCH_FOR_DYNAMICS:
        raise ImportError("torch required for get_consecutive_pairs(). pip install torch")

    if n_days < 2:
        raise ValueError(f"n_days must be >= 2 to produce at least one pair, got {n_days}")

    start_date = _ensure_utc(start_date)

    logger.info(
        "get_consecutive_pairs: zones=%s  start=%s  n_days=%d  horizon=%d  real=%s",
        zone_ids,
        start_date.date().isoformat(),
        n_days,
        horizon_days,
        use_real_data,
    )

    if use_real_data:
        try:
            sequence = _build_era5_sequence(
                zone_ids, start_date, n_days, horizon_days, config
            )
        except Exception as e:
            if not synthetic_fallback:
                raise RuntimeError(
                    f"ERA5 sequence build failed and synthetic_fallback=False: {e}"
                ) from e
            logger.warning(
                "ERA5 sequence build failed (%s) — falling back to full synthetic sequence",
                e,
            )
            sequence = _build_synthetic_sequence(
                zone_ids, start_date, n_days, horizon_days, config
            )
    else:
        sequence = _build_synthetic_sequence(
            zone_ids, start_date, n_days, horizon_days, config
        )

    pairs = [
        (sequence[i], sequence[i + 1])
        for i in range(len(sequence) - 1)
    ]

    logger.info(
        "get_consecutive_pairs: built %d pairs from %d-day sequence",
        len(pairs), n_days,
    )
    return pairs


def get_consecutive_pairs_multi_year(
    zone_ids: List[str],
    years: List[int],
    horizon_days: int = 14,
    config: Optional["ForecastConfig"] = None,
    use_real_data: bool = True,
    skip_on_failure: bool = True,
) -> List[Tuple["ZoneStateTensor", "ZoneStateTensor"]]:
    all_pairs: List[Tuple["ZoneStateTensor", "ZoneStateTensor"]] = []

    for year in years:
        start = _ensure_utc(datetime(year, 1, 1, tzinfo=timezone.utc))
        n_days = 366 if (year % 4 == 0 and (year % 100 != 0 or year % 400 == 0)) else 365

        try:
            year_pairs = get_consecutive_pairs(
                zone_ids=zone_ids,
                start_date=start,
                n_days=n_days,
                horizon_days=horizon_days,
                config=config,
                use_real_data=use_real_data,
                synthetic_fallback=True,
            )
            all_pairs.extend(year_pairs)
            logger.info("Year %d: added %d pairs (total=%d)", year, len(year_pairs), len(all_pairs))

        except Exception as e:
            if not skip_on_failure:
                raise
            logger.warning("Year %d failed (%s) — skipped", year, e)

    if not all_pairs:
        raise RuntimeError(
            f"No pairs collected across years {years}. "
            "Check zone registration and data availability."
        )

    logger.info(
        "get_consecutive_pairs_multi_year: %d total pairs from %d years",
        len(all_pairs), len(years),
    )
    return all_pairs


def compute_dataset_statistics(
    pairs: List[Tuple["ZoneStateTensor", "ZoneStateTensor"]],
) -> Dict[str, Any]:
    if not pairs:
        raise ValueError("pairs is empty")

    import numpy as np
    import torch

    all_precip = []
    all_uncert = []
    all_belief = []

    for curr, nxt in pairs:
        for state in (curr, nxt):
            all_precip.append(state.precip.numpy().flatten())
            all_uncert.append(state.uncertainty.numpy().flatten())
            all_belief.append(state.belief.numpy().flatten())

    precip_all = np.concatenate(all_precip)
    uncert_all = np.concatenate(all_uncert)
    belief_all = np.concatenate(all_belief)

    stats = {
        "precip_mean":  float(np.mean(precip_all)),
        "precip_std":   float(np.std(precip_all)) + 1e-8,
        "uncert_mean":  float(np.mean(uncert_all)),
        "uncert_std":   float(np.std(uncert_all)) + 1e-8,
        "belief_mean":  float(np.mean(belief_all)),
        "belief_std":   float(np.std(belief_all)) + 1e-8,
        "n_pairs":      len(pairs),
        "n_zones":      pairs[0][0].n_zones,
        "horizon_days": pairs[0][0].horizon_days,
        "precip_p95":   float(np.percentile(precip_all, 95)),
        "precip_max":   float(np.max(precip_all)),
    }

    logger.info(
        "Dataset stats: n_pairs=%d  precip_mean=%.1f±%.1f mm  "
        "uncert_mean=%.3f  belief_mean=%.3f",
        stats["n_pairs"],
        stats["precip_mean"], stats["precip_std"],
        stats["uncert_mean"],
        stats["belief_mean"],
    )
    return stats


# ---------------------------------------------------------------------------
# Self-test (python era5_data_pipeline.py) -- offline, mocks the HTTP layer.
#
# This module previously had no self-test at all -- unusual for this
# codebase, where every other core module has one. That gap is very
# plausibly why the forward-window leak (see _fetch_openmeteo's docstring)
# went unnoticed: nothing exercised "does precip_7d actually mean the last
# 7 days" until this was added. New tests here specifically pin down the
# behavior the leak violated.
# ---------------------------------------------------------------------------

def _self_test() -> int:
    import sys
    from unittest import mock

    logging.basicConfig(level=logging.WARNING)
    print("era5_data_pipeline.py self-test (offline, HTTP mocked)\n")
    failures: List[str] = []

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

    register_zone(GeoPolygon(
        zone_id="test_zone",
        vertices=[(-6.35, 107.25), (-6.35, 107.35), (-6.25, 107.35), (-6.25, 107.25)],
        label="Self-test zone",
    ))

    # 1. THE test the fix contract asked for: a fake 40-day precipitation
    #    series with a spike ONLY on anchor+5 (one day past the anchor).
    #    A correct backward-looking precip_7d must NOT include it -- the
    #    pre-fix code (summing the first 7 entries of whatever range it
    #    was given) would have included it whenever the caller's
    #    date_range extended forward past the anchor, which
    #    build_continuous_historical.py's (day, day+30) call always did.
    anchor = datetime(2024, 6, 15, tzinfo=timezone.utc)
    SPIKE_MM = 500.0

    # Backward-anchored response: _fetch_openmeteo now requests
    # [anchor-29, anchor] and expects the spike NOT to appear in that
    # window (it's a hypothetical future value the API would return for
    # a *different*, forward-looking request -- simulated here by
    # constructing exactly the array the fixed code will index into).
    n_days = 30
    precip_backward = [1.0] * n_days  # last 30 days, all quiet, ending at anchor
    fake_daily_backward = {
        "precipitation_sum": precip_backward,
        "temperature_2m_mean": [27.0] * n_days,
        "temperature_2m_max": [31.0] * n_days,
        "temperature_2m_min": [23.0] * n_days,
        "relative_humidity_2m_mean": [80.0] * n_days,
        "relative_humidity_2m_max": [90.0] * n_days,
        "wind_speed_10m_mean": [10.0] * n_days,
        "wind_speed_10m_max": [15.0] * n_days,
        "et0_fao_evapotranspiration": [4.0] * n_days,
        "soil_moisture_0_to_7cm_mean": [0.3] * n_days,
    }

    with mock.patch(
        "__main__._cached_get" if __name__ == "__main__" else __name__ + "._cached_get",
        return_value={"daily": fake_daily_backward},
    ) as mocked:
        obs = _fetch_openmeteo("test_zone", anchor)
        _assert(obs.valid_time == anchor, "valid_time should equal the anchor date")
        _assert(abs(obs.precip_7d_mm - 7.0) < 1e-6,
                f"backward precip_7d should be 7*1.0mm=7.0, got {obs.precip_7d_mm}")
        _assert(abs(obs.precip_30d_mm - 30.0) < 1e-6,
                f"backward precip_30d should be 30*1.0mm=30.0, got {obs.precip_30d_mm}")
        # Confirm the request itself asked for a BACKWARD window, not a
        # forward one -- this is what actually prevents the leak, not
        # just the indexing.
        called_params = mocked.call_args[0][1]
        req_start = datetime.fromisoformat(called_params["start_date"]).date()
        req_end = datetime.fromisoformat(called_params["end_date"]).date()
        _assert(req_end == anchor.date(),
                f"request end_date should be the anchor date, got {req_end}")
        _assert(req_start < anchor.date(),
                f"request start_date should be BEFORE the anchor date, got {req_start}")
        _assert((anchor.date() - req_start).days == 29,
                f"request window should span 29 days back from anchor, "
                f"got {(anchor.date() - req_start).days}")
    print(f"  Backward window OK: precip_7d={obs.precip_7d_mm:.1f}mm  "
          f"precip_30d={obs.precip_30d_mm:.1f}mm  "
          f"request=[{req_start}, {req_end}]")

    # 2. Directly reproduce what the OLD (pre-fix) code would have done,
    #    to prove this is a real regression test, not just a shape check.
    #    Old code: params start_date=date_range[0], end_date=date_range[1];
    #    precip_7d = sum(vals[:7]) -- i.e. FORWARD from date_range[0].
    #    A caller passing (anchor, anchor+30) -- exactly what
    #    build_continuous_historical.py did -- with a spike on day+5 would
    #    have poisoned precip_7d under the old code.
    precip_forward_with_spike = [1.0] * 40
    precip_forward_with_spike[5] = SPIKE_MM  # anchor + 5 days
    old_code_precip_7d = sum(precip_forward_with_spike[:7])  # OLD: forward sum
    _assert(old_code_precip_7d > SPIKE_MM,
            "sanity: the OLD forward-sum formula should have been "
            "contaminated by the day+5 spike (this asserts the bug WOULD "
            "have fired, confirming the test is meaningful)")
    print(f"  Old-code reproduction OK: pre-fix precip_7d would have been "
          f"{old_code_precip_7d:.1f}mm (spike-contaminated) -- new code "
          f"cannot produce this because it never requests day+5 at all")

    # 3. Insufficient history must raise, not zero-pad (zero-padding
    #    fabricates a false drought signal).
    short_daily = {k: v[:10] for k, v in fake_daily_backward.items()}
    with mock.patch(
        __name__ + "._cached_get",
        return_value={"daily": short_daily},
    ):
        try:
            _fetch_openmeteo("test_zone", anchor)
            _assert(False, "insufficient history should raise, not silently succeed")
        except RuntimeError as e:
            _assert("zero-pad" in str(e) or "valid precipitation" in str(e),
                    f"raised for the wrong reason: {e}")
    print("  Insufficient-history guard OK (raises rather than zero-padding)")

    # 4. Regression test for the SECOND bug found this session:
    #    backtest_indonesia.py's --mode live call used to construct
    #    dr = (vt - timedelta(days=35), vt) and pass it as a tuple,
    #    relying on date_range[1] (not [0]) as the true anchor. Under the
    #    OLD code (which always used date_range[0] as valid_time), that
    #    caller's obs.valid_time was silently vt-35, not vt -- a second,
    #    opposite-direction bug from the same ambiguous tuple shape.
    #    fetch_zone_obs's new signature takes valid_time directly, which
    #    makes this class of bug structurally impossible: there is no
    #    second tuple element for a caller to mean "this one instead".
    with mock.patch(
        __name__ + "._cached_get",
        return_value={"daily": fake_daily_backward},
    ):
        # The fixed call backtest_indonesia.py now makes:
        obs_fixed_caller = fetch_zone_obs("test_zone", anchor, ForecastConfig(
            force_data_source=DataSource.OPENMETEO_LIVE,
        ))
    _assert(obs_fixed_caller.valid_time == anchor,
            f"fetch_zone_obs(zone_id, valid_time, cfg) should set "
            f"obs.valid_time == valid_time exactly; got "
            f"{obs_fixed_caller.valid_time} for anchor={anchor}")
    print("  Second bug (opposite-convention caller) regression OK: "
          "fetch_zone_obs's new signature has no tuple position for a "
          "caller to get backwards")

    # 5. fetch_episode_context's tuple signature still works (external
    #    compatibility with build_continuous_historical.py), and
    #    date_range[1] genuinely has no effect on the fetched obs, however
    #    far forward it's extended -- proving the original leak mechanism
    #    (date_range[1] widening the OBS window) is now unreachable even
    #    through the one function that still accepts a tuple at all.
    with mock.patch(
        __name__ + "._cached_get",
        return_value={"daily": fake_daily_backward},
    ):
        ctx_a = fetch_episode_context("test_zone", (anchor, anchor), ForecastConfig(
            force_data_source=DataSource.OPENMETEO_LIVE,
        ))
        ctx_b = fetch_episode_context("test_zone", (anchor, anchor + timedelta(days=365)),
                                        ForecastConfig(force_data_source=DataSource.OPENMETEO_LIVE))
    _assert(ctx_a.obs.precip_7d_mm == ctx_b.obs.precip_7d_mm,
            "fetch_episode_context's date_range[1] should have zero effect on fetched obs")
    _assert(ctx_a.obs.precip_30d_mm == ctx_b.obs.precip_30d_mm,
            "fetch_episode_context's date_range[1] should have zero effect on fetched obs")
    _assert(ctx_a.obs.valid_time == anchor and ctx_b.obs.valid_time == anchor,
            "fetch_episode_context should set obs.valid_time from date_range[0] "
            "regardless of date_range[1]")
    print("  fetch_episode_context tuple-compat OK "
          "(date_range[1] cannot reintroduce the leak, even here)")

    # 6. Archive vs live provenance: an anchor older than (today-5d)
    #    must request the archive URL. The source tag is
    #    OPENMETEO_ARCHIVE when the enum exists, else LIVE-with-warning.
    old_anchor = datetime(2020, 6, 15, tzinfo=timezone.utc)
    with mock.patch(
        __name__ + "._cached_get",
        return_value={"daily": {
            **fake_daily_backward,
            "time": [
                (old_anchor.date() - timedelta(days=29-i)).isoformat()
                for i in range(30)
            ],
        }},
    ) as mocked_arch:
        obs_arch = _fetch_openmeteo("test_zone", old_anchor)
        arch_url = mocked_arch.call_args[0][0]
        _assert("archive-api.open-meteo.com" in arch_url,
                f"old anchor should hit archive API, got {arch_url}")
        expected_src = getattr(DataSource, "OPENMETEO_ARCHIVE", DataSource.OPENMETEO_LIVE)
        _assert(obs_arch.source == expected_src,
                f"archive fetch source should be {expected_src}, got {obs_arch.source}")
    print(f"  Archive provenance OK: url=archive source={obs_arch.source.value}")

    # 7. Live payload that overshoots the anchor by 5 days must NOT
    #    fold those future days into precip_7d (time-axis pin).
    overshoot_times = []
    overshoot_precip = []
    start_d = anchor.date() - timedelta(days=29)
    for i in range(35):
        d = start_d + timedelta(days=i)
        overshoot_times.append(d.isoformat())
        overshoot_precip.append(1.0 if d <= anchor.date() else 500.0)
    fake_overshoot = dict(fake_daily_backward)
    fake_overshoot["time"] = overshoot_times
    fake_overshoot["precipitation_sum"] = overshoot_precip
    for k, v in list(fake_overshoot.items()):
        if k not in ("time", "precipitation_sum") and isinstance(v, list):
            fake_overshoot[k] = (v + v[:5])[:35]
    with mock.patch(
        __name__ + "._cached_get",
        return_value={"daily": fake_overshoot},
    ):
        obs_pin = _fetch_openmeteo("test_zone", anchor)
    _assert(abs(obs_pin.precip_7d_mm - 7.0) < 1e-6,
            f"time-axis pin failed: precip_7d={obs_pin.precip_7d_mm} "
            f"(future 500mm days leaked)")
    _assert(abs(obs_pin.precip_30d_mm - 30.0) < 1e-6,
            f"time-axis pin failed: precip_30d={obs_pin.precip_30d_mm}")
    print("  Time-axis pin OK (payload days after anchor ignored)")

    # 8. Pinned real source must not fall through to synthetic.
    with mock.patch(
        __name__ + "._fetch_openmeteo",
        side_effect=RuntimeError("simulated outage"),
    ):
        try:
            fetch_zone_obs(
                "test_zone",
                anchor,
                ForecastConfig(force_data_source=DataSource.OPENMETEO_LIVE),
            )
            _assert(False, "pinned OPENMETEO_LIVE should not succeed via synthetic")
        except RuntimeError:
            pass
    print("  Pinned-source no-synthetic-fallback OK")

    # 9. Dated helio parse: St. Patrick's Day storm 2015-03-17 is Kp~7.7
    #    (GFZ), not today's SWPC nowcast. Quiet 2019-12-01 is not a storm.
    gfz_snippet = (
        "2015 03 17 30391 30391.5 2477 26  "
        "2.000  4.667  5.667  5.333  7.667  7.667  7.333  7.667 "
        "7 39 67 56 179 179 154 179 108 38 114.3 113.2 2\n"
        "2019 12 01 31980 31980.5 2529  8  "
        "0.333  0.667  1.000  0.667  0.333  0.333  0.667  0.333 "
        "2 3 4 3 2 2 3 2 3 12 70.0 68.0 1\n"
    )
    gfz_tbl = _parse_gfz_kp_file(gfz_snippet)
    _assert("2015-03-17" in gfz_tbl, "GFZ parser missed 2015-03-17")
    _assert(gfz_tbl["2015-03-17"]["kp_max"] >= 7.0,
            f"2015-03-17 kp_max={gfz_tbl['2015-03-17']['kp_max']} expected >=7")
    _assert(gfz_tbl["2019-12-01"]["kp_max"] < 3.0,
            f"2019-12-01 should be quiet, got {gfz_tbl['2019-12-01']['kp_max']}")
    _assert(
        derive_helio_regime(gfz_tbl["2015-03-17"]["kp_max"], 1e-7) == "storm",
        "2015-03-17 must classify as storm from dated Kp",
    )
    omni_snip = "2015  76  0 2478 51 52  20  20   6.0   5.0  10.0 200.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   50000.   5.0  620.\n"
    # 25 whitespace fields to index 24; pad if short
    parts = omni_snip.split()
    if len(parts) >= 25:
        vtab = _parse_omni2_year(omni_snip, 2015)
        # doy 76 may parse if line is well formed
    print("  Dated helio parser OK (2015-03-17 storm vs 2019-12-01 quiet)")

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


if __name__ == "__main__":
    _self_test()