field_indices.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize Vegetation Indices in Timeseries\n",
    "\n",
    "*Author: Alessandro Joshua Pierro*\n",
    "\n",
    "*Affiliation: Agrscope*\n",
    "\n",
    "*Date: December 2023*\n",
    "\n",
    "## Overview\n",
    "\n",
    "This Python notebook is focused on visualizing various vegetation indices, specifically NDVI, NDRE, and EVI, over time using Sentinel-2 satellite imagery. The indices are vital for assessing the health and growth patterns of agricultural crops. This visualization helps in understanding the temporal changes in crop conditions and the effectiveness of agricultural practices.\n",
    "## Objective\n",
    "The primary goal is to provide a clear visual representation of changes in vegetation indices over time for specific agricultural parcels. This helps in identifying patterns, anomalies, and trends in crop growth and health, facilitating informed decision-making in agriculture management.\n",
    "## Methodology\n",
    "The process involves preprocessing Sentinel-2 imagery, extracting essential spectral data, and utilizing shapefiles for precise geographic targeting. \n",
    "## Key functions\n",
    "preprocess_sentinel2_scenes: Handles the resampling and masking of clouds, shadows, and snow in Sentinel-2 scenes.\n",
    "\n",
    "extract_s2_data: Extracts and filters Sentinel-2 data based on cloud cover and time period for given parcels.\n",
    "\n",
    "save_scene_as_geotiff: Saves processed Sentinel-2 scenes as GeoTIFF files for further analysis or archiving.\n",
    "\n",
    "plot_scenes: Generates plots for each vegetation index per scene, offering a visual representation of the data.\n",
    "\n",
    "process_shapefiles_in_subfolders: Automates the processing of shapefiles in subdirectories, facilitating batch processing of data.\n",
    "## Execution\n",
    "Run the code to process shapefiles in specified subfolders within a base directory. The script will iterate through each shapefile, extract relevant Sentinel-2 data, and generate plots for the NDVI, NDRE, and EVI indices for each parcel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Libaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import numpy as np\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import datetime\n",
    "import warnings\n",
    "\n",
    "from datetime import datetime\n",
    "from eodal.core.scene import SceneCollection\n",
    "from eodal.core.sensors.sentinel2 import Sentinel2\n",
    "from eodal.mapper.feature import Feature\n",
    "from eodal.mapper.filter import Filter\n",
    "from eodal.mapper.mapper import Mapper, MapperConfigs\n",
    "from eodal.config import get_settings\n",
    "\n",
    "# Set environment variable to enable working with STAC from Planetary Computer\n",
    "settings = get_settings()\n",
    "settings.USE_STAC = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define eodal functions and plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def preprocess_sentinel2_scenes(\n",
    "    ds: Sentinel2,\n",
    "    target_resolution: int,\n",
    ") -> Sentinel2:\n",
    "    \"\"\"\n",
    "    Resample Sentinel-2 scenes and mask clouds, shadows, and snow\n",
    "    based on the Scene Classification Layer (SCL).\n",
    "\n",
    "    NOTE:\n",
    "        Depending on your needs, the pre-processing function can be\n",
    "        fully customized using the full power of EOdal and its\n",
    "    interfacing libraries!\n",
    "\n",
    "    :param target_resolution:\n",
    "        spatial target resolution to resample all bands to.\n",
    "    :returns:\n",
    "        resampled, cloud-masked Sentinel-2 scene.\n",
    "    \"\"\"\n",
    "    # resample scene\n",
    "    ds.resample(inplace=True, target_resolution=target_resolution)\n",
    "    # mask clouds, shadows, and snow\n",
    "    ds.mask_clouds_and_shadows(inplace=True, cloud_classes=[3, 8, 9, 10, 11])\n",
    "    return ds\n",
    "\n",
    "def extract_s2_data(\n",
    "        parcel: gpd.GeoDataFrame,\n",
    "        time_start: datetime,\n",
    "        time_end: datetime,\n",
    "        scene_cloud_cover_threshold: float = 20,\n",
    "        feature_cloud_cover_threshold: float = 20,\n",
    "        spatial_resolution: int = 10\n",
    "    ) -> SceneCollection:\n",
    "    \"\"\"\n",
    "    Extracts Sentinel-2 data from the STAC SAT archive for a given area and time period.\n",
    "    Scenes that are too cloudy or contain nodata (blackfill), only, are discarded.\n",
    "\n",
    "    The processing level of the Sentinel-2 data is L2A (surface reflectance factors).\n",
    "\n",
    "    :param parcel:\n",
    "        field parcel geometry (defines the spatial extent to extract)\n",
    "    :param time_start:\n",
    "        start of the time period to extract\n",
    "    :param end_time:\n",
    "        end of the time period to extract\n",
    "    :param scene_cloud_cover_threshold:\n",
    "        scene-wide cloudy pixel percentage (from Sentinel-2 metadata) to filter out scenes\n",
    "        with too high cloud coverage values [0-100%]\n",
    "    :param feature_cloud_cover_threshold:\n",
    "        cloudy pixel percentage [0-100%] on the parcel level. Only if the parcel has a\n",
    "        lower percentual share of cloudy pixels (based on the scene classification layer) than\n",
    "        the threshold specified, the Sentinel-2 observation is kept\n",
    "    :param spatial_resolution:\n",
    "        spatial resolution of the Sentinel-2 data in meters (Def: 10m)\n",
    "    :param resampling_method:\n",
    "        spatial resampling method for those Sentinel-2 bands not available in the target\n",
    "        resolution. Nearest Neighbor by default\n",
    "    :returns:\n",
    "        dictionary with the list of scenes for the field parcel (`feature_scenes`), the\n",
    "        DataFrame of (un)used scenes and the reason for not using plus some basic scene\n",
    "        metadata (`scene_properties`)\n",
    "    \"\"\"\n",
    "    # setup the metadata filters (cloud cover and processing level)\n",
    "    metadata_filters = [\n",
    "        Filter('cloudy_pixel_percentage','<', scene_cloud_cover_threshold),\n",
    "        Filter('processing_level', '==', 'Level-2A')\n",
    "    ]\n",
    "    # setup the spatial feature for extracting data\n",
    "    feature = Feature.from_geoseries(parcel.geometry)\n",
    "    \n",
    "    # set up mapping configs\n",
    "    mapper_configs = MapperConfigs(\n",
    "        collection='sentinel2-msi',\n",
    "        time_start=time_start,\n",
    "        time_end=time_end,\n",
    "        feature=feature,\n",
    "        metadata_filters=metadata_filters\n",
    "    )\n",
    "\n",
    "    # get a new mapper instance. Set sensor to `sentinel2`\n",
    "    mapper = Mapper(mapper_configs)\n",
    "\n",
    "    # query the STAC (looks for available scenes in the selected spatio-temporal range)\n",
    "    mapper.query_scenes()\n",
    "\n",
    "    # get observations (loads the actual Sentinel2 scenes)\n",
    "    # the data is extract for the extent of the parcel\n",
    "    scene_kwargs = {\n",
    "        'scene_constructor': Sentinel2.from_safe,            # this tells the mapper how to read and load the data (i.e., Sentinel-2 scenes)\n",
    "        'scene_constructor_kwargs': {},                      # here you could specify which bands to read\n",
    "        'scene_modifier': preprocess_sentinel2_scenes,       # this tells the mapper about (optional) pre-processing of the loaded scenes (must be a callable)\n",
    "        'scene_modifier_kwargs': {'target_resolution': 10}   # here, you have to specify the value of the arguments the `scene_modifier` requires\n",
    "    }\n",
    "    mapper.load_scenes(scene_kwargs=scene_kwargs)\n",
    "\n",
    "    # loop over available Sentinel-2 scenes stored in mapper.data as a EOdal SceneCollection and check\n",
    "    # for no-data. These scenes will be removed from the SceneCollection\n",
    "    scenes_to_del = []\n",
    "    mapper.metadata['scene_used'] = 'yes'\n",
    "    for scene_id, scene in mapper.data:\n",
    "\n",
    "        # check if scene is blackfilled (nodata); if yes continue\n",
    "        if scene.is_blackfilled:\n",
    "            scenes_to_del.append(scene_id)\n",
    "            mapper.metadata.loc[mapper.metadata.sensing_time.dt.strftime('%Y-%m-%d %H:%M') == scene_id.strftime('%Y-%m-%d %H:%M')[0:16], 'scene_used'] = 'No [blackfill]'\n",
    "            continue\n",
    "\n",
    "        # check cloud coverage (including shadows and snow) of the field parcel\n",
    "        feature_cloud_cover = scene.get_cloudy_pixel_percentage(cloud_classes=[3, 8, 9, 10, 11])\n",
    "\n",
    "        # if the scene is too cloudy, we skip it\n",
    "        if feature_cloud_cover > feature_cloud_cover_threshold:\n",
    "            scenes_to_del.append(scene_id)\n",
    "            mapper.metadata.loc[mapper.metadata.sensing_time.dt.strftime('%Y-%m-%d %H:%M') == scene_id.strftime('%Y-%m-%d %H:%M')[0:16], 'scene_used'] = 'No [clouds]'\n",
    "            continue\n",
    "\n",
    "        # calculate the NDVI, NDRE and EVI\n",
    "        scene.calc_si('NDVI', inplace=True)\n",
    "        scene.calc_si('NDRE', inplace=True)\n",
    "        scene.calc_si('EVI', inplace=True)\n",
    "\n",
    "    # delete scenes too cloudy or containing only no-data\n",
    "    for scene_id in scenes_to_del:\n",
    "        del mapper.data[scene_id]\n",
    "    \n",
    "    return mapper\n",
    "\n",
    "def ensure_directory_exists(directory):\n",
    "    if not os.path.exists(directory):\n",
    "        os.makedirs(directory)\n",
    "        \n",
    "def save_scene_as_geotiff(scene, shapefile_name, metadata, scene_id, res):\n",
    "    indices = ['ndvi', 'ndre', 'evi']\n",
    "    base_paths = {\n",
    "        'ndvi': \"C:\\\\Users\\\\aless\\\\Desktop\\\\Agroscope\\\\FerN\\\\results\\\\parcel_visualisation\\\\field_indices\\\\raster\\\\ndvi\",\n",
    "        'ndre': \"C:\\\\Users\\\\aless\\\\Desktop\\\\Agroscope\\\\FerN\\\\results\\\\parcel_visualisation\\\\field_indices\\\\raster\\\\ndre\",\n",
    "        'evi': \"C:\\\\Users\\\\aless\\\\Desktop\\\\Agroscope\\\\FerN\\\\results\\\\parcel_visualisation\\\\field_indices\\\\raster\\\\evi\"\n",
    "    }\n",
    "\n",
    "    # Get the list of sensing dates marked with 'yes' in 'scene_used' column\n",
    "    sensing_dates_with_scene_used_list = metadata[metadata['scene_used'] == 'yes']['sensing_date'].tolist()\n",
    "\n",
    "    for idx, (scene_id, scene) in enumerate(res):\n",
    "        sensing_date_string = sensing_dates_with_scene_used_list[idx].strftime('%Y-%m-%d') if idx < len(sensing_dates_with_scene_used_list) else \"No Date\"\n",
    "        \n",
    "        # Create subfolders for each index and shapefile name\n",
    "        for index in indices:\n",
    "            subfolder_path = os.path.join(base_paths[index], shapefile_name)\n",
    "            ensure_directory_exists(subfolder_path)\n",
    "            \n",
    "            file_path = os.path.join(subfolder_path, f\"{shapefile_name}_{index}_{sensing_date_string}.tiff\")\n",
    "            scene.to_rasterio(file_path, band_selection=[index], as_cog=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def suppress_warnings():\n",
    "    warnings.filterwarnings('ignore')\n",
    "    \n",
    "def plot_scenes(res: SceneCollection, shapefile_name: str, metadata, save_path: str, scene_id) -> None:\n",
    "    f, axes = plt.subplots(ncols=len(res), nrows=3, figsize=(26, 14))  # Only 3 rows\n",
    "    # Check if axes is one-dimensional and convert to two-dimensional if necessary\n",
    "    \n",
    "    # Get list of sensing dates marked with 'yes' in 'scene_used' column\n",
    "    sensing_dates_with_scene_used_list = metadata[metadata['scene_used'] == 'yes']['sensing_date'].tolist()\n",
    "    \n",
    "    if len(res) == 1:\n",
    "        axes = np.array([[ax] for ax in axes]).reshape(3, -1)\n",
    "    \n",
    "    # Set the main title with the shapefile name\n",
    "    f.suptitle(shapefile_name, fontsize=16)\n",
    "    \n",
    "    for idx, (scene_id, scene) in enumerate(res):\n",
    "        # Get the sensing date for the current column and format it as a string\n",
    "        sensing_date_string = sensing_dates_with_scene_used_list[idx].strftime('%Y-%m-%d') if idx < len(sensing_dates_with_scene_used_list) else \"No Date\"\n",
    "        \n",
    "        # Sensing date titles\n",
    "        y_position = 0.92  # Decrease this value to move the text lower\n",
    "        plt.figtext(0.5/len(res) + (idx * (1/len(res))), y_position, sensing_date_string, ha='center', fontsize=8)\n",
    "        \n",
    "        # Plot NDVI on the first row\n",
    "        scene.plot_band(\n",
    "            'NDVI',\n",
    "            colormap='summer',\n",
    "            colorbar_label='NDVI',\n",
    "            vmin=0,\n",
    "            vmax=0.8,\n",
    "            ax=axes[0, idx]\n",
    "        )\n",
    "        axes[0, idx].set_title(f\"NDVI\")\n",
    "\n",
    "        # Plot NDRE on the second row\n",
    "        scene.plot_band(\n",
    "            'NDRE',\n",
    "            colormap='winter',\n",
    "            colorbar_label='NDRE',\n",
    "            vmin=0,\n",
    "            vmax=1.0,\n",
    "            ax=axes[1, idx]\n",
    "        )\n",
    "        axes[1, idx].set_title(f\"NDRE\")\n",
    "\n",
    "        # Plot EVI on the third row\n",
    "        scene.plot_band(\n",
    "            'EVI',\n",
    "            colormap='autumn',\n",
    "            colorbar_label='EVI',\n",
    "            vmin=0,\n",
    "            vmax=1.0,\n",
    "            ax=axes[2, idx]\n",
    "        )\n",
    "        axes[2, idx].set_title(f\"EVI\")\n",
    "\n",
    "        if idx > 0:\n",
    "            for jdx in range(3):\n",
    "                axes[jdx, idx].get_yaxis().set_ticks([])\n",
    "                axes[jdx, idx].set_ylabel('')\n",
    "                \n",
    "\n",
    "    plt.tight_layout(rect=[0, 0.03, 1, 0.95])  # Adjust the rect to make space for the main title\n",
    "\n",
    "    # Save plots as PNG and PDF\n",
    "    png_file = os.path.join(save_path, 'png', f\"{shapefile_name}.png\")\n",
    "    pdf_file = os.path.join(save_path, 'pdf', f\"{shapefile_name}.pdf\")\n",
    "    plt.savefig(png_file, format='png')\n",
    "    plt.savefig(pdf_file, format='pdf')\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_shapefiles_in_subfolders(base_folder: str, time_start: datetime, time_end: datetime):\n",
    "    suppress_warnings()\n",
    "    # Define the paths for saving the plots\n",
    "    save_path = \"C:\\\\Users\\\\aless\\\\Desktop\\\\Agroscope\\\\FerN\\\\results\\\\parcel_visualisation\\\\field_indices\"\n",
    "    # Iterate over all subdirectories in the base_folder\n",
    "    for subdir in os.listdir(base_folder):\n",
    "        subdir_path = os.path.join(base_folder, subdir)\n",
    "        # Check if it's a directory\n",
    "        if os.path.isdir(subdir_path):\n",
    "            # Use glob to find all .shp files in the current subdirectory\n",
    "            shapefiles = glob.glob(os.path.join(subdir_path, \"*.shp\"))\n",
    "            for shp_file in shapefiles:\n",
    "                # Extract the name without the '.shp' extension\n",
    "                shapefile_name = os.path.splitext(os.path.basename(shp_file))[0]\n",
    "                print(f\"Processing {shapefile_name}\")\n",
    "                # Read the shapefile\n",
    "                parcel = gpd.read_file(shp_file)\n",
    "                # Extract Sentinel-2 data for the shapefile\n",
    "                res = extract_s2_data(\n",
    "                    parcel=parcel,\n",
    "                    time_start=time_start,\n",
    "                    time_end=time_end\n",
    "                )\n",
    "                metadata = res.metadata[['product_uri', 'sensing_date', 'scene_used']]\n",
    "                for scene_id, scene in res.data:\n",
    "                   save_scene_as_geotiff(scene,shapefile_name,metadata,scene_id,res.data)   \n",
    "                # Plot the scenes if data is available\n",
    "                if res.data:\n",
    "                    plot_scenes(res.data,shapefile_name,metadata,save_path,scene_id)\n",
    "                else:\n",
    "                    print(f\"No data available for {shapefile_name}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Execute the functions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Bellechasse_Colza\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-12-21 11:32:19,045 eodal        INFO     Starting extraction of sentinel2 scenes\n",
      "2023-12-21 11:32:47,808 eodal        INFO     Finished extraction of sentinel2 scenes\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Blé_Grangeneuve\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-12-21 11:33:04,032 eodal        INFO     Starting extraction of sentinel2 scenes\n",
      "2023-12-21 11:33:39,875 eodal        INFO     Finished extraction of sentinel2 scenes\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Blé_Sorens\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-12-21 11:33:58,935 eodal        INFO     Starting extraction of sentinel2 scenes\n",
      "2023-12-21 11:34:32,219 eodal        INFO     Finished extraction of sentinel2 scenes\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Colza_Grangeneuve\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-12-21 11:34:48,816 eodal        INFO     Starting extraction of sentinel2 scenes\n",
      "2023-12-21 11:35:21,590 eodal        INFO     Finished extraction of sentinel2 scenes\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Epeautre_Bellechasse\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-12-21 11:35:38,799 eodal        INFO     Starting extraction of sentinel2 scenes\n",
      "2023-12-21 11:36:04,398 eodal        INFO     Finished extraction of sentinel2 scenes\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Mais_Combremont\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-12-21 11:36:20,705 eodal        INFO     Starting extraction of sentinel2 scenes\n",
      "2023-12-21 11:37:28,722 eodal        INFO     Finished extraction of sentinel2 scenes\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Mais_Payerne\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-12-21 11:37:55,037 eodal        INFO     Starting extraction of sentinel2 scenes\n",
      "2023-12-21 11:39:30,428 eodal        INFO     Finished extraction of sentinel2 scenes\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Pdt_Flamatt\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-12-21 11:39:55,992 eodal        INFO     Starting extraction of sentinel2 scenes\n",
      "2023-12-21 11:40:56,652 eodal        INFO     Finished extraction of sentinel2 scenes\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Pdt_Wuennewil\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-12-21 11:41:24,521 eodal        INFO     Starting extraction of sentinel2 scenes\n",
      "2023-12-21 11:42:23,922 eodal        INFO     Finished extraction of sentinel2 scenes\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Triticale_Grangeneuve\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-12-21 11:42:52,916 eodal        INFO     Starting extraction of sentinel2 scenes\n",
      "2023-12-21 11:43:26,571 eodal        INFO     Finished extraction of sentinel2 scenes\n"
     ]
    }
   ],
   "source": [
    "# Define the base folder path and the time period for data extraction\n",
    "base_folder_path = \"C:\\\\Users\\\\aless\\\\Desktop\\\\Agroscope\\\\FerN\\\\field_information\\\\Shapefiles_Fields\"\n",
    "time_start = datetime(2023, 1, 1)\n",
    "time_end = datetime(2023, 12, 21)\n",
    "\n",
    "# Call the function to process all shapefiles in subfolders\n",
    "process_shapefiles_in_subfolders(base_folder_path, time_start, time_end)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}