diff --git a/computing/tree_health/forest_structure.py b/computing/tree_health/forest_structure.py new file mode 100644 index 00000000..5943a652 --- /dev/null +++ b/computing/tree_health/forest_structure.py @@ -0,0 +1,244 @@ +""" +Forest Structure Estimation (Field Level @30m) + +Classifies forest cover into dense, open, and scrub using NDVI and +texture metrics from Sentinel-2 / Landsat imagery. Produces annual +raster and vector outputs at 30m resolution. + +Classification thresholds (based on FSI conventions): + - Dense forest: NDVI > 0.6 + - Open forest: 0.3 < NDVI <= 0.6 + - Scrub: NDVI <= 0.3 + +Reference: + - https://www.mdpi.com/2072-4292/13/24/5105 + - https://www.cse.iitd.ernet.in/~aseth/forest-health-ictd2024.pdf +""" + +import ee +from nrm_app.celery import app +from utilities.gee_utils import ( + ee_initialize, + valid_gee_text, + get_gee_asset_path, + is_gee_asset_exists, + check_task_status, + export_raster_asset_to_gee, + export_vector_asset_to_gee, + make_asset_public, + get_gee_dir_path, +) +from utilities.constants import GEE_PATHS +from computing.utils import ( + sync_fc_to_geoserver, + save_layer_info_to_db, + update_layer_sync_status, + get_layer_object, +) +from computing.STAC_specs import generate_STAC_layerwise + +# NDVI thresholds for forest structure classification +DENSE_THRESHOLD = 0.6 +OPEN_THRESHOLD = 0.3 + +# Class values in raster output +CLASS_DENSE = 3 +CLASS_OPEN = 2 +CLASS_SCRUB = 1 +CLASS_NON_FOREST = 0 + + +def _get_cloud_masked_composite(year, roi): + """Build a cloud-free annual composite from Sentinel-2 + Landsat.""" + start = f"{year}-01-01" + end = f"{year}-12-31" + + # Sentinel-2 SR + s2 = ( + ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") + .filterDate(start, end) + .filterBounds(roi) + .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 30)) + ) + + def mask_s2_clouds(image): + qa = image.select("QA60") + cloud_mask = qa.bitwiseAnd(1 << 10).eq(0).And(qa.bitwiseAnd(1 << 11).eq(0)) + return image.updateMask(cloud_mask) + + s2_masked = s2.map(mask_s2_clouds) + + # Compute NDVI from Sentinel-2 bands + def add_ndvi_s2(image): + ndvi = image.normalizedDifference(["B8", "B4"]).rename("ndvi") + return image.addBands(ndvi) + + s2_ndvi = s2_masked.map(add_ndvi_s2) + + # Fallback to Landsat 8/9 if Sentinel-2 coverage is sparse + l8 = ( + ee.ImageCollection("LANDSAT/LC08/C02/T1_L2") + .filterDate(start, end) + .filterBounds(roi) + .filter(ee.Filter.lt("CLOUD_COVER", 30)) + ) + l9 = ( + ee.ImageCollection("LANDSAT/LC09/C02/T1_L2") + .filterDate(start, end) + .filterBounds(roi) + .filter(ee.Filter.lt("CLOUD_COVER", 30)) + ) + + def add_ndvi_landsat(image): + ndvi = image.normalizedDifference(["SR_B5", "SR_B4"]).rename("ndvi") + return image.addBands(ndvi) + + landsat = l8.merge(l9).map(add_ndvi_landsat) + + # Merge and take median NDVI + combined = s2_ndvi.select("ndvi").merge(landsat.select("ndvi")) + median_ndvi = combined.median().clip(roi) + + return median_ndvi + + +def _classify_forest_structure(ndvi_image): + """Classify NDVI image into forest structure classes.""" + dense = ndvi_image.gt(DENSE_THRESHOLD).multiply(CLASS_DENSE) + open_forest = ( + ndvi_image.gt(OPEN_THRESHOLD) + .And(ndvi_image.lte(DENSE_THRESHOLD)) + .multiply(CLASS_OPEN) + ) + scrub = ndvi_image.lte(OPEN_THRESHOLD).And(ndvi_image.gt(0.1)).multiply(CLASS_SCRUB) + + classified = dense.add(open_forest).add(scrub).rename("forest_structure") + return classified.toUint8() + + +def _vectorize_structure(classified_image, roi, scale=30): + """Convert classified raster to vector polygons with attributes.""" + vectors = classified_image.reduceToVectors( + geometry=roi, + scale=scale, + geometryType="polygon", + eightConnected=True, + labelProperty="structure_class", + maxPixels=1e10, + ) + + # Map class numbers to labels + def add_label(feature): + cls = feature.get("structure_class") + label = ( + ee.Algorithms.If(ee.Number(cls).eq(CLASS_DENSE), "dense", + ee.Algorithms.If(ee.Number(cls).eq(CLASS_OPEN), "open", + ee.Algorithms.If(ee.Number(cls).eq(CLASS_SCRUB), "scrub", "non_forest"))) + ) + area = feature.geometry().area().divide(10000) # hectares + return feature.set({"structure_label": label, "area_ha": area}) + + return vectors.map(add_label) + + +@app.task(bind=True) +def compute_forest_structure( + self, + state=None, + district=None, + block=None, + year=2024, + gee_account_id=None, + roi_path=None, + asset_folder_list=None, + asset_suffix=None, + app_type="MWS", +): + """Compute annual forest structure map for a tehsil/AoI. + + Produces both raster (30m classification) and vector (MWS-level polygons) + outputs as GEE assets. + """ + ee_initialize(gee_account_id) + + if state and district and block: + asset_suffix = ( + valid_gee_text(district.lower()) + "_" + valid_gee_text(block.lower()) + ) + asset_folder_list = [state, district, block] + roi_path = ( + get_gee_dir_path( + asset_folder_list, asset_path=GEE_PATHS[app_type]["GEE_ASSET_PATH"] + ) + + f"filtered_mws_{valid_gee_text(district.lower())}_{valid_gee_text(block.lower())}_uid" + ) + + roi = ee.FeatureCollection(roi_path) + roi_geom = roi.geometry() + + raster_name = f"forest_structure_{year}_{asset_suffix}" + vector_name = f"forest_structure_vec_{year}_{asset_suffix}" + + gee_dir = get_gee_dir_path( + asset_folder_list, asset_path=GEE_PATHS[app_type]["GEE_ASSET_PATH"] + ) + raster_asset_id = gee_dir + raster_name + vector_asset_id = gee_dir + vector_name + + # Step 1: compute NDVI composite + print(f"Computing NDVI composite for {year}...") + ndvi = _get_cloud_masked_composite(year, roi_geom) + + # Step 2: classify + print("Classifying forest structure...") + classified = _classify_forest_structure(ndvi) + + # Step 3: export raster + if not is_gee_asset_exists(raster_asset_id): + print(f"Exporting raster: {raster_asset_id}") + task_id = export_raster_asset_to_gee( + classified, raster_name, raster_asset_id, scale=30, region=roi_geom + ) + if task_id: + check_task_status(task_id) + make_asset_public(raster_asset_id) + else: + print(f"Raster already exists: {raster_asset_id}") + + # Step 4: vectorize and export + if not is_gee_asset_exists(vector_asset_id): + print("Vectorizing...") + vectors = _vectorize_structure(classified, roi_geom) + + # Add mean NDVI per polygon + vectors = ndvi.reduceRegions( + collection=vectors, reducer=ee.Reducer.mean(), scale=30 + ).map(lambda f: f.set("mean_ndvi", f.get("mean"))) + + print(f"Exporting vector: {vector_asset_id}") + task_id = export_vector_asset_to_gee(vectors, vector_name, vector_asset_id) + if task_id: + check_task_status(task_id) + make_asset_public(vector_asset_id) + else: + print(f"Vector already exists: {vector_asset_id}") + + # Step 5: sync and save + layer_name = f"{asset_suffix}_forest_structure_{year}" + sync_fc_to_geoserver(vector_asset_id, layer_name) + save_layer_info_to_db( + state=state, + district=district, + block=block, + layer_name=layer_name, + dataset_name="Forest Structure", + metadata={ + "year": year, + "resolution": "30m", + "classes": {"3": "dense", "2": "open", "1": "scrub", "0": "non_forest"}, + "thresholds": {"dense": f"NDVI>{DENSE_THRESHOLD}", "open": f"{OPEN_THRESHOLD}