Mask_Extract_Aggregate_NDVI(RME)_p1.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "This code was used to mask quality NDVI, extracting NDVI time series, aggregating pixel counts, and saving them to Google Drive\n",
        "\n",
        "Created on Mon July 22 2024\n",
        "\n",
        "# @author: Pius N.Nwachukwu"
      ],
      "metadata": {
        "id": "dip5jkxTqryg"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# # Install necessary libraries\n",
        "\n",
        "# !pip install geemap\n",
        "# !pip install earthengine-api\n",
        "# !pip install pycrs==1.0.2"
      ],
      "metadata": {
        "id": "atSm5T5-MyJT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "igUcRSbgEs7N"
      },
      "outputs": [],
      "source": [
        "import os\n",
        "import shutil\n",
        "import ee\n",
        "import geemap\n",
        "import pandas as pd\n",
        "from google.colab import drive"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Authenticate and initialize the Earth Engine module.\n",
        "ee.Authenticate()\n",
        "ee.Initialize(project='ee-')\n",
        "\n",
        "# Mount Google Drive to Colab\n",
        "drive.mount('/content/drive')"
      ],
      "metadata": {
        "id": "Yvq8131EFCGt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Ensure that your shape files are prepared and stored in GDrive\n",
        "\n",
        "# Path to the folder in Google Drive where the shapefiles are stored\n",
        "shapefile_folder = '/content/drive/My Drive/New_Cont2'\n",
        "\n",
        "# Temporary directory to store extracted files\n",
        "temp_dir = '/content/drive/My Drive/New_Cont2'\n",
        "os.makedirs(temp_dir, exist_ok=True)\n"
      ],
      "metadata": {
        "id": "hpEMTknSFHzL"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "# Extract shapefiles from Google Drive folder\n",
        "if os.path.exists(shapefile_folder):\n",
        "    for file_name in os.listdir(shapefile_folder):\n",
        "        if file_name.endswith('.zip'):\n",
        "            shutil.unpack_archive(os.path.join(shapefile_folder, file_name), temp_dir)\n",
        "else:\n",
        "    print(f\"Error: Folder not found - {shapefile_folder}\")\n",
        "\n",
        "# Initialize a dictionary to hold the geometries and their names\n",
        "rois = {}\n",
        "# Load shapefiles as ROIs\n",
        "for file_name in os.listdir(temp_dir):\n",
        "    if file_name.endswith('.shp'):\n",
        "        shapefile_path = os.path.join(temp_dir, file_name)\n",
        "        roi = geemap.shp_to_ee(shapefile_path)\n",
        "        if roi is not None:\n",
        "            roi_name = os.path.splitext(file_name)[0]  # Use the shapefile name (without extension) as the ROI name\n",
        "            rois[roi_name] = roi  # Add to the dictionary\n",
        "        else:\n",
        "            print(f\"Warning: Failed to load shapefile: {shapefile_path}\")\n",
        "\n",
        "# Confirm: Check that the ROIs are loaded\n",
        "if not rois:\n",
        "    print(\"No shapefiles were loaded as ROIs. Please check the shapefile folder.\")\n",
        "else:\n",
        "    print(f\"{len(rois)} shapefiles loaded as ROIs.\")\n",
        "    for name in rois:\n",
        "        print(name)\n",
        "\n",
        "# Function to aggregate pixel count\n",
        "def aggregate_pixel_count(product, geometry, start_date, end_date):\n",
        "    # Load MODIS ImageCollection\n",
        "    collection = ee.ImageCollection(product) \\\n",
        "        .filterBounds(geometry) \\\n",
        "        .filterDate(start_date, end_date)\n",
        "\n",
        "    # Reduce the collection to get the total pixel count in the region of interest\n",
        "    count_image = collection.reduce(ee.Reducer.count()).select([0], ['pixel_count'])\n",
        "\n",
        "    # Clip to the geometry\n",
        "    count_image = count_image.clip(geometry)\n",
        "\n",
        "    return count_image\n",
        "\n",
        "# Function to mask invalid NDVI pixels based on the quality layer\n",
        "def mask_invalid_pixels(image):\n",
        "    ndvi = image.select('NDVI')\n",
        "    quality = image.select('DetailedQA')\n",
        "    valid_mask = quality.bitwiseAnd(0b11).eq(0)\n",
        "    masked_ndvi = ndvi.updateMask(valid_mask)\n",
        "    return masked_ndvi\n",
        "\n",
        "# Function to calculate the percentage of valid NDVI pixels\n",
        "def calculate_valid_percentage(image, geometry):\n",
        "    # Mask invalid NDVI values (e.g., -3000) and consider valid values (e.g., greater than -3000)\n",
        "    valid_pixels = image.select('NDVI').gt(-2000)\n",
        "    total_pixels = image.unmask().reduceRegion(\n",
        "        reducer=ee.Reducer.count(),\n",
        "        geometry=geometry,\n",
        "        scale=1000,\n",
        "        maxPixels=1e9\n",
        "    ).values().get(0)\n",
        "\n",
        "    valid_pixel_count = valid_pixels.reduceRegion(\n",
        "        reducer=ee.Reducer.count(),\n",
        "        geometry=geometry,\n",
        "        scale=1000,\n",
        "        maxPixels=1e9\n",
        "    ).values().get(0)\n",
        "\n",
        "    percentage_valid = ee.Number(valid_pixel_count).divide(ee.Number(total_pixels)).multiply(100)\n",
        "    return image.set('percentage_valid', percentage_valid).set('date', image.date().format('YYYY-MM-dd'))\n",
        "\n",
        "# Function to extract lat/lon values for each pixel\n",
        "def extract_lat_lon(image, geometry):\n",
        "    coords = ee.Image.pixelCoordinates(image.projection())\n",
        "    latlon = image.addBands(coords)\n",
        "    latlon = latlon.sample(region=geometry, scale=1000, geometries=True)\n",
        "    return latlon\n",
        "\n",
        "# Example usage:\n",
        "product = 'MODIS/061/MOD13A3'\n",
        "start_date = '2000-02-01'\n",
        "end_date = '2024-04-30'\n",
        "\n",
        "pixel_counts = {}\n",
        "for roi_name, roi in rois.items():\n",
        "    geometry = roi.geometry()\n",
        "    pixel_count_image = aggregate_pixel_count(product, geometry, start_date, end_date)\n",
        "    # Get the total pixel count for the ROI\n",
        "    pixel_count = pixel_count_image.reduceRegion(\n",
        "        reducer=ee.Reducer.sum(),\n",
        "        geometry=geometry,\n",
        "        scale=1000,  # Adjust scale as needed\n",
        "        maxPixels=1e13  # Increase the maximum number of pixels allowed\n",
        "    ).get('pixel_count').getInfo()\n",
        "    pixel_counts[roi_name] = pixel_count\n",
        "\n",
        "# Print the pixel counts for each ROI\n",
        "for roi_name, count in pixel_counts.items():\n",
        "    print(f'Pixel count for {roi_name}: {count}')\n",
        "\n",
        "# Extract NDVI\n",
        "# Load the MODIS NDVI dataset\n",
        "dataset = ee.ImageCollection(\"MODIS/061/MOD13A3\").select(['NDVI', 'DetailedQA'])\n",
        "\n",
        "# Initialize an empty list to store data\n",
        "all_data = []\n",
        "\n",
        "# Iterate over each ROI and extract data\n",
        "for roi_name, roi in rois.items():\n",
        "    geometry = roi.geometry()\n",
        "\n",
        "    # Get the images for the specified date range and apply the mask\n",
        "    images = dataset.filterDate(start_date, end_date).map(mask_invalid_pixels)\n",
        "\n",
        "    # Calculate the percentage of valid pixels for each image\n",
        "    images_with_percentage = images.map(lambda img: calculate_valid_percentage(img, geometry))\n",
        "\n",
        "    # Extract data for each image\n",
        "    for img in images_with_percentage.toList(images_with_percentage.size()).getInfo():\n",
        "        image = ee.Image(img['id'])\n",
        "        date = img['properties']['date']\n",
        "        percentage_valid = img['properties']['percentage_valid']\n",
        "\n",
        "        # Extract latitude and longitude\n",
        "        latlon = extract_lat_lon(image, geometry)  # Pass the geometry to the function\n",
        "\n",
        "        # Convert to a list of dictionaries\n",
        "        latlon_list = latlon.getInfo()['features']\n",
        "\n",
        "        # Extract the data and add ROI name and date\n",
        "        data = [{'ROI': roi_name,  # Add ROI name\n",
        "                 'date': date,\n",
        "                 'percentage_valid': percentage_valid,\n",
        "                 'latitude': feature['geometry']['coordinates'][1],\n",
        "                 'longitude': feature['geometry']['coordinates'][0],\n",
        "                 'NDVI': feature['properties']['NDVI']} for feature in latlon_list]\n",
        "        all_data.extend(data)  # Append data to the list\n",
        "\n",
        "# Convert to a pandas DataFrame\n",
        "df4 = pd.DataFrame(all_data)\n",
        "print(df4.head())\n"
      ],
      "metadata": {
        "id": "pQW76hxVFLWf"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "# # Define the file path for the new CSV file\n",
        "time_series_csv_file_path = '/content/drive/My Drive/FoRes/df4ndvi_timeseries_data.csv' #Store in my Gdrive\n",
        "\n",
        "# # Export the DataFrame to a CSV file\n",
        "df4.to_csv(time_series_csv_file_path, index=False)\n",
        "\n",
        "# # Confirm the file has been saved\n",
        "print(f\"Time series data has been saved to {time_series_csv_file_path}\")"
      ],
      "metadata": {
        "id": "0dCNgJV0FPxi"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}