""" climatology.py ============== Per-zone day-of-year climatology and anomaly (z-score) computation. This module turns absolute real-world readings into the z-score anomalies the rest of the pipeline was designed around. DATA SOURCES (two tiers, matching the codebase's degrade-safely philosophy) ---------------------------------------------------------------------------- 1. Real: Open-Meteo historical archive API (free, no API key -- the same endpoint era5_data_pipeline._fetch_openmeteo already uses). Daily precipitation_sum + temperature_2m_mean are well-established archive variables. soil_moisture_0_to_7cm_mean is requested as documented in the Open-Meteo archive docs at write time; if the API rejects it or returns nothing, soil climatology degrades to the precip-tracked model below, logged -- never silently zeroed. 2. Synthetic: a deterministic, latitude-aware maritime-continent monsoon model (see _synthetic_climatology). It is a HEURISTIC, calibrated to the broad shape of the Indonesian wet/dry season (SH monsoon: wet Dec-Mar, dry Jun-Sep; weaker/bimodal near the equator; shifted peak for northern Sumatra). It exists so the pipeline keeps producing sensible anomalies offline and in tests. It is NOT a retrieval -- treat its absolute values as plausible shapes, not measurements. """ from __future__ import annotations import json import logging import math import os from dataclasses import dataclass, field from datetime import datetime, timedelta, timezone from pathlib import Path from typing import Any, Dict, List, Optional, Tuple import zone_observation as _zo assert _zo.SCHEMA_VERSION == 3, ( f"climatology: zone_observation schema mismatch " f"(expected 3, got {_zo.SCHEMA_VERSION})" ) from zone_observation import DataSource, ZoneObs, _clip, _stable_seed logger = logging.getLogger(__name__) try: import requests _REQUESTS_AVAILABLE = True except ImportError: _REQUESTS_AVAILABLE = False logger.info("requests not installed -- climatology will use the synthetic model") # --------------------------------------------------------------------------- # Constants # --------------------------------------------------------------------------- _OPENMETEO_ARCHIVE_URL = "https://archive-api.open-meteo.com/v1/archive" _TIMEOUT_S = int(os.environ.get("WEATHER_HTTP_TIMEOUT", "60")) _CACHE_DIR = Path(os.environ.get("WEATHER_CACHE_DIR", ".cache/era5")) / "climatology" _CACHE_DIR.mkdir(parents=True, exist_ok=True) _CACHE_TTL_DAYS = 90 # climatology drifts slowly; quarterly refresh is ample _DAYS_PER_YEAR = 365.25 _TABLE_LEN = 366 # DOY table indexed doy-1; DOY 60 = Feb 29 (leap mapping below) # Std floors: prevent division blow-ups in convectively uniform seasons. _PRECIP_STD_FLOOR = 1.5 # mm/day _TEMP_STD_FLOOR = 0.4 # deg C _SOIL_STD_FLOOR = 2.0 # percent _RH_STD_FLOOR = 3.0 # percent -- RH is bounded [0,100] and often near- # saturated in the tropics, so day-to-day variance # is naturally small; floor prevents z-score blowup # Trailing window for the precipitation anomaly. Matches ZoneObs.precip_30d_mm, # the longest aggregate every real fetcher populates. _PRECIP_WINDOW_DAYS = 30 # --------------------------------------------------------------------------- # Day-of-year helpers (fixed 366-entry table regardless of leap years) # --------------------------------------------------------------------------- def _is_leap(year: int) -> bool: return year % 4 == 0 and (year % 100 != 0 or year % 400 == 0) def doy_index(dt: datetime) -> int: doy = dt.timetuple().tm_yday if not _is_leap(dt.year) and doy >= 60: doy += 1 return min(doy, _TABLE_LEN) def _circular_smooth(values: List[float], half_window: int = 7) -> List[float]: n = len(values) out = [] for i in range(n): acc = 0.0 cnt = 0 for j in range(-half_window, half_window + 1): acc += values[(i + j) % n] cnt += 1 out.append(acc / cnt) return out # --------------------------------------------------------------------------- # ZoneClimatology # --------------------------------------------------------------------------- @dataclass class ZoneClimatology: zone_id: str source: str # 'openmeteo_archive' | 'synthetic_model' | 'mixed' n_years: int period_start_year: int period_end_year: int precip_mean_mm: List[float] = field(default_factory=list) # mm/day precip_std_mm: List[float] = field(default_factory=list) temp_mean_c: List[float] = field(default_factory=list) temp_std_c: List[float] = field(default_factory=list) soil_mean_pct: List[float] = field(default_factory=list) soil_std_pct: List[float] = field(default_factory=list) rh_mean_pct: List[float] = field(default_factory=list) rh_std_pct: List[float] = field(default_factory=list) def __post_init__(self) -> None: for name in ("precip_mean_mm", "precip_std_mm", "temp_mean_c", "temp_std_c", "soil_mean_pct", "soil_std_pct", "rh_mean_pct", "rh_std_pct"): v = getattr(self, name) if len(v) != _TABLE_LEN: raise ValueError( f"ZoneClimatology('{self.zone_id}'): {name} has length " f"{len(v)}, expected {_TABLE_LEN}" ) # --- Window statistics ------------------------------------------------ def window_precip_stats(self, dt: datetime, window_days: int) -> Tuple[float, float]: idx0 = doy_index(dt) - 1 mean = 0.0 var = 0.0 for k in range(window_days): i = (idx0 - k) % _TABLE_LEN mean += self.precip_mean_mm[i] var += self.precip_std_mm[i] ** 2 return mean, max(math.sqrt(var), _PRECIP_STD_FLOOR) def daily_temp_stats(self, dt: datetime) -> Tuple[float, float]: i = doy_index(dt) - 1 return self.temp_mean_c[i], max(self.temp_std_c[i], _TEMP_STD_FLOOR) def daily_soil_stats(self, dt: datetime) -> Tuple[float, float]: i = doy_index(dt) - 1 return self.soil_mean_pct[i], max(self.soil_std_pct[i], _SOIL_STD_FLOOR) def daily_rh_stats(self, dt: datetime) -> Tuple[float, float]: i = doy_index(dt) - 1 return self.rh_mean_pct[i], max(self.rh_std_pct[i], _RH_STD_FLOOR) # --- Serialisation (JSON cache) --------------------------------------- def to_dict(self) -> Dict[str, Any]: return { "zone_id": self.zone_id, "source": self.source, "n_years": self.n_years, "period_start_year": self.period_start_year, "period_end_year": self.period_end_year, "precip_mean_mm": self.precip_mean_mm, "precip_std_mm": self.precip_std_mm, "temp_mean_c": self.temp_mean_c, "temp_std_c": self.temp_std_c, "soil_mean_pct": self.soil_mean_pct, "soil_std_pct": self.soil_std_pct, "rh_mean_pct": self.rh_mean_pct, "rh_std_pct": self.rh_std_pct, } @classmethod def from_dict(cls, d: Dict[str, Any]) -> "ZoneClimatology": rh_mean = d.get("rh_mean_pct") rh_std = d.get("rh_std_pct") if rh_mean is None or len(rh_mean) != _TABLE_LEN: rh_mean = [85.0] * _TABLE_LEN if rh_std is None or len(rh_std) != _TABLE_LEN: rh_std = [_RH_STD_FLOOR] * _TABLE_LEN return cls( zone_id=d["zone_id"], source=d.get("source", "unknown"), n_years=int(d.get("n_years", 0)), period_start_year=int(d.get("period_start_year", 0)), period_end_year=int(d.get("period_end_year", 0)), precip_mean_mm=[float(v) for v in d["precip_mean_mm"]], precip_std_mm=[float(v) for v in d["precip_std_mm"]], temp_mean_c=[float(v) for v in d["temp_mean_c"]], temp_std_c=[float(v) for v in d["temp_std_c"]], soil_mean_pct=[float(v) for v in d["soil_mean_pct"]], soil_std_pct=[float(v) for v in d["soil_std_pct"]], rh_mean_pct=[float(v) for v in rh_mean], rh_std_pct=[float(v) for v in rh_std], ) # --------------------------------------------------------------------------- # Tier 2 -- deterministic synthetic maritime-continent climatology # --------------------------------------------------------------------------- def _synthetic_climatology(zone_id: str, lat: float, n_years: int = 0) -> ZoneClimatology: seed = _stable_seed(f"clim_{zone_id}") # Deterministic jitter in [0.9, 1.1] from the seed's low bits. jitter = 0.9 + 0.2 * ((seed % 1000) / 1000.0) abs_lat = abs(lat) # Wet-season peak: late Jan in the south, shifting earlier north of ~1N. peak_doy = 30.0 if lat <= 1.0 else max(300.0, 30.0 - 12.0 * lat) amp_scale = _clip(abs_lat / 8.0, 0.35, 1.0) precip_base = max(3.0, (7.0 - 0.25 * abs_lat) * jitter) # mm/day precip_amp = 4.5 * amp_scale * jitter # mm/day temp_base = 27.0 - 0.30 * abs_lat # SH zones: coolest around DOY ~200 (mid-Jul). NH: weaker, reversed. temp_amp = 1.3 if lat < 0.0 else -0.5 precip_mean, precip_std = [], [] temp_mean, temp_std = [], [] soil_mean, soil_std = [], [] rh_mean, rh_std = [], [] daily_precip_for_soil: List[float] = [] for doy in range(1, _TABLE_LEN + 1): phase = 2.0 * math.pi * (doy - peak_doy) / _DAYS_PER_YEAR p = precip_base + precip_amp * math.cos(phase) p = max(0.8, p) daily_precip_for_soil.append(p) precip_mean.append(p) precip_std.append(max(_PRECIP_STD_FLOOR, 0.9 * p)) t_phase = 2.0 * math.pi * (doy - 200.0) / _DAYS_PER_YEAR t = temp_base - temp_amp * math.cos(t_phase) temp_mean.append(t) temp_std.append(0.7) r = 86.0 + amp_scale * 6.0 * math.cos(phase) rh_mean.append(_clip(r, 78.0, 94.0)) rh_std.append(max(_RH_STD_FLOOR, 3.5)) # Soil tracks precip with a 20-day lag. for doy in range(1, _TABLE_LEN + 1): lagged = daily_precip_for_soil[(doy - 1 - 20) % _TABLE_LEN] s = _clip(16.0 + 2.4 * lagged, 8.0, 52.0) soil_mean.append(s) soil_std.append(max(_SOIL_STD_FLOOR, 4.0)) return ZoneClimatology( zone_id=zone_id, source="synthetic_model", n_years=n_years, period_start_year=0, period_end_year=0, precip_mean_mm=_circular_smooth(precip_mean), precip_std_mm=precip_std, temp_mean_c=_circular_smooth(temp_mean), temp_std_c=temp_std, soil_mean_pct=_circular_smooth(soil_mean), soil_std_pct=soil_std, rh_mean_pct=_circular_smooth(rh_mean), rh_std_pct=rh_std, ) # --------------------------------------------------------------------------- # Tier 1 -- real climatology from the Open-Meteo archive # --------------------------------------------------------------------------- def _fetch_openmeteo_climatology( zone_id: str, lat: float, lon: float, years: int, end_year: Optional[int] = None, ) -> ZoneClimatology: if not _REQUESTS_AVAILABLE: raise RuntimeError("requests not installed") last_full_year = (end_year if end_year is not None else datetime.now(timezone.utc).year - 1) start_year = last_full_year - years + 1 params = { "latitude": lat, "longitude": lon, "start_date": f"{start_year}-01-01", "end_date": f"{last_full_year}-12-31", "daily": ",".join([ "precipitation_sum", "temperature_2m_mean", "soil_moisture_0_to_7cm_mean", "relative_humidity_2m_mean", ]), "timezone": "UTC", } resp = requests.get(_OPENMETEO_ARCHIVE_URL, params=params, timeout=_TIMEOUT_S) resp.raise_for_status() data = resp.json() daily = data.get("daily", {}) dates = daily.get("time", []) if not dates: raise RuntimeError(f"Open-Meteo archive returned no daily rows for {zone_id}") precip_series = daily.get("precipitation_sum", []) temp_series = daily.get("temperature_2m_mean", []) soil_series = daily.get("soil_moisture_0_to_7cm_mean", []) rh_series = daily.get("relative_humidity_2m_mean", []) p_sum = [0.0] * _TABLE_LEN p_sq = [0.0] * _TABLE_LEN p_n = [0] * _TABLE_LEN t_sum = [0.0] * _TABLE_LEN t_sq = [0.0] * _TABLE_LEN t_n = [0] * _TABLE_LEN s_sum = [0.0] * _TABLE_LEN s_sq = [0.0] * _TABLE_LEN s_n = [0] * _TABLE_LEN r_sum = [0.0] * _TABLE_LEN r_sq = [0.0] * _TABLE_LEN r_n = [0] * _TABLE_LEN def _val(series: List[Any], i: int) -> Optional[float]: if i >= len(series): return None v = series[i] if v is None: return None try: return float(v) except (TypeError, ValueError): return None for i, date_str in enumerate(dates): try: dt = datetime.fromisoformat(date_str).replace(tzinfo=timezone.utc) except ValueError: continue k = doy_index(dt) - 1 p = _val(precip_series, i) if p is not None: p_sum[k] += p p_sq[k] += p * p p_n[k] += 1 t = _val(temp_series, i) if t is not None: t_sum[k] += t t_sq[k] += t * t t_n[k] += 1 s = _val(soil_series, i) if s is not None: s_pct = s * 100.0 # m3/m3 -> % (matches ZoneObs.soil_moisture_pct) s_sum[k] += s_pct s_sq[k] += s_pct * s_pct s_n[k] += 1 r = _val(rh_series, i) if r is not None: r_sum[k] += r r_sq[k] += r * r r_n[k] += 1 if sum(p_n) < 300 * years or sum(t_n) < 300 * years: raise RuntimeError( f"Open-Meteo archive coverage too thin for {zone_id}: " f"precip_days={sum(p_n)} temp_days={sum(t_n)} over {years}y" ) def _mean_std(sums, sqs, ns, floor): means, stds = [], [] for k in range(_TABLE_LEN): n = ns[k] if n == 0: # Should not happen with full-year coverage; guard anyway. means.append(0.0) stds.append(floor) continue m = sums[k] / n var = max(0.0, sqs[k] / n - m * m) means.append(m) stds.append(max(floor, math.sqrt(var))) return means, stds precip_mean, precip_std = _mean_std(p_sum, p_sq, p_n, _PRECIP_STD_FLOOR) temp_mean, temp_std = _mean_std(t_sum, t_sq, t_n, _TEMP_STD_FLOOR) soil_days = sum(s_n) soil_thin = soil_days < 300 * years if not soil_thin: soil_mean, soil_std = _mean_std(s_sum, s_sq, s_n, _SOIL_STD_FLOOR) else: logger.warning( "climatology: soil_moisture_0_to_7cm_mean coverage thin for %s " "(%d days over %dy) -- deriving soil tables from the real precip " "series (lagged mapping). Provenance marked 'mixed'.", zone_id, soil_days, years, ) soil_mean, soil_std = [], [] for doy in range(1, _TABLE_LEN + 1): lagged = precip_mean[(doy - 1 - 20) % _TABLE_LEN] soil_mean.append(_clip(16.0 + 2.4 * lagged, 8.0, 52.0)) soil_std.append(max(_SOIL_STD_FLOOR, 4.0)) rh_days = sum(r_n) rh_thin = rh_days < 300 * years if not rh_thin: rh_mean, rh_std = _mean_std(r_sum, r_sq, r_n, _RH_STD_FLOOR) else: logger.warning( "climatology: relative_humidity_2m_mean coverage thin for %s " "(%d days over %dy) -- deriving RH tables from the real precip " "series (wet-season correlation, same phase). Provenance marked " "'mixed'.", zone_id, rh_days, years, ) rh_mean, rh_std = [], [] p_min, p_max = min(precip_mean), max(precip_mean) p_span = max(p_max - p_min, 1e-6) for doy in range(1, _TABLE_LEN + 1): p_frac = (precip_mean[doy - 1] - p_min) / p_span # 0..1 rh_mean.append(_clip(78.0 + 12.0 * p_frac, 78.0, 94.0)) rh_std.append(max(_RH_STD_FLOOR, 3.5)) source = "openmeteo_archive" if not (soil_thin or rh_thin) else "mixed" return ZoneClimatology( zone_id=zone_id, source=source, n_years=years, period_start_year=start_year, period_end_year=last_full_year, precip_mean_mm=_circular_smooth(precip_mean), precip_std_mm=precip_std, temp_mean_c=_circular_smooth(temp_mean), temp_std_c=temp_std, soil_mean_pct=_circular_smooth(soil_mean), soil_std_pct=soil_std, rh_mean_pct=_circular_smooth(rh_mean), rh_std_pct=rh_std, ) # --------------------------------------------------------------------------- # Public API: cached climatology + anomaly application # --------------------------------------------------------------------------- def _cache_path(zone_id: str, years: int, end_year: Optional[int]) -> Path: key = _stable_seed(f"{zone_id}|{years}|{end_year}") return _CACHE_DIR / f"{zone_id}_{key}.json" def get_zone_climatology( zone_id: str, lat: float, lon: float, years: int = 10, end_year: Optional[int] = None, prefer_real: bool = True, use_cache: bool = True, ) -> ZoneClimatology: path = _cache_path(zone_id, years, end_year) if use_cache and path.exists(): age_days = ( datetime.now(timezone.utc) - datetime.fromtimestamp(path.stat().st_mtime, tz=timezone.utc) ).days if age_days < _CACHE_TTL_DAYS: try: with open(path) as f: return ZoneClimatology.from_dict(json.load(f)) except Exception as e: logger.warning("climatology: cache read failed (%s) -- rebuilding", e) clim: Optional[ZoneClimatology] = None if prefer_real and _REQUESTS_AVAILABLE: try: clim = _fetch_openmeteo_climatology(zone_id, lat, lon, years, end_year) logger.info( "climatology: built real %dy climatology for %s (%d-%d)", years, zone_id, clim.period_start_year, clim.period_end_year, ) except Exception as e: logger.warning( "climatology: real fetch failed for %s (%s) -- synthetic model", zone_id, e, ) clim = None if clim is None: clim = _synthetic_climatology(zone_id, lat, n_years=years) if use_cache: try: with open(path, "w") as f: json.dump(clim.to_dict(), f) except Exception as e: logger.warning("climatology: cache write failed (%s) -- continuing", e) return clim def apply_climatology_anomalies(obs: ZoneObs, clim: ZoneClimatology) -> ZoneObs: if obs.source == DataSource.SYNTHETIC: logger.debug( "climatology: %s source is SYNTHETIC -- anomalies left as injected", obs.zone_id, ) return obs d = obs.to_dict() d.pop("_schema_version", None) has_precip = (obs.precip_30d_mm > 0.0) or (obs.precip_14d_mm > 0.0) \ or (obs.precip_7d_mm > 0.0) or (obs.precip_24h_mm > 0.0) if has_precip: mean_w, std_w = clim.window_precip_stats(obs.valid_time, _PRECIP_WINDOW_DAYS) d["precip_anomaly_idx"] = _clip( (obs.precip_30d_mm - mean_w) / std_w, -5.0, 5.0 ) if obs.temp_mean_c != 0.0: t_mean, t_std = clim.daily_temp_stats(obs.valid_time) d["temp_anomaly_idx"] = _clip((obs.temp_mean_c - t_mean) / t_std, -5.0, 5.0) if obs.soil_moisture_pct > 0.0: s_mean, s_std = clim.daily_soil_stats(obs.valid_time) d["soil_moisture_anom"] = _clip( (obs.soil_moisture_pct - s_mean) / s_std, -5.0, 5.0 ) if obs.rh_mean_pct > 0.0: r_mean, r_std = clim.daily_rh_stats(obs.valid_time) d["rh_anomaly_idx"] = _clip( (obs.rh_mean_pct - r_mean) / r_std, -5.0, 5.0 ) return ZoneObs.from_dict(d) def apply_anomalies_by_zone_id( obs: ZoneObs, lat: float, lon: float, years: int = 10, prefer_real: bool = True, ) -> ZoneObs: clim = get_zone_climatology(obs.zone_id, lat, lon, years=years, prefer_real=prefer_real) return apply_climatology_anomalies(obs, clim) # --------------------------------------------------------------------------- # Self-test (python climatology.py) -- fully offline # --------------------------------------------------------------------------- if __name__ == "__main__": import sys from datetime import timezone as _tz from zone_observation import make_synthetic_zone_obs logging.basicConfig(level=logging.WARNING) print("climatology.py self-test (offline: prefer_real=False)\n") failures: List[str] = [] def _assert(cond: bool, msg: str) -> None: if not cond: failures.append(msg) print(f" FAIL: {msg}") LAT, LON = -6.3, 107.3 # Karawang, West Java # 1. Synthetic climatology: shape, length, round-trip clim = _synthetic_climatology("test_zone", LAT) _assert(len(clim.precip_mean_mm) == _TABLE_LEN, "precip table length") d = clim.to_dict() clim2 = ZoneClimatology.from_dict(d) _assert(clim2.zone_id == clim.zone_id, "ZoneClimatology round-trip zone_id") _assert(abs(clim2.precip_mean_mm[100] - clim.precip_mean_mm[100]) < 1e-12, "ZoneClimatology round-trip values") # 2. Seasonality: Java should be much wetter in Jan than in Aug jan_mean = sum(clim.precip_mean_mm[0:31]) / 31.0 aug_mean = sum(clim.precip_mean_mm[212:243]) / 31.0 _assert(jan_mean > aug_mean * 1.3, f"monsoon shape wrong: Jan={jan_mean:.1f} vs Aug={aug_mean:.1f} mm/day") print(f" Seasonality OK: Jan {jan_mean:.1f} mm/day vs Aug {aug_mean:.1f} mm/day") # 3. doy_index leap mapping: Mar 1 maps to the same entry in every year d1 = doy_index(datetime(2023, 3, 1, tzinfo=_tz.utc)) d2 = doy_index(datetime(2024, 3, 1, tzinfo=_tz.utc)) _assert(d1 == d2 == 61, f"Mar 1 mapping inconsistent: {d1} vs {d2}") _assert(doy_index(datetime(2024, 2, 29, tzinfo=_tz.utc)) == 60, "Feb 29 mapping") print(f" doy_index OK (Mar 1 -> {d1}, Feb 29 -> 60)") # 4. Window stats: 30-day wet-season aggregate exceeds dry-season wet_dt = datetime(2024, 1, 31, tzinfo=_tz.utc) dry_dt = datetime(2024, 8, 31, tzinfo=_tz.utc) wet_mean, wet_std = clim.window_precip_stats(wet_dt, 30) dry_mean, _ = clim.window_precip_stats(dry_dt, 30) _assert(wet_mean > dry_mean, "window aggregate seasonality wrong") _assert(wet_std >= _PRECIP_STD_FLOOR, "window std floor violated") print(f" Window stats OK: wet30={wet_mean:.0f}mm dry30={dry_mean:.0f}mm") # 5. apply_climatology_anomalies: real-source obs gets anomalies obs = make_synthetic_zone_obs("realish_zone", seed=1) od = obs.to_dict() od.pop("_schema_version", None) od["source"] = DataSource.OPENMETEO_LIVE.value # pretend real # A real fetcher leaves anomaly fields at 0.0 ("unset") -- mirror that so # this test measures exactly what this module adds. od["precip_anomaly_idx"] = od["temp_anomaly_idx"] = od["soil_moisture_anom"] = 0.0 real_obs = ZoneObs.from_dict(od) pre_precip_z = real_obs.precip_anomaly_idx out = apply_climatology_anomalies(real_obs, clim) _assert(out is not real_obs, "apply should return a NEW object") _assert(real_obs.precip_anomaly_idx == pre_precip_z, "input obs was mutated!") # The obs built by make_synthetic_zone_obs has neutral-ish aggregates; # anomaly must be finite and within clip range. _assert(-5.0 <= out.precip_anomaly_idx <= 5.0, "anomaly outside clip") _assert(-5.0 <= out.temp_anomaly_idx <= 5.0, "temp anomaly outside clip") _assert(-5.0 <= out.soil_moisture_anom <= 5.0, "soil anomaly outside clip") print(f" Anomaly application OK: precip_z={out.precip_anomaly_idx:+.2f} " f"temp_z={out.temp_anomaly_idx:+.2f} soil_z={out.soil_moisture_anom:+.2f}") # 6. Drought/wet extremes produce correctly-signed anomalies # (keep aggregates monotonic: 24h <= 7d <= 14d <= 30d) dry_obs_d = dict(od) dry_obs_d["precip_24h_mm"] = 0.0 dry_obs_d["precip_7d_mm"] = 0.01 * wet_mean / 4.0 dry_obs_d["precip_14d_mm"] = 0.02 * wet_mean / 2.0 dry_obs_d["precip_30d_mm"] = 0.05 * wet_mean # 5% of wet climatology dry_obs_d["valid_time"] = wet_dt.isoformat() dry_out = apply_climatology_anomalies(ZoneObs.from_dict(dry_obs_d), clim) _assert(dry_out.precip_anomaly_idx < -1.0, f"dry obs should get negative anomaly, got {dry_out.precip_anomaly_idx}") wet_obs_d = dict(od) wet_obs_d["precip_24h_mm"] = 2.5 * dry_mean / 30.0 wet_obs_d["precip_7d_mm"] = 2.5 * dry_mean / 4.0 wet_obs_d["precip_14d_mm"] = 2.5 * dry_mean / 2.0 wet_obs_d["precip_30d_mm"] = 2.5 * dry_mean wet_obs_d["valid_time"] = dry_dt.isoformat() wet_out = apply_climatology_anomalies(ZoneObs.from_dict(wet_obs_d), clim) _assert(wet_out.precip_anomaly_idx > 1.0, f"wet obs should get positive anomaly, got {wet_out.precip_anomaly_idx}") print(f" Sign check OK: dry_z={dry_out.precip_anomaly_idx:+.2f} " f"wet_z={wet_out.precip_anomaly_idx:+.2f}") # 7. SYNTHETIC-source obs is skipped unchanged syn = make_synthetic_zone_obs("syn_zone", drought=True, seed=2) syn_out = apply_climatology_anomalies(syn, clim) _assert(syn_out.precip_anomaly_idx == syn.precip_anomaly_idx, "SYNTHETIC obs anomaly was modified (should be skipped)") print(" SYNTHETIC skip OK") # 8. Zero-precip obs (smap-style real fetch: anomaly fields unset at 0.0) # must NOT read as a catastrophic false drought. zero_d = dict(od) zero_d["precip_24h_mm"] = zero_d["precip_7d_mm"] = 0.0 zero_d["precip_14d_mm"] = zero_d["precip_30d_mm"] = 0.0 zero_d["precip_anomaly_idx"] = 0.0 # real fetchers leave this unset zero_out = apply_climatology_anomalies(ZoneObs.from_dict(zero_d), clim) _assert(zero_out.precip_anomaly_idx == 0.0, "precip-less fetch should keep anomaly 0.0 (no false drought)") print(" Zero-precip guard OK") # 9. Cache round-trip through get_zone_climatology (offline path) clim_c = get_zone_climatology("cache_zone", LAT, LON, years=5, prefer_real=False, use_cache=True) clim_c2 = get_zone_climatology("cache_zone", LAT, LON, years=5, prefer_real=False, use_cache=True) _assert(clim_c.precip_mean_mm == clim_c2.precip_mean_mm, "cached climatology not identical") print(" Cache round-trip OK") print() if failures: print(f"FAILED {len(failures)} test(s):") for f in failures: print(f" - {f}") sys.exit(1) else: print("All 9 test groups passed.")