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