Enhanced Vegetation Index (EVI) from MODIS (MOD13Q1).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "from google.colab import auth\n",
    "import geemap\n",
    "from ipyleaflet import Marker, Popup, WidgetControl\n",
    "from ipywidgets import HTML\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors\n",
    "\n",
    "# Authenticate to access Earth Engine\n",
    "auth.authenticate_user()\n",
    "\n",
    "# Initialize the Earth Engine module\n",
    "ee.Initialize(project='ee-2june2024')\n",
    "\n",
    "# Load the uploaded shapefile containing the updated boundaries for both districts and upazilas\n",
    "boundaries2020 = ee.FeatureCollection(\"projects/ee-2june2024/assets/bgd_adm_03\")\n",
    "\n",
    "# Filter for a specific district using the shapefile\n",
    "district_name = 'Narayanganj'  # Specify the district name here\n",
    "specific_district = boundaries2020.filter(ee.Filter.eq('district', district_name))\n",
    "\n",
    "# Load MODIS EVI data\n",
    "modis = ee.ImageCollection('MODIS/006/MOD13Q1') \\\n",
    "    .filterBounds(specific_district.geometry()) \\\n",
    "    .filterDate('2023-01-01', '2023-12-31') \\\n",
    "    .select('EVI')\n",
    "\n",
    "# Create a composite image of the median EVI\n",
    "evi_composite = modis.median().clip(specific_district.geometry())\n",
    "\n",
    "# Calculate the minimum and maximum EVI values in the clipped region\n",
    "min_max = evi_composite.reduceRegion(\n",
    "    reducer=ee.Reducer.minMax(),\n",
    "    geometry=specific_district.geometry(),\n",
    "    scale=500,\n",
    "    maxPixels=1e9\n",
    ")\n",
    "min_value = min_max.get('EVI_min').getInfo() if min_max.get('EVI_min') else 0\n",
    "max_value = min_max.get('EVI_max').getInfo() if min_max.get('EVI_max') else 1\n",
    "\n",
    "# Define legend keys based on dynamic min and max values\n",
    "num_legend_items = 11  # Number of legend items\n",
    "step = (max_value - min_value) / (num_legend_items - 1)\n",
    "legend_keys = [f\"{min_value + i * step:.2f}\" for i in range(num_legend_items)]\n",
    "\n",
    "# Generate gradient colors\n",
    "def generate_gradient_colors(num_colors):\n",
    "    cmap = plt.get_cmap('RdYlGn')  # You can choose other colormaps like 'plasma', 'inferno', 'magma', 'viridis' etc.\n",
    "    gradient_colors = [mcolors.rgb2hex(cmap(i / num_colors)) for i in range(num_colors)]\n",
    "    return gradient_colors\n",
    "\n",
    "gradient_colors = generate_gradient_colors(num_legend_items)\n",
    "\n",
    "# Reverse the gradient colors\n",
    "gradient_colors.reverse()\n",
    "\n",
    "# Function to get the highest intensity (EVI) and coordinates for each upazila\n",
    "def get_max_intensity(feature):\n",
    "    # Find the maximum EVI value within the feature geometry\n",
    "    max_value = evi_composite.reduceRegion(\n",
    "        reducer=ee.Reducer.max(),\n",
    "        geometry=feature.geometry(),\n",
    "        scale=500,\n",
    "        maxPixels=1e9\n",
    "    ).get('EVI')\n",
    "\n",
    "    # Convert the max_value to an image\n",
    "    max_value_image = ee.Image.constant(max_value)\n",
    "\n",
    "    # Find all pixels with the maximum EVI value\n",
    "    max_points = evi_composite.eq(max_value_image).selfMask().reduceToVectors(\n",
    "        geometry=feature.geometry(),\n",
    "        scale=500,\n",
    "        maxPixels=1e9,\n",
    "        geometryType='centroid',\n",
    "        eightConnected=False\n",
    "    )\n",
    "\n",
    "    return feature.set({\n",
    "        'max_intensity': max_value,\n",
    "        'max_points': max_points.geometry().coordinates()\n",
    "    })\n",
    "\n",
    "# Apply the function to each upazila within the specific district\n",
    "district_with_intensity = specific_district.map(get_max_intensity)\n",
    "\n",
    "# Export relevant properties for each upazila to a CSV file\n",
    "export_data = district_with_intensity.select(['upazila', 'max_intensity', 'max_points'])\n",
    "\n",
    "# Export the data to a CSV file in Google Drive\n",
    "task = ee.batch.Export.table.toDrive(\n",
    "    collection=export_data,\n",
    "    description='Narayanganj_Upazila_EVI',\n",
    "    fileFormat='CSV',\n",
    "    selectors=['upazila', 'max_intensity', 'max_points']\n",
    ")\n",
    "\n",
    "# Start the export task\n",
    "task.start()\n",
    "\n",
    "# Add layers to the map\n",
    "Map = geemap.Map()\n",
    "Map.addLayer(specific_district, {}, 'Narayanganj Upazilas')\n",
    "Map.addLayer(evi_composite, {'min': min_value, 'max': max_value, 'palette': gradient_colors}, 'EVI')\n",
    "\n",
    "# Center map on the specified district\n",
    "Map.centerObject(specific_district, 10)\n",
    "\n",
    "# Add markers for each upazila with highest EVI\n",
    "def add_markers(feature):\n",
    "    max_points = feature.get('max_points').getInfo()\n",
    "    intensity = feature.get('max_intensity').getInfo()\n",
    "    upazila_name = feature.get('upazila').getInfo()\n",
    "\n",
    "    # Handle nested structure of coordinates\n",
    "    if isinstance(max_points[0], list):  # Check if it's nested\n",
    "        for point in max_points:\n",
    "            lon, lat = point\n",
    "\n",
    "            # Verify the point is within the specific district\n",
    "            point_geom = ee.Geometry.Point(lon, lat)\n",
    "            if specific_district.geometry().contains(point_geom).getInfo():\n",
    "                # Create a marker\n",
    "                marker = Marker(location=(lat, lon))\n",
    "\n",
    "                # Add a popup\n",
    "                popup = Popup(\n",
    "                    location=(lat, lon),\n",
    "                    child=HTML(f\"<b>Upazila:</b> {upazila_name}<br><b>Max EVI:</b> {intensity}<br><b>Coordinates:</b> ({lat:.6f}, {lon:.6f})\"),\n",
    "                    close_button=False,\n",
    "                    auto_close=False,\n",
    "                    close_on_escape_key=False\n",
    "                )\n",
    "\n",
    "                # Add the marker and popup to the map\n",
    "                Map.add_layer(marker)\n",
    "                Map.add_layer(popup)\n",
    "    else:\n",
    "        lon, lat = max_points\n",
    "\n",
    "        # Verify the point is within the specific district\n",
    "        point_geom = ee.Geometry.Point(lon, lat)\n",
    "        if specific_district.geometry().contains(point_geom).getInfo():\n",
    "            # Create a marker\n",
    "            marker = Marker(location=(lat, lon))\n",
    "\n",
    "            # Add a popup\n",
    "            popup = Popup(\n",
    "                location=(lat, lon),\n",
    "                child=HTML(f\"<b>Upazila:</b> {upazila_name}<br><b>Max EVI:</b> {intensity}<br><b>Coordinates:</b> ({lat:.6f}, {lon:.6f})\"),\n",
    "                close_button=False,\n",
    "                auto_close=False,\n",
    "                close_on_escape_key=False\n",
    "            )\n",
    "\n",
    "            # Add the marker and popup to the map\n",
    "            Map.add_layer(marker)\n",
    "            Map.add_layer(popup)\n",
    "\n",
    "# Get the list of features and add markers for each\n",
    "features = district_with_intensity.getInfo()['features']\n",
    "for feature in features:\n",
    "    add_markers(ee.Feature(feature))\n",
    "\n",
    "# Add dynamic legend to the map\n",
    "legend_dict = {legend_keys[i]: gradient_colors[i] for i in range(num_legend_items)}\n",
    "Map.add_legend(title=\"EVI\", legend_dict=legend_dict, position=\"bottomleft\")\n",
    "\n",
    "# Display the map\n",
    "Map\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}