Update era5_data_pipeline.py

era5_data_pipeline.py

@@ -1,6 +1,23 @@
 """
 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
@@ -37,6 +54,25 @@ from zone_observation import (
 
 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
 # ---------------------------------------------------------------------------
@@ -337,55 +373,61 @@ def _build_era5_obs_xarray(
     start:   datetime,
     nc_path: Path,
 ) -> ZoneObs:
-    """xarray fallback for _build_era5_obs. Same backward-window,
-    end-indexed convention -- see _build_era5_obs's docstring."""
+    """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 _mean(var: str) -> float:
+        def _windowed(var: str, hours: int) -> "np.ndarray":
             if var not in ds:
-                return 0.0
-            return float(ds[var].values.flatten()[np.isfinite(
-                ds[var].values.flatten()
-            )].mean()) if len(ds[var].values.flatten()) > 0 else 0.0
+                return np.array([])
+            vals = ds[var].values.flatten()[-hours:]
+            return vals[np.isfinite(vals)]
 
-        def _max(var: str) -> float:
-            if var not in ds:
-                return 0.0
-            vals = ds[var].values.flatten()
-            valid = 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 _sum(var: str, scale: float = 1.0) -> float:
-            if var not in ds:
-                return 0.0
-            vals = ds[var].values.flatten()
-            valid = vals[np.isfinite(vals)]
+        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
 
-        temp_mean_c    = _mean("t2m") - 273.15
-        temp_max_c     = _max("t2m")  - 273.15
-        if "t2m" in ds:
-            _t2m_flat = ds["t2m"].values.flatten()
-            _t2m_valid = _t2m_flat[np.isfinite(_t2m_flat)]
-            temp_min_c = float(np.min(_t2m_valid)) - 273.15 if len(_t2m_valid) > 0 else temp_mean_c
-        else:
-            temp_min_c = temp_mean_c
+        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 = _mean("d2m") - 273.15
-        dewpoint_max_c  = _max("d2m")  - 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, _sum("tp", scale=1000.0))
-
-        if "tp" in ds:
-            total_hours = len(ds["tp"].values.flatten())
-        else:
-            total_hours = 0
+        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(
@@ -396,21 +438,16 @@ def _build_era5_obs_xarray(
                 f"widen the ERA5 request window."
             )
 
-        def _window_sum_xr(var: str, hours: int, scale: float = 1.0) -> float:
-            if var not in ds:
-                return 0.0
-            vals = ds[var].values.flatten()[-hours:]
-            valid = vals[np.isfinite(vals)]
-            return float(valid.sum() * scale) if len(valid) > 0 else 0.0
         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 = _mean("u10")
-        v = _mean("v10")
+
+        u = _daily_mean("u10")
+        v = _daily_mean("v10")
         wind_mean = math.sqrt(u**2 + v**2)
-        wind_max  = math.sqrt(_max("u10")**2 + _max("v10")**2)
-        soil_pct  = _mean("swvl1") * 100.0
-        et0_mm    = abs(_sum("pev")) * 1000.0
+        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,
@@ -493,8 +530,33 @@ def _fetch_openmeteo(zone_id: str, valid_time: datetime) -> ZoneObs:
     data  = _cached_get(url, params)
     daily = data.get("daily", {})
 
-    # Chronological, ending at `anchor` -- index from the END, not the
-    # start, so "last N days" genuinely means the N most recent days
+    # 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)
@@ -534,7 +596,7 @@ def _fetch_openmeteo(zone_id: str, valid_time: datetime) -> ZoneObs:
     return ZoneObs(
         zone_id=zone_id,
         valid_time=anchor,
-        source=DataSource.OPENMETEO_LIVE,
+        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),
@@ -558,8 +620,12 @@ def _fetch_era5(zone_id: str, valid_time: datetime) -> ZoneObs:
     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:
-        logger.debug("cdsapi unavailable — falling back to synthetic for %s", zone_id)
-        return _fetch_synthetic(zone_id, valid_time)
+        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)
@@ -622,41 +688,47 @@ def _fetch_era5(zone_id: str, valid_time: datetime) -> ZoneObs:
             )
             logger.info("ERA5 download complete: %s", nc_path)
         except Exception as e:
-            logger.warning(
-                "ERA5 CDS request failed for %s: %s — falling back to synthetic",
-                zone_id, e
-            )
             if nc_path.exists():
                 nc_path.unlink()
-            return _fetch_synthetic(zone_id, anchor)
+            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:
-        logger.warning(
-            "ERA5 NetCDF parse failed for %s: %s — falling back to synthetic",
-            zone_id, e
-        )
-        return _fetch_synthetic(zone_id, anchor)
+        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]."""
+    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)
+        window_start = start - timedelta(days=n_days - 1)
         coll = (
             ee.ImageCollection("NASA/GPM_L3/IMERG_V07")
-            .filterDate(window_start.isoformat(), start.isoformat())
+            .filterDate(window_start.isoformat(), inclusive_end.isoformat())
             .filterBounds(region)
             .select("precipitation")
         )
@@ -1102,10 +1174,29 @@ def _select_source(cfg: ForecastConfig, rng: random.Random) -> DataSource:
     )
 
 
-def _fallback_chain(primary: DataSource) -> List[DataSource]:
+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):
-        return [primary, DataSource.ERA5_REANALYSIS, DataSource.SYNTHETIC]
-    return [primary, DataSource.SYNTHETIC]
+        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(
@@ -1133,6 +1224,15 @@ def fetch_zone_obs(
     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,
@@ -1140,9 +1240,12 @@ def fetch_zone_obs(
         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):
+    for source in _fallback_chain(primary, allow_synthetic=allow_synthetic):
         try:
             obs = fetchers[source](zone_id, valid_time)
             ZoneObs.validate(obs, strict=True)
@@ -1764,6 +1867,72 @@ def _self_test() -> int:
     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")
+
     print()
     if failures:
         print(f"FAILED  {len(failures)} test(s):")
@@ -1771,7 +1940,7 @@ def _self_test() -> int:
             print(f"  - {f}")
         sys.exit(1)
     else:
-        print("All 5 test groups passed.")
+        print("All 8 test groups passed.")
 
 
 if __name__ == "__main__":