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..7538ec1f --- /dev/null +++ b/computing/farm_boundaries/et_intersection.py @@ -0,0 +1,797 @@ +""" +Phase 3 — Intersect AET & PET rasters with farm boundary polygons. + +Downloads the AET (Actual Evapotranspiration) and PET (Potential +Evapotranspiration) rasters from Google Earth Engine for the tehsil's +bounding box, runs zonal statistics against each farm polygon, computes +MAI (Moisture Adequacy Index = AET/PET), and produces an enhanced +GeoParquet with per-farm monthly AET, PET, MAI values and water stress +indicators. + +Data sources: + AET asset: projects/corestack-datasets-alpha/assets/datasets/ + et_downscale/aet_aez__ + PET asset: projects/core-stack-dev-3-helper/assets/ + et_downscale/pet_aez__ + 13 bands: b1–b12 (monthly in mm/day), b13 (annual) + Resolution: 30 meters + NoData: -9999 + +Water stress methodology (Drought Manual 2016 / Shivani-Shuvam): + MAI thresholds: 76–100% no stress, 51–75% mild, 26–50% moderate, 0–25% severe + Kharif water stress: MAI in moderate/severe range (MAI ≤ 50%) during Jul–Oct + Frequency: Return period = N / #kharif_water_stress_years + Intensity: Mean MAI over Kharif months in water stress years + Note: MAI-based classification indicates crop water stress, not drought. + Drought requires additional indicators (SPI/SPEI, VCI) per GoI definition. + +Output: + data/farm_boundaries////farm_boundaries_et.parquet + +Usage (standalone): + from computing.farm_boundaries.et_intersection import intersect_et_with_farms + result = intersect_et_with_farms("rajasthan", "jaipur", "sanganer", year=2017) +""" + +import logging +import os +import tempfile +import zipfile + +import ee +import geopandas as gpd +import numpy as np +import rasterio +import rasterio.features +import requests + +from utilities.constants import FARM_BOUNDARIES_PATH + +logger = logging.getLogger(__name__) + +# ── GEE asset path template ────────────────────────────────────────────────── +AET_ASSET_TEMPLATE = ( + "projects/corestack-datasets-alpha/assets/datasets/" + "et_downscale/aet_aez_{aez}_{year}" +) +PET_ASSET_TEMPLATE = ( + "projects/core-stack-dev-3-helper/assets/" + "et_downscale/pet_aez_{aez}_{year}" +) + +# AEZ zone mapping — which zone covers which region +# (extend this as more zones become available) +AEZ_ZONE_MAP = { + "rajasthan": 2, + "gujarat": 2, + "punjab": 2, + "haryana": 2, +} + +# Default GEE project for authentication +DEFAULT_GEE_PROJECT = "stackd-conversion123" + +# Nodata value used in the AET/PET rasters +AET_NODATA = -9999 + +# MAI (Moisture Adequacy Index) thresholds — Drought Manual 2016, Table 3.5 +# MAI = (AET / PET), expressed as ratio (manual uses percentage) +# 76–100% (0.76–1.00) → No drought +# 51–75% (0.51–0.75) → Mild drought +# 26–50% (0.26–0.50) → Moderate drought +# 0–25% (0.00–0.25) → Severe drought +MAI_NO_DROUGHT_THRESHOLD = 0.76 # MAI ≥ 0.76 → no drought +MAI_MILD_THRESHOLD = 0.51 # 0.51 ≤ MAI < 0.76 → mild drought +MAI_MODERATE_THRESHOLD = 0.26 # 0.26 ≤ MAI < 0.51 → moderate drought +# Below 0.26 → severe drought + +# Kharif water stress: MAI ∈ {moderate, severe} means MAI ≤ 0.50 +# Note: This indicates crop water stress, not drought (per professor's clarification) +KHARIF_WATER_STRESS_MAI_THRESHOLD = 0.50 + +# Kharif season months (July–October), 0-indexed for list access +KHARIF_MONTH_INDICES = [6, 7, 8, 9] # Jul=6, Aug=7, Sep=8, Oct=9 +KHARIF_MONTH_NAMES = ["jul", "aug", "sep", "oct"] + +# Month names for column labeling +MONTH_NAMES = [ + "jan", "feb", "mar", "apr", "may", "jun", + "jul", "aug", "sep", "oct", "nov", "dec", +] + +CRS = "EPSG:4326" +OUTPUT_PARQUET_NAME = "farm_boundaries_et.parquet" + + +# ── helpers ─────────────────────────────────────────────────────────────────── + + +def _farm_parquet_path(state: str, district: str, block: str) -> str: + return os.path.join( + FARM_BOUNDARIES_PATH, state, district, block, "farm_boundaries.parquet" + ) + + +def _output_parquet_path(state: str, district: str, block: str) -> str: + return os.path.join( + FARM_BOUNDARIES_PATH, state, district, block, OUTPUT_PARQUET_NAME + ) + + +def _output_tiff_path(state: str, district: str, block: str, year: int, raster_type: str = "aet") -> str: + return os.path.join( + FARM_BOUNDARIES_PATH, state, district, block, f"{raster_type}_{year}.tif" + ) + + +def _get_aez_zone(state: str) -> int: + """Determine the AEZ zone for a given state.""" + zone = AEZ_ZONE_MAP.get(state) + if zone is None: + raise ValueError( + f"No AEZ zone mapping for state '{state}'. " + f"Known states: {list(AEZ_ZONE_MAP.keys())}" + ) + return zone + + +# ── GEE raster download ────────────────────────────────────────────────────── + + +def _initialize_gee(gee_project: str): + """Initialize Google Earth Engine with the given project.""" + try: + ee.Initialize(project=gee_project) + logger.info("GEE initialized with project: %s", gee_project) + except Exception as exc: + raise RuntimeError( + f"Failed to initialize GEE with project '{gee_project}'. " + f"Run 'earthengine authenticate' first. Error: {exc}" + ) from exc + + +def _download_single_band(img, band_name, region, tmp_dir, band_idx): + """Download a single band from GEE as a GeoTIFF file.""" + single = img.select(band_name) + url = single.getDownloadURL({ + "scale": 30, + "crs": CRS, + "region": region, + "format": "GEO_TIFF", + }) + + response = requests.get(url, timeout=120) + response.raise_for_status() + + tmp_file = os.path.join(tmp_dir, f"band_{band_idx}.tif") + + # GEE may return a zip or raw tiff + with tempfile.NamedTemporaryFile(suffix=".zip", delete=False, dir=tmp_dir) as tmp: + tmp.write(response.content) + tmp_path = tmp.name + + if zipfile.is_zipfile(tmp_path): + with zipfile.ZipFile(tmp_path, "r") as zf: + tif_files = [f for f in zf.namelist() if f.endswith(".tif")] + if tif_files: + zf.extract(tif_files[0], tmp_dir) + extracted = os.path.join(tmp_dir, tif_files[0]) + os.rename(extracted, tmp_file) + os.remove(tmp_path) + else: + os.rename(tmp_path, tmp_file) + + return tmp_file + + +def _download_gee_raster( + asset_path: str, + tiff_path: str, + bbox: tuple, + gee_project: str, + label: str = "raster", +) -> str: + """ + Generic: download any GEE image asset to a local multi-band GeoTIFF. + + Downloads each band individually to stay under the 50MB GEE download + limit, then merges them into a single multi-band GeoTIFF locally. + + Parameters + ---------- + asset_path : str + Full GEE asset path (e.g. projects/.../aet_aez_2_2017). + tiff_path : str + Local path to save the merged multi-band GeoTIFF. + bbox : tuple + (minx, miny, maxx, maxy) in EPSG:4326. + gee_project : str + GEE cloud project ID for authentication. + label : str + Human-readable label for log messages (e.g. 'AET', 'PET'). + + Returns + ------- + str + Path to the downloaded multi-band GeoTIFF file. + """ + # Skip download if already cached + if os.path.exists(tiff_path): + logger.info("%s raster already cached at %s — skipping download.", label, tiff_path) + return tiff_path + + _initialize_gee(gee_project) + logger.info("Loading GEE asset: %s", asset_path) + + img = ee.Image(asset_path) + + minx, miny, maxx, maxy = bbox + region = ee.Geometry.Rectangle([minx, miny, maxx, maxy]) + + band_names = img.bandNames().getInfo() + logger.info( + "Downloading %s: %d bands (bbox=%.3f,%.3f,%.3f,%.3f, 30m)", + label, len(band_names), minx, miny, maxx, maxy, + ) + + os.makedirs(os.path.dirname(tiff_path), exist_ok=True) + + with tempfile.TemporaryDirectory() as tmp_dir: + band_files = [] + for idx, band_name in enumerate(band_names, 1): + logger.info(" [%s] Downloading band %d/%d: %s", label, idx, len(band_names), band_name) + band_file = _download_single_band( + img, band_name, region, tmp_dir, idx, + ) + band_files.append(band_file) + + logger.info("Merging %d bands into single GeoTIFF...", len(band_files)) + + with rasterio.open(band_files[0]) as src0: + meta = src0.meta.copy() + meta.update(count=len(band_files), dtype="float32") + + with rasterio.open(tiff_path, "w", **meta) as dst: + for band_idx, band_file in enumerate(band_files, 1): + with rasterio.open(band_file) as src: + dst.write(src.read(1).astype("float32"), band_idx) + + with rasterio.open(tiff_path) as src: + logger.info( + "%s raster saved: %s — %d bands, shape=%s, crs=%s", + label, tiff_path, src.count, src.shape, src.crs, + ) + + return tiff_path + + +def _download_aet_raster(state, district, block, year, bbox, gee_project): + """Download AET raster from GEE for the given tehsil and year.""" + aez = _get_aez_zone(state) + asset_path = AET_ASSET_TEMPLATE.format(aez=aez, year=year) + tiff_path = _output_tiff_path(state, district, block, year, "aet") + return _download_gee_raster(asset_path, tiff_path, bbox, gee_project, label="AET") + + +def _download_pet_raster(state, district, block, year, bbox, gee_project): + """Download PET raster from GEE for the given tehsil and year.""" + aez = _get_aez_zone(state) + asset_path = PET_ASSET_TEMPLATE.format(aez=aez, year=year) + tiff_path = _output_tiff_path(state, district, block, year, "pet") + return _download_gee_raster(asset_path, tiff_path, bbox, gee_project, label="PET") + + +# ── zonal statistics ────────────────────────────────────────────────────────── + + +def _rasterize_farms(gdf, transform, out_shape): + """ + Burn all farm polygons into a single labelled raster. + + Each pixel gets the 1-based index of the farm polygon it falls within. + Pixels that don't overlap any farm get 0. + + Returns + ------- + np.ndarray + Integer array of shape *out_shape* with farm labels (0 = background). + """ + shapes = ( + (geom, idx) + for idx, geom in enumerate(gdf.geometry, start=1) + ) + labels = rasterio.features.rasterize( + shapes, + out_shape=out_shape, + transform=transform, + fill=0, + dtype="int32", + all_touched=True, # include pixels that even partially overlap + ) + return labels + + +def _extract_band_means(labels, band_data, num_farms): + """ + Given a labelled raster and a band's pixel values, compute the + mean value per farm label using vectorised NumPy operations. + + Parameters + ---------- + labels : np.ndarray (int32) + Farm-label raster (0 = background, 1..N = farm index). + band_data : np.ndarray (float32) + AET pixel values for one band. + num_farms : int + Total number of farms (N). + + Returns + ------- + np.ndarray of shape (num_farms,) + Mean AET for each farm. NaN where no valid pixels exist. + """ + # Mask invalid AET values + valid = ( + ~np.isnan(band_data) & + ~np.isinf(band_data) & + (band_data > AET_NODATA) & + (band_data >= 0) + ) + + # Only work with valid pixels + valid_labels = labels[valid] + valid_values = band_data[valid] + + # Accumulate sums and counts per label using np.bincount + # Labels are 1-based (farm 0 doesn't exist), so index 0 = background + sums = np.bincount(valid_labels, weights=valid_values, minlength=num_farms + 1) + counts = np.bincount(valid_labels, minlength=num_farms + 1) + + # Compute mean (farms are at indices 1..N) + with np.errstate(invalid="ignore", divide="ignore"): + means = sums[1:] / counts[1:] + + # Farms with zero valid pixels → NaN + means[counts[1:] == 0] = np.nan + + return means + + +def _run_zonal_stats( + gdf: gpd.GeoDataFrame, aet_tiff_path: str, pet_tiff_path: str = None, +) -> gpd.GeoDataFrame: + """ + Extract per-farm monthly AET (and optionally PET) values using vectorised + rasterisation, then compute MAI and water stress indicators. + + Approach: + 1. Rasterise ALL farm polygons into a single labelled grid at the + same resolution as the raster image (30 m). + 2. For each band, compute the mean of all pixels that fall within + each farm polygon using fast NumPy bincount operations. + 3. If PET is available, compute MAI = AET / PET per month and + classify water stress using Drought Manual 2016 thresholds. + + New columns added to gdf: + aet_jan..aet_dec, aet_annual (AET monthly + annual) + pet_jan..pet_dec, pet_annual (PET monthly + annual, if available) + mai_jan..mai_dec, mai_annual (MAI = AET/PET, if PET available) + water_stress_months (count of stressed months) + kharif_water_stress (bool: does this farm have Kharif water stress?) + kharif_mai (mean MAI over Kharif months) + """ + logger.info("Running vectorised zonal statistics for %d farms...", len(gdf)) + + # ── AET processing ────────────────────────────────────────────────── + with rasterio.open(aet_tiff_path) as src: + num_bands = src.count + transform = src.transform + out_shape = (src.height, src.width) + logger.info( + "AET raster: %d bands, shape=%s, crs=%s", + num_bands, out_shape, src.crs, + ) + + # Step 1: Rasterise all farm polygons at once + logger.info( + "Rasterising %d farm polygons into a labelled grid (%d×%d)...", + len(gdf), out_shape[0], out_shape[1], + ) + labels = _rasterize_farms(gdf, transform, out_shape) + labelled_count = np.unique(labels[labels > 0]).size + logger.info( + "Labelled grid ready: %d/%d farms have at least one pixel.", + labelled_count, len(gdf), + ) + + # Step 2: Extract mean AET per farm for each band + aet_monthly_cols = [] + num_farms = len(gdf) + + for band_idx in range(1, min(num_bands, 12) + 1): + col_name = f"aet_{MONTH_NAMES[band_idx - 1]}" + logger.info(" Processing AET band %d/%d → %s", band_idx, num_bands, col_name) + + band_data = src.read(band_idx).astype("float32") + means = _extract_band_means(labels, band_data, num_farms) + + gdf[col_name] = np.round(means, 4) + aet_monthly_cols.append(col_name) + + # Band 13 = annual average (if present) + if num_bands >= 13: + logger.info(" Processing AET band 13/%d → aet_annual", num_bands) + band_data = src.read(13).astype("float32") + means = _extract_band_means(labels, band_data, num_farms) + gdf["aet_annual"] = np.round(means, 4) + else: + gdf["aet_annual"] = gdf[aet_monthly_cols].mean(axis=1).round(4) + + # ── PET processing (if available) ──────────────────────────────────── + pet_monthly_cols = [] + has_pet = pet_tiff_path is not None and os.path.exists(pet_tiff_path) + + if has_pet: + logger.info("Processing PET raster: %s", pet_tiff_path) + with rasterio.open(pet_tiff_path) as pet_src: + pet_bands = pet_src.count + pet_transform = pet_src.transform + pet_shape = (pet_src.height, pet_src.width) + + # Re-rasterise farms if PET raster has different grid + if pet_shape != out_shape or pet_transform != transform: + logger.info("PET grid differs from AET — re-rasterising farms...") + pet_labels = _rasterize_farms(gdf, pet_transform, pet_shape) + else: + pet_labels = labels + + for band_idx in range(1, min(pet_bands, 12) + 1): + col_name = f"pet_{MONTH_NAMES[band_idx - 1]}" + logger.info(" Processing PET band %d/%d → %s", band_idx, pet_bands, col_name) + + band_data = pet_src.read(band_idx).astype("float32") + means = _extract_band_means(pet_labels, band_data, num_farms) + + gdf[col_name] = np.round(means, 4) + pet_monthly_cols.append(col_name) + + if pet_bands >= 13: + logger.info(" Processing PET band 13/%d → pet_annual", pet_bands) + band_data = pet_src.read(13).astype("float32") + means = _extract_band_means(pet_labels, band_data, num_farms) + gdf["pet_annual"] = np.round(means, 4) + else: + gdf["pet_annual"] = gdf[pet_monthly_cols].mean(axis=1).round(4) + else: + logger.info("No PET raster available — skipping MAI computation.") + + # ── MAI and water stress indicators ────────────────────────────────── + if has_pet and len(pet_monthly_cols) == 12: + logger.info("Computing MAI (AET/PET) and water stress indicators...") + + # Compute monthly MAI = AET / PET + mai_monthly_cols = [] + for i, month in enumerate(MONTH_NAMES): + aet_col = f"aet_{month}" + pet_col = f"pet_{month}" + mai_col = f"mai_{month}" + + # Safe division: NaN where PET is 0 or missing + with np.errstate(invalid="ignore", divide="ignore"): + mai_values = gdf[aet_col].values / gdf[pet_col].values + # Replace inf/nan from division by zero + mai_values = np.where(np.isfinite(mai_values), mai_values, np.nan) + + gdf[mai_col] = np.round(mai_values, 4) + mai_monthly_cols.append(mai_col) + + # MAI annual = mean of monthly MAIs + gdf["mai_annual"] = gdf[mai_monthly_cols].mean(axis=1).round(4) + + # Water stress classification per month using Drought Manual 2016, Table 3.5 + # Count months in moderate or severe crop water stress (MAI ≤ 0.50) + mai_df = gdf[mai_monthly_cols] + gdf["water_stress_months"] = ( + (mai_df <= KHARIF_WATER_STRESS_MAI_THRESHOLD) & (mai_df.notna()) + ).sum(axis=1).astype(int) + + # Kharif water stress indicator + # A farm has Kharif water stress if MAI ∈ {moderate, severe} + # during any Kharif month (Jul, Aug, Sep, Oct) + # Moderate + Severe = MAI ≤ 50% (ratio ≤ 0.50) + # Note: This is crop water stress, not drought (per GoI definition) + kharif_mai_cols = [f"mai_{m}" for m in KHARIF_MONTH_NAMES] + kharif_mai_df = gdf[kharif_mai_cols] + + # Kharif water stress: any Kharif month has MAI ≤ 0.50 (moderate/severe) + gdf["kharif_water_stress"] = ( + (kharif_mai_df <= KHARIF_WATER_STRESS_MAI_THRESHOLD) & (kharif_mai_df.notna()) + ).any(axis=1) + + # Mean MAI over Kharif months (for intensity computation in multi-year) + gdf["kharif_mai"] = kharif_mai_df.mean(axis=1).round(4) + + # Summary logging + valid = gdf["mai_annual"].notna().sum() + stress_count = gdf["kharif_water_stress"].sum() + logger.info( + "MAI stats complete: %d/%d farms with valid data, " + "avg MAI=%.3f, Kharif water stress farms=%d, avg stress months=%.1f", + valid, len(gdf), + gdf["mai_annual"].mean() if valid > 0 else 0, + stress_count, + gdf["water_stress_months"].mean(), + ) + else: + # Fallback: no PET available, use simple AET threshold + STRESS_THRESHOLD = 0.5 # mm/day — placeholder + monthly_df = gdf[aet_monthly_cols] + gdf["water_stress_months"] = ( + (monthly_df < STRESS_THRESHOLD) & (monthly_df.notna()) + ).sum(axis=1).astype(int) + + valid = gdf["aet_annual"].notna().sum() + logger.info( + "Zonal stats complete (AET only, no PET): %d/%d farms with valid data, " + "avg annual AET=%.3f mm/day, avg stress months=%.1f", + valid, len(gdf), + gdf["aet_annual"].mean() if valid > 0 else 0, + gdf["water_stress_months"].mean(), + ) + + return gdf + + +# ── public entry point ──────────────────────────────────────────────────────── + + +def intersect_et_with_farms( + state: str, + district: str, + block: str, + year: int = 2017, + gee_project: str = DEFAULT_GEE_PROJECT, + overwrite: bool = False, +) -> dict: + """ + Phase 3 pipeline: intersect AET & PET raster data with farm polygons. + + Downloads the AET and PET rasters from Google Earth Engine (clipped to + the tehsil bounding box), runs zonal statistics for each farm polygon, + computes MAI (AET/PET), and saves an enhanced GeoParquet with monthly + AET, PET, MAI values and water stress indicators. + + Parameters + ---------- + state, district, block : str + Lower-cased administrative names. + year : int + Year of ET data to use (2017–2024). + gee_project : str + GEE cloud project ID for authentication. + overwrite : bool + If False and output parquet exists, skip processing. + + Returns + ------- + dict + Summary with output path and statistics. + """ + out_path = _output_parquet_path(state, district, block) + + if not overwrite and os.path.exists(out_path): + logger.info("ET parquet already exists at %s — skipping Phase 3.", out_path) + return {"path": out_path, "skipped": True} + + logger.info( + "Phase 3 — ET intersection for %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 from %s", len(gdf), farm_path) + + # 2. Download AET raster from GEE ---------------------------------------- + bbox = gdf.total_bounds # (minx, miny, maxx, maxy) + logger.info("Farm bounding box: %s", bbox) + + aet_tiff_path = _download_aet_raster( + state, district, block, year, bbox, gee_project, + ) + + # 3. Download PET raster from GEE ---------------------------------------- + try: + pet_tiff_path = _download_pet_raster( + state, district, block, year, bbox, gee_project, + ) + logger.info("PET raster downloaded successfully.") + except Exception as exc: + logger.warning( + "PET download failed (%s). Proceeding with AET only.", exc + ) + pet_tiff_path = None + + # 4. Run zonal statistics (AET + PET + MAI) ------------------------------ + gdf = _run_zonal_stats(gdf, aet_tiff_path, pet_tiff_path) + + # 5. Save enhanced parquet ----------------------------------------------- + os.makedirs(os.path.dirname(out_path), exist_ok=True) + gdf.to_parquet(out_path, index=False) + logger.info("Enhanced parquet saved → %s", out_path) + + # 6. Summary statistics -------------------------------------------------- + summary = { + "state": state, + "district": district, + "block": block, + "year": year, + "farm_count": len(gdf), + "columns": list(gdf.columns), + "avg_aet_annual": round(gdf["aet_annual"].mean(), 4) if gdf["aet_annual"].notna().any() else None, + "avg_water_stress_months": round(gdf["water_stress_months"].mean(), 1), + "has_pet": pet_tiff_path is not None, + "path": out_path, + } + + if "mai_annual" in gdf.columns: + summary["avg_mai_annual"] = round(gdf["mai_annual"].mean(), 4) if gdf["mai_annual"].notna().any() else None + summary["kharif_water_stress_farms"] = int(gdf["kharif_water_stress"].sum()) if "kharif_water_stress" in gdf.columns else 0 + + logger.info("Phase 3 complete: %s", summary) + return summary + + +# ── multi-year water stress indicators ──────────────────────────────────────── + + +def compute_multi_year_water_stress( + state: str, + district: str, + block: str, + start_year: int = 2017, + end_year: int = 2024, + gee_project: str = DEFAULT_GEE_PROJECT, +) -> dict: + """ + Compute cross-year water stress indicators (frequency & intensity) + as prescribed by the Drought Manual 2016 and approved by the professor. + Note: MAI-based classification indicates crop water stress, not drought. + Drought requires SPI/SPEI + VCI in addition to MAI (per GoI definition). + + For each year in [start_year, end_year]: + - Downloads AET & PET rasters + - Computes per-farm monthly MAI = AET / PET + - Classifies Kharif water stress years (any Kharif month MAI ≤ 0.50) + + Then across all years: + - Frequency: Return period = N / #kharif_water_stress_years (NA if 0) + - Intensity: Mean MAI over Kharif months in water stress years, + averaged across all water stress years + + Parameters + ---------- + state, district, block : str + Lower-cased administrative names. + start_year, end_year : int + Year range (inclusive) for multi-year analysis. + gee_project : str + GEE cloud project ID for authentication. + + Returns + ------- + dict + Summary with output path, indicators, and statistics. + """ + out_path = os.path.join( + FARM_BOUNDARIES_PATH, state, district, block, + "farm_boundaries_water_stress.parquet", + ) + + logger.info( + "Multi-year water stress analysis for %s/%s/%s (%d–%d)", + state, district, block, start_year, end_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." + ) + + base_gdf = gpd.read_parquet(farm_path) + num_farms = len(base_gdf) + logger.info("Loaded %d farm polygons.", num_farms) + + bbox = base_gdf.total_bounds + years = list(range(start_year, end_year + 1)) + N = len(years) + + # Per-farm accumulators: track Kharif water stress years and their MAI + kharif_stress_count = np.zeros(num_farms, dtype=int) + kharif_mai_sum_stress = np.zeros(num_farms, dtype=float) + years_processed = 0 + + # 2. Process each year --------------------------------------------------- + for year in years: + logger.info("── Processing year %d ──", year) + try: + aet_tiff = _download_aet_raster( + state, district, block, year, bbox, gee_project, + ) + pet_tiff = _download_pet_raster( + state, district, block, year, bbox, gee_project, + ) + except Exception as exc: + logger.warning("Skipping year %d — download failed: %s", year, exc) + continue + + # Run zonal stats on a copy (we only need MAI columns) + year_gdf = base_gdf.copy() + year_gdf = _run_zonal_stats(year_gdf, aet_tiff, pet_tiff) + + if "kharif_water_stress" not in year_gdf.columns: + logger.warning("Year %d: MAI not computed (PET missing?). Skipping.", year) + continue + + years_processed += 1 + + # Accumulate water stress counts and intensity values + is_drought = year_gdf["kharif_water_stress"].values.astype(bool) + kharif_mai_values = year_gdf["kharif_mai"].values + + kharif_stress_count += is_drought.astype(int) + + # Add Kharif MAI only for water stress years (for intensity averaging) + stress_mask = is_drought & np.isfinite(kharif_mai_values) + kharif_mai_sum_stress[stress_mask] += kharif_mai_values[stress_mask] + + # 3. Compute cross-year indicators --------------------------------------- + result_gdf = base_gdf.copy() + + result_gdf["total_years"] = years_processed + result_gdf["kharif_water_stress_years"] = kharif_stress_count + + # Frequency: Return period = N / #kharif_water_stress_years + with np.errstate(invalid="ignore", divide="ignore"): + return_period = years_processed / kharif_stress_count.astype(float) + return_period[kharif_stress_count == 0] = np.nan # NA for zero water stress years + result_gdf["return_period_years"] = np.round(return_period, 2) + + # Intensity: Mean MAI over Kharif months in water stress years + with np.errstate(invalid="ignore", divide="ignore"): + stress_intensity = kharif_mai_sum_stress / kharif_stress_count.astype(float) + stress_intensity[kharif_stress_count == 0] = np.nan + result_gdf["water_stress_intensity_mai"] = np.round(stress_intensity, 4) + + # 4. Save water stress parquet ------------------------------------------- + os.makedirs(os.path.dirname(out_path), exist_ok=True) + result_gdf.to_parquet(out_path, index=False) + logger.info("Multi-year water stress parquet saved → %s", out_path) + + # 5. Summary ------------------------------------------------------------- + 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_water_stress": int((kharif_stress_count > 0).sum()), + "avg_return_period": round(result_gdf["return_period_years"].mean(), 2) if valid > 0 else None, + "avg_water_stress_intensity": round(result_gdf["water_stress_intensity_mai"].mean(), 4) if valid > 0 else None, + "columns": list(result_gdf.columns), + "path": out_path, + } + logger.info("Multi-year water stress 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..22eeac46 --- /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): + """ + 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/utilities/constants.py b/utilities/constants.py index b0bb5dd4..e25025a4 100644 --- a/utilities/constants.py +++ b/utilities/constants.py @@ -7,14 +7,15 @@ NREGA_ASSETS_INPUT_DIR = "data/nrega_assets/input" NREGA_ASSETS_OUTPUT_DIR = "data/nrega_assets/output" - -ANTYODAYA_2020 = "data/antyodaya/output/pan_india_antyodaya_2020.gpkg" -LIVESTOCKS = "data/livestock/pan_india_livestock.gpkg" + +ANTYODAYA_2020 = "data/antyodaya/output/pan_india_antyodaya_2020.gpkg" +LIVESTOCKS = "data/livestock/pan_india_livestock.gpkg" MERGE_MWS_PATH = "data/merge_mws" RASTERS_PATH = "data/rasters" CROP_GRID_PATH = "data/crop_grid" +FARM_BOUNDARIES_PATH = "data/farm_boundaries" KML_PATH = "data/kml/" SHAPEFILE_DIR = "data/kml/shapefiles"