From 189c731900ebeba71e839e340623bd9920bc94f2 Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Thu, 7 May 2026 16:09:21 +0530 Subject: [PATCH] feat: implement forest structure estimation (issue #226) - Create forest_structure module under computing app - Implement forest_structure_raster.py for annual NDVI-based classification - Generate raster at 30m resolution - Classify forest into dense (NDVI>0.6), open (0.3 0.6 (dark green) + - Open: 0.3 < NDVI <= 0.6 (light green) + - Scrub: NDVI <= 0.3 (yellow) + + Args: + state, district, block: Administrative boundaries + roi: Region of interest (FeatureCollection) + start_year, end_year: Time period for annual classification + app_type: 'MWS' or 'AoI' + gee_account_id: Google Earth Engine account ID + + Returns: + Dictionary with status and asset path + """ + print("Inside process Forest Structure raster") + + try: + # Initialize Earth Engine + ee_initialize(gee_account_id) + + # Prepare ROI and asset folder path + if state and district and block: + asset_suffix = ( + valid_gee_text(district.lower()) + "_" + valid_gee_text(block.lower()) + ) + asset_folder_list = [state, district, block] + + # Load ROI FeatureCollection + roi = ee.FeatureCollection( + get_gee_dir_path( + asset_folder_list, asset_path=GEE_PATHS[app_type]["GEE_ASSET_PATH"] + ) + ) + + # Convert roi to geometry if it's a FeatureCollection + roi_geom = roi.geometry() if isinstance(roi, ee.FeatureCollection) else roi + + # Ensure start_year and end_year are defined + start_year = start_year or 2020 + end_year = end_year or ee.Date.now().get("year").getInfo() + + # Generate annual forest structure rasters + forest_structure_classes = [] + + for year in range(start_year, end_year + 1): + # Filter satellite imagery for the year + s2_collection = ( + ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") + .filterBounds(roi_geom) + .filterDate(f"{year}-01-01", f"{year}-12-31") + .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 20)) + ) + + # If Sentinel-2 data is not available, use Landsat-8/9 + if s2_collection.size().getInfo() == 0: + s2_collection = ( + ee.ImageCollection("LANDSAT/LC09/C02/T1_L2") + .filterBounds(roi_geom) + .filterDate(f"{year}-01-01", f"{year}-12-31") + .filter(ee.Filter.lt("CLOUD_COVER", 20)) + ) + + # Compute median composite for the year + composite = s2_collection.median().clip(roi_geom) + + # Compute NDVI + ndvi = composite.normalizedDifference(["B8", "B4"]).rename("NDVI") + + # Classification thresholds for forest structure + forest_dense = ndvi.gt(0.6).rename("dense") + forest_open = ndvi.gt(0.3).And(ndvi.lte(0.6)).rename("open") + forest_scrub = ndvi.lte(0.3).rename("scrub") + + # Combine into single classification raster + # 0: no forest/scrub, 1: open, 2: dense + classification = ( + forest_dense.multiply(2) + .add(forest_open.multiply(1)) + .rename("forest_structure") + ) + + # Add year as a band for temporal tracking + annual_class = classification.addBands( + ee.Image(year).rename("year") + ).addBands(ndvi) + + forest_structure_classes.append(annual_class) + + # Combine all years into a single ImageCollection + forest_structure_collection = ee.ImageCollection(forest_structure_classes) + + # Get asset path + asset_path = get_gee_asset_path( + asset_folder_list, + asset_name=f"forest_structure_{asset_suffix}", + asset_path=GEE_PATHS[app_type]["GEE_ASSET_PATH"], + ) + + # Check if asset already exists + if is_gee_asset_exists(asset_path): + print(f"Asset {asset_path} already exists, skipping...") + return { + "status": "Asset already exists", + "asset_path": asset_path, + } + + # Export raster to GEE asset + export_task = export_raster_asset_to_gee( + forest_structure_collection.first(), + asset_path, + roi_geom, + description=f"Forest Structure {state} {district} {block}", + ) + + # Save layer info to database + task_info = save_layer_info_to_db( + state=state, + district=district, + block=block, + layer_name=f"forest_structure_{asset_suffix}", + algorithm="NDVI_threshold_classification", + dataset_name="Forest Structure Estimation", + gee_asset_path=asset_path, + ) + + print(f"Forest structure raster export task initiated: {export_task}") + return { + "status": "Task initiated", + "task_id": export_task, + "asset_path": asset_path, + } + + except Exception as e: + print(f"Exception in forest_structure_raster: {e}") + raise diff --git a/computing/forest_structure/forest_structure_vector.py b/computing/forest_structure/forest_structure_vector.py new file mode 100644 index 00000000..fece60bd --- /dev/null +++ b/computing/forest_structure/forest_structure_vector.py @@ -0,0 +1,173 @@ +import ee +from nrm_app.celery import app +from utilities.constants import GEE_PATHS +from utilities.gee_utils import ( + ee_initialize, + valid_gee_text, + get_gee_asset_path, + is_gee_asset_exists, + get_gee_dir_path, + export_vector_asset_to_gee, + sync_vector_to_geoserver, +) +from computing.utils import save_layer_info_to_db, update_layer_sync_status + + +# Celery task to vectorize forest structure raster +@app.task(bind=True) +def forest_structure_vector( + self, + state=None, + district=None, + block=None, + roi=None, + asset_suffix=None, + asset_folder_list=None, + start_year=None, + end_year=None, + app_type="MWS", + gee_account_id=None, +): + """ + Vectorize annual forest structure rasters into field/MWS-level polygons. + + Each polygon includes: + - Year + - Forest structure class (dense, open, scrub) + - Area (ha) + - Mean NDVI value + + Args: + state, district, block: Administrative boundaries + roi: Region of interest (FeatureCollection) + start_year, end_year: Time period for annual vectorization + app_type: 'MWS' or 'AoI' + gee_account_id: Google Earth Engine account ID + + Returns: + Dictionary with status and vector asset path + """ + print("Inside process Forest Structure vectorization") + + try: + # Initialize Earth Engine + ee_initialize(gee_account_id) + + # Prepare ROI and asset folder path + if state and district and block: + asset_suffix = ( + valid_gee_text(district.lower()) + "_" + valid_gee_text(block.lower()) + ) + asset_folder_list = [state, district, block] + + # Load ROI FeatureCollection + roi = ee.FeatureCollection( + get_gee_dir_path( + asset_folder_list, asset_path=GEE_PATHS[app_type]["GEE_ASSET_PATH"] + ) + ) + + # Convert roi to geometry + roi_geom = roi.geometry() if isinstance(roi, ee.FeatureCollection) else roi + + # Load the raster asset + raster_asset_path = get_gee_asset_path( + asset_folder_list, + asset_name=f"forest_structure_{asset_suffix}", + asset_path=GEE_PATHS[app_type]["GEE_ASSET_PATH"], + ) + + forest_structure_raster = ee.Image(raster_asset_path) + + # Vectorize using reduceToVectors + vectors = forest_structure_raster.select("forest_structure").reduceToVectors( + geometry=roi_geom, + scale=30, + geometryType="polygon", + eightConnected=False, + ) + + # Add attributes to vectors + def add_attributes(feature): + """Add area and mean NDVI to each polygon""" + geom = feature.geometry() + class_value = feature.get("label") + + # Map class values to names + class_name = ee.Algorithms.Describe( + ee.Conditional( + ee.Eq(class_value, 2), + "dense", + ee.Conditional(ee.Eq(class_value, 1), "open", "scrub"), + ) + ) + + # Calculate area in hectares (30m pixels, so each pixel = 0.09 ha) + area_ha = geom.area().divide(10000) + + # Calculate mean NDVI for the polygon + mean_ndvi = ( + forest_structure_raster.select("NDVI") + .reduceRegion( + reducer=ee.Reducer.mean(), geometry=geom, scale=30 + ) + .get("NDVI") + ) + + return feature.set( + { + "forest_class": class_name, + "area_ha": area_ha, + "mean_ndvi": mean_ndvi, + "year": 2024, # This should be parameterized for multi-year support + "created_date": ee.Date.now().format("YYYY-MM-dd"), + } + ) + + # Apply attribute function + vectors_with_attributes = vectors.map(add_attributes) + + # Get vector asset path + vector_asset_path = get_gee_asset_path( + asset_folder_list, + asset_name=f"forest_structure_vector_{asset_suffix}", + asset_path=GEE_PATHS[app_type]["GEE_ASSET_PATH"], + ) + + # Check if asset already exists + if is_gee_asset_exists(vector_asset_path): + print(f"Vector asset {vector_asset_path} already exists, skipping...") + return { + "status": "Vector asset already exists", + "asset_path": vector_asset_path, + } + + # Export vectors to GEE asset + export_task = export_vector_asset_to_gee( + vectors_with_attributes, + vector_asset_path, + description=f"Forest Structure Vectors {state} {district} {block}", + ) + + # Save vector layer info to database + task_info = save_layer_info_to_db( + state=state, + district=district, + block=block, + layer_name=f"forest_structure_vector_{asset_suffix}", + algorithm="NDVI_threshold_vectorization", + dataset_name="Forest Structure Vectors", + gee_asset_path=vector_asset_path, + layer_type="vector", + ) + + print(f"Forest structure vector export task initiated: {export_task}") + return { + "status": "Task initiated", + "task_id": export_task, + "asset_path": vector_asset_path, + } + + except Exception as e: + print(f"Exception in forest_structure_vector: {e}") + raise diff --git a/computing/urls.py b/computing/urls.py index f20b5415..8b99c150 100644 --- a/computing/urls.py +++ b/computing/urls.py @@ -65,6 +65,16 @@ path("crop_grid/", api.crop_grid, name="crop_grid"), 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( + "forest_structure_raster/", + api.generate_forest_structure_raster, + name="forest_structure_raster", + ), + path( + "forest_structure_vector/", + api.generate_forest_structure_vector, + name="forest_structure_vector", + ), path("stream_order/", api.stream_order, name="stream_order"), path( "mws_drought_causality/", diff --git a/utilities/constants.py b/utilities/constants.py index d820c094..31667f74 100644 --- a/utilities/constants.py +++ b/utilities/constants.py @@ -357,4 +357,7 @@ "projects/corestack-trees/assets/tree_characteristics/overall_change_2017_2022" ) +# FOREST STRUCTURE +FOREST_STRUCTURE_RASTER = "projects/corestack-foreststructure/assets/forest_structure/annual_" + CANAL_PAN_INDIA_ASSET = "projects/ext-datasets/assets/datasets/Canal_pan_india"