DHDRL commited on
Commit
cabae60
·
verified ·
1 Parent(s): aea11ca

Update indonesia_zones.py

Browse files
Files changed (1) hide show
  1. indonesia_zones.py +2 -152
indonesia_zones.py CHANGED
@@ -3,46 +3,6 @@ indonesia_zones.py
3
  ==================
4
  Indonesia grounding layer: real agricultural zone registry, rice crop
5
  calendars, monsoon-onset context, and planting-window planning outputs.
6
-
7
- WHY THIS MODULE EXISTS
8
- ----------------------
9
- The rest of the pipeline (era5_data_pipeline, crop_risk_scorer,
10
- weather_forecast_env, hierarchical_search) is fully generic -- it knows
11
- what a ZoneObs is but nothing about WHERE the zones are or WHAT is growing
12
- there. Concretely, before this module:
13
-
14
- * No Indonesian zone was registered anywhere; register_zone() was only
15
- ever called with ad-hoc polygons (tests, hierarchical_search nodes).
16
- * Real fetchers never populated crop_stage / days_to_harvest, so
17
- crop_risk_scorer._optimal_harvest_window never fired on real data and
18
- _harvest_pressure always used the 0.3 default branch.
19
- * "Planning" outputs did not exist: the system produced alert levels but
20
- no actionable windows (when to plant, whether irrigation is needed).
21
-
22
- This module is deliberately self-contained: it depends only on
23
- zone_observation (+ crop_risk_scorer for the advisory's risk summary) so it
24
- can be imported by era5_data_pipeline callers, the backtester, or a service
25
- layer without pulling in the RL stack.
26
-
27
- DATA PROVENANCE (honest boundaries)
28
- -----------------------------------
29
- * Zone geometries: rectangular approximations around well-known
30
- agricultural-region centroids (+/-0.35 deg). They are NOT administrative
31
- boundaries. For production alerting routed to real officials, replace
32
- to_polygon() geometries with BPS/BIG (Badan Informasi Geospasial)
33
- regency-level shapefiles.
34
- * Crop calendars: heuristic two-season irrigated-rice calendar for Java /
35
- Sumatra / Sulawesi lowlands and a single rainfed season for the
36
- south-eastern islands (NTB/NTT), consistent with the broadly documented
37
- wet (Oct-Apr) / dry (May-Sep) season structure. Real planting dates
38
- vary by regency, variety (e.g. Ciherang vs Inpari), and irrigation
39
- district (bendungan/waduk service area). Treat stages as
40
- climatological-typical, not field-observed.
41
- * Monsoon onset: climatological heuristic with ENSO/IOD sign adjustments
42
- (El Nino / positive IOD delay onset over the southern archipelago;
43
- La Nina / negative IOD advance it). BMKG publishes official onset
44
- forecasts per ZOM (zona musim) each year -- wiring THAT product in is
45
- the production upgrade; see monsoon_onset_estimate() docstring.
46
  """
47
 
48
  from __future__ import annotations
@@ -81,11 +41,6 @@ logger = logging.getLogger(__name__)
81
 
82
  @dataclass(frozen=True)
83
  class IndonesiaZone:
84
- """One agricultural zone in the Indonesian registry.
85
-
86
- half_extent_deg: the rectangle half-width in degrees (applied to both
87
- lat and lon). 0.35 deg ~= 39 km -- roughly regency scale in Java.
88
- """
89
  zone_id: str
90
  label: str
91
  province: str
@@ -99,8 +54,6 @@ class IndonesiaZone:
99
  notes: str = ""
100
 
101
  def to_polygon(self) -> GeoPolygon:
102
- """Rectangular approximation. NOT an administrative boundary --
103
- see module docstring for the BPS/BIG production upgrade path."""
104
  h = self.half_extent_deg
105
  return GeoPolygon(
106
  zone_id=self.zone_id,
@@ -187,8 +140,7 @@ INDONESIA_ZONES: Tuple[IndonesiaZone, ...] = (
187
 
188
  _ZONE_BY_ID: Dict[str, IndonesiaZone] = {z.zone_id: z for z in INDONESIA_ZONES}
189
 
190
- # Indonesia bounding box (matches hierarchical_search.py's self-test value).
191
- INDONESIA_BBOX = (-11.0, 6.0, 95.0, 141.0) # (lat_min, lat_max, lon_min, lon_max)
192
 
193
 
194
  def get_zone(zone_id: str) -> IndonesiaZone:
@@ -201,11 +153,6 @@ def get_zone(zone_id: str) -> IndonesiaZone:
201
 
202
 
203
  def register_indonesia_zones() -> List[str]:
204
- """Register every INDONESIA_ZONES polygon with era5_data_pipeline's zone
205
- registry (idempotent). Returns the registered zone_ids.
206
-
207
- Local import so this module stays importable without the pipeline.
208
- """
209
  from era5_data_pipeline import register_zone
210
  ids = []
211
  for z in INDONESIA_ZONES:
@@ -221,14 +168,6 @@ def register_indonesia_zones() -> List[str]:
221
 
222
  @dataclass(frozen=True)
223
  class RiceSeason:
224
- """One cropping season, in day-of-year windows (1..366).
225
-
226
- plant_*: the typical transplanting/direct-seeding window.
227
- harvest_*: the typical harvest window. Growing duration between them is
228
- derived, not stored. Windows may wrap around year end (wet-season rice
229
- planted in November is harvested the following February) -- handled by
230
- the _in_doy_window helper.
231
- """
232
  name: str
233
  plant_start: int
234
  plant_end: int
@@ -236,42 +175,31 @@ class RiceSeason:
236
  harvest_end: int
237
 
238
 
239
- # Heuristic calendars (see module docstring for provenance and limits).
240
  CROP_CALENDARS: Dict[str, Tuple[RiceSeason, ...]] = {
241
- # Java / Bali / Sulawesi lowland irrigated: wet + dry season rice.
242
  "java_double": (
243
  RiceSeason("wet_rice", plant_start=305, plant_end=365,
244
  harvest_start=46, harvest_end=105), # Nov -> Feb/Mar
245
  RiceSeason("dry_rice", plant_start=105, plant_end=151,
246
  harvest_start=213, harvest_end=258), # Apr/May -> Aug/Sep
247
  ),
248
- # Southern Sumatra: wet season slightly earlier, same double structure.
249
  "sumatra_double": (
250
  RiceSeason("wet_rice", plant_start=290, plant_end=350,
251
  harvest_start=31, harvest_end=90),
252
  RiceSeason("dry_rice", plant_start=100, plant_end=146,
253
  harvest_start=205, harvest_end=250),
254
  ),
255
- # Equatorial Sumatra / Kalimantan: wet season peaks Oct-Dec; less
256
- # pronounced dry season -- second crop is partial (represented as a
257
- # narrower, less reliable window).
258
  "sumatra_equatorial": (
259
  RiceSeason("main_rice", plant_start=274, plant_end=334,
260
  harvest_start=15, harvest_end=75),
261
  RiceSeason("second_rice", plant_start=90, plant_end=135,
262
  harvest_start=195, harvest_end=240),
263
  ),
264
- # NTT single rainfed crop (maize/upland rice): plant with the onset
265
- # rains Dec-Jan, harvest Apr-May.
266
  "ntt_single": (
267
  RiceSeason("rainfed_main", plant_start=335, plant_end=31,
268
  harvest_start=100, harvest_end=140),
269
  ),
270
  }
271
 
272
- # Phase fractions of the growing period (plant_end -> harvest_start):
273
- # vegetative first ~45%, reproductive next ~25%, grain filling next ~20%,
274
- # maturation the final ~10% before harvest opens.
275
  _PHASE_FRACTIONS: Tuple[Tuple[CropStage, float], ...] = (
276
  (CropStage.VEGETATIVE, 0.45),
277
  (CropStage.REPRODUCTIVE, 0.25),
@@ -281,15 +209,12 @@ _PHASE_FRACTIONS: Tuple[Tuple[CropStage, float], ...] = (
281
 
282
 
283
  def _in_doy_window(doy: int, start: int, end: int) -> bool:
284
- """True if doy is inside [start, end], with wrap-around support
285
- (e.g. start=335, end=31 covers Dec -> Jan)."""
286
  if start <= end:
287
  return start <= doy <= end
288
  return doy >= start or doy <= end
289
 
290
 
291
  def _doy_distance_forward(from_doy: int, to_doy: int) -> int:
292
- """Number of days from from_doy forward to to_doy (wrap-aware, 0..365)."""
293
  return (to_doy - from_doy) % 366 if (to_doy - from_doy) % 366 != 0 else 0
294
 
295
 
@@ -301,40 +226,22 @@ def crop_stage_for_date(
301
  zone_id: str,
302
  dt: datetime,
303
  ) -> Tuple[CropStage, Optional[int], Optional[str]]:
304
- """Climatological-typical crop stage for a zone on a date.
305
-
306
- Returns (crop_stage, days_to_harvest, season_name).
307
- * days_to_harvest is None when the zone is FALLOW or the crop is
308
- already inside the harvest window (0 would mean "harvest today";
309
- we return 0 then, not None -- None means "no crop in the ground").
310
- * season_name is None when FALLOW.
311
-
312
- This is a CALENDAR model, not a field observation. It exists so real
313
- pipeline fetches can populate crop_stage / days_to_harvest (which the
314
- scorer's harvest-window and harvest-pressure logic needs) instead of
315
- leaving them UNKNOWN/None forever on real data.
316
- """
317
  z = get_zone(zone_id)
318
  doy = dt.timetuple().tm_yday
319
 
320
  for season in CROP_CALENDARS[z.calendar]:
321
  plant_len = _window_len(season.plant_start, season.plant_end)
322
 
323
- # Inside the planting window -> PLANTING; harvest is one growing
324
- # period ahead of the END of the window (approximation: crop goes in
325
- # mid-window).
326
  if _in_doy_window(doy, season.plant_start, season.plant_end):
327
  mid_plant = (season.plant_start + plant_len // 2) % 366 or 366
328
  grow_len = _doy_distance_forward(mid_plant, season.harvest_start)
329
  dth = _doy_distance_forward(doy, season.harvest_start)
330
  return CropStage.PLANTING, max(0, min(dth, grow_len + plant_len)), season.name
331
 
332
- # Two weeks before the planting window -> LAND_PREP.
333
  prep_start = (season.plant_start - 14) % 366 or 366
334
  if _in_doy_window(doy, prep_start, season.plant_start):
335
  return CropStage.LAND_PREP, None, season.name
336
 
337
- # Growing period: plant_end -> harvest_start, phase by fraction.
338
  grow_len = _doy_distance_forward(season.plant_end, season.harvest_start)
339
  if grow_len > 0 and _in_doy_window(doy, season.plant_end, season.harvest_start):
340
  elapsed = _doy_distance_forward(season.plant_end, doy)
@@ -348,7 +255,6 @@ def crop_stage_for_date(
348
  break
349
  return stage, _doy_distance_forward(doy, season.harvest_start), season.name
350
 
351
- # Harvest window itself.
352
  if _in_doy_window(doy, season.harvest_start, season.harvest_end):
353
  return CropStage.HARVEST, 0, season.name
354
 
@@ -359,14 +265,7 @@ def crop_stage_for_date(
359
  # Monsoon onset
360
  # ---------------------------------------------------------------------------
361
 
362
- # Climatological wet-season onset day-of-year by region group, with the
363
- # typical interannual std. HEURISTIC values consistent with the broadly
364
- # documented SE->NW onset progression over the maritime continent (earliest
365
- # over the south-eastern islands in November, reaching N Sumatra as its
366
- # Oct-Dec rainfall peak). BMKG's official per-ZOM onset forecast is the
367
- # production replacement -- see monsoon_onset_estimate() docstring.
368
  _ONSET_CLIMATOLOGY: Dict[str, Tuple[int, int]] = {
369
- # region_group: (base_onset_doy, std_days)
370
  "ntt": (320, 18), # mid-November
371
  "bali_nt": (325, 16),
372
  "java": (330, 15), # late Nov / early Dec
@@ -379,7 +278,6 @@ _ONSET_CLIMATOLOGY: Dict[str, Tuple[int, int]] = {
379
 
380
  @dataclass(frozen=True)
381
  class MonsoonOnset:
382
- """Wet-season onset estimate for one zone and season year."""
383
  zone_id: str
384
  season_year: int # the year the wet season STARTS in
385
  base_onset_doy: int
@@ -399,28 +297,6 @@ def monsoon_onset_estimate(
399
  year: int,
400
  basin: Optional[BasinContext] = None,
401
  ) -> MonsoonOnset:
402
- """Estimate wet-season onset for a zone.
403
-
404
- Base value is the regional climatological onset (see
405
- _ONSET_CLIMATOLOGY). When a BasinContext is supplied, its ENSO/IOD
406
- phases shift the estimate with the DOCUMENTED sign of the
407
- teleconnection and an explicitly heuristic magnitude:
408
-
409
- * El Nino (ONI > +0.5): onset delayed over the southern archipelago.
410
- +8 days per +1.0 degC ONI anomaly (clipped to [0, +20]).
411
- * La Nina (ONI < -0.5): onset advanced. Same magnitude, negative.
412
- * Positive IOD (DMI > +0.4): suppressed rainfall over the SE islands
413
- -> additional delay of +6 days per +1.0 degC DMI (clipped [0, +15]),
414
- halved for Sumatra (weaker IOD coupling that far west).
415
- * Negative IOD: mirrored advance.
416
-
417
- These magnitudes are order-of-magnitude heuristics, NOT fitted
418
- regressions. The production upgrade is BMKG's official per-ZOM (zona
419
- musim) onset forecast, published ahead of each wet season; this
420
- function's contract (base + teleconnection adjustment) is designed so
421
- a BMKG product can replace `adjusted_onset_doy` without changing
422
- callers.
423
- """
424
  z = get_zone(zone_id)
425
  base_doy, std = _ONSET_CLIMATOLOGY[z.region_group]
426
 
@@ -461,12 +337,6 @@ def monsoon_onset_estimate(
461
 
462
  @dataclass
463
  class PlantingAdvisory:
464
- """Actionable planting-window recommendation for one zone.
465
-
466
- This is the module's planning output: not an alert (crop_risk_scorer
467
- already does those) but a forward-looking operational recommendation
468
- answering "when should planting happen here, and is water a constraint?".
469
- """
470
  zone_id: str
471
  generated_at: datetime
472
 
@@ -514,25 +384,6 @@ def planting_advisory(
514
  config: Optional[ForecastConfig] = None,
515
  basin: Optional[BasinContext] = None,
516
  ) -> PlantingAdvisory:
517
- """Build a PlantingAdvisory for a zone from its current obs + forecast.
518
-
519
- Window-selection rule (documented heuristic, tuned for rice land
520
- preparation at the start of the wet season):
521
-
522
- * Scan the forecast horizon with a 7-day sliding window.
523
- * A window QUALIFIES when:
524
- - cumulative precip inside the window is in [MIN, MAX] mm
525
- (rainfed: [25, 150] -- enough to rewet the profile, not so much
526
- that puddling/transplanting is disrupted;
527
- irrigated: [10, 200] -- canal supply relaxes the rain floor),
528
- - no single day exceeds 80 mm (waterlogging/transplant damage),
529
- - mean prob_drought_day over the window <= 0.6.
530
- * The FIRST qualifying window is returned (earliest acceptable start).
531
-
532
- These are agronomic heuristics, not calibrated values; they are
533
- constants at the top of the function precisely so they can be tuned
534
- against field data later.
535
- """
536
  RAINFED_MIN, RAINFED_MAX = 25.0, 150.0
537
  IRRIGATED_MIN, IRRIGATED_MAX = 10.0, 200.0
538
  DAY_CAP_MM = 80.0
@@ -545,7 +396,6 @@ def planting_advisory(
545
 
546
  stage, dth, season = crop_stage_for_date(zone_id, vt)
547
 
548
- # --- Scan forecast for a qualifying window ---
549
  p_min, p_max = (IRRIGATED_MIN, IRRIGATED_MAX) if z.irrigation == "irrigated" \
550
  else (RAINFED_MIN, RAINFED_MAX)
551
  precip = list(forecast.precip_mm)
@@ -758,4 +608,4 @@ if __name__ == "__main__":
758
  print(f" - {f}")
759
  sys.exit(1)
760
  else:
761
- print("All 8 test groups passed.")
 
3
  ==================
4
  Indonesia grounding layer: real agricultural zone registry, rice crop
5
  calendars, monsoon-onset context, and planting-window planning outputs.
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
6
  """
7
 
8
  from __future__ import annotations
 
41
 
42
  @dataclass(frozen=True)
43
  class IndonesiaZone:
 
 
 
 
 
44
  zone_id: str
45
  label: str
46
  province: str
 
54
  notes: str = ""
55
 
56
  def to_polygon(self) -> GeoPolygon:
 
 
57
  h = self.half_extent_deg
58
  return GeoPolygon(
59
  zone_id=self.zone_id,
 
140
 
141
  _ZONE_BY_ID: Dict[str, IndonesiaZone] = {z.zone_id: z for z in INDONESIA_ZONES}
142
 
143
+ INDONESIA_BBOX = (-11.0, 6.0, 95.0, 141.0)
 
144
 
145
 
146
  def get_zone(zone_id: str) -> IndonesiaZone:
 
153
 
154
 
155
  def register_indonesia_zones() -> List[str]:
 
 
 
 
 
156
  from era5_data_pipeline import register_zone
157
  ids = []
158
  for z in INDONESIA_ZONES:
 
168
 
169
  @dataclass(frozen=True)
170
  class RiceSeason:
 
 
 
 
 
 
 
 
171
  name: str
172
  plant_start: int
173
  plant_end: int
 
175
  harvest_end: int
176
 
177
 
 
178
  CROP_CALENDARS: Dict[str, Tuple[RiceSeason, ...]] = {
 
179
  "java_double": (
180
  RiceSeason("wet_rice", plant_start=305, plant_end=365,
181
  harvest_start=46, harvest_end=105), # Nov -> Feb/Mar
182
  RiceSeason("dry_rice", plant_start=105, plant_end=151,
183
  harvest_start=213, harvest_end=258), # Apr/May -> Aug/Sep
184
  ),
 
185
  "sumatra_double": (
186
  RiceSeason("wet_rice", plant_start=290, plant_end=350,
187
  harvest_start=31, harvest_end=90),
188
  RiceSeason("dry_rice", plant_start=100, plant_end=146,
189
  harvest_start=205, harvest_end=250),
190
  ),
 
 
 
191
  "sumatra_equatorial": (
192
  RiceSeason("main_rice", plant_start=274, plant_end=334,
193
  harvest_start=15, harvest_end=75),
194
  RiceSeason("second_rice", plant_start=90, plant_end=135,
195
  harvest_start=195, harvest_end=240),
196
  ),
 
 
197
  "ntt_single": (
198
  RiceSeason("rainfed_main", plant_start=335, plant_end=31,
199
  harvest_start=100, harvest_end=140),
200
  ),
201
  }
202
 
 
 
 
203
  _PHASE_FRACTIONS: Tuple[Tuple[CropStage, float], ...] = (
204
  (CropStage.VEGETATIVE, 0.45),
205
  (CropStage.REPRODUCTIVE, 0.25),
 
209
 
210
 
211
  def _in_doy_window(doy: int, start: int, end: int) -> bool:
 
 
212
  if start <= end:
213
  return start <= doy <= end
214
  return doy >= start or doy <= end
215
 
216
 
217
  def _doy_distance_forward(from_doy: int, to_doy: int) -> int:
 
218
  return (to_doy - from_doy) % 366 if (to_doy - from_doy) % 366 != 0 else 0
219
 
220
 
 
226
  zone_id: str,
227
  dt: datetime,
228
  ) -> Tuple[CropStage, Optional[int], Optional[str]]:
 
 
 
 
 
 
 
 
 
 
 
 
 
229
  z = get_zone(zone_id)
230
  doy = dt.timetuple().tm_yday
231
 
232
  for season in CROP_CALENDARS[z.calendar]:
233
  plant_len = _window_len(season.plant_start, season.plant_end)
234
 
 
 
 
235
  if _in_doy_window(doy, season.plant_start, season.plant_end):
236
  mid_plant = (season.plant_start + plant_len // 2) % 366 or 366
237
  grow_len = _doy_distance_forward(mid_plant, season.harvest_start)
238
  dth = _doy_distance_forward(doy, season.harvest_start)
239
  return CropStage.PLANTING, max(0, min(dth, grow_len + plant_len)), season.name
240
 
 
241
  prep_start = (season.plant_start - 14) % 366 or 366
242
  if _in_doy_window(doy, prep_start, season.plant_start):
243
  return CropStage.LAND_PREP, None, season.name
244
 
 
245
  grow_len = _doy_distance_forward(season.plant_end, season.harvest_start)
246
  if grow_len > 0 and _in_doy_window(doy, season.plant_end, season.harvest_start):
247
  elapsed = _doy_distance_forward(season.plant_end, doy)
 
255
  break
256
  return stage, _doy_distance_forward(doy, season.harvest_start), season.name
257
 
 
258
  if _in_doy_window(doy, season.harvest_start, season.harvest_end):
259
  return CropStage.HARVEST, 0, season.name
260
 
 
265
  # Monsoon onset
266
  # ---------------------------------------------------------------------------
267
 
 
 
 
 
 
 
268
  _ONSET_CLIMATOLOGY: Dict[str, Tuple[int, int]] = {
 
269
  "ntt": (320, 18), # mid-November
270
  "bali_nt": (325, 16),
271
  "java": (330, 15), # late Nov / early Dec
 
278
 
279
  @dataclass(frozen=True)
280
  class MonsoonOnset:
 
281
  zone_id: str
282
  season_year: int # the year the wet season STARTS in
283
  base_onset_doy: int
 
297
  year: int,
298
  basin: Optional[BasinContext] = None,
299
  ) -> MonsoonOnset:
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
300
  z = get_zone(zone_id)
301
  base_doy, std = _ONSET_CLIMATOLOGY[z.region_group]
302
 
 
337
 
338
  @dataclass
339
  class PlantingAdvisory:
 
 
 
 
 
 
340
  zone_id: str
341
  generated_at: datetime
342
 
 
384
  config: Optional[ForecastConfig] = None,
385
  basin: Optional[BasinContext] = None,
386
  ) -> PlantingAdvisory:
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
387
  RAINFED_MIN, RAINFED_MAX = 25.0, 150.0
388
  IRRIGATED_MIN, IRRIGATED_MAX = 10.0, 200.0
389
  DAY_CAP_MM = 80.0
 
396
 
397
  stage, dth, season = crop_stage_for_date(zone_id, vt)
398
 
 
399
  p_min, p_max = (IRRIGATED_MIN, IRRIGATED_MAX) if z.irrigation == "irrigated" \
400
  else (RAINFED_MIN, RAINFED_MAX)
401
  precip = list(forecast.precip_mm)
 
608
  print(f" - {f}")
609
  sys.exit(1)
610
  else:
611
+ print("All 8 test groups passed.")