""" indonesia_zones.py ================== Indonesia grounding layer: real agricultural zone registry, rice crop calendars, monsoon-onset context, and planting-window planning outputs. """ from __future__ import annotations import logging import math from dataclasses import dataclass, field, asdict from datetime import datetime, timedelta, timezone from typing import Any, Dict, List, Optional, Tuple import zone_observation as _zo assert _zo.SCHEMA_VERSION == 3, ( f"indonesia_zones: zone_observation schema mismatch " f"(expected 3, got {_zo.SCHEMA_VERSION})" ) from zone_observation import ( AlertLevel, BasinContext, CropStage, ForecastConfig, ForecastResult, GeoPolygon, RiskScore, ZoneObs, _clip, ) logger = logging.getLogger(__name__) # --------------------------------------------------------------------------- # Zone registry # --------------------------------------------------------------------------- @dataclass(frozen=True) class IndonesiaZone: zone_id: str label: str province: str lat: float lon: float primary_crop: str # 'rice_lowland' | 'rice_upland' | 'maize' | 'mixed_horticulture' irrigation: str # 'irrigated' | 'rainfed' | 'supplemental' calendar: str # key into CROP_CALENDARS region_group: str # key into _ONSET_CLIMATOLOGY half_extent_deg: float = 0.35 notes: str = "" def to_polygon(self) -> GeoPolygon: h = self.half_extent_deg return GeoPolygon( zone_id=self.zone_id, label=self.label, vertices=[ (self.lat - h, self.lon - h), (self.lat - h, self.lon + h), (self.lat + h, self.lon + h), (self.lat + h, self.lon - h), ], ) # Major food-crop producing regions. Centroids are well-known geographic # centres of the named production belts. INDONESIA_ZONES: Tuple[IndonesiaZone, ...] = ( IndonesiaZone( "karawang_rice", "Karawang rice belt", "West Java", -6.30, 107.30, "rice_lowland", "irrigated", "java_double", "java", notes="Jatiluhur-irrigated plain; one of the highest-yield lowland rice areas in Indonesia.", ), IndonesiaZone( "indramayu_rice", "Indramayu rice belt", "West Java", -6.45, 108.10, "rice_lowland", "irrigated", "java_double", "java", notes="North-coast Java plain; double cropping dominant.", ), IndonesiaZone( "central_java_rice", "Central Java lowlands (Semarang-Demak)", "Central Java", -6.95, 110.40, "rice_lowland", "irrigated", "java_double", "java", ), IndonesiaZone( "east_java_rice", "East Java lowlands (Ngawi-Bojonegoro)", "East Java", -7.45, 111.60, "rice_lowland", "irrigated", "java_double", "java", notes="Bengawan Solo irrigation command area.", ), IndonesiaZone( "banten_rice", "Banten lowlands", "Banten", -6.35, 106.10, "rice_lowland", "irrigated", "java_double", "java", ), IndonesiaZone( "lampung_rice", "Lampung lowlands", "Lampung", -5.00, 105.30, "rice_lowland", "irrigated", "sumatra_double", "sumatra_south", notes="Way Sekampung / Way Rarem irrigation; also major cassava/maize area.", ), IndonesiaZone( "south_sumatra_rice", "South Sumatra lowlands (Palembang)", "South Sumatra", -3.20, 104.70, "rice_lowland", "rainfed", "sumatra_double", "sumatra_south", notes="Large tidal-swamp (pasang surut) rice area; drainage/flood risk is structural.", ), IndonesiaZone( "west_sumatra_rice", "West Sumatra valleys", "West Sumatra", -0.60, 100.60, "rice_lowland", "irrigated", "sumatra_equatorial", "sumatra_north", ), IndonesiaZone( "north_sumatra_rice", "North Sumatra lowlands (Deli Serdang)", "North Sumatra", 3.30, 98.90, "rice_lowland", "irrigated", "sumatra_equatorial", "sumatra_north", notes="Wet season peaks Oct-Dec here -- calendar shifted vs Java.", ), IndonesiaZone( "south_sulawesi_rice", "South Sulawesi (Bone-Wajo)", "South Sulawesi", -4.60, 120.20, "rice_lowland", "irrigated", "java_double", "sulawesi", notes="Sidrap-Wajo-Bone is Sulawesi's main rice surplus area.", ), IndonesiaZone( "central_kalimantan_rice", "Central Kalimantan lowlands", "Central Kalimantan", -2.40, 113.90, "rice_lowland", "rainfed", "sumatra_equatorial", "kalimantan", notes="Peatland-adjacent; former mega-rice area. Drainage and fire risk both relevant.", ), IndonesiaZone( "bali_rice", "Bali subak lowlands", "Bali", -8.50, 115.20, "rice_lowland", "irrigated", "java_double", "bali_nt", notes="Subak cooperative irrigation (UNESCO); strong dry season Jun-Sep.", ), IndonesiaZone( "lombok_rice", "Lombok lowlands", "West Nusa Tenggara", -8.65, 116.30, "rice_lowland", "irrigated", "java_double", "bali_nt", ), IndonesiaZone( "kupang_dryland", "Kupang dryland (maize)", "East Nusa Tenggara", -10.20, 123.60, "maize", "rainfed", "ntt_single", "ntt", notes="Single rainfed crop; strongest monsoon seasonality in Indonesia; chronic dry-season water deficit.", ), ) _ZONE_BY_ID: Dict[str, IndonesiaZone] = {z.zone_id: z for z in INDONESIA_ZONES} INDONESIA_BBOX = (-11.0, 6.0, 95.0, 141.0) def get_zone(zone_id: str) -> IndonesiaZone: if zone_id not in _ZONE_BY_ID: raise KeyError( f"Unknown Indonesian zone '{zone_id}'. " f"Registered: {sorted(_ZONE_BY_ID)}" ) return _ZONE_BY_ID[zone_id] def register_indonesia_zones() -> List[str]: from era5_data_pipeline import register_zone ids = [] for z in INDONESIA_ZONES: register_zone(z.to_polygon()) ids.append(z.zone_id) logger.info("register_indonesia_zones: %d zones registered", len(ids)) return ids # --------------------------------------------------------------------------- # Crop calendars # --------------------------------------------------------------------------- @dataclass(frozen=True) class RiceSeason: name: str plant_start: int plant_end: int harvest_start: int harvest_end: int CROP_CALENDARS: Dict[str, Tuple[RiceSeason, ...]] = { "java_double": ( RiceSeason("wet_rice", plant_start=305, plant_end=365, harvest_start=46, harvest_end=105), # Nov -> Feb/Mar RiceSeason("dry_rice", plant_start=105, plant_end=151, harvest_start=213, harvest_end=258), # Apr/May -> Aug/Sep ), "sumatra_double": ( RiceSeason("wet_rice", plant_start=290, plant_end=350, harvest_start=31, harvest_end=90), RiceSeason("dry_rice", plant_start=100, plant_end=146, harvest_start=205, harvest_end=250), ), "sumatra_equatorial": ( RiceSeason("main_rice", plant_start=274, plant_end=334, harvest_start=15, harvest_end=75), RiceSeason("second_rice", plant_start=90, plant_end=135, harvest_start=195, harvest_end=240), ), "ntt_single": ( RiceSeason("rainfed_main", plant_start=335, plant_end=31, harvest_start=100, harvest_end=140), ), } _PHASE_FRACTIONS: Tuple[Tuple[CropStage, float], ...] = ( (CropStage.VEGETATIVE, 0.45), (CropStage.REPRODUCTIVE, 0.25), (CropStage.GRAIN_FILLING, 0.20), (CropStage.MATURATION, 0.10), ) def _in_doy_window(doy: int, start: int, end: int) -> bool: if start <= end: return start <= doy <= end return doy >= start or doy <= end def _doy_distance_forward(from_doy: int, to_doy: int) -> int: return (to_doy - from_doy) % 366 if (to_doy - from_doy) % 366 != 0 else 0 def _window_len(start: int, end: int) -> int: return (end - start) % 366 + 1 def crop_stage_for_date( zone_id: str, dt: datetime, ) -> Tuple[CropStage, Optional[int], Optional[str]]: z = get_zone(zone_id) doy = dt.timetuple().tm_yday for season in CROP_CALENDARS[z.calendar]: plant_len = _window_len(season.plant_start, season.plant_end) if _in_doy_window(doy, season.plant_start, season.plant_end): mid_plant = (season.plant_start + plant_len // 2) % 366 or 366 grow_len = _doy_distance_forward(mid_plant, season.harvest_start) dth = _doy_distance_forward(doy, season.harvest_start) return CropStage.PLANTING, max(0, min(dth, grow_len + plant_len)), season.name prep_start = (season.plant_start - 14) % 366 or 366 if _in_doy_window(doy, prep_start, season.plant_start): return CropStage.LAND_PREP, None, season.name grow_len = _doy_distance_forward(season.plant_end, season.harvest_start) if grow_len > 0 and _in_doy_window(doy, season.plant_end, season.harvest_start): elapsed = _doy_distance_forward(season.plant_end, doy) frac = elapsed / max(grow_len, 1) acc = 0.0 stage = CropStage.MATURATION for s, f in _PHASE_FRACTIONS: acc += f if frac <= acc: stage = s break return stage, _doy_distance_forward(doy, season.harvest_start), season.name if _in_doy_window(doy, season.harvest_start, season.harvest_end): return CropStage.HARVEST, 0, season.name return CropStage.FALLOW, None, None # --------------------------------------------------------------------------- # Monsoon onset # --------------------------------------------------------------------------- _ONSET_CLIMATOLOGY: Dict[str, Tuple[int, int]] = { "ntt": (320, 18), # mid-November "bali_nt": (325, 16), "java": (330, 15), # late Nov / early Dec "sulawesi": (330, 16), "kalimantan": (305, 16), "sumatra_south": (300, 15), # late Oct / early Nov "sumatra_north": (285, 15), # mid-Oct } @dataclass(frozen=True) class MonsoonOnset: zone_id: str season_year: int # the year the wet season STARTS in base_onset_doy: int adjusted_onset_doy: float onset_std_days: int enso_adjustment_days: float iod_adjustment_days: float confidence: str # 'climatological-heuristic' notes: str = "" def to_dict(self) -> Dict[str, Any]: return asdict(self) def monsoon_onset_estimate( zone_id: str, year: int, basin: Optional[BasinContext] = None, ) -> MonsoonOnset: z = get_zone(zone_id) base_doy, std = _ONSET_CLIMATOLOGY[z.region_group] enso_adj = 0.0 iod_adj = 0.0 if basin is not None: if basin.enso_oni > 0.5: enso_adj = _clip(8.0 * basin.enso_oni, 0.0, 20.0) elif basin.enso_oni < -0.5: enso_adj = _clip(8.0 * basin.enso_oni, -20.0, 0.0) iod_scale = 3.0 if z.region_group.startswith("sumatra") else 6.0 if basin.iod_dmi > 0.4: iod_adj = _clip(iod_scale * basin.iod_dmi, 0.0, 15.0) elif basin.iod_dmi < -0.4: iod_adj = _clip(iod_scale * basin.iod_dmi, -15.0, 0.0) adjusted = base_doy + enso_adj + iod_adj return MonsoonOnset( zone_id=zone_id, season_year=year, base_onset_doy=base_doy, adjusted_onset_doy=adjusted, onset_std_days=std, enso_adjustment_days=enso_adj, iod_adjustment_days=iod_adj, confidence="climatological-heuristic", notes=( f"base={base_doy} ({z.region_group}), " f"ENSO {enso_adj:+.0f}d, IOD {iod_adj:+.0f}d" + (" (neutral basin context)" if basin is None else "") ), ) # --------------------------------------------------------------------------- # Planting advisory (the "planning" output) # --------------------------------------------------------------------------- @dataclass class PlantingAdvisory: zone_id: str generated_at: datetime # Current crop-calendar position crop_stage_now: str days_to_harvest_now: Optional[int] season_now: Optional[str] # Recommended planting window (None when none found in horizon) window_found: bool recommended_window_start: Optional[datetime] = None recommended_window_end: Optional[datetime] = None window_precip_total_mm: float = 0.0 # Water constraint soil_ready: Optional[bool] = None # None = soil state unknown irrigation_needed: bool = False water_note: str = "" # Context monsoon_onset: Optional[Dict[str, Any]] = None risk_alert_level: str = "none" risk_action_notes: str = "" method_notes: str = "" def to_dict(self) -> Dict[str, Any]: d = asdict(self) d["generated_at"] = self.generated_at.isoformat() d["recommended_window_start"] = ( self.recommended_window_start.isoformat() if self.recommended_window_start else None ) d["recommended_window_end"] = ( self.recommended_window_end.isoformat() if self.recommended_window_end else None ) return d def planting_advisory( zone_id: str, obs: ZoneObs, forecast: ForecastResult, config: Optional[ForecastConfig] = None, basin: Optional[BasinContext] = None, ) -> PlantingAdvisory: RAINFED_MIN, RAINFED_MAX = 25.0, 150.0 IRRIGATED_MIN, IRRIGATED_MAX = 10.0, 200.0 DAY_CAP_MM = 80.0 DROUGHT_PROB_CAP = 0.6 WINDOW = 7 cfg = config or ForecastConfig() z = get_zone(zone_id) vt = obs.valid_time if obs.valid_time.tzinfo else obs.valid_time.replace(tzinfo=timezone.utc) stage, dth, season = crop_stage_for_date(zone_id, vt) p_min, p_max = (IRRIGATED_MIN, IRRIGATED_MAX) if z.irrigation == "irrigated" \ else (RAINFED_MIN, RAINFED_MAX) precip = list(forecast.precip_mm) drought_probs = list(forecast.prob_drought_day) if forecast.prob_drought_day else [0.0] * len(precip) win_start = win_end = None win_total = 0.0 for i in range(0, max(0, len(precip) - WINDOW + 1)): chunk = precip[i:i + WINDOW] total = sum(chunk) if not (p_min <= total <= p_max): continue if max(chunk) > DAY_CAP_MM: continue mean_drought = sum(drought_probs[i:i + WINDOW]) / WINDOW if mean_drought > DROUGHT_PROB_CAP: continue win_start = vt + timedelta(days=i) win_end = vt + timedelta(days=i + WINDOW) win_total = total break # --- Water constraint assessment --- soil_ready: Optional[bool] = None if obs.soil_moisture_pct > 0.0: soil_ready = obs.soil_moisture_pct >= 25.0 irrigation_needed = (z.irrigation == "rainfed") and (win_start is None) water_bits = [] if soil_ready is None: water_bits.append("soil moisture unknown (no observation)") else: water_bits.append( f"soil {'adequate' if soil_ready else 'dry'} ({obs.soil_moisture_pct:.0f}%)" ) if irrigation_needed: water_bits.append( "no qualifying wet window in forecast horizon -- rainfed planting " "should wait or plan supplemental irrigation" ) elif win_start is not None and z.irrigation == "rainfed": water_bits.append("forecast window meets rainfed land-prep minimum") # --- Monsoon context (when a planting window is upcoming) --- monsoon_dict = None for s in CROP_CALENDARS[z.calendar]: days_to_plant = _doy_distance_forward(vt.timetuple().tm_yday, s.plant_start) if days_to_plant <= 120: onset = monsoon_onset_estimate(zone_id, vt.year, basin) monsoon_dict = onset.to_dict() break # --- Risk summary (reuse the scorer; do not re-derive) --- from crop_risk_scorer import compute_risk_score risk = compute_risk_score(obs, forecast, cfg) return PlantingAdvisory( zone_id=zone_id, generated_at=datetime.now(timezone.utc), crop_stage_now=stage.value, days_to_harvest_now=dth, season_now=season, window_found=win_start is not None, recommended_window_start=win_start, recommended_window_end=win_end, window_precip_total_mm=round(win_total, 1), soil_ready=soil_ready, irrigation_needed=irrigation_needed, water_note="; ".join(water_bits), monsoon_onset=monsoon_dict, risk_alert_level=risk.alert_level.value, risk_action_notes=risk.action_notes, method_notes=( f"window rule: 7d total in [{p_min:.0f},{p_max:.0f}]mm, " f"daily cap {DAY_CAP_MM:.0f}mm, mean drought prob <= {DROUGHT_PROB_CAP}; " f"zone irrigation class: {z.irrigation}" ), ) # --------------------------------------------------------------------------- # Self-test (python indonesia_zones.py) -- fully offline # --------------------------------------------------------------------------- if __name__ == "__main__": import json import sys from zone_observation import ( make_synthetic_basin_context, make_synthetic_forecast_result, make_synthetic_zone_obs, ) logging.basicConfig(level=logging.WARNING) print("indonesia_zones.py self-test (offline)\n") failures: List[str] = [] def _assert(cond: bool, msg: str) -> None: if not cond: failures.append(msg) print(f" FAIL: {msg}") # 1. Registry integrity: zones valid, centroids inside Indonesia bbox _assert(len(INDONESIA_ZONES) >= 12, "registry should have >= 12 zones") for z in INDONESIA_ZONES: p = z.to_polygon() lat_min, lat_max, lon_min, lon_max = INDONESIA_BBOX _assert(lat_min - 1 <= z.lat <= lat_max + 1, f"{z.zone_id} lat outside bbox") _assert(lon_min - 1 <= z.lon <= lon_max + 1, f"{z.zone_id} lon outside bbox") _assert(p.approx_area_km2 > 100.0, f"{z.zone_id} polygon area implausible") print(f" Registry OK: {len(INDONESIA_ZONES)} zones, polygons valid") # 2. Calendar progression: a year of months hits the key stages stages_seen = set() for month in range(1, 13): dt = datetime(2025, month, 15, tzinfo=timezone.utc) stage, dth, season = crop_stage_for_date("karawang_rice", dt) stages_seen.add(stage) if stage == CropStage.VEGETATIVE: _assert(dth is not None and dth > 0, "growing stage needs days_to_harvest") if stage == CropStage.FALLOW: _assert(dth is None, "FALLOW must have days_to_harvest=None") for needed in (CropStage.PLANTING, CropStage.VEGETATIVE, CropStage.HARVEST): _assert(needed in stages_seen, f"stage {needed.value} never reached in a year") print(f" Calendar OK: stages seen = {sorted(s.value for s in stages_seen)}") # 3. days_to_harvest decreases as the season advances d1, dth1, _ = crop_stage_for_date("karawang_rice", datetime(2025, 1, 10, tzinfo=timezone.utc)) d2, dth2, _ = crop_stage_for_date("karawang_rice", datetime(2025, 1, 31, tzinfo=timezone.utc)) if dth1 and dth2: _assert(dth2 < dth1, f"days_to_harvest should decrease: {dth1} -> {dth2}") print(f" days_to_harvest monotonic OK: {dth1} -> {dth2}") # 4. Monsoon onset: NTT base + El Nino delay onset_neutral = monsoon_onset_estimate("kupang_dryland", 2025) _assert(onset_neutral.base_onset_doy == 320, "NTT base onset should be DOY 320") el_nino = BasinContext( valid_date=datetime(2025, 9, 1, tzinfo=timezone.utc), enso_oni=1.5, iod_dmi=1.0, ) onset_nino = monsoon_onset_estimate("kupang_dryland", 2025, el_nino) _assert(onset_nino.adjusted_onset_doy > onset_neutral.adjusted_onset_doy, "El Nino + positive IOD should delay onset") la_nina = BasinContext( valid_date=datetime(2025, 9, 1, tzinfo=timezone.utc), enso_oni=-1.2, iod_dmi=-0.8, ) onset_nina = monsoon_onset_estimate("kupang_dryland", 2025, la_nina) _assert(onset_nina.adjusted_onset_doy < onset_neutral.adjusted_onset_doy, "La Nina + negative IOD should advance onset") print(f" Monsoon OK: neutral={onset_neutral.adjusted_onset_doy:.0f} " f"nino={onset_nino.adjusted_onset_doy:.0f} nina={onset_nina.adjusted_onset_doy:.0f}") # 5. Advisory: normal (typical onset-season) forecast -> window found # for irrigated Java zone; a full FLOOD forecast must NOT qualify # (you do not transplant into a flood -- the daily cap exists for # exactly this). obs = make_synthetic_zone_obs("karawang_rice", seed=7) wet_fc = make_synthetic_forecast_result("karawang_rice", valid_time=obs.valid_time, seed=7) adv = planting_advisory("karawang_rice", obs, wet_fc) _assert(adv.window_found, "normal forecast should yield a planting window (irrigated)") if adv.window_found: _assert(adv.recommended_window_start is not None, "window start missing") _assert(adv.recommended_window_end > adv.recommended_window_start, "window order") d = adv.to_dict() json.dumps(d) # must be JSON-serialisable print(f" Advisory (normal) OK: window {d['recommended_window_start'][:10]} " f"-> {d['recommended_window_end'][:10]} ({adv.window_precip_total_mm}mm)") flood_fc = make_synthetic_forecast_result("karawang_rice", valid_time=obs.valid_time, flood=True, seed=7) adv_flood = planting_advisory("karawang_rice", obs, flood_fc) _assert(not adv_flood.window_found, "flood forecast should NOT yield a planting window (daily cap)") print(f" Advisory (flood-rejected) OK: window_found={adv_flood.window_found}") # 6. Advisory: dry forecast on rainfed NTT zone -> irrigation flag obs_d = make_synthetic_zone_obs("kupang_dryland", drought=True, seed=8) dry_fc = make_synthetic_forecast_result("kupang_dryland", valid_time=obs_d.valid_time, drought=True, seed=8) adv_dry = planting_advisory("kupang_dryland", obs_d, dry_fc) _assert(adv_dry.irrigation_needed, "rainfed zone + dry forecast should flag irrigation") _assert(not adv_dry.window_found, "dry forecast should NOT yield a rainfed window") _assert(adv_dry.soil_ready is False, "drought obs soil should read as not ready") print(f" Advisory (dry) OK: irrigation_needed={adv_dry.irrigation_needed} " f"note='{adv_dry.water_note[:60]}...'") # 7. Advisory embeds scorer alert level _assert(adv_dry.risk_alert_level in ("none", "watch", "advisory", "warning", "critical"), "advisory should embed a valid alert level") print(f" Advisory risk embedding OK (alert={adv_dry.risk_alert_level})") # 8. Registration with the pipeline (import-time optional) try: ids = register_indonesia_zones() _assert(len(ids) == len(INDONESIA_ZONES), "registered id count mismatch") from era5_data_pipeline import _resolve_latlon lat, lon = _resolve_latlon("karawang_rice") _assert(abs(lat - (-6.30)) < 0.01 and abs(lon - 107.30) < 0.01, "pipeline centroid resolution wrong") print(f" Pipeline registration OK ({len(ids)} zones)") except ImportError as e: print(f" Pipeline registration SKIPPED ({e})") print() if failures: print(f"FAILED {len(failures)} test(s):") for f in failures: print(f" - {f}") sys.exit(1) else: print("All 8 test groups passed.")