SustainabilityLabIITGN/NDVI_PERG_API
0
1import os2from datetime import datetime3import ee4import json5import numpy as np6import geemap.foliumap as gee_folium7import leafmap.foliumap as leaf_folium8import gradio as gr9import pandas as pd10import geopandas as gpd11import plotly.express as px12import branca.colormap as cm13from shapely.ops import transform14import pyproj15from io import BytesIO16import requests17import kml2geojson18import folium19import xml.etree.ElementTree as ET20from fastapi import FastAPI, HTTPException, Response21from urllib.parse import unquote22from pydantic import BaseModel, HttpUrl23 24 25app = FastAPI()26 27# --- Helper Functions ---28 29def one_time_setup():30 """Initializes the Earth Engine API."""31 try:32 # Attempt to initialize with default credentials33 ee.Initialize()34 except Exception:35 try:36 # Fallback to service account credentials if default init fails37 credentials_path = os.path.expanduser("~/.config/earthengine/credentials.json")38 ee_credentials = os.environ.get("EE")39 if ee_credentials:40 os.makedirs(os.path.dirname(credentials_path), exist_ok=True)41 with open(credentials_path, "w") as f:42 f.write(ee_credentials)43 credentials = ee.ServiceAccountCredentials('ujjwal@ee-ujjwaliitd.iam.gserviceaccount.com', credentials_path)44 ee.Initialize(credentials, project='ee-ujjwaliitd')45 except Exception as inner_e:46 # If the fallback also fails, print the error47 print(f"Earth Engine initialization failed: {inner_e}")48 49 50def _process_spatial_data(data_bytes: BytesIO) -> gpd.GeoDataFrame:51 """Core function to process bytes of a KML or GeoJSON file."""52 # Read the first few bytes to determine file type without consuming the stream53 start_of_file = data_bytes.read(100)54 data_bytes.seek(0) # Reset stream position55 56 # Check if the file is KML (XML-based)57 if start_of_file.strip().lower().startswith(b'<?xml'):58 try:59 geojson_data = kml2geojson.convert(data_bytes)60 if not geojson_data or not geojson_data[0].get("features"):61 raise ValueError("KML file is empty or has no features.")62 features = geojson_data[0]["features"]63 input_gdf = gpd.GeoDataFrame.from_features(features, crs="EPSG:4326")64 except Exception as e:65 raise ValueError(f"Failed to process KML data: {e}")66 # Otherwise, assume it's a format geopandas can read from databytes67 else:68 try:69 geojson_str = data_bytes.read().decode('utf-8')70 input_gdf = gpd.read_file(geojson_str)71 except Exception as e:72 raise ValueError(f"Failed to read GeoJSON or other vector data: {e}")73 return input_gdf74 75def get_gdf_from_file(file_obj):76 """Reads a KML or GeoJSON file from a Gradio file object and returns a GeoDataFrame."""77 if file_obj is None:78 return None79 with open(file_obj.name, 'rb') as f:80 data_bytes = BytesIO(f.read())81 return _process_spatial_data(data_bytes)82 83def get_gdf_from_url(url: str) -> gpd.GeoDataFrame:84 """Downloads and reads a KML/GeoJSON from a URL."""85 if not url or not url.strip():86 return None87 88 # Handle Google Drive URLs89 if "drive.google.com" in url:90 if "/file/d/" in url:91 file_id = url.split('/d/')[1].split('/')[0]92 elif "open?id=" in url:93 file_id = url.split('open?id=')[1].split('&')[0]94 else:95 raise ValueError("Unsupported Google Drive URL format. Please provide a direct link or a shareable link with 'open?id=' or '/file/d/'.")96 download_url = f"https://drive.google.com/uc?export=download&id={file_id}"97 else:98 download_url = url99 100 try:101 response = requests.get(download_url, timeout=30)102 response.raise_for_status()103 data_bytes = BytesIO(response.content)104 return _process_spatial_data(data_bytes)105 except requests.exceptions.RequestException as e:106 raise ValueError(f"Failed to download file from URL: {e}")107 108 109def find_best_epsg(geometry) -> int:110 """Finds the most suitable EPSG code for a given geometry based on its centroid."""111 if geometry.geom_type == "Polygon":112 centroid = geometry.centroid113 else:114 raise ValueError("Geometry is not a Polygon.")115 116 common_epsg_codes = [117 7761, # Gujarat118 7774, # Rajasthan119 7766, # MadhyaPradesh120 7767, # Maharastra121 7755, # India122 # Add other relevant state/country EPSG codes here123 ]124 125 for epsg in common_epsg_codes:126 try:127 crs = pyproj.CRS.from_epsg(epsg)128 area_of_use = crs.area_of_use.bounds129 if (area_of_use[0] <= centroid.x <= area_of_use[2]) and \130 (area_of_use[1] <= centroid.y <= area_of_use[3]):131 return epsg132 except pyproj.exceptions.CRSError:133 continue134 return 4326 # Default to WGS84 if no suitable projection is found135 136def shape_3d_to_2d(shape):137 """Converts a 3D geometry to 2D."""138 if shape.has_z:139 return transform(lambda x, y, z: (x, y), shape)140 return shape141 142def preprocess_gdf(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:143 """Preprocesses a GeoDataFrame by converting geometries to 2D and fixing invalid ones."""144 gdf["geometry"] = gdf["geometry"].apply(shape_3d_to_2d)145 gdf["geometry"] = gdf.buffer(0)146 return gdf147 148def to_best_crs(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:149 """Converts a GeoDataFrame to the most suitable CRS."""150 if not gdf.empty and gdf["geometry"].iloc[0] is not None:151 best_epsg_code = find_best_epsg(gdf.to_crs(epsg=4326)["geometry"].iloc[0])152 return gdf.to_crs(epsg=best_epsg_code)153 return gdf154 155def is_valid_polygon(geometry_gdf):156 """Checks if the geometry in a GeoDataFrame is a valid, non-empty Polygon."""157 if geometry_gdf.empty:158 return False159 geometry = geometry_gdf.geometry.item()160 return (geometry.type == 'Polygon') and (not geometry.is_empty)161 162def add_geometry_to_map(m, geometry_gdf, buffer_geometry_gdf, opacity=0.3):163 """Adds geometry and its buffer to a folium map."""164 if buffer_geometry_gdf is not None and not buffer_geometry_gdf.empty:165 folium.GeoJson(166 buffer_geometry_gdf.to_crs(epsg=4326),167 name="Geometry Buffer",168 style_function=lambda x: {"color": "red", "fillOpacity": opacity, "fillColor": "red"}169 ).add_to(m)170 if geometry_gdf is not None and not geometry_gdf.empty:171 folium.GeoJson(172 geometry_gdf.to_crs(epsg=4326),173 name="Geometry",174 style_function=lambda x: {"color": "blue", "fillOpacity": opacity, "fillColor": "blue"}175 ).add_to(m)176 return m177 178def get_wayback_data():179 """Fetches and parses Wayback imagery data from ArcGIS."""180 try:181 url = "https://wayback.maptiles.arcgis.com/arcgis/rest/services/World_Imagery/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"182 response = requests.get(url)183 response.raise_for_status() # Ensure request was successful184 185 # Parse XML186 root = ET.fromstring(response.content)187 188 ns = {189 "wmts": "https://www.opengis.net/wmts/1.0",190 "ows": "https://www.opengis.net/ows/1.1",191 "xlink": "https://www.w3.org/1999/xlink",192 }193 194 # Use a robust XPath to find all 'Layer' elements anywhere in the document.195 # This is less brittle than specifying the full path.196 layers = root.findall(".//wmts:Contents/wmts:Layer", ns)197 198 layer_data = []199 for layer in layers:200 title = layer.find("ows:Title", ns)201 identifier = layer.find("ows:Identifier", ns)202 resource = layer.find("wmts:ResourceURL", ns) # Tile URL template203 204 title_text = title.text if title is not None else "N/A"205 identifier_text = identifier.text if identifier is not None else "N/A"206 url_template = resource.get("template") if resource is not None else "N/A"207 208 layer_data.append({"Title": title_text, "ResourceURL_Template": url_template})209 210 wayback_df = pd.DataFrame(layer_data)211 wayback_df["date"] = pd.to_datetime(wayback_df["Title"].str.extract(r"(\d{4}-\d{2}-\d{2})").squeeze(), errors="coerce")212 wayback_df.set_index("date", inplace=True)213 return wayback_df.sort_index(ascending=False)214 215 except Exception as e:216 print(f"Could not fetch or parse Wayback data: {e}")217 return pd.DataFrame()218 219 220def get_dem_slope_maps(ee_geometry, map_bounds, wayback_url, wayback_title, zoom=12,):221 """Creates DEM and Slope maps from SRTM data, using wayback tiles as a basemap if available."""222 223 print(wayback_url, wayback_title)224 225 # --- DEM Map ---226 dem_map = gee_folium.Map(zoom_start=zoom)227 if wayback_url:228 dem_map.add_tile_layer(url=wayback_url, name=wayback_title, attribution="Esri")229 230 dem_map_html = "<div>No DEM data available for this area.</div>"231 try:232 dem_layer = ee.Image("USGS/SRTMGL1_003").resample("bilinear").reproject(crs="EPSG:4326", scale=30).clip(ee_geometry)233 stats = dem_layer.reduceRegion(reducer=ee.Reducer.minMax(), geometry=ee_geometry, scale=30, maxPixels=1e9).getInfo()234 print(stats)235 236 if stats and stats.get('elevation_min') is not None:237 min_val, max_val = stats['elevation_min'], stats['elevation_max']238 vis_params = {"min": min_val, "max": max_val, "palette": ['#0000FF', '#00FF00', '#FFFF00', '#FF0000']}239 dem_map.addLayer(dem_layer, vis_params, "Elevation")240 dem_map.add_colorbar(vis_params=vis_params, label="Elevation (m)")241 242 dem_map.addLayerControl()243 dem_map.fit_bounds(map_bounds, padding=(10, 10))244 dem_map_html = dem_map._repr_html_()245 246 except Exception as e:247 print(f"Error creating DEM map: {e}")248 dem_map_html = f"<div>Error creating DEM map: {e}</div>"249 250 # ---Slope Map --- #251 slope_map = gee_folium.Map(zoom_start=zoom)252 if wayback_url:253 slope_map.add_tile_layer(url=wayback_url, name=wayback_title, attribution="Esri")254 255 slope_map_html = "<div>No Slope data available for this area.</div>"256 try:257 dem_for_slope = ee.Image("USGS/SRTMGL1_003")258 # Calculate slope. The result is an image with slope values in degrees.259 slope_layer = ee.Terrain.slope(dem_for_slope).clip(ee_geometry)260 261 stats = slope_layer.reduceRegion(reducer=ee.Reducer.minMax(), geometry=ee_geometry, scale=30, maxPixels=1e9).getInfo()262 print(stats)263 264 if stats and stats.get('slope_min') is not None:265 min_val, max_val = stats['slope_min'], stats['slope_max']266 vis_params = {"min": min_val, "max": max_val, "palette": ['#0000FF', '#00FF00', '#FFFF00', '#FF0000']}267 slope_map.addLayer(slope_layer, vis_params, "Slope")268 slope_map.add_colorbar(vis_params=vis_params, label="Slope (degrees)")269 270 slope_map.addLayerControl()271 slope_map.fit_bounds(map_bounds, padding=(10, 10))272 slope_map_html = slope_map._repr_html_()273 274 except Exception as e:275 print(f"Error creating Slope map: {e}")276 slope_map_html = f"<div>Error creating Slope map: {e}</div>"277 278 return dem_map_html, slope_map_html279 280def add_indices(image, nir_band, red_band, blue_band, green_band, swir_band, swir2_band, evi_vars):281 """Calculates and adds multiple vegetation indices to an Earth Engine image."""282 nir = image.select(nir_band).divide(10000)283 red = image.select(red_band).divide(10000)284 blue = image.select(blue_band).divide(10000)285 green = image.select(green_band).divide(10000)286 swir = image.select(swir_band).divide(10000)287 swir2 = image.select(swir2_band).divide(10000)288 289 # Previously existing indices290 ndvi = image.normalizedDifference([nir_band, red_band]).rename('NDVI')291 evi = image.expression(292 'G * ((NIR - RED) / (NIR + C1 * RED - C2 * BLUE + L))', {293 'NIR': nir, 'RED': red, 'BLUE': blue,294 'G': evi_vars['G'], 'C1': evi_vars['C1'], 'C2': evi_vars['C2'], 'L': evi_vars['L']295 }).rename('EVI')296 evi2 = image.expression(297 'G * (NIR - RED) / (NIR + L + C * RED)', {298 'NIR': nir, 'RED': red,299 'G': evi_vars['G'], 'L': evi_vars['L'], 'C': evi_vars['C']300 }).rename('EVI2')301 try:302 table = ee.FeatureCollection('projects/in793-aq-nb-24330048/assets/cleanedVDI').select(303 ["B2", "B4", "B8", "cVDI"], ["Blue", "Red", "NIR", 'cVDI'])304 classifier = ee.Classifier.smileRandomForest(50).train(305 features=table, classProperty='cVDI', inputProperties=['Blue', 'Red', 'NIR'])306 rf = image.classify(classifier).multiply(ee.Number(0.2)).add(ee.Number(0.1)).rename('RandomForest')307 except Exception as e:308 print(f"Random Forest calculation failed: {e}")309 rf = ee.Image.constant(0).rename('RandomForest')310 ci = image.expression(311 '(-3.98 * (BLUE/NIR) + 12.54 * (GREEN/NIR) - 5.49 * (RED/NIR) - 0.19) / ' +312 '(-21.87 * (BLUE/NIR) + 12.4 * (GREEN/NIR) + 19.98 * (RED/NIR) + 1) * 2.29', {313 'NIR': nir, 'RED': red, 'BLUE': blue, 'GREEN': green314 }).clamp(0, 1).rename('CI')315 gujvdi = image.expression(316 '0.5 * (NIR - RED) / (NIR + 6 * RED - 8.25 * BLUE - 0.01)', {317 'NIR': nir, 'RED': red, 'BLUE': blue318 }).rename('GujVDI')319 mndwi = image.normalizedDifference([green_band, swir_band]).rename('MNDWI')320 321 # Newly added indices322 savi = image.expression('(1 + L) * (NIR - RED) / (NIR + RED + L)', {323 'NIR': nir, 'RED': red, 'L': 0.5324 }).rename('SAVI')325 mvi = image.expression('(NIR - (GREEN + SWIR)) / (NIR + (GREEN + SWIR))', {326 'NIR': nir, 'GREEN': green, 'SWIR': swir327 }).rename('MVI')328 nbr = image.normalizedDifference([nir_band, swir2_band]).rename('NBR')329 gci = image.expression('(NIR - GREEN) / GREEN', {330 'NIR': nir, 'GREEN': green331 }).rename('GCI')332 333 return image.addBands([ndvi, evi, evi2, rf, ci, gujvdi, mndwi, savi, mvi, nbr, gci])334 335 336# --- Gradio App Logic ---337 338# Initialize GEE and fetch wayback data once at the start339one_time_setup()340WAYBACK_DF = get_wayback_data()341 342def process_and_display(file_obj, url_str, buffer_m, progress=gr.Progress()):343 """Main function to process the uploaded file or URL and generate initial outputs."""344 if file_obj is None and not (url_str and url_str.strip()):345 return None, "Please upload a file or provide a URL.", None, None, None, None, None346 347 348 progress(0, desc="Reading and processing geometry...")349 try:350 input_gdf = get_gdf_from_file(file_obj) if file_obj is not None else get_gdf_from_url(url_str)351 input_gdf = preprocess_gdf(input_gdf)352 geometry_gdf = next((input_gdf.iloc[[i]] for i in range(len(input_gdf)) if is_valid_polygon(input_gdf.iloc[[i]])), None)353 if geometry_gdf is None:354 return None, "No valid polygon found in the provided file.", None, None, None, None, None355 geometry_gdf = to_best_crs(geometry_gdf)356 outer_geometry_gdf = geometry_gdf.copy()357 outer_geometry_gdf["geometry"] = outer_geometry_gdf["geometry"].buffer(buffer_m)358 buffer_geometry_gdf = gpd.GeoDataFrame(359 geometry=[outer_geometry_gdf.unary_union.difference(geometry_gdf.unary_union)],360 crs=geometry_gdf.crs361 )362 except Exception as e:363 return None, f"Error processing file: {e}", None, None, None, None, None364 365 progress(0.5, desc="Generating maps and stats...")366 m = folium.Map()367 wayback_url = None368 wayback_title = "Esri Satellite"369 if not WAYBACK_DF.empty:370 latest_item = WAYBACK_DF.iloc[0]371 wayback_title = f"Esri Wayback ({latest_item.name.strftime('%Y-%m-%d')})"372 wayback_url = (373 latest_item["ResourceURL_Template"]374 .replace("{TileMatrixSet}", "GoogleMapsCompatible")375 .replace("{TileMatrix}", "{z}")376 .replace("{TileRow}", "{y}")377 .replace("{TileCol}", "{x}")378 )379 folium.TileLayer(tiles=wayback_url, attr="Esri", name=wayback_title).add_to(m)380 381 m = add_geometry_to_map(m, geometry_gdf, buffer_geometry_gdf, opacity=0.3)382 383 bounds = geometry_gdf.to_crs(epsg=4326).total_bounds384 map_bounds = [[bounds[1], bounds[0]], [bounds[3], bounds[2]]]385 m.fit_bounds(map_bounds, padding=(10, 10))386 folium.LayerControl().add_to(m)387 388 ee_geometry = ee.Geometry(json.loads(geometry_gdf.to_crs(4326).to_json())['features'][0]['geometry'])389 dem_html, slope_html = get_dem_slope_maps(ee_geometry, map_bounds, wayback_url=wayback_url, wayback_title=wayback_title)390 391 stats_df = pd.DataFrame({392 "Area (ha)": [f"{geometry_gdf.area.item() / 10000:.2f}"],393 "Perimeter (m)": [f"{geometry_gdf.length.item():.2f}"],394 "Centroid (Lat, Lon)": [f"({geometry_gdf.to_crs(4326).centroid.y.iloc[0]:.6f}, {geometry_gdf.to_crs(4326).centroid.x.iloc[0]:.6f})"]395 })396 geometry_json = geometry_gdf.to_json()397 buffer_geometry_json = buffer_geometry_gdf.to_json()398 progress(1, desc="Done!")399 return m._repr_html_(), None, stats_df, dem_html, slope_html, geometry_json, buffer_geometry_json400 401@app.get("/api/geometry")402def calculate_geometry_metrics(file_url: str):403 """404 Accepts a URL to a KML/GeoJSON file, calculates the area and405 perimeter of the first valid polygon, and returns the results406 in a CSV format compatible with Google Sheets' IMPORTDATA.407 """408 try:409 decoded_url = unquote(file_url)410 input_gdf = get_gdf_from_url(decoded_url)411 412 if input_gdf is None or input_gdf.empty:413 raise ValueError("Could not read geometry from the provided URL.")414 415 input_gdf = preprocess_gdf(input_gdf)416 geometry_gdf = next((input_gdf.iloc[[i]] for i in range(len(input_gdf)) if is_valid_polygon(input_gdf.iloc[[i]])), None)417 418 if geometry_gdf is None:419 raise ValueError("No valid polygon found in the provided file.")420 421 projected_gdf = to_best_crs(geometry_gdf)422 area_hectares = projected_gdf.area.item() / 10000423 perimeter_meters = projected_gdf.length.item()424 425 centroid_gdf = projected_gdf.to_crs(epsg=4326)426 centroid_point = centroid_gdf.centroid.item()427 428 data_row = (429 f"area_hectares, {round(area_hectares, 4)},"430 f"perimeter_meters, {round(perimeter_meters, 4)},"431 f"latitude, {round(centroid_point.y, 4)},"432 f"longitude, {round(centroid_point.x, 4)}"433 )434 csv_output = f"{data_row}"435 return Response(content=csv_output, media_type="text/csv")436 437 except ValueError as e:438 # Handle specific errors with a 400 Bad Request439 raise HTTPException(status_code=400, detail=str(e))440 except Exception as e:441 # Handle any other unexpected errors with a 500442 print(f"An unexpected error occurred in /api/geometry: {e}")443 raise HTTPException(status_code=500, detail="An unexpected server error occurred.")444 445 446def calculate_indices(447 geometry_json, buffer_geometry_json, veg_indices, evi_vars, date_range,448 min_year, max_year, progress=gr.Progress()449):450 """Calculates vegetation indices based on user inputs."""451 one_time_setup()452 453 if not all([geometry_json, buffer_geometry_json, veg_indices]):454 return "Please process a file and select at least one index first.", None, None, None455 456 try:457 geometry_gdf = gpd.read_file(geometry_json)458 buffer_geometry_gdf = gpd.read_file(buffer_geometry_json)459 ee_geometry = ee.Geometry(json.loads(geometry_gdf.to_crs(4326).to_json())['features'][0]['geometry'])460 buffer_ee_geometry = ee.Geometry(json.loads(buffer_geometry_gdf.to_crs(4326).to_json())['features'][0]['geometry'])461 462 start_day, start_month = date_range[0].day, date_range[0].month463 end_day, end_month = date_range[1].day, date_range[1].month464 dates = [465 (f"{year}-{start_month:02d}-{start_day:02d}", f"{year}-{end_month:02d}-{end_day:02d}")466 for year in range(min_year, max_year + 1)467 ]468 469 collection = (470 ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")471 .select(472 ["B2", "B3", "B4", "B8", "B11", "B12", "MSK_CLDPRB"],473 ["Blue", "Green", "Red", "NIR", "SWIR", "SWIR2", "MSK_CLDPRB"]474 )475 .map(lambda img: add_indices(img, 'NIR', 'Red', 'Blue', 'Green', 'SWIR', 'SWIR2', evi_vars))476 )477 478 result_rows = []479 for i, (start_date, end_date) in enumerate(dates):480 progress((i + 1) / len(dates), desc=f"Processing {start_date} to {end_date}")481 filtered_collection = collection.filterDate(start_date, end_date).filterBounds(ee_geometry)482 if filtered_collection.size().getInfo() == 0:483 continue484 485 year_val = int(start_date.split('-')[0])486 row = {'Year': year_val, 'Date Range': f"{start_date}_to_{end_date}"}487 488 for veg_index in veg_indices:489 mosaic = filtered_collection.qualityMosaic(veg_index)490 mean_val = mosaic.reduceRegion(reducer=ee.Reducer.mean(), geometry=ee_geometry, scale=10, maxPixels=1e9).get(veg_index).getInfo()491 buffer_mean_val = mosaic.reduceRegion(reducer=ee.Reducer.mean(), geometry=buffer_ee_geometry, scale=10, maxPixels=1e9).get(veg_index).getInfo()492 493 row[veg_index] = mean_val494 row[f"{veg_index}_buffer"] = buffer_mean_val495 row[f"{veg_index}_ratio"] = (mean_val / buffer_mean_val) if buffer_mean_val and buffer_mean_val != 0 else np.nan496 result_rows.append(row)497 498 if not result_rows:499 return "No satellite imagery found for the selected dates.", None, None, None500 501 result_df = pd.DataFrame(result_rows)502 result_df = result_df.round(3)503 504 plots = []505 for veg_index in veg_indices:506 plot_cols = [veg_index, f"{veg_index}_buffer", f"{veg_index}_ratio"]507 existing_plot_cols = [col for col in plot_cols if col in result_df.columns]508 509 plot_df = result_df[['Year'] + existing_plot_cols].dropna()510 511 if not plot_df.empty:512 fig = px.line(plot_df, x='Year', y=existing_plot_cols, markers=True, title=f"{veg_index} Time Series")513 fig.update_layout(xaxis_title="Year", yaxis_title="Index Value")514 fig.update_xaxes(dtick=1)515 plots.append(fig)516 517 return None, result_df, plots, "Calculation complete."518 519 except Exception as e:520 import traceback521 traceback.print_exc()522 return f"An error occurred during calculation: {e}", None, None, None523 524def get_histogram(index_name, image, geometry, bins):525 """Calculates the histogram for an image within a given geometry using GEE."""526 try:527 # Request histogram data from Earth Engine528 hist_info = image.reduceRegion(529 reducer=ee.Reducer.fixedHistogram(min=bins[0], max=bins[-1], steps=len(bins)-1),530 geometry=geometry,531 scale=10, # Scale in meters appropriate for Sentinel-2532 maxPixels=1e9533 ).get(index_name).getInfo()534 535 # Extract histogram counts536 if hist_info:537 histogram = [item[1] for item in hist_info]538 return np.array(histogram), bins539 else:540 # Return empty histogram if no data541 return np.array([0] * (len(bins) - 1)), bins542 except Exception as e:543 print(f"Could not compute histogram for {index_name}: {e}")544 return np.array([0] * (len(bins) - 1)), bins545 546def generate_comparison_maps(geometry_json, selected_index, selected_years, evi_vars, date_start_str, date_end_str, progress=gr.Progress()):547 """Generates side-by-side maps for a selected index and two selected years with a custom HTML legend."""548 if not geometry_json or not selected_index or not selected_years:549 return "Please process a file and select an index and years first.", "", ""550 if len(selected_years) != 2:551 return "Please select exactly two years to compare.", "", ""552 553 one_time_setup()554 geometry_gdf = gpd.read_file(geometry_json).to_crs(4326)555 ee_geometry = ee.Geometry(json.loads(geometry_gdf.to_json())['features'][0]['geometry'])556 bounds = geometry_gdf.total_bounds557 map_bounds = [[bounds[1], bounds[0]], [bounds[3], bounds[2]]]558 559 start_month, start_day = map(int, date_start_str.split('-'))560 end_month, end_day = map(int, date_end_str.split('-'))561 562 maps_html = []563 for i, year in enumerate(selected_years):564 progress((i + 1) / 2, desc=f"Generating map for {year}")565 start_date = f"{year}-{start_month:02d}-{start_day:02d}"566 end_date = f"{year}-{end_month:02d}-{end_day:02d}"567 568 wayback_url = None569 wayback_title = "Default Satellite"570 if not WAYBACK_DF.empty:571 try:572 target_date = datetime(int(year), start_month, 15)573 nearest_idx = WAYBACK_DF.index.get_indexer([target_date], method='nearest')[0]574 wayback_item = WAYBACK_DF.iloc[nearest_idx]575 wayback_title = f"Esri Wayback ({wayback_item.name.strftime('%Y-%m-%d')})"576 wayback_url = (577 wayback_item["ResourceURL_Template"]578 .replace("{TileMatrixSet}", "GoogleMapsCompatible")579 .replace("{TileMatrix}", "{z}")580 .replace("{TileRow}", "{y}")581 .replace("{TileCol}", "{x}")582 )583 except Exception as e:584 print(f"Could not find a suitable Wayback basemap for {year}: {e}")585 wayback_url = None586 wayback_title = "Default Satellite"587 588 collection = (589 ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")590 .filterDate(start_date, end_date)591 .filterBounds(ee_geometry)592 .select(593 ["B2", "B3", "B4", "B8", "B11", "B12", "MSK_CLDPRB"],594 ["Blue", "Green", "Red", "NIR", "SWIR", "SWIR2", "MSK_CLDPRB"]595 )596 .map(lambda img: add_indices(img, 'NIR', 'Red', 'Blue', 'Green', 'SWIR', 'SWIR2', evi_vars))597 )598 599 if collection.size().getInfo() == 0:600 maps_html.append(f"<div style='text-align:center; padding-top: 50px;'>No data found for {year}.</div>")601 continue602 603 mosaic = collection.qualityMosaic(selected_index)604 m = gee_folium.Map(zoom_start=14)605 if wayback_url:606 m.add_tile_layer(wayback_url, name=wayback_title, attribution="Esri")607 else:608 m.add_basemap("SATELLITE")609 610 if selected_index in ["NDVI", "RandomForest", "GujVDI", "CI", "EVI", "EVI2"]:611 bins = [0, 0.2, 0.4, 0.6, 0.8, 1]612 histogram, _ = get_histogram(selected_index, mosaic.select(selected_index), ee_geometry, bins)613 614 total_pix = np.sum(histogram)615 formatted_histogram = ["0.00"] * len(histogram)616 if total_pix > 0:617 formatted_histogram = [f"{h*100/total_pix:.2f}" for h in histogram]618 619 # Define visualization parameters for the classified layer620 ind_vis_params = {621 "min": 0, "max": 1,622 "palette": ["#FF0000", "#FFFF00", "#FFA500", "#00FE00", "#00A400"],623 }624 m.addLayer(mosaic.select(selected_index).clip(ee_geometry), ind_vis_params, f"{selected_index} Classified Layer ({year})")625 626 legend_items = {627 f"0-0.2: Open/Sparse Vegetation ({formatted_histogram[0]}%)": "#FF0000",628 f"0.2-0.4: Low Vegetation ({formatted_histogram[1]}%)": "#FFFF00",629 f"0.4-0.6: Moderate Vegetation ({formatted_histogram[2]}%)": "#FFA500",630 f"0.6-0.8: Dense Vegetation ({formatted_histogram[3]}%)": "#00FE00",631 f"0.8-1: Very Dense Vegetation ({formatted_histogram[4]}%)": "#00A400",632 }633 634 legend_html = f'''635 <div style="636 position: fixed; 637 bottom: 20px; 638 right: 10px; 639 z-index:9999; 640 font-size:14px;641 background-color: rgba(255, 255, 255, 0.85);642 border:2px solid grey; 643 padding: 10px;644 border-radius: 6px;645 ">646 <b>{selected_index} Classification ({year})</b>647 <ul style="list-style-type: none; padding-left: 0; margin: 5px 0 0 0;">648 '''649 for text, color in legend_items.items():650 legend_html += f'<li><i style="background:{color}; width:20px; height:15px; display:inline-block; margin-right:5px; vertical-align:middle; border: 1px solid black;"></i> {text}</li>'651 legend_html += "</ul></div>"652 653 # Add the raw HTML to the map654 m.get_root().html.add_child(folium.Element(legend_html))655 656 else: # Fallback for indices without a color bar657 vis_params = {"min": 0.0, "max": 1.0, "palette": ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901', '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01', '012E01', '011D01', '011301']}658 clipped_image = mosaic.select(selected_index).clip(ee_geometry)659 m.addLayer(clipped_image, vis_params, f"{selected_index} {year}")660 m.add_colorbar(vis_params=vis_params, label=f"{selected_index} Value")661 662 folium.GeoJson(geometry_gdf, name="Geometry", style_function=lambda x: {"color": "yellow", "fillOpacity": 0, "weight": 2.5}).add_to(m)663 m.fit_bounds(map_bounds, padding=(10, 10))664 m.addLayerControl()665 maps_html.append(m._repr_html_())666 667 while len(maps_html) < 2:668 maps_html.append("")669 670 return f"Comparison generated for {selected_years[0]} and {selected_years[1]}.", maps_html[0], maps_html[1]671 672# --- Gradio Interface ---673theme = gr.themes.Soft(primary_hue="teal", secondary_hue="orange").set(674 background_fill_primary="white"675)676 677with gr.Blocks(theme=theme, title="Kamlan: KML Analyzer") as demo:678 # Hidden state to store data between steps679 geometry_data = gr.State()680 buffer_geometry_data = gr.State()681 timeseries_df_state = gr.State()682 683 gr.HTML("""684 <div style="display: flex; justify-content: space-between; align-items: center;">685 <img src="https://huggingface.co/spaces/SustainabilityLabIITGN/NDVI_PERG/resolve/main/Final_IITGN-Logo-symmetric-Color.png" style="width: 10%; margin-right: auto;">686 <h1 style="text-align: center;">Kamlan: KML Analyzer</h1>687 <img src="https://huggingface.co/spaces/SustainabilityLabIITGN/NDVI_PERG/resolve/main/IFS.jpg" style="width: 10%; margin-left: auto;">688 </div>689 """)690 691 with gr.Row():692 with gr.Column(scale=1):693 gr.Markdown("## 1. Provide Input Geometry")694 gr.Markdown("Use either file upload OR a URL.")695 file_input = gr.File(label="Upload KML/GeoJSON File", file_types=[".kml", ".geojson"])696 url_input = gr.Textbox(label="Or Provide File URL", placeholder="e.g., https://.../my_file.kml")697 buffer_input = gr.Number(label="Buffer (meters)", value=50)698 process_button = gr.Button("Process Input", variant="primary")699 700 with gr.Accordion("Advanced Settings", open=False):701 gr.Markdown("### Select Vegetation Indices")702 all_veg_indices = ["GujVDI", "NDVI", "EVI", "EVI2", "RandomForest", "CI", "MNDWI", "SAVI", "MVI", "NBR", "GCI"]703 veg_indices_checkboxes = gr.CheckboxGroup(all_veg_indices, label="Indices", value=["NDVI"])704 705 gr.Markdown("### EVI/EVI2 Parameters")706 with gr.Row():707 evi_g = gr.Number(label="G", value=2.5)708 evi_c1 = gr.Number(label="C1", value=6.0)709 evi_c2 = gr.Number(label="C2", value=7.5)710 with gr.Row():711 evi_l = gr.Number(label="L", value=1.0)712 evi_c = gr.Number(label="C", value=2.4)713 714 gr.Markdown("### Date Range")715 today = datetime.now()716 date_start_input = gr.Textbox(label="Start Date (MM-DD)", value="11-15")717 date_end_input = gr.Textbox(label="End Date (MM-DD)", value="12-15")718 719 with gr.Row():720 min_year_input = gr.Number(label="Start Year", value=2019, precision=0)721 max_year_input = gr.Number(label="End Year", value=today.year, precision=0)722 723 calculate_button = gr.Button("Calculate Vegetation Indices", variant="primary")724 725 726 with gr.Column(scale=2):727 gr.Markdown("## 2. Results")728 info_box = gr.Textbox(label="Status", interactive=False, placeholder="Status messages will appear here...")729 map_output = gr.HTML(label="Map View")730 stats_output = gr.DataFrame(label="Geometry Metrics")731 732 gr.Markdown("### Digital Elevation Model (DEM) and Slope")733 with gr.Row():734 dem_map_output = gr.HTML(label="DEM Map")735 slope_map_output = gr.HTML(label="Slope Map")736 737 with gr.Tabs():738 with gr.TabItem("Time Series Plot"):739 plot_output = gr.Plot(label="Time Series Plot")740 with gr.TabItem("Time Series Data"):741 timeseries_table = gr.DataFrame(label="Time Series Data")742 743 gr.Markdown("---") # Visual separator744 gr.Markdown("## 3. Year-over-Year Index Comparison")745 with gr.Row():746 comparison_index_select = gr.Radio(all_veg_indices, label="Select Index for Comparison", value="NDVI")747 comparison_years_select = gr.CheckboxGroup(label="Select Two Years to Compare", choices=[])748 749 compare_button = gr.Button("Generate Comparison Maps", variant="secondary")750 751 with gr.Row():752 map_year_1_output = gr.HTML(label="Comparison Map 1")753 map_year_2_output = gr.HTML(label="Comparison Map 2")754 755 756 # --- Event Handlers ---757 def process_on_load(request: gr.Request):758 """Checks for a 'file_url' query parameter when the app loads."""759 return request.query_params.get("file_url", "")760 761 demo.load(process_on_load, None, url_input)762 763 process_button.click(764 fn=process_and_display,765 inputs=[file_input, url_input, buffer_input],766 outputs=[map_output, info_box, stats_output, dem_map_output, slope_map_output, geometry_data, buffer_geometry_data]767 )768 769 def calculate_wrapper(geometry_json, buffer_json, veg_indices,770 g, c1, c2, l, c, start_date_str, end_date_str,771 min_year, max_year, progress=gr.Progress()):772 """Wrapper to parse inputs and handle outputs for the main calculation function."""773 try:774 evi_vars = {'G': g, 'C1': c1, 'C2': c2, 'L': l, 'C': c}775 start_month, start_day = map(int, start_date_str.split('-'))776 end_month, end_day = map(int, end_date_str.split('-'))777 date_range = (datetime(2000, start_month, start_day), datetime(2000, end_month, end_day))778 779 error_msg, df, plots, success_msg = calculate_indices(780 geometry_json, buffer_json, veg_indices,781 evi_vars, date_range, int(min_year), int(max_year), progress782 )783 784 status_message = error_msg or success_msg785 first_plot = plots[0] if plots else None786 df_display = df.round(3) if df is not None else None787 788 available_years = []789 if df is not None and 'Year' in df.columns:790 available_years = sorted(df['Year'].unique().tolist())791 792 return status_message, df_display, df, first_plot, gr.update(choices=available_years, value=[])793 794 except Exception as e:795 return f"An error occurred in the wrapper: {e}", None, None, None, gr.update(choices=[], value=[])796 797 calculate_button.click(798 fn=calculate_wrapper,799 inputs=[800 geometry_data, buffer_geometry_data, veg_indices_checkboxes,801 evi_g, evi_c1, evi_c2, evi_l, evi_c,802 date_start_input, date_end_input,803 min_year_input, max_year_input804 ],805 outputs=[info_box, timeseries_table, timeseries_df_state, plot_output, comparison_years_select]806 )807 808 def comparison_wrapper(geometry_json, selected_index, selected_years, g, c1, c2, l, c, start_date_str, end_date_str, progress=gr.Progress()):809 """Wrapper for the comparison map generation."""810 try:811 evi_vars = {'G': g, 'C1': c1, 'C2': c2, 'L': l, 'C': c}812 status, map1, map2 = generate_comparison_maps(813 geometry_json, selected_index, selected_years, evi_vars,814 start_date_str, end_date_str, progress815 )816 return status, map1, map2817 except Exception as e:818 return f"Error during comparison: {e}", "", ""819 820 compare_button.click(821 fn=comparison_wrapper,822 inputs=[823 geometry_data, comparison_index_select, comparison_years_select,824 evi_g, evi_c1, evi_c2, evi_l, evi_c,825 date_start_input, date_end_input826 ],827 outputs=[info_box, map_year_1_output, map_year_2_output]828 )829 830 gr.HTML("""831 <div style="text-align: center; margin-top: 20px;">832 <p>Developed by <a href="https://sustainability-lab.github.io/">Sustainability Lab</a>, <a href="https://www.iitgn.ac.in/">IIT Gandhinagar</a></p>833 <p>Supported by <a href="https://forests.gujarat.gov.in/">Gujarat Forest Department</a></p>834 </div>835 """)836 837def compute_index_histogram(file_url: str, index_name: str, start_date: str, end_date: str):838 """Computes histogram percentage for a vegetation index within a polygon and date range."""839 one_time_setup()840 decoded_url = unquote(file_url)841 gdf = get_gdf_from_url(decoded_url)842 gdf = preprocess_gdf(gdf)843 844 geometry_gdf = next((gdf.iloc[[i]] for i in range(len(gdf)) if is_valid_polygon(gdf.iloc[[i]])), None)845 if geometry_gdf is None:846 raise HTTPException(status_code=400, detail="No valid polygon found.")847 848 ee_geometry = ee.Geometry(json.loads(geometry_gdf.to_crs(4326).to_json())['features'][0]['geometry'])849 850 evi_vars = {'G': 2.5, 'C1': 6.0, 'C2': 7.5, 'L': 1.0, 'C': 2.4}851 collection = (852 ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")853 .filterDate(start_date, end_date)854 .select(855 ["B2", "B3", "B4", "B8", "B11", "B12", "MSK_CLDPRB"],856 ["Blue", "Green", "Red", "NIR", "SWIR", "SWIR2", "MSK_CLDPRB"]857 )858 .map(lambda img: add_indices(img, 'NIR', 'Red', 'Blue', 'Green', 'SWIR', 'SWIR2', evi_vars))859 .filterBounds(ee_geometry)860 )861 862 if collection.size().getInfo() == 0:863 raise HTTPException(status_code=404, detail="No imagery found for the polygon in given date range.")864 865 # Use median composite to avoid cloud spikes866 # mosaic = collection.median()867 mosaic = collection.qualityMosaic(index_name)868 mean_val = mosaic.reduceRegion(869 reducer=ee.Reducer.mean(),870 geometry=ee_geometry,871 scale=10,872 maxPixels=1e9873 ).get(index_name).getInfo()874 875 if mean_val is None:876 raise HTTPException(status_code=500, detail=f"Failed to compute mean for {index_name}.")877 878 bins = [0, 0.2, 0.4, 0.6, 0.8, 1.0]879 histogram, _ = get_histogram(index_name, mosaic.select(index_name), ee_geometry, bins)880 881 total = histogram.sum()882 883 # result = {f'mean': mean_val}884 # for i in range(len(bins) - 1):885 # label = f"{bins[i]}-{bins[i+1]}"886 # pct = (histogram[i] / total * 100) if total > 0 else 0887 # result[label] = round(pct, 2)888 889 890 output_parts = [f'mean, {round(mean_val, 4)}']891 892 for i in range(len(bins) - 1):893 label = f"{bins[i]}-{bins[i+1]}"894 pct = (histogram[i] / total * 100) if total > 0 else 0895 # Append each key-value pair string to the list896 output_parts.append(f'{label}, {round(pct, 2)}')897 898 # Join all parts into a single comma-separated string899 final_string = ", ".join(output_parts)900 901 # Return the string as a plain text response902 return Response(content=final_string, media_type="text/csv")903 904# Register endpoints for multiple indices with date range905for idx in ["NDVI", "EVI", "EVI2", "RandomForest", "CI", "GujVDI", "MNDWI", "SAVI", "MVI", "NBR", "GCI"]:906 endpoint_path = f"/api/{idx}"907 908 async def index_hist_endpoint(file_url: str, start_date: str, end_date: str, index_name=idx):909 return compute_index_histogram(file_url, index_name, start_date, end_date)910 911 app.get(endpoint_path)(index_hist_endpoint)912 913# Mount the Gradio app onto the FastAPI app914app = gr.mount_gradio_app(app, demo, path="/")915 916 917if __name__ == "__main__":918 import uvicorn919 uvicorn.run(app, host="0.0.0.0", port=7860)920