diff --git a/KT_DOCUMENT.md b/KT_DOCUMENT.md new file mode 100644 index 00000000..2fc1a580 --- /dev/null +++ b/KT_DOCUMENT.md @@ -0,0 +1,363 @@ +# Farm ET Pipeline — Knowledge Transfer Document + +## Overview + +This document covers the **Farm-level ET (Evapotranspiration) Pipeline** — the system that intersects AET/PET rasters with farm boundary polygons to compute **Moisture Adequacy Index (MAI)** and kharif water stress indicators per farm. + +The pipeline has been validated against Shuvam Chakraborty's GEE-based drought analysis rasters. This document explains the architecture, data flow, recent bug fixes, and how to run/extend the pipeline. + +--- + +## Table of Contents + +1. [Architecture & Data Flow](#1-architecture--data-flow) +2. [Directory Structure](#2-directory-structure) +3. [Raster Band Ordering (Critical)](#3-raster-band-ordering-critical) +4. [Pipeline Execution](#4-pipeline-execution) +5. [Key Functions in `et_intersection.py`](#5-key-functions-in-et_intersectionpy) +6. [Missing Data Handling & Gap-Fill](#6-missing-data-handling--gap-fill) +7. [Output Schema](#7-output-schema) +8. [Validation Against Shuvam's Rasters](#8-validation-against-shuvams-rasters) +9. [Recent Bug Fixes & Changes](#9-recent-bug-fixes--changes) +10. [How to Extend to New Regions/Years](#10-how-to-extend-to-new-regionsyears) +11. [Known Limitations](#11-known-limitations) + +--- + +## 1. Architecture & Data Flow + +``` +Input Rasters (COG GeoTIFF) Farm Boundaries (GeoParquet) + merge_AET___cog.tif farm_boundaries.parquet + merge_PET___cog.tif ├── farm_id + ├── 13 bands (crop-year order) ├── geometry (polygon) + ├── Band 1-12: monthly mm/day └── area_m2 + ├── Band 13: annual + └── NoData: -9999 + │ │ + └───────────────┬───────────────────────────┘ + │ + ┌───────────▼───────────┐ + │ et_intersection.py │ + │ │ + │ 1. Read rasters │ + │ 2. Mask nodata→NaN │ + │ 3. Zonal stats │ + │ 4. Gap-fill AET/PET │ + │ 5. Compute MAI │ + │ 6. Kharif stress │ + └───────────┬───────────┘ + │ + ┌─────────────────┼─────────────────┐ + ▼ ▼ ▼ + farm_static.parquet farm_annual.parquet farm_monthly.parquet + (geometry, area) (MAI, stress/yr) (AET, PET, MAI/month) +``` + +--- + +## 2. Directory Structure + +``` +core-stack-backend/ +├── computing/farm_boundaries/ +│ ├── et_intersection.py ← Main pipeline (Phase 3) +│ ├── fetch_raw.py ← Phase 1: download farm boundaries +│ ├── convert.py ← Phase 2: convert to GeoParquet +│ └── farm_boundary.py ← API / orchestrator +│ +├── data/ +│ ├── et_rasters/ +│ │ ├── merge_AET_4_2018_cog.tif (17.6 GB, AEZ zone 4, year 2018) +│ │ ├── merge_PET_4_2018_cog.tif (12.7 GB) +│ │ ├── shuvam_mai_sanganer_2018.tif (validation raster) +│ │ └── shuvam_mai_dudu_2018.tif (validation raster) +│ │ +│ └── farm_boundaries/rajasthan/jaipur/ +│ ├── sanganer/ +│ │ ├── farm_boundaries.parquet (94,214 farms) +│ │ ├── farm_static.parquet (geometry + area) +│ │ ├── farm_annual.parquet (annual MAI, kharif stress) +│ │ ├── farm_monthly.parquet (monthly AET, PET, MAI) +│ │ └── mai_tally_2018.parquet (validation output) +│ └── dudu/ +│ └── ... (same structure, 92,394 farms) +│ +├── compare_mai_tally.py ← Validation script vs Shuvam +├── plot_tally_timeseries.py ← Time-series comparison plots +└── regen_all_phase3.py ← Re-run Phase 3 for all tehsils +``` + +--- + +## 3. Raster Band Ordering (Critical) + +**The AET/PET rasters use crop-year band ordering, NOT calendar-year.** + +``` +Band 1 = July (year Y) +Band 2 = August (year Y) +Band 3 = September (year Y) +Band 4 = October (year Y) +Band 5 = November (year Y) +Band 6 = December (year Y) +Band 7 = January (year Y+1) +Band 8 = February (year Y+1) +Band 9 = March (year Y+1) +Band 10 = April (year Y+1) +Band 11 = May (year Y+1) +Band 12 = June (year Y+1) +Band 13 = Annual mean +``` + +This is encoded in the code as: +```python +CROP_YEAR_BAND_TO_MONTH = [7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6] +``` + +Where index `i` maps raster band `i` (0-indexed) to calendar month number. For example, band index 0 → month 7 (July). + +**Why this matters:** If you assume band 1 = January (calendar order), every monthly column will have the wrong month's data, and kharif stress will be computed from Oct–Jan instead of Jul–Oct. + +--- + +## 4. Pipeline Execution + +### Running for a specific tehsil + +```python +from computing.farm_boundaries.et_intersection import intersect_et_with_farms + +result = intersect_et_with_farms( + state="rajasthan", + district="jaipur", + block="sanganer", # tehsil name + year=2018, + overwrite=True # set True to regenerate +) +``` + +### Re-running for all tehsils + +```bash +python regen_all_phase3.py +``` + +### Prerequisites +- Conda environment: `corestackenv` +- `.env` file in project root with `LOCAL_ET_RASTERS_PATH` set +- Farm boundary parquets already generated (Phases 1 & 2) +- AET/PET rasters downloaded to `data/et_rasters/` + +--- + +## 5. Key Functions in `et_intersection.py` + +### `_read_raster_clipped(raster_path, bbox)` +- Reads a COG raster using windowed reading (only loads the bbox region, not the full 17 GB file) +- Converts nodata sentinel values to NaN using the raster's metadata (`src.nodata`), with a fallback to -9999 +- Returns `(data, transform)` where `data` is shape `(bands, height, width)` as float32 + +### `_extract_band_means(labels, band_data, num_farms)` +- Computes per-farm mean of one raster band using vectorised `np.bincount` +- Ignores NaN, Inf, and negative values +- Farms with zero valid pixels → NaN (not imputed) + +### `_gap_fill_monthly_farms(monthly_matrix)` +- Temporal interpolation on a `(n_farms, 12)` matrix in calendar order (col 0 = Jan, col 11 = Dec) +- Rules mirror Shuvam's `fill_monthly_collection()` — see Section 6 + +### `_run_zonal_stats(gdf, aet_data, aet_transform, pet_data, pet_transform)` +- Core function that: + 1. Rasterizes farm polygons onto the raster grid + 2. Extracts per-farm monthly AET using `CROP_YEAR_BAND_TO_MONTH` mapping + 3. Gap-fills AET, then PET + 4. Computes MAI = AET/PET (capped to [0, 1]) + 5. Computes kharif stress flags from Jul–Oct MAI values + +### `_save_monthly_parquet(gdf, state, district, block, year)` +- Converts wide-format columns to long format (one row per farm per month) +- Assigns correct calendar dates: Jul–Dec → year Y, Jan–Jun → year Y+1 + +--- + +## 6. Missing Data Handling & Gap-Fill + +### Nodata at Read Time +- Raster nodata value (typically -9999) is converted to NaN immediately when the raster is read +- Any residual values ≤ -9999 are also masked (belt-and-suspenders) + +### Temporal Gap-Fill Rules +Applied at the **farm level** on both AET and PET separately, before computing MAI. + +The rules respect the crop-year boundary (July → June) and mirror Shuvam's GEE logic: + +| Month | Fill source | Rationale | +|---|---|---| +| **July** | August only | Crop-year start — no backward crossing to previous June | +| **June** | May only | Crop-year end — no forward crossing to next July | +| **All others** | Mean of ±1 neighbouring month | Standard temporal interpolation | +| **No valid neighbour** | Stays NaN | Cannot impute without any reference | + +This is implemented in `_GAP_FILL_NEIGHBOURS` dict and `_gap_fill_monthly_farms()`. + +### MAI Capping +After computing MAI = AET/PET: +- Values > 1.0 are capped to 1.0 (AET cannot physically exceed PET; values > 1 are raster artifacts) +- Non-finite values → NaN +- A log warning is emitted when farms have MAI > 1 + +--- + +## 7. Output Schema + +### `farm_static.parquet` +| Column | Type | Description | +|---|---|---| +| farm_id | string | Unique farm identifier | +| geometry | geometry | Farm polygon | +| area_m2 | float | Farm area in square metres | +| bbox | dict | Bounding box {xmin, ymin, xmax, ymax} | + +### `farm_annual.parquet` +| Column | Type | Description | +|---|---|---| +| farm_id | string | Unique farm identifier | +| tehsil, district, state | string | Administrative hierarchy | +| area_in_ha | float | Farm area in hectares | +| year | int | Crop year (e.g. 2018 = Jul 2018 – Jun 2019) | +| aet_annual | float | Annual mean AET (mm/day) | +| pet_annual | float | Annual mean PET (mm/day) | +| mai_annual | float | Annual mean MAI [0, 1] | +| kharif_mai | float | Mean MAI for Jul–Oct | +| kharif_water_stress | bool | True if any kharif month MAI ≤ 0.50 | +| kharif_severe_stress | bool | True if any kharif month MAI ≤ 0.25 | + +### `farm_monthly.parquet` +| Column | Type | Description | +|---|---|---| +| farm_id | string | Unique farm identifier | +| tehsil, district, state | string | Administrative hierarchy | +| year | int | Crop year | +| date | datetime | Calendar date (e.g. 2018-07-01 for July) | +| aet | float | Monthly AET (mm/day) | +| pet | float | Monthly PET (mm/day) | +| mai | float | Monthly MAI [0, 1] | + +**Note:** For crop year 2018, `date` spans 2018-07-01 to 2019-06-01. + +--- + +## 8. Validation Against Shuvam's Rasters + +### Running the tally +```bash +python compare_mai_tally.py +``` + +This reads Shuvam's exported MAI GeoTIFFs (`shuvam_mai_sanganer_2018.tif`, `shuvam_mai_dudu_2018.tif`) and computes per-farm zonal means, then compares against our pipeline's output. + +### Current tally results (2018) + +| Metric | Sanganer (63,934 farms) | Dudu (80,273 farms) | +|---|---|---| +| Annual MAI — mean abs diff | 0.030 | 0.016 | +| Kharif MAI — mean abs diff | 0.036 | 0.019 | +| Farms within 0.01 (annual) | 30.8% | 47.7% | +| Farms within 0.01 (kharif) | 27.4% | 37.8% | + +### Sources of remaining difference +1. **Annual definition:** Shuvam's annual = crop-year mean (Jul–Jun); ours = calendar-year mean (Jan–Dec) +2. **Gap-fill granularity:** Shuvam applies gap-fill at the pixel level in GEE; we apply it at the farm level after zonal stats +3. **Spatial coverage:** Shuvam's raster clips to tehsil boundary in GEE; our rasters extend slightly beyond, so NaN counts differ + +### Generating time-series plots +```bash +python plot_tally_timeseries.py +``` +Output: `data/tally_plots/mai_timeseries_tally_sanganer_2018.png` and `gap_fill_and_tally_demo_2018.png` + +--- + +## 9. Recent Bug Fixes & Changes + +### Bug Fix 1: Band-to-Month Mapping (Critical) +**Problem:** The pipeline assumed raster band 1 = January (calendar-year order). The actual rasters use crop-year order where band 1 = July. + +**Impact:** Every monthly column (`aet_jan`, `mai_jul`, etc.) had the wrong month's data. Kharif stress was being computed from Oct–Jan instead of Jul–Oct. + +**Fix:** Added `CROP_YEAR_BAND_TO_MONTH` mapping array. Each raster band index is now explicitly mapped to its correct calendar month before assigning to the column. + +**File:** `et_intersection.py`, lines 78–87, 267–274, 299–302 + +--- + +### Bug Fix 2: Monthly Parquet Date Assignment +**Problem:** Monthly parquet assigned `date = 2018-01-01` for band 1, even though band 1 is actually July. All 12 months were dated within the same calendar year. + +**Impact:** Any downstream consumer of `farm_monthly.parquet` filtering by date would get wrong data. + +**Fix:** Jan–Jun months are now assigned to `year + 1`. For crop year 2018: Jul 2018 → Dec 2018 + Jan 2019 → Jun 2019. + +**File:** `et_intersection.py`, lines 453–467 + +--- + +### Bug Fix 3: Nodata Masking at Read Time +**Problem:** `_read_raster_clipped()` was reading raw raster values without checking `src.nodata`. Nodata pixels (value -9999) were treated as valid data, contaminating zonal statistics (e.g., pulling farm means to large negative values). + +**Impact:** Farm-level AET/PET means were incorrect for farms near raster boundaries. + +**Fix:** Now reads `src.nodata` from raster metadata and converts all matching pixels to NaN. Also applies a belt-and-suspenders check for any value ≤ -9999. + +**File:** `et_intersection.py`, lines 126–154 + +--- + +### New Feature: Temporal Gap-Fill +**What:** Added `_gap_fill_monthly_farms()` function that fills missing monthly values using neighbouring months, following Shuvam's crop-year boundary rules. + +**Why:** Some farms have NaN for individual months due to cloud cover or missing MODIS composites. Without gap-fill, these NaNs propagate to annual and kharif means, artificially reducing coverage. + +**File:** `et_intersection.py`, lines 195–248 + +--- + +### New Feature: MAI Capping to [0, 1] +**What:** MAI values > 1.0 are capped to 1.0 with a log warning. + +**Why:** AET cannot physically exceed PET, but raster misalignment or model artifacts can produce MAI > 1. Previously these were left as-is, distorting farm-level statistics. + +**File:** `et_intersection.py`, lines 326–336 + +--- + +## 10. How to Extend to New Regions/Years + +### Adding a new tehsil +1. Run Phase 1 (`fetch_raw.py`) and Phase 2 (`convert.py`) to generate `farm_boundaries.parquet` +2. Ensure the AET/PET rasters for the correct AEZ zone exist in `data/et_rasters/` +3. Check `AEZ_ZONE_MAP` in `et_intersection.py` — add the state if it's not mapped +4. Run `intersect_et_with_farms(state, district, tehsil, year)` + +### Adding a new year +1. Obtain the AET/PET rasters: `merge_AET___cog.tif` and `merge_PET___cog.tif` +2. Place them in `data/et_rasters/` +3. Run `intersect_et_with_farms(state, district, tehsil, year=)` +4. The pipeline appends to existing parquets (idempotent — re-running the same year replaces old data) + +### Adding a new AEZ zone +1. Add the state→zone mapping in `AEZ_ZONE_MAP` +2. Ensure rasters follow the naming convention `merge_AET___cog.tif` + +--- + +## 11. Known Limitations + +1. **Gap-fill granularity:** We gap-fill at the farm level (after zonal stats), not at the pixel level like Shuvam's GEE pipeline. This means if an entire farm is NaN for a month, we interpolate using neighbouring months' farm means, not neighbouring months' pixel values. + +2. **Annual MAI definition:** Our `mai_annual` is the mean of 12 monthly MAI values (Jan–Dec calendar order). Shuvam's is crop-year (Jul–Jun). For apples-to-apples comparison, use `kharif_mai` which covers the same Jul–Oct period in both. + +3. **Raster coverage at edges:** Our rasters may not cover all farms at tehsil edges. Farms outside the raster extent get NaN. Shuvam's GEE pipeline clips exactly to the tehsil boundary, so he has different NaN counts. + +4. **Large raster files:** The AET/PET rasters are 12–17 GB each. Windowed COG reading keeps memory manageable, but initial processing for a tehsil takes 30–60 seconds. diff --git a/compare_mai_tally.py b/compare_mai_tally.py new file mode 100644 index 00000000..16c02802 --- /dev/null +++ b/compare_mai_tally.py @@ -0,0 +1,159 @@ +""" +Fast MAI tally — compares our farm-level annual & kharif MAI against Shuvam's GEE raster. + +Shuvam's raster bands (crop year Jul 2018 – Jun 2019): + b1=Jul, b2=Aug, b3=Sep, b4=Oct, b5=Nov, b6=Dec, + b7=Jan, b8=Feb, b9=Mar, b10=Apr, b11=May, b12=Jun, + b13=Annual (crop-year mean) + +Kharif months Jul-Oct (bands 1-4) are the SAME calendar period in both +pipelines, so kharif comparison is the fairest apples-to-apples tally. +""" +import dotenv +dotenv.load_dotenv('.env') + +import geopandas as gpd +import pandas as pd +import numpy as np +import rasterio +import rasterio.features +import rasterio.windows +import time + +TEHSILS = [ + { + "name": "sanganer", "state": "rajasthan", + "district": "jaipur", "block": "sanganer", + "shuvam_raster": "data/et_rasters/shuvam_mai_sanganer_2018.tif", + }, + { + "name": "dudu", "state": "rajasthan", + "district": "jaipur", "block": "dudu", + "shuvam_raster": "data/et_rasters/shuvam_mai_dudu_2018.tif", + }, +] + +# Shuvam crop-year band mapping (1-indexed for rasterio) +# Jul=1, Aug=2, Sep=3, Oct=4 -> kharif bands +SHUVAM_KHARIF_BANDS = [1, 2, 3, 4] # Jul, Aug, Sep, Oct +SHUVAM_ANNUAL_BAND = 13 + + +def _zonal_mean_band(gdf_reprojected, bbox, raster_path, band_number): + """Read one band and compute per-farm zonal mean.""" + minx, miny, maxx, maxy = bbox + with rasterio.open(raster_path) as src: + window = rasterio.windows.from_bounds(minx, miny, maxx, maxy, src.transform) + transform = rasterio.windows.transform(window, src.transform) + data = src.read(band_number, window=window).astype("float32") + nodata_val = float(src.nodata) if src.nodata is not None else -9999.0 + + data[data == nodata_val] = np.nan + data[data <= -9999] = np.nan + + out_shape = data.shape + n_farms = len(gdf_reprojected) + shapes = ((geom, idx) for idx, geom in enumerate(gdf_reprojected.geometry, start=1)) + labels = rasterio.features.rasterize( + shapes, out_shape=out_shape, transform=transform, + fill=0, dtype="int32", all_touched=True, + ) + + valid = np.isfinite(data) + v_labels = labels[valid] + v_vals = data[valid] + + sums = np.bincount(v_labels, weights=v_vals, minlength=n_farms + 1) + counts = np.bincount(v_labels, minlength=n_farms + 1) + with np.errstate(invalid="ignore", divide="ignore"): + means = sums[1:] / counts[1:] + means[counts[1:] == 0] = np.nan + return means + + +def fast_zonal_means(gdf, raster_path, bands): + """ + Compute per-farm mean for each band in `bands`. + Returns a dict: band_number -> np.ndarray of length n_farms. + """ + with rasterio.open(raster_path) as src: + raster_crs = src.crs + + farms = gdf.to_crs(raster_crs) + bbox = farms.total_bounds + + results = {} + for band in bands: + results[band] = _zonal_mean_band(farms, bbox, raster_path, band) + return results + + +def print_diff(label, ours, theirs): + both = ~(np.isnan(ours) | np.isnan(theirs)) + d = np.abs(ours[both] - theirs[both]) + n = both.sum() + if n == 0: + print(f" {label}: no overlapping data") + return + print(f" {label}: n={n:,} | mean_diff={d.mean():.4f} | median={np.median(d):.4f} " + f"| max={d.max():.4f} | <0.01: {(d<0.01).sum():,} ({(d<0.01).mean()*100:.1f}%) " + f"| >0.05: {(d>0.05).sum():,} ({(d>0.05).mean()*100:.1f}%)") + + +for t in TEHSILS: + t0 = time.time() + print(f"\n{'='*65}") + print(f"TALLY — {t['name'].upper()} 2018") + print(f"{'='*65}") + + base = f"data/farm_boundaries/{t['state']}/{t['district']}/{t['block']}" + static = gpd.read_parquet(f"{base}/farm_static.parquet") + annual = pd.read_parquet(f"{base}/farm_annual.parquet") + annual = annual[annual['year'] == 2018] + monthly = pd.read_parquet(f"{base}/farm_monthly.parquet") + monthly = monthly[monthly['year'] == 2018] + + # Our kharif MAI = mean of Jul, Aug, Sep, Oct monthly MAI + monthly['month_str'] = monthly['date'].astype(str).str[5:7] + kharif_monthly = monthly[monthly['month_str'].isin(['07','08','09','10'])] + kharif_mai_ours = kharif_monthly.groupby('farm_id')['mai'].mean().reset_index() + kharif_mai_ours.columns = ['farm_id', 'kharif_mai_ours'] + + gdf = static[['farm_id','geometry']].merge( + annual[['farm_id','mai_annual']], on='farm_id', how='inner' + ).merge(kharif_mai_ours, on='farm_id', how='left') + print(f"Loaded {len(gdf):,} farm polygons") + + # Compute zonal means for annual band + kharif bands + bands_to_read = [SHUVAM_ANNUAL_BAND] + SHUVAM_KHARIF_BANDS + print(f"Computing zonal means (bands {bands_to_read})...") + band_means = fast_zonal_means(gdf, t['shuvam_raster'], bands_to_read) + + gdf['shuvam_annual'] = band_means[SHUVAM_ANNUAL_BAND] + # Shuvam kharif MAI = mean of Jul(b1)+Aug(b2)+Sep(b3)+Oct(b4) per farm + kharif_stack = np.stack([band_means[b] for b in SHUVAM_KHARIF_BANDS], axis=1) + kharif_mean = np.nanmean(kharif_stack, axis=1) + kharif_mean[np.all(np.isnan(kharif_stack), axis=1)] = np.nan + gdf['shuvam_kharif'] = kharif_mean + + print(f"Done in {time.time()-t0:.0f}s\n") + + print("--- ANNUAL MAI (our Jan-Dec vs Shuvam crop-year Jul18-Jun19) ---") + print_diff("annual", gdf['mai_annual'].values, gdf['shuvam_annual'].values) + + print("\n--- KHARIF MAI Jul-Oct 2018 (same period in both pipelines) ---") + print_diff("kharif", gdf['kharif_mai_ours'].values, gdf['shuvam_kharif'].values) + + print("\n--- Coverage ---") + print(f" NaN in our annual: {gdf['mai_annual'].isna().sum():,}") + print(f" NaN in Shuvam annual: {gdf['shuvam_annual'].isna().sum():,}") + print(f" NaN in our kharif: {gdf['kharif_mai_ours'].isna().sum():,}") + print(f" NaN in Shuvam kharif: {gdf['shuvam_kharif'].isna().sum():,}") + + # Save full tally + out = f"{base}/mai_tally_2018.parquet" + gdf[['farm_id','mai_annual','shuvam_annual','kharif_mai_ours','shuvam_kharif']].to_parquet(out, index=False) + print(f"\nFull tally saved -> {out}") + print(f"Total time: {time.time()-t0:.0f}s") + +print("\nAll done.") diff --git a/computing/api.py b/computing/api.py index 450b16ed..89ff284c 100644 --- a/computing/api.py +++ b/computing/api.py @@ -94,6 +94,48 @@ from .misc.digital_elevation_model import generate_dem_layer from .misc.canal_layer import canal_vector from .STAC_specs.stac_collection import generate_stac_collection_task +from .farm_boundaries.farm_boundary import build_farm_boundary_map + + +@api_security_check(allowed_methods="POST") +@schema(None) +def generate_farm_boundaries(request): + print("Inside generate_farm_boundaries API.") + try: + state = request.data.get("state", "").lower().strip() + district = request.data.get("district", "").lower().strip() + block = request.data.get("block", "").lower().strip() + api_key = request.data.get("api_key", "").strip() + year = request.data.get("year", None) + overwrite = request.data.get("overwrite", False) + + if not all([state, district, block, api_key]): + return Response( + {"Error": "state, district, block, and api_key are required."}, + status=status.HTTP_400_BAD_REQUEST, + ) + + if year is not None: + year = int(year) + if year < 2017 or year > 2024: + return Response( + {"Error": "year must be between 2017 and 2024."}, + status=status.HTTP_400_BAD_REQUEST, + ) + + build_farm_boundary_map.apply_async( + args=[state, district, block, api_key, year, overwrite], + queue="nrm", + ) + + msg = "Farm boundary pipeline initiated." + if year: + msg += f" ET intersection enabled for year {year}." + + return Response({"Success": msg}, status=status.HTTP_200_OK) + except Exception as e: + print("Exception in generate_farm_boundaries api :: ", e) + return Response({"Exception": str(e)}, status=status.HTTP_500_INTERNAL_SERVER_ERROR) @api_security_check(allowed_methods="POST") diff --git a/computing/farm_boundaries/__init__.py b/computing/farm_boundaries/__init__.py new file mode 100644 index 00000000..e555d280 --- /dev/null +++ b/computing/farm_boundaries/__init__.py @@ -0,0 +1 @@ +# Farm Boundaries pipeline module diff --git a/computing/farm_boundaries/convert.py b/computing/farm_boundaries/convert.py new file mode 100644 index 00000000..430fc7fc --- /dev/null +++ b/computing/farm_boundaries/convert.py @@ -0,0 +1,368 @@ +""" +Phase 2 — Convert raw per-cell JSON files into clipped GeoParquets. + +Strategy: + 1. Use DuckDB to read all raw JSON files rapidly and extract ALL landscape + features (field, trees, dug_well, farm_pond, other_water) into memory. + 2. Hand off to GeoPandas for spatial operations: WKB geometry parsing, + polygon clipping to the tehsil boundary, geometry validation. + 3. Write one GeoParquet file per structure type — allowing downstream + pipelines to consume each layer independently. + +Outputs: + data/farm_boundaries////farm_boundaries.parquet + data/farm_boundaries////trees.parquet + data/farm_boundaries////dug_wells.parquet + data/farm_boundaries////farm_ponds.parquet + data/farm_boundaries////other_water.parquet + +Usage (standalone / debug): + from computing.farm_boundaries.convert import convert_to_geoparquet + result = convert_to_geoparquet("rajasthan", "jaipur", "sanganer") + print(result) # {"farm_boundaries": {...}, "trees": {...}, ...} +""" + +import json +import logging +import os + +import geopandas as gpd +import pandas as pd +from shapely.geometry import shape +from shapely.validation import make_valid + +from utilities.constants import FARM_BOUNDARIES_PATH, SOI_TEHSIL + +logger = logging.getLogger(__name__) + +CRS = "EPSG:4326" + +# Map from alu_type value in the API response to output parquet filename +ALU_TYPE_TO_PARQUET = { + "field": "farm_boundaries.parquet", + "trees": "trees.parquet", + "dug_well": "dug_wells.parquet", + "farm_pond": "farm_ponds.parquet", + "other_water": "other_water.parquet", +} + + +# ── helpers ─────────────────────────────────────────────────────────────────── + + +def _get_tehsil_polygon(state: str, district: str, block: str): + """Re-load the tehsil polygon (shared with fetch_raw.py).""" + soi = gpd.read_file(SOI_TEHSIL) + mask = ( + (soi["STATE"].str.lower() == state) + & (soi["District"].str.lower() == district) + & (soi["TEHSIL"].str.lower() == block) + ) + subset = soi[mask] + if subset.empty: + raise ValueError( + f"Tehsil not found: state={state}, district={district}, block={block}" + ) + return subset.dissolve().geometry.iloc[0] + + +def _raw_dir(state: str, district: str, block: str) -> str: + return os.path.join(FARM_BOUNDARIES_PATH, state, district, block, "raw") + + +def _manifest_path(state: str, district: str, block: str) -> str: + return os.path.join(FARM_BOUNDARIES_PATH, state, district, block, "manifest.json") + + +def _output_dir(state: str, district: str, block: str) -> str: + return os.path.join(FARM_BOUNDARIES_PATH, state, district, block) + + +def _load_fetched_tokens(manifest_file: str) -> list: + """Return only the tokens that have actual landscape data.""" + if not os.path.exists(manifest_file): + raise FileNotFoundError( + f"Manifest not found at {manifest_file}. " + "Run Phase 1 (fetch_raw_boundaries) first." + ) + with open(manifest_file) as f: + manifest = json.load(f) + return manifest.get("fetched", []) + + +# ── DuckDB extraction ───────────────────────────────────────────────────────── + + +def _extract_all_features_with_duckdb(raw_dir: str, tokens: list) -> list: + """ + Extract ALL landscape features from raw cell JSON files. + + Uses the Python parser directly (DuckDB UNNEST-on-JSON is not supported + in the installed DuckDB version). The Python path is fast enough for + typical tehsil sizes (< 1000 cells). + + Returns a list of dicts, each with keys: + cell_token, farm_uid, alu_type, geometry_geojson, properties_json + """ + return _extract_all_features_python_fallback(raw_dir, tokens) + + +def _extract_all_features_python_fallback(raw_dir: str, tokens: list) -> list: + """ + Pure-Python fallback: reads every cell JSON file and collects ALL + landscape features. Used when DuckDB JSON parsing fails. + """ + features = [] + for token in tokens: + path = os.path.join(raw_dir, f"{token}.json") + if not os.path.exists(path): + continue + try: + with open(path) as f: + data = json.load(f) + except Exception as exc: + logger.warning("Could not parse %s: %s", path, exc) + continue + + landscape = data.get("landscape", {}) + geojson_raw = landscape.get("geojson", "") + if not geojson_raw: + continue + + try: + fc = json.loads(geojson_raw) if isinstance(geojson_raw, str) else geojson_raw + except json.JSONDecodeError as exc: + logger.warning("Invalid GeoJSON in cell %s: %s", token, exc) + continue + + for feat in fc.get("features", []): + props = feat.get("properties", {}) + alu_type = props.get("alu_type", "") + if not alu_type: + continue + features.append( + { + "cell_token": token, + "plus_code": feat.get("id", ""), + "farm_uid": feat.get("id", ""), + "alu_type": alu_type, + "geometry_geojson": json.dumps(feat.get("geometry", {})), + "properties_json": json.dumps(props), + } + ) + return features + + +# ── GeoPandas spatial processing ────────────────────────────────────────────── + + +def _build_geodataframe(records: list) -> gpd.GeoDataFrame: + """ + Convert the flat list of feature dicts into a GeoDataFrame. + Parses the embedded GeoJSON geometry string into Shapely geometries. + """ + if not records: + return gpd.GeoDataFrame( + columns=["farm_uid", "cell_token", "alu_type", "geometry"], + geometry="geometry", + crs=CRS, + ) + + rows = [] + for rec in records: + try: + geom_dict = ( + json.loads(rec["geometry_geojson"]) + if isinstance(rec["geometry_geojson"], str) + else rec["geometry_geojson"] + ) + geom = shape(geom_dict) + if not geom.is_valid: + geom = make_valid(geom) + except Exception as exc: + logger.debug("Skipping invalid geometry: %s", exc) + continue + + # Parse properties blob for any extra attributes we want to keep. + try: + props = ( + json.loads(rec["properties_json"]) + if rec.get("properties_json") + else {} + ) + except Exception: + props = {} + + rows.append( + { + "farm_uid": rec.get("farm_uid", "") or rec.get("plus_code", ""), + "cell_token": rec.get("cell_token", ""), + "alu_type": rec.get("alu_type") or props.get("alu_type", "field"), + "plus_code": rec.get("plus_code", ""), # from feature-level id + "area_m2": props.get("area_sq_m", None), + "class_confidence": props.get("class_confidence", None), + "capture_date": props.get("capture_timestamp_sec", None), + "geometry": geom, + } + ) + + if not rows: + return gpd.GeoDataFrame(geometry=[], crs=CRS) + + gdf = gpd.GeoDataFrame(rows, geometry="geometry", crs=CRS) + return gdf + + +def _assign_farm_ids( + gdf: gpd.GeoDataFrame, state: str, district: str, block: str +) -> gpd.GeoDataFrame: + """ + Assign a unique, human-readable farm_id to every row. + Format: ___ + """ + prefix = f"{state}_{district}_{block}" + gdf = gdf.reset_index(drop=True) + gdf["farm_id"] = [f"{prefix}_{i:06d}" for i in gdf.index] + return gdf + + +# ── public entry point ─────────────────────────────────────────────────────── + + +def convert_to_geoparquet( + state: str, + district: str, + block: str, + overwrite: bool = False, +) -> dict: + """ + Phase 2 pipeline: read raw cell JSON files, extract ALL structure types, + clip them to the tehsil boundary, and write one GeoParquet per type. + + Structure types: + field -> farm_boundaries.parquet + trees -> trees.parquet + dug_well -> dug_wells.parquet + farm_pond -> farm_ponds.parquet + other_water -> other_water.parquet + + Parameters + ---------- + state, district, block : str + Lower-cased administrative names. + overwrite : bool + If False (default), skip any structure whose parquet already exists. + If True, regenerate all parquets. + + Returns + ------- + dict + Top-level keys: each alu_type name -> {path, count, skipped}. + Also includes 'path' pointing to farm_boundaries.parquet for + backward compatibility with the Celery task. + """ + out_dir = _output_dir(state, district, block) + os.makedirs(out_dir, exist_ok=True) + + logger.info( + "Phase 2 — converting raw JSON to GeoParquets for %s/%s/%s", + state, district, block, + ) + + # 1. Load manifest -------------------------------------------------------- + manifest_file = _manifest_path(state, district, block) + fetched_tokens = _load_fetched_tokens(manifest_file) + logger.info("%d cells with landscape data to process.", len(fetched_tokens)) + + if not fetched_tokens: + logger.warning("No data cells found in manifest. Returning empty result.") + return {"path": None, "farm_count": 0, "state": state, "district": district, "block": block} + + # 2. Check which structure types still need processing -------------------- + structures_to_process = {} + all_results = {} + for alu_type, parquet_name in ALU_TYPE_TO_PARQUET.items(): + out_path = os.path.join(out_dir, parquet_name) + if not overwrite and os.path.exists(out_path): + logger.info(" [%s] Already exists at %s — skipping.", alu_type, out_path) + all_results[alu_type] = {"path": out_path, "skipped": True} + else: + structures_to_process[alu_type] = out_path + + if not structures_to_process: + logger.info("All structure parquets already exist. Use overwrite=True to regenerate.") + # Backward compatibility: return path to farm_boundaries.parquet + farm_path = os.path.join(out_dir, ALU_TYPE_TO_PARQUET["field"]) + return {"path": farm_path, "skipped": True, "all_structures": all_results} + + # 3. Extract ALL features from raw JSON ----------------------------------- + raw_dir = _raw_dir(state, district, block) + logger.info("Extracting all features from %d cells...", len(fetched_tokens)) + all_records = _extract_all_features_with_duckdb(raw_dir, fetched_tokens) + logger.info("Total features extracted across all types: %d", len(all_records)) + + # 4. Load tehsil boundary once (shared for all structure clips) ----------- + tehsil_geom = _get_tehsil_polygon(state, district, block) + tehsil_gdf = gpd.GeoDataFrame(geometry=[tehsil_geom], crs=CRS) + + # 5. Build GeoDataFrame for ALL records ----------------------------------- + full_gdf = _build_geodataframe(all_records) + logger.info("Built GeoDataFrame: %d total features.", len(full_gdf)) + + # 6. Per-structure-type: filter, clip, assign IDs, save ------------------- + for alu_type, out_path in structures_to_process.items(): + logger.info(" [%s] Processing...", alu_type) + + # Filter to this structure type + subset = full_gdf[full_gdf["alu_type"] == alu_type].copy() + if subset.empty: + logger.info(" [%s] No features found — writing empty parquet.", alu_type) + subset.to_parquet(out_path, index=False) + all_results[alu_type] = {"path": out_path, "count": 0, "skipped": False} + continue + + logger.info(" [%s] %d raw features before clipping.", alu_type, len(subset)) + + # Clip to tehsil boundary + subset = gpd.clip(subset, tehsil_gdf) + logger.info(" [%s] %d features after clipping.", alu_type, len(subset)) + + # Assign unique IDs (prefix differs per type) + prefix_map = { + "field": f"{state}_{district}_{block}", + "trees": f"{state}_{district}_{block}_tree", + "dug_well": f"{state}_{district}_{block}_well", + "farm_pond": f"{state}_{district}_{block}_pond", + "other_water": f"{state}_{district}_{block}_water", + } + prefix = prefix_map.get(alu_type, f"{state}_{district}_{block}_{alu_type}") + subset = subset.reset_index(drop=True) + subset["feature_id"] = [f"{prefix}_{i:06d}" for i in subset.index] + # Keep farm_id alias for fields (backward compatibility) + if alu_type == "field": + subset["farm_id"] = subset["feature_id"] + + # Reorder columns + priority_cols = ["feature_id", "farm_id", "farm_uid", "cell_token", + "alu_type", "plus_code", "area_m2", "class_confidence", + "capture_date", "geometry"] + existing = [c for c in priority_cols if c in subset.columns] + subset = subset[existing] + + # Save + subset.to_parquet(out_path, index=False) + logger.info(" [%s] Saved %d features -> %s", alu_type, len(subset), out_path) + all_results[alu_type] = {"path": out_path, "count": len(subset), "skipped": False} + + # Backward compatibility: top-level 'path' points to farm_boundaries.parquet + farm_info = all_results.get("field", {}) + summary = { + "state": state, + "district": district, + "block": block, + "path": farm_info.get("path"), + "farm_count": farm_info.get("count", 0), + "all_structures": all_results, + } + logger.info("Phase 2 complete: %s", summary) + return summary diff --git a/computing/farm_boundaries/et_intersection.py b/computing/farm_boundaries/et_intersection.py new file mode 100644 index 00000000..1d40c739 --- /dev/null +++ b/computing/farm_boundaries/et_intersection.py @@ -0,0 +1,935 @@ +""" +Phase 3 — Intersect AET & PET rasters with farm boundary polygons. + +Reads locally stored COG (Cloud Optimized GeoTIFF) rasters for AET and PET, +runs zonal statistics against each farm polygon, computes MAI (Moisture +Adequacy Index = AET/PET), and produces three parquets per the core-lens schema: + + farm_static.parquet — one row per farm (geometry + static properties) + farm_annual.parquet — one row per farm per year (annual ET metrics) + farm_monthly.parquet — one row per farm per month (date, AET, PET, MAI) + +Data sources: + Local COG rasters at LOCAL_ET_RASTERS_PATH: + merge_AET___cog.tif (13 bands: b1-b12 monthly mm/day, b13 annual) + merge_PET___cog.tif (same structure) + Resolution: 30 metres | NoData: -9999 | CRS: EPSG:4326 + +Water stress methodology (aligned with Shuvam Chakraborty / ET Applications): + MAI = AET / PET (ratio, per pixel, only where both AET & PET are valid and PET > 0) + Moderate kharif stress : mean kharif MAI <= 0.50 + Severe kharif stress : mean kharif MAI <= 0.25 + Kharif months : July, August, September, October + +Missing data protocol (mirrors Shuvam's divide_where_valid approach): + - Pixel-level : MAI = NaN if AET is NaN, PET is NaN, or PET = 0 + - Farm-level : column = NaN if the farm has zero valid pixels for that band + - Annual MAI : mean of all valid monthly MAI values (NaN months excluded) + - No imputation is performed on missing farms or missing months. +""" + +import logging +import os +import warnings +from datetime import date + +import geopandas as gpd +import numpy as np +import pandas as pd +import rasterio +import rasterio.features +import rasterio.merge +import rasterio.windows +from shapely.geometry import box + +from utilities.constants import ( + AEZ_GEOJSON, + FARM_BOUNDARIES_PATH, + LOCAL_ET_RASTERS_PATH, + SOI_TEHSIL, +) + +logger = logging.getLogger(__name__) + +AET_NODATA = -9999 + +# MAI thresholds — aligned with Shuvam Chakraborty / ET Applications (GEE pipeline) +# Moderate stress : MAI <= 0.50 (farm is water-stressed but not severely) +# Severe stress : MAI <= 0.25 (farm is severely water-stressed) +MAI_MODERATE_THRESHOLD = 0.50 +MAI_SEVERE_THRESHOLD = 0.25 + +# Kept for backward compatibility +KHARIF_WATER_STRESS_MAI_THRESHOLD = MAI_MODERATE_THRESHOLD + +KHARIF_MONTH_NAMES = ["jul", "aug", "sep", "oct"] + +# Calendar month names (used for column naming) +MONTH_NAMES = [ + "jan", "feb", "mar", "apr", "may", "jun", + "jul", "aug", "sep", "oct", "nov", "dec", +] + +# Rasters use crop-year band ordering: band 1 = July, band 2 = August, +# ..., band 6 = December (year Y), band 7 = January, ..., band 12 = June (year Y+1). +# This list maps band index 0..11 to the correct calendar month number 1..12. +CROP_YEAR_BAND_TO_MONTH = [7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6] + +# Reverse: calendar month number -> column index in MONTH_NAMES +_MONTH_NUM_TO_NAME = { + 1: "jan", 2: "feb", 3: "mar", 4: "apr", 5: "may", 6: "jun", + 7: "jul", 8: "aug", 9: "sep", 10: "oct", 11: "nov", 12: "dec", +} + +CRS = "EPSG:4326" + + +# ── path helpers ─────────────────────────────────────────────────────────────── + +def _block_dir(state, district, block): + return os.path.join(FARM_BOUNDARIES_PATH, state, district, block) + +def _farm_parquet_path(state, district, block): + return os.path.join(_block_dir(state, district, block), "farm_boundaries.parquet") + +def _static_parquet_path(state, district, block): + return os.path.join(_block_dir(state, district, block), "farm_static.parquet") + +def _annual_parquet_path(state, district, block): + return os.path.join(_block_dir(state, district, block), "farm_annual_vectorize.parquet") + +def _monthly_parquet_path(state, district, block): + return os.path.join(_block_dir(state, district, block), "farm_monthly_vectorize.parquet") + +def _local_aet_path(aez, year): + return os.path.join(LOCAL_ET_RASTERS_PATH, f"merge_AET_{aez}_{year}_cog.tif") + +def _local_pet_path(aez, year): + return os.path.join(LOCAL_ET_RASTERS_PATH, f"merge_PET_{aez}_{year}_cog.tif") + +def _get_tehsil_polygon(state, district, block): + """ + Load the tehsil polygon from the shared SOI tehsil boundaries GeoJSON + (SOI_TEHSIL), matched case-insensitively on state/district/tehsil name. + Dissolves to a single geometry in case the tehsil has multiple rows. + """ + soi = gpd.read_file(SOI_TEHSIL) + mask = ( + (soi["STATE"].str.lower() == state) + & (soi["District"].str.lower() == district) + & (soi["TEHSIL"].str.lower() == block) + ) + subset = soi[mask] + if subset.empty: + raise ValueError( + f"Tehsil not found in {SOI_TEHSIL}: state={state}, district={district}, block={block}" + ) + return subset.dissolve().geometry.iloc[0] + + +AEZ_MIN_OVERLAP_FRAC = 0.001 # 0.1% of tehsil area — filters boundary-snapping slivers + + +def _get_aez_zones(state, district, block, min_overlap_frac=AEZ_MIN_OVERLAP_FRAC): + """ + Determine ALL AEZ zones a tehsil's boundary genuinely overlaps (not just + the single largest one) by intersecting it against the pan-India AEZ + polygons (AEZ_GEOJSON) — AEZ boundaries don't follow tehsil/state lines, + so a tehsil near a zone boundary can straddle more than one zone, and + farms sitting in the minority zone need that zone's own raster, not + whichever zone covers the most of the tehsil. + + `min_overlap_frac` filters out negligible sliver overlaps caused by + boundary-snapping/precision mismatches between the independently + digitized SOI tehsil and AEZ datasets — not genuine multi-zone tehsils. + + Returns + ------- + list[int] AEZ zone codes (ae_regcode), ordered by overlap area + descending (largest first). + """ + tehsil_geom = _get_tehsil_polygon(state, district, block) + tehsil_area = tehsil_geom.area + if tehsil_area <= 0: + raise ValueError(f"Tehsil {state}/{district}/{block} has zero-area geometry.") + + aez_gdf = gpd.read_file(AEZ_GEOJSON) + if aez_gdf.crs is not None and str(aez_gdf.crs) != CRS: + aez_gdf = aez_gdf.to_crs(CRS) + + # Only used to compare overlap areas *within* one tehsil, not as an + # absolute measurement — geographic-CRS area distortion is negligible + # at this scale, so the degree-based warning geopandas raises is safe + # to suppress here. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + overlap_areas = aez_gdf.geometry.intersection(tehsil_geom).area + + significant = overlap_areas[(overlap_areas / tehsil_area) >= min_overlap_frac] + if significant.empty: + raise ValueError( + f"Tehsil {state}/{district}/{block} does not intersect any AEZ zone in {AEZ_GEOJSON}." + ) + + ordered = significant.sort_values(ascending=False) + zones = [int(aez_gdf.loc[i, "ae_regcode"]) for i in ordered.index] + + if len(zones) > 1: + pct = [round(100 * overlap_areas.loc[i] / tehsil_area, 1) for i in ordered.index] + logger.info( + "Tehsil %s/%s/%s straddles %d AEZ zones: %s (overlap%% of tehsil area: %s)", + state, district, block, len(zones), zones, pct, + ) + else: + logger.info("AEZ zone for %s/%s/%s: %d", state, district, block, zones[0]) + + return zones + + +# ── local raster reading ─────────────────────────────────────────────────────── + +def _read_raster_clipped(raster_paths, bbox): + """ + Read one or more local COG rasters — one per AEZ zone the tehsil + straddles — and mosaic them into a single array windowed to `bbox` via + rasterio.merge. When a tehsil sits entirely inside one zone, this is + just `raster_paths` of length 1 and behaves like a plain windowed read. + When it straddles multiple zones, farms are matched against whichever + zone's raster actually covers their location, without needing to split + farms into per-zone groups beforehand. + + Converts the rasters' nodata value (-9999 per spec, confirmed via + src.nodata) to NaN immediately so all downstream code works cleanly + with NaN semantics. + + Parameters + ---------- + raster_paths : str | list[str] + One or more raster file paths (all same resolution/CRS/band layout). + bbox : tuple + (minx, miny, maxx, maxy) in EPSG:4326. + + Returns + ------- + data : np.ndarray shape (bands, height, width), float32 + transform : affine transform for the merged/clipped window + """ + if isinstance(raster_paths, str): + raster_paths = [raster_paths] + + existing_paths = [p for p in raster_paths if os.path.exists(p)] + if not existing_paths: + raise FileNotFoundError(f"None of the expected rasters exist: {raster_paths}") + missing_paths = set(raster_paths) - set(existing_paths) + if missing_paths: + logger.warning( + "%d/%d AEZ-zone raster(s) missing, proceeding with the rest: %s", + len(missing_paths), len(raster_paths), sorted(missing_paths), + ) + + minx, miny, maxx, maxy = bbox + srcs = [rasterio.open(p) for p in existing_paths] + try: + nodata_val = next( + (float(s.nodata) for s in srcs if s.nodata is not None), float(AET_NODATA) + ) + merged, transform = rasterio.merge.merge( + srcs, bounds=(minx, miny, maxx, maxy), nodata=nodata_val, + ) + finally: + for s in srcs: + s.close() + + data = merged.astype("float32") + n_nodata = int(np.sum(data == nodata_val)) + if n_nodata > 0: + logger.debug("Merged raster (%d source file(s)): masking %d nodata pixels (value=%.0f)", + len(existing_paths), n_nodata, nodata_val) + + # Convert nodata sentinel AND any residual AET_NODATA values to NaN + data[data == nodata_val] = np.nan + data[data <= float(AET_NODATA)] = np.nan # belt-and-suspenders for -9999 variants + + return data, transform + + +# ── zonal statistics (true polygon-raster geometric intersection) ────────────── +# No rasterize/paint step anywhere below: each farm's own polygon is intersected +# directly against the raster's pixel grid using shapely, and each pixel's +# contribution is weighted by its exact overlap area with that farm. This is +# slower than a single vectorised rasterize+bincount pass over all farms at +# once, but avoids both of that approach's inaccuracies: (a) two farms sharing +# a boundary pixel can no longer "steal" it from each other — each farm is +# evaluated independently against its own geometry — and (b) a pixel that's +# only partially inside a farm contributes proportionally, not all-or-nothing. + +def _farm_pixel_window(geom, transform, arr_height, arr_width): + """ + Compute the row/col index range of the raster array that could possibly + contain pixels overlapping `geom`'s bounding box. This is only a cheap + pre-filter (plain arithmetic on the affine transform, no rasterization) + so we don't test every pixel in the array against every farm — the real + geometric intersection test happens per-candidate-pixel afterwards. + + Returns + ------- + row_start, row_stop, col_start, col_stop : int + Half-open index ranges into the (height, width) array, clipped to + the array's own bounds, with a 1-pixel buffer on each side. + """ + minx, miny, maxx, maxy = geom.bounds + inv = ~transform # map (x, y) -> fractional (col, row) pixel coords + + col_a, row_a = inv * (minx, maxy) + col_b, row_b = inv * (maxx, miny) + + row_start = max(int(np.floor(min(row_a, row_b))) - 1, 0) + row_stop = min(int(np.ceil(max(row_a, row_b))) + 1, arr_height) + col_start = max(int(np.floor(min(col_a, col_b))) - 1, 0) + col_stop = min(int(np.ceil(max(col_a, col_b))) + 1, arr_width) + + return row_start, row_stop, col_start, col_stop + + +def _farm_pixel_weights(geom, transform, row_start, row_stop, col_start, col_stop): + """ + For every candidate pixel in the window, build its exact map-space + rectangle from the raster's affine transform and compute how much of + that rectangle's area truly overlaps `geom` via a direct shapely + intersection — an exact vector-on-vector overlap, not an approximation. + + Returns + ------- + list of (row, col, weight) tuples, one per pixel that genuinely overlaps + the farm polygon (weight > 0). `weight` is the overlap area, in the + raster's native coordinate units (deg² for EPSG:4326). + """ + weighted_pixels = [] + for row in range(row_start, row_stop): + for col in range(col_start, col_stop): + x0, y0 = transform * (col, row) + x1, y1 = transform * (col + 1, row + 1) + pixel_box = box(min(x0, x1), min(y0, y1), max(x0, x1), max(y0, y1)) + + if not geom.intersects(pixel_box): + continue + + overlap = geom.intersection(pixel_box) + weight = overlap.area + if weight > 0: + weighted_pixels.append((row, col, weight)) + + return weighted_pixels + + +def _weighted_band_means(band_stack, weighted_pixels, num_bands): + """ + Given a raster band stack (bands, height, width) and the list of + (row, col, weight) overlap pixels for one farm, compute the + area-weighted mean value per band for that farm. + + A pixel is excluded from a given band's average only if that band's + value at that pixel is itself NaN/invalid (nodata) — validity is + checked per band, independently, since nodata coverage can differ + month to month. + + Returns + ------- + np.ndarray of shape (num_bands,) — area-weighted mean per band. + NaN where the farm has zero valid weight for that band. + """ + if not weighted_pixels: + return np.full(num_bands, np.nan) + + rows = np.array([p[0] for p in weighted_pixels]) + cols = np.array([p[1] for p in weighted_pixels]) + weights = np.array([p[2] for p in weighted_pixels]) + + values = band_stack[:, rows, cols] # shape (num_bands, n_pixels) + valid = np.isfinite(values) & (values >= 0) # per-band validity mask + + weighted_vals = np.where(valid, values * weights, 0.0) + weighted_wts = np.where(valid, weights, 0.0) + + sum_vals = weighted_vals.sum(axis=1) + sum_wts = weighted_wts.sum(axis=1) + + with np.errstate(invalid="ignore", divide="ignore"): + means = sum_vals / sum_wts + means[sum_wts == 0] = np.nan # farm has no valid weighted pixels → NaN + return means + + +def _extract_farm_zonal_means(band_stack, geom, transform): + """ + Compute area-weighted per-band means for one farm polygon against a + raster band stack, via true polygon-pixel geometric intersection. + + Combines the three steps above: find the candidate pixel window, compute + exact overlap weights against `geom`, then take the area-weighted mean + per band. + """ + num_bands, arr_height, arr_width = band_stack.shape + + if geom is None or geom.is_empty: + return np.full(num_bands, np.nan) + + row_start, row_stop, col_start, col_stop = _farm_pixel_window( + geom, transform, arr_height, arr_width + ) + if row_start >= row_stop or col_start >= col_stop: + return np.full(num_bands, np.nan) + + weighted_pixels = _farm_pixel_weights( + geom, transform, row_start, row_stop, col_start, col_stop + ) + return _weighted_band_means(band_stack, weighted_pixels, num_bands) + + +# ── temporal gap-filling ─────────────────────────────────────────────────────── +# Mirrors Shuvam Chakraborty's fill_monthly_collection() in ET_Applications/helper.py. +# Crop year: July (agri-month 1) → June (agri-month 12). +# Rules: +# July (agri_month 1) → neighbour: August only (no backward crossing crop-year start) +# June (agri_month 12) → neighbour: May only (no forward crossing crop-year end) +# All others → previous and next calendar month (±1 month) +# NaN farm-months with no valid neighbour remain NaN. + +# Calendar-month index (0=Jan … 11=Dec) → list of neighbour indices +_GAP_FILL_NEIGHBOURS: dict = { + 0: [11, 1], # Jan: Dec, Feb + 1: [0, 2], # Feb: Jan, Mar + 2: [1, 3], # Mar: Feb, Apr + 3: [2, 4], # Apr: Mar, May + 4: [3, 5], # May: Apr, Jun + 5: [4], # Jun: May only (crop-year end — no forward crossing) + 6: [7], # Jul: Aug only (crop-year start — no backward crossing) + 7: [6, 8], # Aug: Jul, Sep + 8: [7, 9], # Sep: Aug, Oct + 9: [8, 10], # Oct: Sep, Nov + 10: [9, 11], # Nov: Oct, Dec + 11: [10, 0], # Dec: Nov, Jan +} + + +def _gap_fill_monthly_farms(monthly_matrix: np.ndarray) -> np.ndarray: + """ + Gap-fill a (n_farms × 12) monthly matrix following Shuvam's crop-year rules. + + Parameters + ---------- + monthly_matrix : np.ndarray, shape (n_farms, 12) + Columns are calendar months Jan–Dec (indices 0–11). + NaN = missing / nodata. + + Returns + ------- + filled : np.ndarray, same shape. + NaN cells replaced with the nanmean of valid neighbours. + Cells with no valid neighbour remain NaN. + """ + filled = monthly_matrix.copy() + + for m, neighbours in _GAP_FILL_NEIGHBOURS.items(): + missing = np.isnan(filled[:, m]) + if not missing.any(): + continue + neighbour_vals = np.stack([filled[:, n] for n in neighbours], axis=1) + fill_vals = np.nanmean(neighbour_vals, axis=1) + can_fill = missing & np.isfinite(fill_vals) + filled[can_fill, m] = fill_vals[can_fill] + + return filled + + +def _extract_all_farms(gdf, band_data, transform, label): + """ + Loop over every farm polygon and compute its area-weighted per-band + means via true geometric intersection (_extract_farm_zonal_means). + + No rasterize/label-grid step: each farm is evaluated independently + against its own geometry, so shared boundary pixels are never + exclusively "won" by one neighbour over another. + + Returns + ------- + np.ndarray of shape (num_farms, num_bands). + """ + num_farms = len(gdf) + num_bands = band_data.shape[0] + means = np.full((num_farms, num_bands), np.nan) + log_every = 5000 + + for i, geom in enumerate(gdf.geometry): + means[i] = _extract_farm_zonal_means(band_data, geom, transform) + if (i + 1) % log_every == 0: + logger.info(" %s: %d/%d farms processed", label, i + 1, num_farms) + + n_with_data = int(np.any(np.isfinite(means), axis=1).sum()) + logger.info("%s: %d/%d farms have at least one valid pixel.", label, n_with_data, num_farms) + return means + + +def _run_zonal_stats(gdf, aet_data, aet_transform, pet_data=None, pet_transform=None): + """ + Compute per-farm monthly AET, PET, MAI from pre-loaded raster arrays. + Adds wide-format columns (aet_jan..aet_dec, pet_jan..pet_dec, mai_jan..mai_dec, + aet_annual, pet_annual, mai_annual, kharif_mai, kharif_water_stress) to gdf. + + Per-farm values come from true polygon-raster geometric intersection + (see _extract_farm_zonal_means) rather than a shared rasterized grid. + """ + num_farms = len(gdf) + + logger.info("Extracting AET for %d farms via geometric intersection...", num_farms) + aet_means = _extract_all_farms(gdf, aet_data, aet_transform, label="AET") + + # AET monthly — map each crop-year band to its correct calendar month column + # Band 0=Jul, 1=Aug, ..., 5=Dec, 6=Jan, ..., 11=Jun (CROP_YEAR_BAND_TO_MONTH) + aet_monthly_cols = [f"aet_{m}" for m in MONTH_NAMES] # ordered jan..dec + num_aet_bands = aet_data.shape[0] + for band_idx in range(min(num_aet_bands, 12)): + cal_month = CROP_YEAR_BAND_TO_MONTH[band_idx] + col = f"aet_{_MONTH_NUM_TO_NAME[cal_month]}" + gdf[col] = np.round(aet_means[:, band_idx], 4) + + # ── Temporal gap-fill AET — matrix in calendar order (jan=col0..dec=col11) ── + aet_matrix = gdf[aet_monthly_cols].values.astype("float64") + aet_filled = _gap_fill_monthly_farms(aet_matrix) + n_filled_aet = int(np.sum(np.isnan(aet_matrix) & np.isfinite(aet_filled))) + logger.info("Gap-fill AET: filled %d farm-month NaN values.", n_filled_aet) + for i, col in enumerate(aet_monthly_cols): + gdf[col] = np.round(aet_filled[:, i], 4) + + if num_aet_bands >= 13: + gdf["aet_annual"] = np.round(aet_means[:, 12], 4) + else: + gdf["aet_annual"] = gdf[aet_monthly_cols].mean(axis=1).round(4) + + # PET monthly + pet_monthly_cols = [f"pet_{m}" for m in MONTH_NAMES] + if pet_data is not None: + logger.info("Extracting PET for %d farms via geometric intersection...", num_farms) + pet_means = _extract_all_farms(gdf, pet_data, pet_transform, label="PET") + + num_pet_bands = pet_data.shape[0] + for band_idx in range(min(num_pet_bands, 12)): + cal_month = CROP_YEAR_BAND_TO_MONTH[band_idx] + col = f"pet_{_MONTH_NUM_TO_NAME[cal_month]}" + gdf[col] = np.round(pet_means[:, band_idx], 4) + + # ── Temporal gap-fill PET ───────────────────────────────────────────── + pet_matrix = gdf[pet_monthly_cols].values.astype("float64") + pet_filled = _gap_fill_monthly_farms(pet_matrix) + n_filled_pet = int(np.sum(np.isnan(pet_matrix) & np.isfinite(pet_filled))) + logger.info("Gap-fill PET: filled %d farm-month NaN values.", n_filled_pet) + for i, col in enumerate(pet_monthly_cols): + gdf[col] = np.round(pet_filled[:, i], 4) + + if num_pet_bands >= 13: + gdf["pet_annual"] = np.round(pet_means[:, 12], 4) + else: + gdf["pet_annual"] = gdf[pet_monthly_cols].mean(axis=1).round(4) + + # MAI + water stress + if len(pet_monthly_cols) == 12: + mai_monthly_cols = [] + for month in MONTH_NAMES: + col = f"mai_{month}" + aet_vals = gdf[f"aet_{month}"].values + pet_vals = gdf[f"pet_{month}"].values + with np.errstate(invalid="ignore", divide="ignore"): + v = aet_vals / pet_vals + # MAI is physically bounded to [0, 1]: AET cannot exceed PET. + # Values > 1 indicate raster misalignment or model artifacts → cap at 1. + n_invalid = int(np.sum(np.isfinite(v) & (v > 1))) + if n_invalid > 0: + logger.warning( + "MAI[%s]: %d farms have MAI > 1 (raster artifact) — capped at 1.0", + month, n_invalid, + ) + # Set to NaN where not finite, cap valid values to [0, 1] + v = np.where(np.isfinite(v), np.clip(v, 0.0, 1.0), np.nan) + gdf[col] = np.round(v, 4) + mai_monthly_cols.append(col) + + gdf["mai_annual"] = gdf[mai_monthly_cols].mean(axis=1).round(4) + + kharif_cols = [f"mai_{m}" for m in KHARIF_MONTH_NAMES] + kharif_df = gdf[kharif_cols] + gdf["kharif_mai"] = kharif_df.mean(axis=1).round(4) + + # Moderate stress: any kharif month with MAI <= 0.50 + gdf["kharif_water_stress"] = ( + (kharif_df <= MAI_MODERATE_THRESHOLD) & kharif_df.notna() + ).any(axis=1) + + # Severe stress: any kharif month with MAI <= 0.25 + gdf["kharif_severe_stress"] = ( + (kharif_df <= MAI_SEVERE_THRESHOLD) & kharif_df.notna() + ).any(axis=1) + + n_nan = int(gdf["mai_annual"].isna().sum()) + n_valid = int(gdf["mai_annual"].notna().sum()) + logger.info( + "MAI complete: avg_annual=%.4f | kharif_stress=%d | severe=%d | nan_farms=%d | valid_farms=%d", + gdf["mai_annual"].mean() if n_valid > 0 else float("nan"), + int(gdf["kharif_water_stress"].sum()), + int(gdf["kharif_severe_stress"].sum()), + n_nan, n_valid, + ) + else: + logger.warning("PET not available — MAI not computed.") + + return gdf + + +# ── output writers ───────────────────────────────────────────────────────────── + +def _save_static_parquet(gdf, state, district, block): + """ + Write farm_static.parquet — one row per farm, geometry + static properties. + Skips if file already exists (static data never changes). + """ + out_path = _static_parquet_path(state, district, block) + if os.path.exists(out_path): + logger.info("farm_static.parquet already exists — skipping.") + return out_path + + keep = ["farm_id", "farm_uid", "cell_token", "alu_type", + "class_confidence", "capture_date", "geometry"] + static = gdf[[c for c in keep if c in gdf.columns]].copy() + static.insert(0, "state", state) + static.insert(0, "district", district) + static.insert(0, "tehsil", block) + if "area_m2" in gdf.columns: + static["area_in_ha"] = (gdf["area_m2"] / 10_000).round(4) + + # Bounding box struct + static["bbox"] = static.geometry.apply( + lambda g: { + "xmin": round(g.bounds[0], 6), "ymin": round(g.bounds[1], 6), + "xmax": round(g.bounds[2], 6), "ymax": round(g.bounds[3], 6), + } + ) + + os.makedirs(os.path.dirname(out_path), exist_ok=True) + static.to_parquet(out_path, index=False) + logger.info("farm_static.parquet saved → %s (%d farms)", out_path, len(static)) + return out_path + + +def _save_annual_parquet(gdf, state, district, block, year): + """ + Append one year of annual ET metrics to farm_annual.parquet. + Replaces any existing rows for the same year (idempotent). + """ + out_path = _annual_parquet_path(state, district, block) + + keep = ["farm_id", "aet_annual", "pet_annual", + "mai_annual", "kharif_mai", "kharif_water_stress", "kharif_severe_stress"] + annual = gdf[[c for c in keep if c in gdf.columns]].copy() + annual["tehsil"] = block + annual["district"] = district + annual["state"] = state + annual["year"] = int(year) + if "area_m2" in gdf.columns: + annual["area_in_ha"] = (gdf["area_m2"] / 10_000).round(4) + + col_order = ["farm_id", "tehsil", "district", "state", "area_in_ha", "year", + "aet_annual", "pet_annual", "mai_annual", "kharif_mai", + "kharif_water_stress", "kharif_severe_stress"] + annual = annual[[c for c in col_order if c in annual.columns]] + + if os.path.exists(out_path): + existing = pd.read_parquet(out_path) + existing = existing[existing["year"] != int(year)] + combined = pd.concat([existing, annual], ignore_index=True) + else: + combined = annual + + os.makedirs(os.path.dirname(out_path), exist_ok=True) + combined.to_parquet(out_path, index=False) + logger.info( + "farm_annual.parquet updated → %s (%d total rows)", out_path, len(combined) + ) + return out_path + + +def _save_monthly_parquet(gdf, state, district, block, year): + """ + Melt monthly AET/PET/MAI wide columns into long format and + append to farm_monthly.parquet. One row per farm per month. + """ + out_path = _monthly_parquet_path(state, district, block) + + farm_ids = gdf["farm_id"].values if "farm_id" in gdf.columns else np.arange(len(gdf)) + area_vals = (gdf["area_m2"] / 10_000).round(4).values if "area_m2" in gdf.columns else np.full(len(gdf), np.nan) + + rows = [] + for month in MONTH_NAMES: # month = 'jan','feb',...'dec' (calendar order) + month_num = list(_MONTH_NUM_TO_NAME.keys())[ + list(_MONTH_NUM_TO_NAME.values()).index(month) + ] # calendar month number 1-12 + # Months Jul-Dec belong to `year`; Jan-Jun belong to the next calendar year + # (crop year starting July spans two calendar years) + cal_year = int(year) if month_num >= 7 else int(year) + 1 + rows.append(pd.DataFrame({ + "farm_id": farm_ids, + "tehsil": block, + "district": district, + "state": state, + "area_in_ha": area_vals, + "year": int(year), + "date": date(cal_year, month_num, 1), + "aet": gdf[f"aet_{month}"].values if f"aet_{month}" in gdf.columns else np.nan, + "pet": gdf[f"pet_{month}"].values if f"pet_{month}" in gdf.columns else np.nan, + "mai": gdf[f"mai_{month}"].values if f"mai_{month}" in gdf.columns else np.nan, + })) + + monthly = pd.concat(rows, ignore_index=True) + monthly["date"] = pd.to_datetime(monthly["date"]) + + if os.path.exists(out_path): + existing = pd.read_parquet(out_path) + existing = existing[existing["year"] != int(year)] + combined = pd.concat([existing, monthly], ignore_index=True) + else: + combined = monthly + + os.makedirs(os.path.dirname(out_path), exist_ok=True) + combined.to_parquet(out_path, index=False) + logger.info( + "farm_monthly.parquet updated → %s (%d total rows)", out_path, len(combined) + ) + return out_path + + +# ── main entry point ─────────────────────────────────────────────────────────── + +def intersect_et_with_farms( + state: str, + district: str, + block: str, + year: int = 2018, + overwrite: bool = False, +) -> dict: + """ + Phase 3: intersect local AET/PET COG rasters with farm polygons. + + Reads local rasters from LOCAL_ET_RASTERS_PATH, runs per-farm geometric + zonal statistics, computes MAI, and writes/updates the three core-lens parquets: + farm_static.parquet, farm_annual.parquet, farm_monthly.parquet + + Parameters + ---------- + state, district, block : str + Lower-cased administrative names. + year : int + Year of ET data to process (e.g. 2018). + overwrite : bool + Re-process even if this year's data already exists. + + Returns + ------- + dict summary with paths and key statistics. + """ + # Skip if year already in annual parquet + annual_path = _annual_parquet_path(state, district, block) + if not overwrite and os.path.exists(annual_path): + existing = pd.read_parquet(annual_path) + if "year" in existing.columns and int(year) in existing["year"].values: + logger.info("Year %d already processed — skipping Phase 3.", year) + return {"skipped": True, "year": year, "path": annual_path} + + logger.info( + "Phase 3 — ET intersection: %s/%s/%s year=%d", state, district, block, year + ) + + # 1. Load farm boundaries + farm_path = _farm_parquet_path(state, district, block) + if not os.path.exists(farm_path): + raise FileNotFoundError( + f"Farm boundaries parquet not found at {farm_path}. " + "Run Phases 1 & 2 first." + ) + gdf = gpd.read_parquet(farm_path) + logger.info("Loaded %d farm polygons.", len(gdf)) + + bbox = gdf.total_bounds # (minx, miny, maxx, maxy) + aez_zones = _get_aez_zones(state, district, block) + print(aez_zones) + + # 2. Load local AET raster(s) — one per AEZ zone the tehsil straddles, + # mosaicked together by _read_raster_clipped via rasterio.merge + # aet_paths = [_local_aet_path(z, year) for z in aez_zones] + # if not any(os.path.exists(p) for p in aet_paths): + # raise FileNotFoundError( + # f"No local AET raster found for zones {aez_zones}: {aet_paths}\n" + # f"Place the file(s) at: {LOCAL_ET_RASTERS_PATH}/merge_AET__{year}_cog.tif" + # ) + # logger.info("Reading local AET raster(s): %s", aet_paths) + # aet_data, aet_transform = _read_raster_clipped(aet_paths, bbox) + # logger.info("AET loaded: %d bands, shape=%s", aet_data.shape[0], aet_data.shape[1:]) + + # # 3. Load local PET raster(s) (optional) + # pet_data, pet_transform = None, None + # pet_paths = [_local_pet_path(z, year) for z in aez_zones] + # if any(os.path.exists(p) for p in pet_paths): + # logger.info("Reading local PET raster(s): %s", pet_paths) + # pet_data, pet_transform = _read_raster_clipped(pet_paths, bbox) + # logger.info("PET loaded: %d bands, shape=%s", pet_data.shape[0], pet_data.shape[1:]) + # else: + # logger.warning( + # "No local PET raster found for zones %s — MAI will not be computed.", aez_zones + # ) + + # # 4. Zonal statistics + # gdf = _run_zonal_stats(gdf, aet_data, aet_transform, pet_data, pet_transform) + + # # 5. Write 3-file schema + # static_path = _save_static_parquet(gdf, state, district, block) + # annual_path = _save_annual_parquet(gdf, state, district, block, year) + # monthly_path = _save_monthly_parquet(gdf, state, district, block, year) + + # summary = { + # "state": state, "district": district, "block": block, "year": year, + # "farm_count": len(gdf), + # "paths": { + # "static": static_path, + # "annual": annual_path, + # "monthly": monthly_path, + # }, + # } + # if "aet_annual" in gdf.columns and gdf["aet_annual"].notna().any(): + # summary["avg_aet_annual"] = round(float(gdf["aet_annual"].mean()), 4) + # if "mai_annual" in gdf.columns and gdf["mai_annual"].notna().any(): + # summary["avg_mai_annual"] = round(float(gdf["mai_annual"].mean()), 4) + # if "kharif_water_stress" in gdf.columns: + # summary["kharif_stress_farms"] = int(gdf["kharif_water_stress"].sum()) + + # logger.info("Phase 3 complete: %s", summary) + # return summary + return {"task" : "done"} + + +# ── multi-year analysis ──────────────────────────────────────────────────────── + +def compute_multi_year_water_stress( + state: str, + district: str, + block: str, + start_year: int = 2017, + end_year: int = 2024, +) -> dict: + """ + Run Phase 3 for each year in [start_year, end_year] and compute + cross-year frequency and intensity indicators: + + kharif_water_stress_years — number of years with kharif stress + return_period_years — N / stress_years (NaN if 0 stress years) + water_stress_intensity_mai — mean kharif MAI over stress years + + The annual parquet is updated with per-year rows. + A separate farm_water_stress_summary.parquet is written with the + cross-year indicators appended to the static columns. + + Returns + ------- + dict summary with paths and aggregate statistics. + """ + logger.info( + "Multi-year water stress: %s/%s/%s %d–%d", + state, district, block, start_year, end_year, + ) + + farm_path = _farm_parquet_path(state, district, block) + if not os.path.exists(farm_path): + raise FileNotFoundError(f"Farm parquet not found: {farm_path}") + + base_gdf = gpd.read_parquet(farm_path) + num_farms = len(base_gdf) + bbox = base_gdf.total_bounds + aez_zones = _get_aez_zones(state, district, block) + years = list(range(start_year, end_year + 1)) + + kharif_stress_count = np.zeros(num_farms, dtype=int) + kharif_mai_sum_stress = np.zeros(num_farms, dtype=float) + years_processed = 0 + + for year in years: + logger.info("── Processing year %d ──", year) + + aet_paths = [_local_aet_path(z, year) for z in aez_zones] + pet_paths = [_local_pet_path(z, year) for z in aez_zones] + + if not any(os.path.exists(p) for p in aet_paths): + logger.warning("AET raster(s) missing for year %d (zones %s) — skipping.", year, aez_zones) + continue + + try: + aet_data, aet_transform = _read_raster_clipped(aet_paths, bbox) + pet_data, pet_transform = ( + _read_raster_clipped(pet_paths, bbox) + if any(os.path.exists(p) for p in pet_paths) + else (None, None) + ) + except Exception as exc: + logger.warning("Year %d: raster read failed — %s", year, exc) + continue + + year_gdf = base_gdf.copy() + year_gdf = _run_zonal_stats(year_gdf, aet_data, aet_transform, pet_data, pet_transform) + + if "kharif_water_stress" not in year_gdf.columns: + logger.warning("Year %d: MAI not computed (PET missing?). Skipping.", year) + continue + + # Save this year into annual + monthly parquets + _save_annual_parquet(year_gdf, state, district, block, year) + _save_monthly_parquet(year_gdf, state, district, block, year) + + years_processed += 1 + is_stress = year_gdf["kharif_water_stress"].values.astype(bool) + kharif_mai_values = year_gdf["kharif_mai"].values + kharif_stress_count += is_stress.astype(int) + stress_mask = is_stress & np.isfinite(kharif_mai_values) + kharif_mai_sum_stress[stress_mask] += kharif_mai_values[stress_mask] + + # Cross-year indicators + result_gdf = base_gdf.copy() + result_gdf["total_years"] = years_processed + result_gdf["kharif_water_stress_years"] = kharif_stress_count + + with np.errstate(invalid="ignore", divide="ignore"): + rp = years_processed / kharif_stress_count.astype(float) + rp[kharif_stress_count == 0] = np.nan + result_gdf["return_period_years"] = np.round(rp, 2) + + with np.errstate(invalid="ignore", divide="ignore"): + intensity = kharif_mai_sum_stress / kharif_stress_count.astype(float) + intensity[kharif_stress_count == 0] = np.nan + result_gdf["water_stress_intensity_mai"] = np.round(intensity, 4) + + # Save summary parquet + out_path = os.path.join( + _block_dir(state, district, block), "farm_water_stress_summary.parquet" + ) + # keep static cols + cross-year indicators, no geometry duplication + summary_cols = [c for c in result_gdf.columns if not c.startswith(("aet_", "pet_", "mai_"))] + result_gdf[summary_cols].to_parquet(out_path, index=False) + logger.info("Water stress summary parquet saved → %s", out_path) + + valid = result_gdf["return_period_years"].notna().sum() + summary = { + "state": state, "district": district, "block": block, + "years_range": f"{start_year}–{end_year}", + "years_processed": years_processed, + "farm_count": num_farms, + "farms_with_any_stress": int((kharif_stress_count > 0).sum()), + "avg_return_period": round(float(result_gdf["return_period_years"].mean()), 2) if valid > 0 else None, + "avg_stress_intensity_mai": round(float(result_gdf["water_stress_intensity_mai"].mean()), 4) if valid > 0 else None, + "path": out_path, + } + logger.info("Multi-year analysis complete: %s", summary) + return summary diff --git a/computing/farm_boundaries/farm_boundary.py b/computing/farm_boundaries/farm_boundary.py new file mode 100644 index 00000000..e9a0030c --- /dev/null +++ b/computing/farm_boundaries/farm_boundary.py @@ -0,0 +1,112 @@ +""" +Celery task that orchestrates the three-phase farm boundary pipeline: + + Phase 1 — fetch_raw.fetch_raw_boundaries() + Queries the AnthroKrishi API per S2 cell and saves raw JSON files + to disk with a crash-safe manifest. + + Phase 2 — convert.convert_to_geoparquet() + Reads raw JSON via DuckDB, filters farm field polygons, clips to + the tehsil boundary with GeoPandas, and writes a GeoParquet file. + + Phase 3 — et_intersection.intersect_et_with_farms() [OPTIONAL] + Downloads AET raster from Google Earth Engine, computes per-farm + monthly ET via zonal statistics, and writes an enhanced parquet. + +The task is wired to the "nrm" Celery queue (same as all other CoRE Stack +pipelines) and supports automatic retries on transient failures. + +Triggered via the Django API: + POST /api/v1/generate_farm_boundaries/ + { + "state": "rajasthan", + "district": "jaipur", + "block": "sanganer", + "api_key": "AIzaSy...", + "year": 2017 ← optional, enables Phase 3 + } +""" + +import logging + +from nrm_app.celery import app + +from .convert import convert_to_geoparquet +from .fetch_raw import fetch_raw_boundaries + +logger = logging.getLogger(__name__) + + +@app.task(bind=True, max_retries=3, default_retry_delay=60) +def build_farm_boundary_map(self, state: str, district: str, block: str, api_key: str, year: int = None, overwrite = False): + """ + Celery task: runs Phase 1, Phase 2, and optionally Phase 3. + + Parameters + ---------- + state, district, block : str + Lower-cased administrative names. + api_key : str + AnthroKrishi / Agricultural Understanding API key. + year : int, optional + If provided, runs Phase 3 (ET intersection) for the given year. + Valid range: 2017–2024. + + Returns + ------- + dict + Combined summary from all phases. + """ + logger.info( + "Farm boundary pipeline started — state=%s district=%s block=%s", + state, district, block, + ) + + try: + # ── Phase 1: Fetch ────────────────────────────────────────────────── + phase1_summary = fetch_raw_boundaries( + state=state, + district=district, + block=block, + api_key=api_key, + resume=True, # safe to retry; already-fetched cells are skipped + ) + logger.info("Phase 1 done: %s", phase1_summary) + + # ── Phase 2: Convert ──────────────────────────────────────────────── + phase2_summary = convert_to_geoparquet( + state=state, + district=district, + block=block, + overwrite=False, # skip if parquet already exists + ) + logger.info("Phase 2 done: %s", phase2_summary) + + # ── Phase 3: ET Intersection (optional) ───────────────────────────── + phase3_summary = None + if year is not None: + from .et_intersection import intersect_et_with_farms + + logger.info("Phase 3 — ET intersection for year %d", year) + phase3_summary = intersect_et_with_farms( + state=state, + district=district, + block=block, + year=year, + ) + logger.info("Phase 3 done: %s", phase3_summary) + + except Exception as exc: + logger.exception( + "Farm boundary pipeline failed for %s/%s/%s: %s", + state, district, block, exc, + ) + raise self.retry(exc=exc) + + result = { + "phase1": phase1_summary, + "phase2": phase2_summary, + "phase3": phase3_summary, + } + logger.info("Farm boundary pipeline completed successfully: %s", result) + return result diff --git a/computing/farm_boundaries/fetch_raw.py b/computing/farm_boundaries/fetch_raw.py new file mode 100644 index 00000000..fb698f93 --- /dev/null +++ b/computing/farm_boundaries/fetch_raw.py @@ -0,0 +1,340 @@ +""" +Phase 1 — Fetch raw farm boundary data from the Google AnthroKrishi +(Agricultural Understanding) API and persist each S2-cell response as a +JSON file on disk. + +**Optimised with async I/O**: Uses ``aiohttp`` to fire concurrent +API requests (controlled by a semaphore) instead of the original +sequential loop. For Sanganer (810 cells) this reduces Phase 1 from +~6 minutes to under 30 seconds. + +Directory layout after a successful run: + data/farm_boundaries////raw/.json + data/farm_boundaries////manifest.json + +The manifest records which cells were successfully fetched so that +Phase 2 (convert.py) and any future resume run know exactly what to +process. Cells that returned an error or an empty landscape are +recorded separately so you can inspect them without re-querying. + +Usage (standalone / debug): + from computing.farm_boundaries.fetch_raw import fetch_raw_boundaries + fetch_raw_boundaries("rajasthan", "jaipur", "sanganer", + api_key="AIzaSy...") +""" + +import asyncio +import json +import logging +import os +import time + +import aiohttp +import geopandas as gpd +import s2sphere +from shapely.geometry import box + +from utilities.constants import FARM_BOUNDARIES_PATH, SOI_TEHSIL + +logger = logging.getLogger(__name__) + +# ── AnthroKrishi REST endpoint ──────────────────────────────────────────────── +ANTHROKRISHI_API_URL = ( + "https://agriculturalunderstanding.googleapis.com/v1:lookupLandscape" +) + +# S2 level 13 ≈ 1 km × 1 km tiles +S2_LEVEL = 13 + +# ── Concurrency tuning ─────────────────────────────────────────────────────── +# Maximum number of API requests in flight at the same time. +# 20 is a conservative default that stays within Google API quota limits. +MAX_CONCURRENT_REQUESTS = 20 + +# How many cells to process before flushing the manifest to disk. +# Lower = more crash-safe but more I/O; higher = faster but riskier. +MANIFEST_FLUSH_INTERVAL = 25 + + +# ── helpers ─────────────────────────────────────────────────────────────────── + + +def _get_tehsil_polygon(state: str, district: str, block: str): + """ + Load the tehsil polygon from the shared SOI shapefile. + Returns a single Shapely geometry in EPSG:4326. + Raises ValueError if the tehsil is not found. + """ + soi = gpd.read_file(SOI_TEHSIL) + mask = ( + (soi["STATE"].str.lower() == state) + & (soi["District"].str.lower() == district) + & (soi["TEHSIL"].str.lower() == block) + ) + subset = soi[mask] + if subset.empty: + raise ValueError( + f"Tehsil not found in SOI shapefile: state={state}, " + f"district={district}, block={block}" + ) + # Dissolve to a single geometry in case the tehsil has multiple rows. + return subset.dissolve().geometry.iloc[0] + + +def _get_s2_cells_for_bbox(tehsil_geom) -> list: + """ + Return all Level-13 S2 cells whose centres fall inside the tehsil's + bounding box. Using the bounding box (not the exact polygon) is + simpler and faster; Phase 2 clips the result to the exact boundary. + + Returns a list of s2sphere.CellId objects. + """ + minx, miny, maxx, maxy = tehsil_geom.bounds + + ll_lo = s2sphere.LatLng.from_degrees(miny, minx) + ll_hi = s2sphere.LatLng.from_degrees(maxy, maxx) + rect = s2sphere.LatLngRect.from_point_pair(ll_lo, ll_hi) + + coverer = s2sphere.RegionCoverer() + coverer.min_level = S2_LEVEL + coverer.max_level = S2_LEVEL + # Allow a generous upper bound so large tehsils are not under-covered. + coverer.max_cells = 1_000_000 + + covering = coverer.get_covering(rect) + return list(covering) + + +def _output_dir(state: str, district: str, block: str) -> str: + """Return (and create) the raw-data directory for this tehsil.""" + path = os.path.join(FARM_BOUNDARIES_PATH, state, district, block, "raw") + os.makedirs(path, exist_ok=True) + return path + + +def _manifest_path(state: str, district: str, block: str) -> str: + return os.path.join(FARM_BOUNDARIES_PATH, state, district, block, "manifest.json") + + +def _load_manifest(manifest_file: str) -> dict: + if os.path.exists(manifest_file): + with open(manifest_file) as f: + return json.load(f) + return {"fetched": [], "empty": [], "errors": []} + + +def _save_manifest(manifest_file: str, manifest: dict): + with open(manifest_file, "w") as f: + json.dump(manifest, f, indent=2) + + +# ── async fetching engine ──────────────────────────────────────────────────── + + +async def _fetch_one_cell_async( + session: aiohttp.ClientSession, + semaphore: asyncio.Semaphore, + cell_id: s2sphere.CellId, + api_key: str, + raw_dir: str, +) -> dict: + """ + Fetch a single S2 cell asynchronously. + + Returns a dict: {"token": str, "status": "fetched"|"empty"|"error"} + """ + token = cell_id.to_token() + payload = { + "locationSpecifier": { + "s2CellId": str(cell_id.id()) + } + } + + async with semaphore: + try: + async with session.post( + ANTHROKRISHI_API_URL, + params={"key": api_key}, + json=payload, + timeout=aiohttp.ClientTimeout(total=30), + ) as response: + response.raise_for_status() + + if not await response.read(): + return {"token": token, "status": "empty"} + + data = await response.json() + + except Exception as exc: + logger.warning("Error fetching cell %s: %s", token, exc) + return {"token": token, "status": "error"} + + # Save raw response to disk + out_path = os.path.join(raw_dir, f"{token}.json") + with open(out_path, "w") as f: + json.dump(data, f) + + has_landscape = bool(data.get("landscape")) + return { + "token": token, + "status": "fetched" if has_landscape else "empty", + } + + +async def _fetch_all_cells_async( + cells_to_fetch: list, + api_key: str, + raw_dir: str, + manifest_file: str, + manifest: dict, + max_concurrent: int = MAX_CONCURRENT_REQUESTS, +) -> dict: + """ + Fetch all cells concurrently using aiohttp with a semaphore + to cap the number of in-flight requests. + + Cells are processed in batches; the manifest is flushed after + each batch for crash safety. + """ + semaphore = asyncio.Semaphore(max_concurrent) + total = len(cells_to_fetch) + + # Split into batches for manifest flush points + batch_size = MANIFEST_FLUSH_INTERVAL + results_summary = {"fetched": 0, "empty": 0, "errors": 0} + + async with aiohttp.ClientSession() as session: + for batch_start in range(0, total, batch_size): + batch = cells_to_fetch[batch_start : batch_start + batch_size] + batch_end = min(batch_start + len(batch), total) + + logger.info( + "Fetching cells %d–%d of %d (concurrency=%d)", + batch_start + 1, batch_end, total, max_concurrent, + ) + + # Fire all requests in this batch concurrently + tasks = [ + _fetch_one_cell_async(session, semaphore, cell_id, api_key, raw_dir) + for cell_id in batch + ] + results = await asyncio.gather(*tasks) + + # Update manifest with batch results + for result in results: + status = result["status"] + token = result["token"] + + if status == "fetched": + manifest["fetched"].append(token) + results_summary["fetched"] += 1 + elif status == "empty": + manifest["empty"].append(token) + results_summary["empty"] += 1 + else: + manifest["errors"].append(token) + results_summary["errors"] += 1 + + # Flush manifest after each batch (crash-safe checkpoint) + _save_manifest(manifest_file, manifest) + + return results_summary + + +# ── public entry point ──────────────────────────────────────────────────────── + + +def fetch_raw_boundaries( + state: str, + district: str, + block: str, + api_key: str, + max_concurrent: int = MAX_CONCURRENT_REQUESTS, + resume: bool = True, +) -> dict: + """ + Phase 1 pipeline: fetch raw AnthroKrishi data for every S2 cell + covering the tehsil bounding box and save each response to disk. + + Uses async I/O to fetch up to ``max_concurrent`` cells in parallel, + reducing wall-clock time by 10-20× compared to the sequential version. + + Parameters + ---------- + state, district, block : str + Lower-cased administrative names (must match SOI shapefile columns). + api_key : str + AnthroKrishi / Agricultural Understanding API key. + max_concurrent : int + Maximum number of simultaneous HTTP requests (default 20). + resume : bool + If True, skip cells that are already recorded in the manifest so + an interrupted run can be safely restarted without re-fetching. + + Returns + ------- + dict + Summary with counts of fetched / empty / error cells and the + path to the raw data directory. + """ + t0 = time.time() + logger.info("Phase 1 — fetching farm boundaries for %s/%s/%s", state, district, block) + + # 1. Load tehsil boundary ------------------------------------------------ + tehsil_geom = _get_tehsil_polygon(state, district, block) + logger.info("Tehsil polygon loaded. Bounding box: %s", tehsil_geom.bounds) + + # 2. Enumerate S2 cells -------------------------------------------------- + cells = _get_s2_cells_for_bbox(tehsil_geom) + logger.info("S2 level-%d cells to query: %d", S2_LEVEL, len(cells)) + + # 3. Set up output paths ------------------------------------------------- + raw_dir = _output_dir(state, district, block) + manifest_file = _manifest_path(state, district, block) + manifest = _load_manifest(manifest_file) + + already_done = set(manifest["fetched"] + manifest["empty"] + manifest["errors"]) + + # 4. Determine which cells still need fetching --------------------------- + if resume: + cells_to_fetch = [c for c in cells if c.to_token() not in already_done] + skipped = len(cells) - len(cells_to_fetch) + if skipped: + logger.info("Resuming: %d cells already done, %d remaining.", skipped, len(cells_to_fetch)) + else: + cells_to_fetch = cells + + # 5. Fetch concurrently -------------------------------------------------- + if cells_to_fetch: + logger.info( + "Starting async fetch: %d cells, max_concurrent=%d", + len(cells_to_fetch), max_concurrent, + ) + results = asyncio.run( + _fetch_all_cells_async( + cells_to_fetch, api_key, raw_dir, manifest_file, manifest, max_concurrent, + ) + ) + logger.info( + "Async fetch complete: fetched=%d, empty=%d, errors=%d", + results["fetched"], results["empty"], results["errors"], + ) + else: + logger.info("All cells already fetched — nothing to do.") + + elapsed = time.time() - t0 + + summary = { + "state": state, + "district": district, + "block": block, + "total_cells": len(cells), + "fetched": len(manifest["fetched"]), + "empty": len(manifest["empty"]), + "errors": len(manifest["errors"]), + "elapsed_seconds": round(elapsed, 1), + "raw_dir": raw_dir, + "manifest": manifest_file, + } + logger.info("Phase 1 complete in %.1fs: %s", elapsed, summary) + return summary diff --git a/computing/urls.py b/computing/urls.py index 203ed0aa..56346b70 100644 --- a/computing/urls.py +++ b/computing/urls.py @@ -63,6 +63,11 @@ name="change_detection_vector", ), path("crop_grid/", api.crop_grid, name="crop_grid"), + path( + "generate_farm_boundaries/", + api.generate_farm_boundaries, + name="generate_farm_boundaries", + ), path("tree_health_raster/", api.tree_health_raster, name="tree_health_raster"), path("tree_health_vector/", api.tree_health_vector, name="tree_health_vector"), path("stream_order/", api.stream_order, name="stream_order"), diff --git a/plot_tally_timeseries.py b/plot_tally_timeseries.py new file mode 100644 index 00000000..d4dd229b --- /dev/null +++ b/plot_tally_timeseries.py @@ -0,0 +1,289 @@ +""" +Time-series tally plots — corrected version. +Selects 4 representative farms that Sir can visually verify. + +Plot 1 (2x2): Four cases showing our vs Shuvam side-by-side +Plot 2: Gap-fill before/after demonstration + pipeline comparison + +Saves to: data/tally_plots/ +""" +import dotenv; dotenv.load_dotenv('.env') + +import geopandas as gpd +import pandas as pd +import numpy as np +import rasterio +import rasterio.windows +import matplotlib.pyplot as plt +import os, sys +sys.path.insert(0, '.') +from computing.farm_boundaries.et_intersection import _gap_fill_monthly_farms + +TEHSIL = "sanganer" +BASE = f"data/farm_boundaries/rajasthan/jaipur/{TEHSIL}" +SHUVAM = f"data/et_rasters/shuvam_mai_{TEHSIL}_2018.tif" +OUT_DIR = "data/tally_plots" +os.makedirs(OUT_DIR, exist_ok=True) + +MONTH_LABELS = ['Jul-18','Aug-18','Sep-18','Oct-18','Nov-18','Dec-18', + 'Jan-19','Feb-19','Mar-19','Apr-19','May-19','Jun-19'] +CROP_DATES = ['2018-07-01','2018-08-01','2018-09-01','2018-10-01', + '2018-11-01','2018-12-01','2019-01-01','2019-02-01', + '2019-03-01','2019-04-01','2019-05-01','2019-06-01'] + +# ── load ───────────────────────────────────────────────────────────────────── +print("Loading...") +static = gpd.read_parquet(f"{BASE}/farm_static.parquet") +monthly = pd.read_parquet(f"{BASE}/farm_monthly.parquet") +tally = pd.read_parquet(f"{BASE}/mai_tally_2018.parquet") + +monthly['ds'] = monthly['date'].astype(str).str[:10] +monthly_cy = monthly[monthly['ds'].isin(CROP_DATES)].copy() +monthly_cy['sk'] = monthly_cy['ds'].map({d:i for i,d in enumerate(CROP_DATES)}) + +# Load Shuvam raster into memory +with rasterio.open(SHUVAM) as src: + rcrs = src.crs + rdata = src.read().astype("float32") + rtf = src.transform + nd = float(src.nodata) if src.nodata else -9999.0 +rdata[rdata == nd] = np.nan +rdata[rdata <= -9999] = np.nan +farms_rp = static.to_crs(rcrs) + + +def shuvam_series(fid): + """12-month MAI from Shuvam's raster for one farm.""" + geom = farms_rp.loc[farms_rp['farm_id']==fid, 'geometry'].values[0] + b = geom.bounds + w = rasterio.windows.from_bounds(*b, rtf) + r0, c0 = max(0, int(w.row_off)), max(0, int(w.col_off)) + r1 = min(rdata.shape[1], int(w.row_off + w.height)+1) + c1 = min(rdata.shape[2], int(w.col_off + w.width)+1) + if r1 <= r0 or c1 <= c0: + return np.full(12, np.nan) + clip = rdata[:12, r0:r1, c0:c1] + return np.array([np.nanmean(clip[i]) if np.isfinite(clip[i]).any() + else np.nan for i in range(12)]) + + +def our_series(fid): + """12-month MAI from our parquet, crop-year order.""" + rows = monthly_cy[monthly_cy['farm_id']==fid].sort_values('sk') + if len(rows) != 12: + return np.full(12, np.nan) + return rows['mai'].values.astype(float) + + +# ── select 4 farms wisely ──────────────────────────────────────────────────── +# CRITICAL: only from farms with BOTH annual values valid +both_ok = tally.dropna(subset=['mai_annual','shuvam_annual']).copy() +both_ok['diff'] = (both_ok['mai_annual'] - both_ok['shuvam_annual']).abs() + +# 1. Best match (smallest diff, > 0 to avoid trivial) +best = both_ok[both_ok['diff'] > 0].nsmallest(10, 'diff')['farm_id'].iloc[0] + +# 2. Typical farm (diff near the median) +med_diff = both_ok['diff'].median() +typical = both_ok.iloc[(both_ok['diff'] - med_diff).abs().argsort()[:1]]['farm_id'].iloc[0] + +# 3. Largest difference +worst = both_ok.nlargest(1, 'diff')['farm_id'].iloc[0] + +# 4. A "rajasthan_jaipur_sanganer_000000" — the well-known reference farm +ref_farm = "rajasthan_jaipur_sanganer_000000" + +FARMS = [ + (best, "Farm A: best match"), + (ref_farm, "Farm B: reference farm"), + (typical, "Farm C: median-diff farm"), + (worst, "Farm D: largest difference"), +] + +print("Selected farms:") +for fid, label in FARMS: + row = both_ok[both_ok['farm_id']==fid] + d = row['diff'].values[0] if len(row) else float('nan') + print(f" {label}: {fid} diff={d:.4f}") + +# ═══════════════════════════════════════════════════════════════════════════════ +# PLOT 1: 2×2 time-series comparison +# ═══════════════════════════════════════════════════════════════════════════════ +fig, axes = plt.subplots(2, 2, figsize=(16, 10)) +fig.suptitle( + "MAI Time-Series Tally — Sanganer 2018\n" + "Our Pipeline vs Shuvam's GEE Raster (crop-year Jul 2018 – Jun 2019)", + fontsize=14, fontweight='bold', y=1.01 +) + +x = np.arange(12) + +for ax, (fid, label) in zip(axes.flat, FARMS): + ours = our_series(fid) + shuv = shuvam_series(fid) + + ax.plot(x, shuv, 'o-', color='#E85D04', linewidth=2, markersize=6, label="Shuvam (GEE)") + ax.plot(x, ours, 's--', color='#0077B6', linewidth=2, markersize=6, label="Ours (local)") + + # Difference shading + valid = np.isfinite(ours) & np.isfinite(shuv) + if valid.any(): + ax.fill_between(x, ours, shuv, where=valid, alpha=0.12, color='purple') + + # Kharif + ax.axvspan(-0.5, 3.5, alpha=0.06, color='green') + ax.text(1.5, 0.97, 'Kharif', transform=ax.get_xaxis_transform(), + ha='center', va='top', fontsize=8, color='green', alpha=0.7) + + # Stats + if valid.any(): + d = np.abs(ours[valid] - shuv[valid]) + ax.set_title(f"{label}\n{fid}\n" + f"mean|diff|={d.mean():.4f} max|diff|={d.max():.4f}", + fontsize=9) + else: + ax.set_title(f"{label}\n{fid}\n(no overlapping data)", fontsize=9) + + ax.axhline(0.5, color='orange', ls=':', lw=1, alpha=0.7, label='Moderate stress (0.5)') + ax.axhline(0.25, color='red', ls=':', lw=1, alpha=0.7, label='Severe stress (0.25)') + ax.set_xticks(x) + ax.set_xticklabels(MONTH_LABELS, rotation=45, ha='right', fontsize=8) + ax.set_ylabel("MAI (AET/PET)", fontsize=9) + ax.set_ylim(-0.05, 1.05) + ax.legend(fontsize=7, loc='upper right') + ax.grid(True, alpha=0.3) + +plt.tight_layout() +p1 = f"{OUT_DIR}/mai_timeseries_tally_sanganer_2018.png" +plt.savefig(p1, dpi=150, bbox_inches='tight') +print(f"Saved: {p1}") +plt.close() + + +# ═══════════════════════════════════════════════════════════════════════════════ +# PLOT 2: Gap-fill demo (left) + pipeline comparison (right) +# ═══════════════════════════════════════════════════════════════════════════════ +print("Generating gap-fill demo...") + +# Use the reference farm (known to have all 12 months) +actual_mai = our_series(ref_farm) + +# Simulate 2 missing months: Aug (idx 1) and Sep (idx 2) in crop order +# But gap-fill works on calendar order. We need to convert. +# Crop order: Jul(0) Aug(1) Sep(2) Oct(3) Nov(4) Dec(5) Jan(6) Feb(7) Mar(8) Apr(9) May(10) Jun(11) +# Calendar: Jan=0 Feb=1 Mar=2 Apr=3 May=4 Jun=5 Jul=6 Aug=7 Sep=8 Oct=9 Nov=10 Dec=11 +# Mapping: crop_idx -> cal_month: {0->6, 1->7, 2->8, 3->9, 4->10, 5->11, 6->0, 7->1, 8->2, 9->3, 10->4, 11->5} +crop_to_cal = [6, 7, 8, 9, 10, 11, 0, 1, 2, 3, 4, 5] + +# Build calendar-order matrix from actual values +cal_complete = np.full(12, np.nan) +for ci, cal_i in enumerate(crop_to_cal): + cal_complete[cal_i] = actual_mai[ci] + +# Now set Aug (cal idx 7) and Sep (cal idx 8) to NaN +cal_missing = cal_complete.copy() +cal_missing[7] = np.nan # Aug +cal_missing[8] = np.nan # Sep + +# Apply gap-fill +cal_filled = _gap_fill_monthly_farms(cal_missing.reshape(1, 12))[0] + +# Convert back to crop order for plotting +actual_crop = actual_mai.copy() +missing_crop = np.array([cal_missing[crop_to_cal[i]] for i in range(12)]) +filled_crop = np.array([cal_filled[crop_to_cal[i]] for i in range(12)]) + +# Shuvam for same farm +shuv_ref = shuvam_series(ref_farm) + +fig2, axes2 = plt.subplots(1, 2, figsize=(18, 6)) +fig2.suptitle( + "Missing Data Handling & Gap-Fill Demonstration\n" + f"Farm: {ref_farm} | Sanganer 2018", + fontsize=13, fontweight='bold' +) + +# ── Left: Before vs After gap-fill ─────────────────────────────────────────── +ax = axes2[0] +ax.set_title("Before vs After Gap-Fill\n(Aug & Sep artificially set to NaN to demonstrate)", fontsize=10) + +ax.plot(x, actual_crop, 'o-', color='green', linewidth=2, markersize=6, + label='Actual (complete data)', zorder=4, alpha=0.5) +ax.plot(x, missing_crop, 'x--', color='#999', linewidth=1.5, markersize=10, + label='Before gap-fill (Aug, Sep = NaN)', zorder=3, markeredgewidth=2) +ax.plot(x, filled_crop, 's-', color='#0077B6', linewidth=2.5, markersize=8, + label='After gap-fill', zorder=5) + +# Annotate the two filled points +for ci in range(12): + if np.isnan(missing_crop[ci]) and np.isfinite(filled_crop[ci]): + cal_m = crop_to_cal[ci] + nbrs = [crop_to_cal.index(n) for n in [6, 8] if cal_m == 7] or \ + [crop_to_cal.index(n) for n in [7, 9] if cal_m == 8] or [] + # Get label text + if cal_m == 7: # Aug filled from Jul & Sep + lbl = f"Filled: {filled_crop[ci]:.3f}\nmean(Jul, Sep)" + elif cal_m == 8: # Sep filled from Aug & Oct + lbl = f"Filled: {filled_crop[ci]:.3f}\nmean(Aug, Oct)" + else: + lbl = f"Filled: {filled_crop[ci]:.3f}" + + yoff = 0.15 if ci == 1 else 0.12 + ax.annotate(lbl, xy=(ci, filled_crop[ci]), + xytext=(ci + 0.7, filled_crop[ci] + yoff), + fontsize=8, color='#0077B6', + arrowprops=dict(arrowstyle='->', color='#0077B6', lw=1.2), + bbox=dict(boxstyle='round,pad=0.2', facecolor='#dbeafe', alpha=0.8)) + +ax.axvspan(-0.5, 3.5, alpha=0.07, color='green') +ax.text(1.5, 0.96, 'Kharif', transform=ax.get_xaxis_transform(), + ha='center', va='top', fontsize=9, color='green', alpha=0.8) +ax.axhline(0.5, color='orange', ls=':', lw=1.2, label='Moderate stress (0.5)') +ax.axhline(0.25, color='red', ls=':', lw=1.2, label='Severe stress (0.25)') + +rule_text = ("Gap-fill rules (mirrors Shuvam's ET_Applications/helper.py):\n" + " Aug → mean(Jul, Sep)\n" + " Sep → mean(Aug, Oct)\n" + " Jul → Aug only (crop-year start boundary)\n" + " Jun → May only (crop-year end boundary)\n" + " If no valid neighbour → stays NaN") +ax.text(0.55, 0.98, rule_text, transform=ax.transAxes, fontsize=8, + va='top', bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.85)) + +ax.set_xticks(x) +ax.set_xticklabels(MONTH_LABELS, rotation=45, ha='right', fontsize=9) +ax.set_ylabel("MAI (AET/PET)", fontsize=10) +ax.set_ylim(-0.05, 1.15) +ax.legend(fontsize=8.5, loc='upper left') +ax.grid(True, alpha=0.3) + +# ── Right: Our pipeline vs Shuvam ──────────────────────────────────────────── +ax2 = axes2[1] +valid = np.isfinite(actual_crop) & np.isfinite(shuv_ref) +d = np.abs(actual_crop[valid] - shuv_ref[valid]) +ax2.set_title(f"Our Pipeline vs Shuvam's GEE Raster (same farm)\n" + f"mean|diff|={d.mean():.4f} max|diff|={d.max():.4f}", fontsize=10) + +ax2.plot(x, shuv_ref, 'o-', color='#E85D04', linewidth=2.5, markersize=8, label="Shuvam (GEE)") +ax2.plot(x, actual_crop, 's--', color='#0077B6', linewidth=2.5, markersize=8, label="Our pipeline") +ax2.fill_between(x, actual_crop, shuv_ref, where=valid, alpha=0.15, color='purple', label='Difference') + +ax2.axvspan(-0.5, 3.5, alpha=0.07, color='green') +ax2.text(1.5, 0.96, 'Kharif', transform=ax2.get_xaxis_transform(), + ha='center', va='top', fontsize=9, color='green', alpha=0.8) +ax2.axhline(0.5, color='orange', ls=':', lw=1.2, label='Moderate stress (0.5)') +ax2.axhline(0.25, color='red', ls=':', lw=1.2, label='Severe stress (0.25)') +ax2.set_xticks(x) +ax2.set_xticklabels(MONTH_LABELS, rotation=45, ha='right', fontsize=9) +ax2.set_ylabel("MAI (AET/PET)", fontsize=10) +ax2.set_ylim(-0.05, 1.15) +ax2.legend(fontsize=8.5) +ax2.grid(True, alpha=0.3) + +plt.tight_layout() +p2 = f"{OUT_DIR}/gap_fill_and_tally_demo_2018.png" +plt.savefig(p2, dpi=150, bbox_inches='tight') +print(f"Saved: {p2}") +plt.close() + +print(f"\nAll plots saved to {OUT_DIR}/") diff --git a/regen_all_phase3.py b/regen_all_phase3.py new file mode 100644 index 00000000..ecc8409e --- /dev/null +++ b/regen_all_phase3.py @@ -0,0 +1,28 @@ +"""Re-run Phase 3 for Sanganer and Dudu with all three fixes applied: +1. -9999 nodata → NaN at read time (using src.nodata from raster metadata) +2. MAI capped at [0,1] — values > 1 are raster artifacts, logged and capped +3. Improved logging: nan_farms, valid_farms, severe count all reported +""" +import os +import dotenv +dotenv.load_dotenv('.env') + +import logging +logging.basicConfig( + level=logging.INFO, + format='%(asctime)s %(levelname)s %(name)s: %(message)s' +) + +from computing.farm_boundaries.et_intersection import intersect_et_with_farms + +print("=" * 60) +print("Re-running Phase 3 — Sanganer, Jaipur 2018") +print("=" * 60) +intersect_et_with_farms('rajasthan', 'jaipur', 'sanganer', year=2018, overwrite=True) + +print("\n" + "=" * 60) +print("Re-running Phase 3 — Dudu, Jaipur 2018") +print("=" * 60) +intersect_et_with_farms('rajasthan', 'jaipur', 'dudu', year=2018, overwrite=True) + +print("\nDone. All parquets regenerated with bug fixes applied.") diff --git a/show_parquets.py b/show_parquets.py new file mode 100644 index 00000000..d3528a45 --- /dev/null +++ b/show_parquets.py @@ -0,0 +1,50 @@ +"""Print sample rows from all three output parquets for KT.""" +import pandas as pd +import geopandas as gpd + +BASE = "data/farm_boundaries/rajasthan/jaipur/sanganer" + +print("=" * 80) +print("1. farm_static.parquet") +print("=" * 80) +static = gpd.read_parquet(f"{BASE}/farm_static.parquet") +print(f"Shape: {static.shape[0]:,} rows × {static.shape[1]} columns") +print(f"Columns: {list(static.columns)}") +print(f"CRS: {static.crs}") +print(f"\nSample (first 5 rows, geometry truncated):") +display = static.head(5).copy() +display['geometry'] = display['geometry'].apply(lambda g: str(g)[:60] + "...") +print(display.to_string(index=False)) +print(f"\nMemory: {static.memory_usage(deep=True).sum()/1e6:.1f} MB") + +print("\n" + "=" * 80) +print("2. farm_annual.parquet") +print("=" * 80) +annual = pd.read_parquet(f"{BASE}/farm_annual.parquet") +print(f"Shape: {annual.shape[0]:,} rows × {annual.shape[1]} columns") +print(f"Columns: {list(annual.columns)}") +print(f"\nSample (first 10 rows):") +print(annual.head(10).to_string(index=False)) +print(f"\nYear distribution: {annual['year'].value_counts().to_dict()}") +print(f"Kharif stress farms: {annual['kharif_water_stress'].sum():,} / {len(annual):,}") +print(f"Severe stress farms: {annual['kharif_severe_stress'].sum():,} / {len(annual):,}") +print(f"MAI annual stats: mean={annual['mai_annual'].mean():.4f}, " + f"median={annual['mai_annual'].median():.4f}, " + f"min={annual['mai_annual'].min():.4f}, max={annual['mai_annual'].max():.4f}") +print(f"NaN count: {annual['mai_annual'].isna().sum():,}") + +print("\n" + "=" * 80) +print("3. farm_monthly.parquet") +print("=" * 80) +monthly = pd.read_parquet(f"{BASE}/farm_monthly.parquet") +print(f"Shape: {monthly.shape[0]:,} rows × {monthly.shape[1]} columns") +print(f"Columns: {list(monthly.columns)}") +print(f"\nDate range: {monthly['date'].min()} to {monthly['date'].max()}") +print(f"\nSample (first 12 rows = one farm, all months):") +one_farm = monthly[monthly['farm_id'] == monthly['farm_id'].iloc[0]].sort_values('date') +print(one_farm.to_string(index=False)) +print(f"\nMonthly MAI stats:") +for d in sorted(monthly['date'].unique()): + subset = monthly[monthly['date'] == d] + m = subset['mai'] + print(f" {str(d)[:7]}: mean={m.mean():.4f} median={m.median():.4f} NaN={m.isna().sum():,}/{len(m):,}") diff --git a/utilities/constants.py b/utilities/constants.py index b0bb5dd4..f2d2026b 100644 --- a/utilities/constants.py +++ b/utilities/constants.py @@ -15,6 +15,8 @@ RASTERS_PATH = "data/rasters" CROP_GRID_PATH = "data/crop_grid" +FARM_BOUNDARIES_PATH = "data/farm_boundaries" +LOCAL_ET_RASTERS_PATH = "data/et_rasters" KML_PATH = "data/kml/" SHAPEFILE_DIR = "data/kml/shapefiles" @@ -306,6 +308,7 @@ WWF_HYDROSHEDS_DRAINAGE_DIRECTION = "WWF/HydroSHEDS/03DIR" PAN_INDIA_RASTER_FABDEM = "projects/corestack-datasets/assets/datasets/terrain/pan_india_terrain_raster_fabdem" SOI_TEHSIL = "data/admin-boundary/input/soi_tehsil.geojson" +AEZ_GEOJSON = "data/AEZ_GeoJSON.geojson" FABDEM = "projects/sat-io/open-datasets/FABDEM" WATERREJUVENATION = "projects/ee-corestackdev/assets/apps/waterrej/proj1" WATERREJ_LULCFORM = "projects/ee-corestackdev/assets/apps/waterrej/lulcfrom"