{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "64a8cca0",
   "metadata": {},
   "source": [
    "# Exercise 2: Analyzing vegetation indices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4aab8777",
   "metadata": {},
   "source": [
    "## 1. Install dependencies if needed\n",
    "\n",
    "Run the commands below to install the dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cab6644",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -q pillow==10.0.1 matplotlib==3.8.0 scikit-image==0.21.0 scikit-learn==1.3.1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9547a9b",
   "metadata": {},
   "source": [
    "## 2. Compute the Normalized Difference Vegetation Index (NDVI) of sentinel images\n",
    "- Download the Sentinel-2 images from [here](https://drive.google.com/file/d/1NGfrldM4r7uC04nRgevD4DUH48Bnm72Y/view?usp=sharing).\n",
    "- Place the downloaded .zip file to the same folder as this notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "811b0ed6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import zipfile\n",
    "import os\n",
    "\n",
    "with zipfile.ZipFile(\"images_week2.zip\", 'r') as zip_ref:\n",
    "    zip_ref.extractall()\n",
    "    \n",
    "print(os.listdir(\"images_week2\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f8e5c01",
   "metadata": {},
   "source": [
    "- The zip file has 5 folders that contain files of 5 different sentinel images acquired in different timestamps. The folder names have the following structure: ***image_year_month***. Read the images using the function ***imread*** of the ***skimage*** library\n",
    "- Compute NDVI of each image and save it on disk\n",
    "- Visualize the NDVI image of the sentinel image taken on Aug 2021"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2afa2ef2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Fill the directory path where the images were uncompressed\n",
    "images_path = \"\"\n",
    "sub_directories = [\"image_2020_10\", \"image_2021_02\", \"image_2021_04\", \"image_2021_06\", \"image_2021_08\"]\n",
    "directories = [os.path.join(images_path, subdir) for subdir in sub_directories]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21379658",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from skimage.io import imsave, imread\n",
    "\n",
    "# Compute NDVI of several sentinel images\n",
    "for directory in directories:\n",
    "    # Define the full image paths to read the Red and NIR band of each sentinel image\n",
    "    band4_path = directory + \"/band4.jp2\" # Red\n",
    "    band8_path = directory + \"/band8.jp2\" # Near-infrared (NIR)\n",
    "    \n",
    "    # Read images\n",
    "    band4_array = imread(band4_path).astype(np.float32)\n",
    "    band8_array = imread(band8_path).astype(np.float32)\n",
    "    \n",
    "    # TODO: compute NDVI indexes\n",
    "    ndvi_array = \n",
    "    \n",
    "    # Save NDVI image\n",
    "    ndvi_path = directory + \"/ndvi.tif\".format(directory)\n",
    "    imsave(ndvi_path, ndvi_array)\n",
    "    print(\"Image saved in: \" + ndvi_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "338607c7",
   "metadata": {},
   "source": [
    "Now, let's visualize the NDVI images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "404a3385",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from  matplotlib.colors import LinearSegmentedColormap\n",
    "\n",
    "def plot_vegetation_index(image_arr):\n",
    "    # Remove outliers\n",
    "    p_min, p_max = np.percentile(image_arr[~np.isnan(image_arr)], (2, 98))\n",
    "    image_arr_clipped = image_arr.clip(p_min, p_max)\n",
    "    # Create a color map\n",
    "    cmap_rg=LinearSegmentedColormap.from_list('rg',\n",
    "                        [\"red\", \"yellow\", \"green\", \"darkgreen\"], N=256) \n",
    "    # Show image\n",
    "    plt.imshow(image_arr_clipped, cmap=cmap_rg)\n",
    "    plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58ecefe0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Visualize the NDVI image (that you saved on disk in previous code blocks) \n",
    "#       of the sentinel image taken on Aug 2021,\n",
    "#       using the function plot_vegetation_index (defined above)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ec42866",
   "metadata": {},
   "source": [
    "## 3. Compute and visualize the average NDVI values of 4 regions in the 5 images provided\n",
    "\n",
    "- Compute the average NDVI of 4 regions (defined by the bounding boxes) in the 5 images provided.\n",
    "- The sentinel images were taken on Oct 2020, Feb 2021,  Apr 2021, Jun 2021, and Aug 2021.\n",
    "- The bounding boxes are defined in the pixel coordinate system\n",
    "- Use the function ***crop_image*** defined below to crop the regions\n",
    "- Visualize the average NDVI values of the 4 regions as temporal series, using matplotlib. It can help to see the code example presented in slide 4 of the exercise 2 (available on moodle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e94ec1b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def crop_image_array(image_array, y_min, x_min, y_max, x_max):\n",
    "    \"\"\"\n",
    "    Crop image array, delimited by a bounding box\n",
    "    \n",
    "    Parameters\n",
    "    ==================\n",
    "    \n",
    "    y_min: int\n",
    "        Minimum vertical coordinate of the bounding box\n",
    "    x_min: int\n",
    "        Minimum horizontal coordinate of the bounding box\n",
    "    y_max: int\n",
    "        Maximum vertical coordinate of the bounding box\n",
    "    x_max: int\n",
    "        Maximum horizontal coordinate of the bounding box\n",
    "    \n",
    "    \"\"\"\n",
    "    return image_array[y_min:y_max+1, x_min:x_max+1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8261dce6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bounding boxes of the regions to be analyzed\n",
    "bounding_boxes = {'vineyard':    {'y_min': 814, 'x_min': 1509, 'y_max': 818, 'x_max': 1513}, \n",
    "                  'trees':       {'y_min': 3173, 'x_min': 3239, 'y_max': 3177, 'x_max': 3243}, \n",
    "                  'sport_field': {'y_min': 1278, 'x_min': 2628, 'y_max': 1282, 'x_max': 2632}, \n",
    "                  'buildings':   {'y_min': 1511, 'x_min': 2794, 'y_max': 1515, 'x_max': 2798}\n",
    "                 }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1d708cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Compute the average NDVI values \n",
    "#       of the 4 regions in the 5 images (sorted in chronological order),\n",
    "#       and stored it in a dictionary with the following structure. \n",
    "#\n",
    "# ndvi_values = {\n",
    "#   \"vineyard\" : [avg_ndvi_of_region_vineyard_in_image1, \n",
    "#                 avg_ndvi_of_region_vineyard_in_image2, \n",
    "#                ..., \n",
    "#                avg_ndvi_of_region_vineyard_in_image5]\n",
    "#\n",
    "#   \"trees\" :    [avg_ndvi_of_region_trees_in_image1, \n",
    "#                 avg_ndvi_of_region_trees_in_image2, \n",
    "#                ..., \n",
    "#                avg_ndvi_of_region_trees_in_image5]\n",
    "#\n",
    "#   \"buildings\" : [avg_ndvi_on_region_buildings_in_image1, \n",
    "#                 avg_ndvi_on_region_buildings_in_image2, \n",
    "#                ..., \n",
    "#                avg_ndvi_on_region_buildings_in_image5]\n",
    "#\n",
    "#   \"sport_field\" : [avg_ndvi_on_region_sport_field_in_image1, \n",
    "#                 avg_ndvi_on_region_sport_field_in_image2, \n",
    "#                ..., \n",
    "#                avg_ndvi_on_region_sport_field_in_image5]\n",
    "# }\n",
    "\n",
    "# Create an empty dictionary\n",
    "ndvi_values = {}\n",
    "# Iterate over all the regions. \n",
    "# The function \".keys()\" retrieves the keys of a dictionary\n",
    "# In the code below bounding_boxes.keys() = [\"vineyard\", \"trees\", \"buildings\", \"sport_field\"]\n",
    "for region in bounding_boxes.keys():\n",
    "    ndvi_values[region] = []\n",
    "    for directory in directories:\n",
    "        # TODO: Read the NDVI images you saved\n",
    "        \n",
    "        \n",
    "        # TODO: Crop image using the bounding box of each \"region\" \n",
    "        #       Use the function \"crop_image_array\" defined above\n",
    "        \n",
    "        \n",
    "        # TODO: Compute average NDVI value for the cropped region (using np.mean)\n",
    "        \n",
    "        \n",
    "        # TODO: Add the computed average NDVI value to the list ndvi_values[region]\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90d586df",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Plot the computed average NDVI values of the 4 regions in the 5 images\n",
    "#       follow the example of the slides of the exercise 2 (available on moodle)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c48a9aa4",
   "metadata": {},
   "source": [
    "## 4. Answer the remaining questions of the PDF instructions file"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}