Update climatology.py

climatology.py

@@ -3,25 +3,8 @@ climatology.py
 ==============
 Per-zone day-of-year climatology and anomaly (z-score) computation.
 
-WHY THIS MODULE EXISTS (the gap it closes)
-------------------------------------------
-Every real fetcher in era5_data_pipeline.py (_build_era5_obs,
-_fetch_openmeteo, _fetch_imerg, _fetch_smap) sets precip_anomaly_idx /
-temp_anomaly_idx / soil_moisture_anom to 0.0 with the comment "requires
-climatology -- set in scorer". Nothing ever supplied that climatology, and
-the scorer never set the fields either. The consequences are structural:
-
-  * ZoneObs.drought_signal() = 0.6 * clip(-precip_anomaly/3)
-                             + 0.4 * clip(-soil_anom/3)
-  * ZoneObs.flood_signal()   = 0.4 * clip(+precip_anomaly/3) + ...
-
-With all anomalies pinned at 0.0, drought_signal() is identically 0.0 and
-flood_signal() loses its primary term on every real observation. The whole
-risk-scoring stack (crop_risk_scorer, the env's belief initialisation, the
-hierarchical search gating) was effectively only sensitive on SYNTHETIC
-data, where anomalies are injected directly by the event flags. This module
-is the missing climatology layer: it turns absolute real-world readings
-into the z-score anomalies the rest of the pipeline was designed around.
+ 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)
 ----------------------------------------------------------------------------
@@ -39,14 +22,6 @@ DATA SOURCES (two tiers, matching the codebase's degrade-safely philosophy)
      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.
-
-IMPORTANT USAGE NOTE
---------------------
-apply_climatology_anomalies() SKIPS observations whose source is
-DataSource.SYNTHETIC. The synthetic generator injects its own meaningful
-anomalies via the drought/flood event flags; re-scoring those against a
-climatology would double-transform a deliberately-constructed signal.
-Climatology anomalies are for REAL observations only.
 """
 
 from __future__ import annotations
@@ -83,15 +58,12 @@ except ImportError:
 # Constants
 # ---------------------------------------------------------------------------
 
-# Same archive endpoint era5_data_pipeline.py uses; duplicated here (rather
-# than imported) so era5_data_pipeline can import THIS module lazily without
-# a circular import at module load time.
 _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
+_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)
@@ -100,7 +72,9 @@ _TABLE_LEN = 366  # DOY table indexed doy-1; DOY 60 = Feb 29 (leap mapping below
 _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; prevents z-score blowup when tropical RH variance is naturally small
+_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.
@@ -116,12 +90,6 @@ def _is_leap(year: int) -> bool:
 
 
 def doy_index(dt: datetime) -> int:
-    """Map a date to a 1..366 index in the fixed climatology table.
-
-    In non-leap years, dates after Feb 28 are shifted up by one so that e.g.
-    Mar 1 always maps to the same table entry (61) in every year. Table entry
-    60 (Feb 29) is only ever hit by leap-year dates.
-    """
     doy = dt.timetuple().tm_yday
     if not _is_leap(dt.year) and doy >= 60:
         doy += 1
@@ -129,7 +97,6 @@ def doy_index(dt: datetime) -> int:
 
 
 def _circular_smooth(values: List[float], half_window: int = 7) -> List[float]:
-    """Circular moving average over the DOY table (Dec wraps to Jan)."""
     n = len(values)
     out = []
     for i in range(n):
@@ -148,12 +115,6 @@ def _circular_smooth(values: List[float], half_window: int = 7) -> List[float]:
 
 @dataclass
 class ZoneClimatology:
-    """Day-of-year climatology for one zone.
-
-    All lists have length 366 and are indexed by (doy_index(dt) - 1).
-    precip is a DAILY mean rate (mm/day); window aggregates are computed by
-    summing daily means over the window (see window_precip_stats).
-    """
     zone_id: str
     source: str                       # 'openmeteo_archive' | 'synthetic_model' | 'mixed'
     n_years: int
@@ -165,9 +126,6 @@ class ZoneClimatology:
     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)
-    # FIX: baseline uses rh_mean (not rh_max -- unreliable in Open-Meteo);
-    # correlates well enough to correct fungi_risk_signal()'s unchanged
-    # max-based threshold.
     rh_mean_pct:    List[float] = field(default_factory=list)
     rh_std_pct:     List[float] = field(default_factory=list)
 
@@ -184,15 +142,6 @@ class ZoneClimatology:
 
     # --- Window statistics ------------------------------------------------
     def window_precip_stats(self, dt: datetime, window_days: int) -> Tuple[float, float]:
-        """Climatological mean and std of a TRAILING `window_days` precip total
-        ending at dt's day-of-year.
-
-        Mean: sum of daily means (exact under the daily model).
-        Std: sqrt(sum of daily variances) -- assumes day-to-day independence,
-             so it UNDERSTATES true variance during correlated multi-day
-             spells, inflating anomaly magnitude for persistent events.
-             Bounded by the floors and the [-5, 5] clip in ZoneObs.
-        """
         idx0 = doy_index(dt) - 1
         mean = 0.0
         var = 0.0
@@ -234,10 +183,6 @@ class ZoneClimatology:
 
     @classmethod
     def from_dict(cls, d: Dict[str, Any]) -> "ZoneClimatology":
-        # FIX: rh_mean_pct/rh_std_pct are new fields; old cached files
-        # predate them. Default to flat 85% (synthetic-model range) instead
-        # of crashing -- self-heals within _CACHE_TTL_DAYS as real RH is
-        # fetched.
         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:
@@ -266,22 +211,6 @@ class ZoneClimatology:
 # ---------------------------------------------------------------------------
 
 def _synthetic_climatology(zone_id: str, lat: float, n_years: int = 0) -> ZoneClimatology:
-    """Deterministic heuristic climatology for the Indonesian maritime continent.
-    A documented heuristic, NOT a retrieval -- see module docstring.
-
-      * Precip: single-harmonic wet season, peak ~DOY 30 (late Jan) for the
-        southern archipelago (Java, Bali, Nusa Tenggara, Sulawesi, S. Sumatra);
-        peak shifts earlier (Oct-Dec) moving north past ~1 deg N (N. Sumatra).
-        Amplitude grows with distance from equator; equatorial belt stays wet
-        year-round.
-      * Temp: weak annual cycle (~2.6 deg C peak-to-peak), coolest Jul-Aug in
-        the south (SH dry season), weaker and phase-reversed north of equator.
-      * Soil: precip-tracked with ~20-day lag, scaled to ERA5 swvl1's typical
-        volumetric-% range for the region.
-
-    Per-zone jitter (+/-10% on base/amp, via _stable_seed) keeps neighbouring
-    zones numerically distinct without changing the seasonal shape.
-    """
     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)
@@ -316,10 +245,6 @@ def _synthetic_climatology(zone_id: str, lat: float, n_years: int = 0) -> ZoneCl
         temp_mean.append(t)
         temp_std.append(0.7)
 
-        # RH baseline tracks the wet season (in phase with precip), range
-        # 78-94%. This is a baseline for ANOMALY detection only --
-        # fungi_risk_signal() still applies its own absolute threshold to
-        # rh_max_pct separately.
         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))
@@ -359,29 +284,6 @@ def _fetch_openmeteo_climatology(
     years: int,
     end_year: Optional[int] = None,
 ) -> ZoneClimatology:
-    """Build a DOY climatology from the Open-Meteo historical archive.
-
-    Downloads `years` full calendar years of daily data in one request and
-    pools by day-of-year. Raises on any failure -- the caller
-    (get_zone_climatology) falls back to the synthetic model, matching the
-    pipeline-wide degrade-safely pattern.
-
-    soil_moisture_0_to_7cm_mean is requested per the Open-Meteo archive
-    documentation at write time (VERIFIED LIVE against the archive API:
-    the variable exists and returns daily means). UNITS: Open-Meteo returns
-    soil moisture in m3/m3; this function converts to percent (x100) so the
-    table matches ZoneObs.soil_moisture_pct and ERA5's swvl1 x 100 handling
-    in era5_data_pipeline._build_era5_obs. (Found via a z-score clipped at
-    +5.0 against an 18% observation -- the raw 0.2-0.4 m3/m3 values were
-    being read as ~0.3%.)
-
-    If the key is absent/empty in the response (API change, or variable not
-    in the daily list for this endpoint), the soil tables are derived from
-    the REAL precip series via the same lagged mapping the synthetic model
-    uses -- so a soil-variable outage degrades one field's provenance, not
-    the whole fetch. The result's `source` is then 'mixed' rather than
-    'openmeteo_archive' so downstream auditing can tell.
-    """
     if not _REQUESTS_AVAILABLE:
         raise RuntimeError("requests not installed")
 
@@ -567,20 +469,6 @@ def get_zone_climatology(
     prefer_real: bool = True,
     use_cache: bool = True,
 ) -> ZoneClimatology:
-    """Return the day-of-year climatology for a zone, cached on disk.
-
-    Resolution order:
-      1. Fresh cache hit (same zone/years/end_year, < _CACHE_TTL_DAYS old).
-      2. Real Open-Meteo archive fetch (if prefer_real and requests present).
-      3. Deterministic synthetic monsoon model (never fails).
-
-    Args:
-        end_year:   Last calendar year included in the climatology period.
-                    Default: the most recent COMPLETE year (now.year - 1).
-                    Pin this explicitly for backtests so the climatology
-                    cannot see the period being backtested (look-ahead).
-        prefer_real: Set False to force the synthetic model (offline tests).
-    """
     path = _cache_path(zone_id, years, end_year)
     if use_cache and path.exists():
         age_days = (
@@ -623,29 +511,6 @@ def get_zone_climatology(
 
 
 def apply_climatology_anomalies(obs: ZoneObs, clim: ZoneClimatology) -> ZoneObs:
-    """Return a NEW ZoneObs with the four anomaly fields populated as
-    z-scores against `clim`. Never mutates the input (to_dict/from_dict
-    round-trip, matching the codebase idiom).
-
-    Skipped/Guarded cases (all deliberate, all logged at debug level):
-      * obs.source == SYNTHETIC: returned unchanged. Synthetic obs carry
-        injected anomalies from the event flags; re-scoring them against a
-        climatology would double-transform the training signal.
-      * precip anomaly only computed when at least one precip aggregate is
-        non-zero (a precip-less fetch like _fetch_smap would otherwise read
-        as a catastrophic false drought: (0 - mean)/std << 0).
-      * temp anomaly only when temp_mean_c != 0.0 (0.0 is the "unset"
-        default, not a real temperature in this pipeline's operating range).
-      * soil anomaly only when soil_moisture_pct > 0.0.
-      * rh anomaly only when rh_mean_pct > 0.0 (0.0 is "unset", not a real
-        humidity reading -- see rh_anomaly_idx's field comment in
-        zone_observation.py for why this exists: fungi_risk_signal()'s pure
-        absolute-RH threshold was flat across ENSO regimes in a tropical
-        climate, so it was masking correctly regime-sensitive drought/flood
-        signals in the actual alert_level output).
-
-    Z-scores are clipped to [-5, 5] by ZoneObs.__post_init__ as usual.
-    """
     if obs.source == DataSource.SYNTHETIC:
         logger.debug(
             "climatology: %s source is SYNTHETIC -- anomalies left as injected",
@@ -690,12 +555,6 @@ def apply_anomalies_by_zone_id(
     years: int = 10,
     prefer_real: bool = True,
 ) -> ZoneObs:
-    """Convenience wrapper: resolve (or build) the cached climatology for
-    obs.zone_id, then apply it. This is the entry point era5_data_pipeline
-    calls; kept separate from apply_climatology_anomalies so callers that
-    already hold a ZoneClimatology (e.g. the backtester looping over days)
-    don't pay the cache lookup per step.
-    """
     clim = get_zone_climatology(obs.zone_id, lat, lon, years=years,
                                 prefer_real=prefer_real)
     return apply_climatology_anomalies(obs, clim)
@@ -831,4 +690,4 @@ if __name__ == "__main__":
             print(f"  - {f}")
         sys.exit(1)
     else:
-        print("All 9 test groups passed.")
+        print("All 9 test groups passed.")