NDVI from COPERNICUSS2_SR_HARMONIZED.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries.\n",
    "# Authenticate to access Earth Engine.\n",
    "# Initialize the Earth Engine module.\n",
    "# Load the shapefile containing boundaries for districts and upazilas.\n",
    "# Filter for a specific district using the shapefile.\n",
    "# Define a function to calculate NDVI.\n",
    "# Load Sentinel-2 image collection and calculate NDVI.\n",
    "# Create a composite image of the median NDVI.\n",
    "# Calculate the minimum and maximum NDVI values in the clipped region.\n",
    "# Define legend keys based on dynamic min and max values.\n",
    "# Generate gradient colors.\n",
    "# Reverse the gradient colors.\n",
    "# Define a function to get the highest NDVI and coordinates for each upazila.\n",
    "# Apply the function to each upazila within the specific district.\n",
    "# Select relevant properties for each upazila to export.\n",
    "# Export the data to a CSV file in Google Drive.\n",
    "# Start the export task.\n",
    "# Add the specific district layer to the map.\n",
    "# Add the NDVI layer to the map.\n",
    "# Center the map on the specified district.\n",
    "# Define a function to add markers for each upazila with the highest NDVI.\n",
    "# Get the list of features and add markers for each.\n",
    "# Add a dynamic legend to the map.\n",
    "# Display the map.\n",
    "\n",
    "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",
    "# Function to calculate NDVI\n",
    "def calculate_ndvi(image):\n",
    "    ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI')\n",
    "    return image.addBands(ndvi)\n",
    "\n",
    "# Load Sentinel-2 image collection and calculate NDVI\n",
    "sentinel2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') \\\n",
    "    .filterBounds(specific_district.geometry()) \\\n",
    "    .filterDate('2023-01-01', '2023-07-31') \\\n",
    "    .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)) \\\n",
    "    .map(calculate_ndvi)\n",
    "\n",
    "# Create a composite image of the median NDVI\n",
    "ndvi_composite = sentinel2.select('NDVI').median().clip(specific_district.geometry())\n",
    "\n",
    "# Calculate the minimum and maximum NDVI values in the clipped region\n",
    "min_max = ndvi_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('NDVI_min').getInfo() if min_max.get('NDVI_min') else 0\n",
    "max_value = min_max.get('NDVI_max').getInfo() if min_max.get('NDVI_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 (NDVI) and coordinates for each upazila\n",
    "def get_max_intensity(feature):\n",
    "    # Find the maximum NDVI value within the feature geometry\n",
    "    max_value = ndvi_composite.reduceRegion(\n",
    "        reducer=ee.Reducer.max(),\n",
    "        geometry=feature.geometry(),\n",
    "        scale=500,\n",
    "        maxPixels=1e9\n",
    "    ).get('NDVI')\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 NDVI value\n",
    "    max_points = ndvi_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_NDVI',\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(ndvi_composite, {'min': min_value, 'max': max_value, 'palette': gradient_colors}, 'NDVI')\n",
    "\n",
    "# Center map on the specified district\n",
    "Map.centerObject(specific_district, 10)\n",
    "\n",
    "# Add markers for each upazila with highest NDVI\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 NDVI:</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 NDVI:</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=\"NDVI\", legend_dict=legend_dict, position=\"bottomleft\")\n",
    "\n",
    "# Display the map\n",
    "Map"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}