drought-monitoring.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://github.com/vhertel/drought-monitoring/blob/master/resources/header.png?raw=1\" width=\"1000\"/>\n",
    "\n",
    "# Drought Monitoring\n",
    "\n",
    "<img src=\"https://github.com/vhertel/drought-monitoring/blob/master/resources/example.png?raw=1\" width=\"1000\"/>\n",
    "\n",
    "***\n",
    "\n",
    "Spectral vegetation indices are among the most commonly used satellite data products for evaluation, monitoring, and measurement of vegetation cover, condition, biophysical processes, and changes. This [Recommended Practice](https://un-spider.org/advisory-support/recommended-practices) shows how to apply a multi-temporal analysis of the MODIS-based Standardized Vegetation Index (SVI) or Vegetation Condition Index (VCI) to support drought monitoring and early warning.\n",
    "\n",
    "This Jupyter Notebook covers the full processing chain from data query and download up to the export of the final data products by utilizing open access MODIS data. The tool's workflow follows the UN-SPIDER Recommended Practice on [Drought Monitoring](https://www.un-spider.org/advisory-support/recommended-practices/recommended-practice-drought-monitoring-using-standard) and is illustrated in the chart below. After entering user specifications, MODIS data can directly be downloaded from [AppEEARS](https://lpdaacsvc.cr.usgs.gov/appeears/) (Application for Extracting and Exploring Analysis Ready Samples). Subsequently, the data is processed and stored in a variety of output formats.\n",
    "\n",
    "<img src=\"https://github.com/vhertel/drought-monitoring/blob/master/resources/charts/chart0.png?raw=1\" width=\"1000\"/>\n",
    "\n",
    "***\n",
    "\n",
    "***File Structure***  \n",
    "The Jupyter Notebook file constitutes the directory of origin. Additional data is contained in subfolders. MODIS data needs to be stored in subfolders called *'input/evi_data'* and *'input/pixel_reliability'*. If no image is provided, the subfolders will automatically be created when accessing and downloading data from [AppEEARS](https://lpdaacsvc.cr.usgs.gov/appeears/) through this tool. The processed data is stored in a subfolder called *'output'*.  \n",
    "\n",
    "***Limitations***  \n",
    "MODIS sensors are mounted on two satellite platforms: Terra (launched in December 1999) and Aqua (launched in May 2002). This means, the time series go back to 2000/2002, which is not very long from a climatologic viewpoint. Since the SVI is using the mean and standard deviation of the time series, it is desirable to use a representative time span. It is, however, possible to extend the time series with AVHRR data provided that a thorough inter-calibration between the two sensor systems is given.  \n",
    "The SVI can only provide a relative comparison of the vegetation condition while the assessed deviation from the mean vegetation condition cannot be translated into an absolute deviation of for example the plant height. Neither can the SVI be interpreted for absolute quantification of agricultural damage.\n",
    "\n",
    "***\n",
    "\n",
    "## Initialization & User Input\n",
    "\n",
    "<img src=\"https://github.com/vhertel/drought-monitoring/blob/master/resources/charts/chart1.png?raw=1\" width=\"1000\"/>\n",
    "\n",
    "Please run the cell below for initialization and specify the working directory as well as whether data shall be downloaded from [AppEEARS](https://lpdaacsvc.cr.usgs.gov/appeears/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     56,
     75,
     162,
     184,
     214,
     237,
     254,
     370,
     436,
     463,
     485,
     520,
     555,
     582
    ]
   },
   "outputs": [],
   "source": [
    "# Click to run\n",
    "\n",
    "#####################################################\n",
    "###################### IMPORTS ######################\n",
    "#####################################################\n",
    "\n",
    "import requests as r\n",
    "import getpass, pprint, time, os, cgi, json\n",
    "import geopandas as gpd\n",
    "import fnmatch\n",
    "import imageio\n",
    "import os\n",
    "import re\n",
    "import gc\n",
    "import builtins\n",
    "import warnings\n",
    "import numpy as np\n",
    "import ipyleaflet\n",
    "import geopandas\n",
    "import json\n",
    "import ipywidgets\n",
    "import functools\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "import datetime\n",
    "from datetime import date\n",
    "from matplotlib import colors\n",
    "from ipyfilechooser import FileChooser\n",
    "from osgeo import gdal\n",
    "gdal.UseExceptions()\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# If tqdm is installed it will be used to display a nice progress bar\n",
    "try:\n",
    "    from tqdm.autonotebook import tqdm\n",
    "    import sys\n",
    "    # Overload built-in enumerate function\n",
    "    def enumerate(iterator, desc=None):\n",
    "        return builtins.enumerate(tqdm(iterator, desc=desc, file=sys.stdout,\n",
    "                                       leave=False))\n",
    "except ImportError:\n",
    "    tqdm = None\n",
    "    # Overload built-in enumerate function\n",
    "    def enumerate(iterator, desc=None):\n",
    "        if desc is not None:\n",
    "            print(desc)\n",
    "        return builtins.enumerate(iterator)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "####################################################\n",
    "############### FUNCTION DEFINITIONS ###############\n",
    "####################################################\n",
    "\n",
    "# Searches the input string for a DOY\n",
    "def get_doy(re_doy, string):\n",
    "    \"\"\"\n",
    "    Searches the input string for a DOY, and if one is found it returns a tuple\n",
    "    containing the DOY; else it returns None.\n",
    "\n",
    "    :param re_doy: compiled re pattern used to match a DOY\n",
    "    :param string: input string that will be searched for a DOY\n",
    "    :return : if a DOY is found a tuple of the form\n",
    "              (<int> year, <str> day of year), else None\n",
    "    \"\"\"\n",
    "    search_doy = re_doy.search(string)\n",
    "    if search_doy is None:  # Case where no doy was found\n",
    "        return None\n",
    "    doy = search_doy.group(1)\n",
    "    year = int(doy[:4])  # Treat year values as integers\n",
    "    day = doy[4:]  # Day values need to remain strings to be zero padded\n",
    "    return year, day\n",
    "\n",
    "# Checks that for each evi file there is a corresponding pixel reliability file\n",
    "def check_prepare_files(evi_files, pr_files):\n",
    "    \"\"\"\n",
    "    Checks that for each evi file there is a corresponding pixel reliability\n",
    "    file. Then a dictionary is returned containing the sorted files lists for\n",
    "    each DOY.\n",
    "\n",
    "    :param evi_files: list of evi files to be processed, containing a tuples of\n",
    "                       the form (<str> file_path, <int> year, <str> day)\n",
    "    :param pr_files: list of pixel reliability files to be processed,\n",
    "                     containing a tuples ofthe form\n",
    "                     (<str> file_path, <int> year, <str> day)\n",
    "    :return: dictionary where keys are the available DOYs and values are\n",
    "             dictionaries for evi_files and pr_files countaining the sorted\n",
    "             evi_files and sorted pr_files.\n",
    "    \"\"\"\n",
    "    # Dictionary keys have properties similar enough to a mathematical set for\n",
    "    # our purposes\n",
    "    doy_dict = dict()\n",
    "\n",
    "    for path, year, day in evi_files:\n",
    "        try:\n",
    "            doy_dict[day]['evi'].append((path, year, day))\n",
    "        except KeyError:\n",
    "            doy_dict[day] = dict()\n",
    "            doy_dict[day]['evi'] = [(path, year, day)]\n",
    "\n",
    "    for path, year, day in pr_files:\n",
    "        try:\n",
    "            doy_dict[day]['pr'].append((path, year, day))\n",
    "        except KeyError:\n",
    "            try:\n",
    "                doy_dict[day]['pr'] = [(path, year, day)]\n",
    "            except KeyError:\n",
    "                doy_dict[day] = dict()\n",
    "                doy_dict[day]['pr'] = [(path, year, day)]\n",
    "\n",
    "    if len(evi_files) != len(pr_files):\n",
    "        mismatch_doys = list()\n",
    "        for day in doy_dict.keys():\n",
    "            try:\n",
    "                if len(doy_dict[day]['evi']) != len(doy_dict[day]['pr']):\n",
    "                    mismatch_doys.append(day)\n",
    "            except KeyError:\n",
    "                mismatch_doys.append(day)\n",
    "        mismatch_files = dict()\n",
    "        for day in mismatch_doys:\n",
    "            if any(i not in doy_dict[day]['pr'] for i in doy_dict[day]['evi']):\n",
    "                evi_years = [i[1] for i in doy_dict[day]['evi']]\n",
    "                pr_years = [i[1] for i in doy_dict[day]['pr']]\n",
    "                mismatch_evi = [i for i in evi_years if i not in pr_years]\n",
    "                mismatch_pr = [i for i in pr_years if i not in evi_years]\n",
    "\n",
    "                # If the lists are empty future list comprehensions will fail\n",
    "                # to correctly exclude all files\n",
    "                if len(mismatch_evi) == 0:\n",
    "                    mismatch_evi.append(-1)\n",
    "                if len(mismatch_pr) == 0:\n",
    "                    mismatch_pr.append(-1)\n",
    "\n",
    "                mismatch_files[day] = dict()\n",
    "                mismatch_files[day]['evi'] = [i[0] for i in doy_dict[day]['evi']\n",
    "                                              if i[1] in mismatch_evi]\n",
    "\n",
    "                mismatch_files[day]['pr'] = [i[0] for i in doy_dict[day]['pr']\n",
    "                                             if i[1] in mismatch_pr]\n",
    "        war_msg = \"There is a mismatch in the number of evi data files and \"\n",
    "        war_msg += \"pixel reliability files.\\n\"\n",
    "        war_msg += \"The the mismatched files for each doy are:\\n\"\n",
    "        war_msg += str(mismatch_files)\n",
    "        war_msg += \"\\n\"\n",
    "        war_msg += \"These files will be ignored.\"\n",
    "        warnings.warn(war_msg)\n",
    "\n",
    "        for day, v in mismatch_files.items():\n",
    "            for k_sub, v_sub in v.items():\n",
    "                doy_dict[day][k_sub] = [i for i in doy_dict[day][k_sub]\n",
    "                                        if i[0] not in v_sub]\n",
    "    else:\n",
    "        for day in doy_dict.keys():\n",
    "            doy_dict[day]['evi'] = sorted(doy_dict[day]['evi'],\n",
    "                                          key=lambda x: x[1])\n",
    "            doy_dict[day]['pr'] = sorted(doy_dict[day]['pr'],\n",
    "                                         key=lambda x: x[1])\n",
    "\n",
    "    return doy_dict\n",
    "\n",
    "# Loads all pixel reliability images into the cloud mask array\n",
    "def load_cloud_mask(files_list, cloud_mask):\n",
    "    \"\"\"\n",
    "    Loads all pixel reliability images into the cloud mask array.\n",
    "\n",
    "    :param files_list: list of files to be processed, containing a tuples of\n",
    "                       the form (<str> file_path, <int> year, <str> day)\n",
    "    :param cloud_mask: matrix containing the processed pixel reliability data\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    for i, image in enumerate(files_list, \"Loading cloud_mask\"):\n",
    "        # Open the file from the pixel reliability file list\n",
    "        dataset = gdal.Open(image[0])\n",
    "        band = dataset.GetRasterBand(1)\n",
    "        # del dataset\n",
    "        data = band.ReadAsArray()\n",
    "\n",
    "        # Write the data from each file to the array\n",
    "        cloud_mask[:, :, i] = data\n",
    "\n",
    "    return cloud_mask\n",
    "\n",
    "# Filters the pixel reliability array\n",
    "def prepare_cloud_mask(cloud_mask):\n",
    "    \"\"\"\n",
    "    Filters the pixel reliability array.\n",
    "\n",
    "    :param cloud_mask: pixel reliability array\n",
    "    :return: filtered and rescaled pixel reliability matrix\n",
    "    \"\"\"\n",
    "    # Exchange value 0 with value 1 (Pixels with value 0 and 1 are used)\n",
    "    cloud_mask[cloud_mask == 0] = 1\n",
    "\n",
    "    # Set values that are of no interest to us to NaN\n",
    "    # cloud_mask[cloud_mask != 1] = np.nan\n",
    "    cloud_mask[cloud_mask != 1] = np.nan\n",
    "\n",
    "    # These are additional filters you may want to use\n",
    "    # Set value 2 to NA (Snow and Ice)\n",
    "    # cloud_mask[cloud_mask == 2] = np.nan\n",
    "\n",
    "    # Set value 3 to NA (Cloud)\n",
    "    # cloud_mask[cloud_mask == 3] = np.nan\n",
    "\n",
    "    # Set no data value (= -1) to NA\n",
    "    # cloud_mask[cloud_mask == -1] = np.nan\n",
    "\n",
    "    # Set all values above 3 to NA\n",
    "    # cloud_mask[cloud_mask > 3] = np.nan\n",
    "\n",
    "    return cloud_mask\n",
    "\n",
    "# Loads a single geotiff into the evi matrix\n",
    "def load_evi(files_list, cloud_mask):\n",
    "    \"\"\"\n",
    "    Loads a single geotiff into the evi matrix.\n",
    "\n",
    "    :param files_list: list of files to be processed, containing a tuples of\n",
    "                       the form (<str> file_path, <int> year, <str> day)\n",
    "    :param cloud_mask: matrix containing the processed pixel reliability data\n",
    "    :return: evi matrix with an additional geotiff of data\n",
    "    \"\"\"\n",
    "    for i, image in enumerate(files_list, \"Loading EVI\"):\n",
    "        # Open the file from the evi file list\n",
    "        dataset = gdal.Open(image[0])\n",
    "        band = dataset.GetRasterBand(1)\n",
    "        data = band.ReadAsArray()\n",
    "\n",
    "        # Apply the cloud mask and write the data from each file to the array\n",
    "        # Note: the evi data is multiplied into the cloud_mask matrix in order\n",
    "        # to save RAM and speed up the computation\n",
    "        cloud_mask[:, :, i] *= data\n",
    "\n",
    "    return cloud_mask\n",
    "\n",
    "# Filters the evi to set all negative values to nan and rescales the data\n",
    "def prepare_evi(evi):\n",
    "    \"\"\"\n",
    "    Filters the evi to set all negative values to nan and rescales the data.\n",
    "\n",
    "    :param evi: numpy array containing the results of multiplying the\n",
    "                cloud_mask data with the evi data\n",
    "    :return: filtered and rescaled input numpy array\n",
    "    \"\"\"\n",
    "    # Set negative values to nan\n",
    "    evi[evi < 0] = np.nan\n",
    "\n",
    "    # Rescale the data\n",
    "    evi *= 0.0001\n",
    "\n",
    "    return evi\n",
    "\n",
    "# Calculates the SVI\n",
    "def calculate_svi(evi):\n",
    "    \"\"\"\n",
    "    Calculates the SVI.\n",
    "\n",
    "    The standard deviation is calculated using the following formula:\n",
    "    std = sqrt((E[X ** 2] - E[X] ** 2) * (N / N - 1))\n",
    "    std = sqrt((mean(x ** 2) - mean(x) ** 2) * (N / N - 1))\n",
    "\n",
    "    Here N is the number of obervations, and the multiplier (N / N - 1) is the\n",
    "    adjustment factor necessary to convert the standard deviation to an\n",
    "    unbiased estimator of the variance of the infinite population.\n",
    "\n",
    "    The SVI itself is calculated using the following formula:\n",
    "    SVI = (EVI - mean(EVI)) / std(EVI)\n",
    "\n",
    "    A lot of the operations are done in place, such as squaring the evi matrix,\n",
    "    using it for something, and then taking the square root of it to get back\n",
    "    the original values. This may seem inefficient however it is much faster\n",
    "    than allocating large amount of memory and doing copy based operations.\n",
    "\n",
    "    When it comes to numpy arrays, abbreviated code like this:\n",
    "    a += b\n",
    "    c /= d\n",
    "    Means that the operations are performed in place.\n",
    "\n",
    "    Whereas operations like this:\n",
    "    a = a + b\n",
    "    c = c / d\n",
    "    First create a copy of the a and c arrays before performing the addition\n",
    "    and division. This leads to additional memory and computational overhead\n",
    "    due to creating a copy of the arrays.\n",
    "\n",
    "    Since the input evi matrix can have only positive values squaring and\n",
    "    taking the square root in place causes no computational errors.\n",
    "\n",
    "\n",
    "    Note: The code for computing the mean and standard deviation is loosely\n",
    "          based on the built-in numpy methods for nanmean and nanstd.\n",
    "\n",
    "    :param evi: numpy array containing the geotiff processing results\n",
    "    :return: evi array which has had its values replaced in-place with the SVI\n",
    "             values\n",
    "    \"\"\"\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.simplefilter('ignore')\n",
    "        # Allocate memory to where intermittent data will be saved\n",
    "        shape = (evi.shape[0], evi.shape[1], 1)\n",
    "        num_non_nan = np.zeros(shape, dtype=np.float64)\n",
    "        sum_ = np.zeros(shape, dtype=np.float64)\n",
    "        sum_squares = np.zeros(shape, dtype=np.float64)\n",
    "\n",
    "        # Create a mask that will say where nan values are, then negate it so\n",
    "        # it shows where the non nan values are\n",
    "        mask = ~np.isnan(evi)\n",
    "\n",
    "        # Sum the non nan values\n",
    "        # Note: setting keepdims to True will allow to do efficient array wide\n",
    "        # broadcasting in future calculations (this greatly increases\n",
    "        # computation speed and reduces memory usage)\n",
    "        np.sum(mask, axis=2, out=num_non_nan, dtype=np.intp, keepdims=True)\n",
    "\n",
    "        # Negate the mask so it shows where the nan values are\n",
    "        np.invert(mask, out=mask)\n",
    "\n",
    "        # Substitute nan values with 0 in place\n",
    "        evi[mask] = 0\n",
    "\n",
    "        # Calculate the sum\n",
    "        np.sum(evi, axis=2, out=sum_, dtype=np.float64, keepdims=True)\n",
    "\n",
    "        # Square all values in place\n",
    "        np.square(evi, out=evi)\n",
    "\n",
    "        # Calculate the sum of squares\n",
    "        np.sum(evi, axis=2, out=sum_squares, dtype=np.float64, keepdims=True)\n",
    "\n",
    "        # Take the square root of evi to get back all the original values\n",
    "        np.sqrt(evi, out=evi)\n",
    "\n",
    "        # Compute the mean by dividing the sum_ values in place\n",
    "        evi_mean = sum_  # Set easy to understand pointer\n",
    "        evi_mean /= num_non_nan  # Divide in place\n",
    "\n",
    "        # Square the mean in place\n",
    "        np.square(evi_mean, out=evi_mean)\n",
    "\n",
    "        # Compute the standard deviation in place (saving it in sum_squares)\n",
    "        evi_std = sum_squares  # Set easy to understand pointer\n",
    "        evi_std /= num_non_nan\n",
    "        evi_std -= evi_mean\n",
    "\n",
    "        # Adjust the standard deviation to be an unbiased estimator of the\n",
    "        # variance of the infinite population\n",
    "        evi_std *= num_non_nan\n",
    "        num_non_nan -= 1\n",
    "        num_non_nan[num_non_nan <= 0] = np.nan  # Avoid zero division errors\n",
    "\n",
    "        # Finish calculating the standard deviation\n",
    "        evi_std /= num_non_nan\n",
    "        np.sqrt(evi_std, out=evi_std)\n",
    "\n",
    "        # Take the square root in place to get back original mean\n",
    "        np.sqrt(evi_mean, out=evi_mean)\n",
    "\n",
    "        # Replace 0 with nan in the evi array\n",
    "        evi[mask] = np.nan\n",
    "\n",
    "        # Calculate the SVI\n",
    "        # Note: the svi is saved into the evi matrix in order to save RAM\n",
    "        svi = evi\n",
    "        svi -= evi_mean\n",
    "        svi /= evi_std\n",
    "\n",
    "        return svi\n",
    "    \n",
    "# Calculates the VCI\n",
    "def calculate_vci(evi):\n",
    "    \"\"\"\n",
    "    Calculates the VCI.\n",
    "\n",
    "    The VCI is calculated using the following formula:\n",
    "    VCI = ((EVI - min(EVI)) / (max(EVI) - min(EVI))) * 100\n",
    "\n",
    "    A lot of the operations are done in place, such as squaring the evi matrix,\n",
    "    using it for something, and then taking the square root of it to get back\n",
    "    the original values. This may seem inefficient however it is much faster\n",
    "    than allocating large amount of memory and doing copy based operations.\n",
    "\n",
    "    When it comes to numpy arrays, abbreviated code like this:\n",
    "    a += b\n",
    "    c /= d\n",
    "    Means that the operations are performed in place.\n",
    "\n",
    "    Whereas operations like this:\n",
    "    a = a + b\n",
    "    c = c / d\n",
    "    First create a copy of the a and c arrays before performing the addition\n",
    "    and division. This leads to additional memory and computational overhead\n",
    "    due to creating a copy of the arrays.\n",
    "\n",
    "    Since the input evi matrix can have only positive values squaring and\n",
    "    taking the square root in place causes no computational errors.\n",
    "\n",
    "\n",
    "    Note: The code for computing the mean and standard deviation is loosely\n",
    "          based on the built-in numpy methods for nanmean and nanstd.\n",
    "\n",
    "    :param evi: numpy array containing the geotiff processing results\n",
    "    :return: evi array which has had its values replaced in-place with the VCI\n",
    "             values\n",
    "    \"\"\"\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.simplefilter('ignore')\n",
    "        # Allocate memory to where intermittent data will be saved\n",
    "        shape = (evi.shape[0], evi.shape[1], 1)\n",
    "        min_ = np.zeros(shape, dtype=np.float64)\n",
    "        max_ = np.zeros(shape, dtype=np.float64)\n",
    "\n",
    "        # Create a mask that will say where nan values are, then negate it so\n",
    "        # it shows where the non nan values are\n",
    "        # mask = np.isnan(evi)\n",
    "\n",
    "        # Substitute nan values with 0 in place\n",
    "        # evi[mask] = 0\n",
    "\n",
    "        # Find the minimum values\n",
    "        np.fmin.reduce(evi, axis=2, out=min_, keepdims=True)\n",
    "\n",
    "        # Find the maximum values\n",
    "        np.fmax.reduce(evi, axis=2, out=max_, keepdims=True)\n",
    "\n",
    "        # Calculate the VCI\n",
    "        # Note: the vci is saved into the evi matrix in order to save RAM\n",
    "        vci = evi\n",
    "        vci -= min_\n",
    "        max_ -= min_\n",
    "        vci /= max_\n",
    "        vci *= 100\n",
    "\n",
    "        return vci\n",
    "\n",
    "# Generate an output png image representing the SVI for the input map area\n",
    "def generate_SVI_png(file_name, data, breaks, extent):\n",
    "    \"\"\"\n",
    "    Generate an output png image representing the SVI for the input map area.\n",
    "\n",
    "    :param file_name: <str> file name to be used when saving the png image\n",
    "    :param data: numpy view onto the svi for a single year\n",
    "    :param breaks: standard deviation values that should be used as breaks in\n",
    "                   the outputted png images\n",
    "    :param extent: tuple containing data for how to modify the output graphic\n",
    "                   in order for it to scale correctly\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    bounds = breaks  # Set easy to understand pointer\n",
    "\n",
    "    # Define the size of the figure (in inches)\n",
    "    fig, ax = plt.subplots(figsize=(5, 5))\n",
    "    plt.title(file_name)\n",
    "    cmap = colors.ListedColormap(['#4C0E0D', '#E72223', '#F19069', '#F9F6C6',\n",
    "                                  '#64B14B', '#04984A', '#00320E'])\n",
    "    norm = colors.BoundaryNorm(bounds, cmap.N)\n",
    "    cax = ax.imshow(data, cmap=cmap, norm=norm, extent=extent)\n",
    "    fig.colorbar(cax, cmap=cmap, norm=norm, boundaries=bounds, ticks=breaks)\n",
    "\n",
    "    plt.savefig(os.path.join(path_png, file_name + \".png\"), dpi=100)\n",
    "    plt.close()\n",
    "\n",
    "# Generate an output png image representing the VCI for the input map area\n",
    "def generate_VCI_png(file_name, data, extent):\n",
    "    \"\"\"\n",
    "    Generate an output png image representing the VCI for the input map area.\n",
    "\n",
    "    :param file_name: <str> file name to be used when saving the png image\n",
    "    :param data: numpy view onto the vci for a single year\n",
    "    :param extent: tuple containing data for how to modify the output graphic\n",
    "                   in order for it to scale correctly\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    # Define the size of the figure (in inches)\n",
    "    fig, ax = plt.subplots(figsize=(5, 5))\n",
    "    plt.title(file_name)\n",
    "    cmap = colors.ListedColormap(['#8B0000', '#FF4500', '#FFFF00', '#9ACD32',\n",
    "                                  '#008000'])\n",
    "    cax = ax.imshow(data, cmap=cmap, extent=extent)\n",
    "    fig.colorbar(cax, cmap=cmap)\n",
    "\n",
    "    plt.savefig(os.path.join(path_png, file_name + \".png\"), dpi=100)\n",
    "    plt.close()\n",
    "\n",
    "# Generate an output geotiff image representing the SVI/VCI for the input map area\n",
    "def generate_geotiff(file_name, data, geo_transform, projection):\n",
    "    \"\"\"\n",
    "    Generate an output geotiff image representing the SVI/VCI for the input map\n",
    "    area.\n",
    "\n",
    "    :param file_name: <str> file name to be used when saving the tif image\n",
    "    :param data: numpy view onto the svi/vci for a single year\n",
    "    :param geo_transform: geo transform data to be used in saving the tif image\n",
    "    :param projection: projection data to be used in saving the tif image\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    # Set geotiff output path\n",
    "    geotiff_path = os.path.join(path_tif, file_name + \".tif\")\n",
    "\n",
    "    # Read columns from data array\n",
    "    cols = data.shape[1]\n",
    "    # Read rows from data array\n",
    "    rows = data.shape[0]\n",
    "\n",
    "    # Set the driver to Geotiff\n",
    "    driver = gdal.GetDriverByName('GTiff')\n",
    "    # Create raster with shape of array as float64\n",
    "    out_raster = driver.Create(geotiff_path, cols, rows, 1, gdal.GDT_Float64)\n",
    "    # Read geo information from input file\n",
    "    out_raster.SetGeoTransform(geo_transform)\n",
    "    # Read band\n",
    "    out_band = out_raster.GetRasterBand(1)\n",
    "    # Set no data value to numpy's nan\n",
    "    out_band.SetNoDataValue(np.nan)\n",
    "    out_band.WriteArray(data)\n",
    "    # Set the projection according to the input file projection\n",
    "    out_raster.SetProjection(projection)\n",
    "    out_band.FlushCache()\n",
    "\n",
    "# Generates all the output png and tif images\n",
    "def generate_SVI_images(files_list, svi, extent, geo_transform, projection):\n",
    "    \"\"\"\n",
    "    Generates all the output png and tif images.\n",
    "\n",
    "    :param files_list: list of files to be processed, containing a tuples of\n",
    "                       the form (<str> file_path, <int> year, <str> day)\n",
    "    :param svi: numpy array containing the svi\n",
    "    :param extent:\n",
    "    :param geo_transform: geo transform data to be used in saving the tif image\n",
    "    :param projection: projection data to be used in saving the tif image\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    day = files_list[0][2]\n",
    "    years = [i[1] for i in files_list]\n",
    "\n",
    "    # Calculate the standard deviation values that will be used to define the\n",
    "    # color scheme\n",
    "    std = np.nanstd(svi[:, :, 0])\n",
    "    std15 = 1.5 * std\n",
    "    std2 = 2 * std\n",
    "    breaks = [-4, -std2, -std15, -std, std, std15, std2, 4]\n",
    "\n",
    "    for i, year in enumerate(years, \"Generating images\"):\n",
    "        file_name = \"SVI_{}_{}_{}\".format(study_area, day, year)\n",
    "        array = svi[:, :, i]\n",
    "        # Generating pngs for every time step\n",
    "        # Use pre-defined function to write pngs\n",
    "        generate_SVI_png(file_name, array, breaks, extent)\n",
    "\n",
    "        # Generating geotiffs for every time step\n",
    "        # Use pre-defined function to write geotiffs\n",
    "        generate_geotiff(file_name, array, geo_transform, projection)\n",
    "        gc.collect()  # Force garbage collection to save RAM\n",
    "\n",
    "# Generates all the output png and tif images\n",
    "def generate_VCI_images(files_list, vci, extent, geo_transform, projection):\n",
    "    \"\"\"\n",
    "    Generates all the output png and tif images.\n",
    "\n",
    "    :param files_list: list of files to be processed, containing a tuples of\n",
    "                       the form (<str> file_path, <int> year, <str> day)\n",
    "    :param vci: numpy array containing the vci\n",
    "    :param extent:\n",
    "    :param geo_transform: geo transform data to be used in saving the tif image\n",
    "    :param projection: projection data to be used in saving the tif image\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    day = files_list[0][2]\n",
    "    years = [i[1] for i in files_list]\n",
    "    for i, year in enumerate(years, \"Generating images\"):\n",
    "        file_name = \"VCI_{}_{}_{}\".format(study_area, day, year)\n",
    "        array = vci[:, :, i]\n",
    "        # Generating pngs for every time step\n",
    "        # Use pre-defined function to write pngs\n",
    "        generate_VCI_png(file_name, array, extent)\n",
    "\n",
    "        # Generating geotiffs for every time step\n",
    "        # Use pre-defined function to write geotiffs\n",
    "        generate_geotiff(file_name, array, geo_transform, projection)\n",
    "        gc.collect()  # Force garbage collection to save RAM\n",
    "   \n",
    "# Checks folder structure\n",
    "def filechooserCallback(chooser):\n",
    "    \n",
    "    global study_area\n",
    "    global input_path\n",
    "    global path_evi\n",
    "    global path_pr\n",
    "    global path_qual\n",
    "    global output_path\n",
    "    global path_png\n",
    "    global path_tif\n",
    "    \n",
    "    # Specify the name of your study area, this will be used for naming the output maps\n",
    "    study_area = chooser.selected_path.split('/')[-1]\n",
    "\n",
    "    # check if input and output folders exist, if note create folders\n",
    "    # input folders\n",
    "    input_path = os.path.join(chooser.selected_path, 'input')\n",
    "    if not os.path.isdir(input_path):\n",
    "        os.mkdir(input_path)\n",
    "    # folder for EVI geotiffs\n",
    "    path_evi = os.path.join(input_path, 'evi_data')\n",
    "    if not os.path.isdir(path_evi):\n",
    "        os.mkdir(path_evi)\n",
    "    # folder for pixel reliability geotiffs\n",
    "    path_pr = os.path.join(input_path, 'pixel_reliability')\n",
    "    if not os.path.isdir(path_pr):\n",
    "        os.mkdir(path_pr)\n",
    "    # folder for quality\n",
    "    path_qual = os.path.join(input_path, 'quality')\n",
    "    if not os.path.isdir(path_qual):\n",
    "        os.mkdir(path_qual)\n",
    "\n",
    "    # output folders\n",
    "    output_path = os.path.join(chooser.selected_path, 'output')\n",
    "    if not os.path.isdir(output_path):\n",
    "        os.mkdir(output_path)\n",
    "    # folder for png files\n",
    "    path_png = os.path.join(output_path, 'SVI_png')\n",
    "    if not os.path.isdir(path_png):\n",
    "        os.mkdir(path_png)\n",
    "    # folder for tif files\n",
    "    path_tif = os.path.join(output_path, 'SVI_tif')\n",
    "    if not os.path.isdir(path_tif):\n",
    "        os.mkdir(path_tif) \n",
    "\n",
    "\n",
    "        \n",
    "        \n",
    "####################################################\n",
    "####################### CODE #######################\n",
    "####################################################  \n",
    "\n",
    "fc = FileChooser(os.getcwd())\n",
    "fc.title = '<b>Please select path to folder:</b>'\n",
    "fc.show_only_dirs = True\n",
    "fc.use_dir_icons = True\n",
    "fc.register_callback(filechooserCallback)\n",
    "display(fc)\n",
    "\n",
    "download = ipywidgets.Dropdown(description = 'Data Download',\n",
    "                               options     = ['Yes', 'No'],\n",
    "                               layout      = ipywidgets.Layout(width='40%'),\n",
    "                               style       = {'description_width': '15ex'})\n",
    "display(download)\n",
    "\n",
    "index = ipywidgets.Dropdown(description = 'Index',\n",
    "                            options     = ['SVI', 'VCI'],\n",
    "                            layout      = ipywidgets.Layout(width='40%'),\n",
    "                            style       = {'description_width': '15ex'})\n",
    "display(index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download Data\n",
    "\n",
    "<img src=\"https://github.com/vhertel/drought-monitoring/blob/master/resources/charts/chart2.png?raw=1\" width=\"1000\"/>\n",
    "\n",
    "This section allows interactive data access and download from [AppEEARS](https://lpdaacsvc.cr.usgs.gov/appeears/). After entering the respective login details and the start date of the analysis, the desired country can be selected on the map. By clicking the *'Request'* button, the corresponding MODIS data will be queried and collected through [AppEEARS](https://lpdaacsvc.cr.usgs.gov/appeears/). Depending on the selected time period and country size, this process may take some minutes to hours. The status can be retrieved by clicking the *'Check Status'* button. The name of the folder selected above will be used as the title of the request. Once the request has been processed, the *'See Products'* button can be used to view a list of all queries on [AppEEARS](https://lpdaacsvc.cr.usgs.gov/appeears/) and to download the data to the folder selected above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     7,
     11,
     19,
     144,
     149,
     176
    ],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Click to run\n",
    "\n",
    "####################################################\n",
    "############### FUNCTION DEFINITIONS ###############\n",
    "####################################################\n",
    "\n",
    "# Shows selected country/area\n",
    "def click_handler(feature, **kwargs):\n",
    "    html.value = '<h5>Selected: <b>{}</b></h5>'.format(feature['properties']['name'])\n",
    "\n",
    "# Date conversion\n",
    "def datestdtojd(stddate):\n",
    "    fmt = '%d-%m-%Y'\n",
    "    sdtdate = datetime.datetime.strptime(stddate, fmt)\n",
    "    sdtdate = sdtdate.timetuple()\n",
    "    jdate = sdtdate.tm_yday\n",
    "    return(jdate)\n",
    "\n",
    "# Requests task at AppEEARS\n",
    "def on_requestButton_clicked(b, user, password):\n",
    "    \n",
    "    # Status update\n",
    "    print('Product was requested. Please check status for data query progress.')\n",
    "    \n",
    "    # Insert API URL, call login service, provide credentials & return json\n",
    "    token_response = r.post('{}login'.format(api), auth=(user.value, password.value)).json()\n",
    "    # Remove user and password information        \n",
    "    del user, password\n",
    "    \n",
    "    # Start a list for products to be requested, beginning with MOD13Q1.006\n",
    "    prods = ['MOD13Q1.006']     \n",
    "    # Request layers for the first product\n",
    "    lst_response = r.get('{}product/{}'.format(api, prods[0])).json()  \n",
    "    # Create tupled list linking desired product with desired layers\n",
    "    layers = [(prods[0],'_250m_16_days_EVI')]  \n",
    "    # Append to tupled list linking desired product with desired layers\n",
    "    layers.append((prods[0],'_250m_16_days_pixel_reliability')) \n",
    "    prodLayer = []\n",
    "    for l in layers:\n",
    "        prodLayer.append({\n",
    "                \"layer\": l[1],\n",
    "                \"product\": l[0]\n",
    "              })\n",
    "        \n",
    "    # Save login token to a variable    \n",
    "    token = token_response['token']   \n",
    "    # Create a header to store token information, needed to submit a request\n",
    "    global head\n",
    "    head = {'Authorization': 'Bearer {}'.format(token)}  \n",
    "    \n",
    "    # Checks whether country/area was selected\n",
    "    country = re.search('<h5>Selected: <b>(.*)</b></h5>', html.value)\n",
    "    if country is None:\n",
    "        raise Exception('Please select country on the map.')\n",
    "    # Extract area of interest and set to variable\n",
    "    selected_country = countries.query('name == \"%s\"' % country.group(1)).to_json()      \n",
    "    aoi_gc = json.loads(selected_country)\n",
    "\n",
    "    # Call to spatial API, return projs as json\n",
    "    projections = r.get('{}spatial/proj'.format(api)).json()  \n",
    "    # Create an empty dictionary\n",
    "    projs = {}                      \n",
    "    # Fill dictionary with `Name` as keys\n",
    "    for p in projections: \n",
    "        projs[p['Name']] = p                    \n",
    "    \n",
    "    # Create relevant dates required below\n",
    "    today = date.today()\n",
    "    # Convert to required format\n",
    "    d1 = today.strftime(\"%m-%d-%Y\")   \n",
    "    #d1_day = datestdtojd(d1)\n",
    "    #d1_year = today.strftime(\"%Y\")\n",
    "    #d1_day_jcd = d1_year + str(d1_day)\n",
    "    \n",
    "    # Create list with existing png files\n",
    "    file_list_folder = [f for f in os.listdir(path_png) if fnmatch.fnmatch(f, '*.png')] \n",
    "    # If files already exist, extract date of late file \n",
    "    if len(file_list_folder) != 0:            \n",
    "        # Sort files in file list according to year and doy.\n",
    "        s = sorted(file_list_folder, key = lambda x: (x.split('_')[-1],x.split('.')[0])) \n",
    "        # Extract date of last png in the folder.\n",
    "        last_jcd = s[-1].split('.')[0].split('_')[-1]+s[-1].split('_')[-2] \n",
    "        # Extract day of last png in folder\n",
    "        last_day = s[-1].split('_')[-2]         \n",
    "        # Extract day of last png in folder\n",
    "        last_year = s[-1].split('.')[0].split('_')[-1]     \n",
    "        # Convert date to normal datetime\n",
    "        d = datetime.datetime.strptime(last_jcd, '%Y%j').date()       \n",
    "        # Add 16 days as next image will be available for that date.\n",
    "        nd = d + datetime.timedelta(days=16)                               \n",
    "        \n",
    "        # Format to meet api date format.\n",
    "        last_image = d.strftime(\"%m-%d-%Y\")                                \n",
    "        next_image = nd.strftime(\"%m-%d-%Y\")\n",
    "\n",
    "        # Status update\n",
    "        print(f'The last image in the folder is from {last_image}. As a result, the start date for the request will be {next_image}.')\n",
    "\n",
    "    else:\n",
    "        # If list does not exist, start from first date for which data is available.\n",
    "        last_jcd = '0'                        \n",
    "        next_image = str(datetime.datetime.strptime(str(start_date_picker.value), '%Y-%m-%d').strftime('%m-%d-%Y'))\n",
    " \n",
    "    # Type of task, area or point\n",
    "    task_type = ['point','area']   \n",
    "    # Set output projection \n",
    "    proj = projs['geographic']['Name']  \n",
    "    # Set output file format type\n",
    "    outFormat = ['geotiff', 'netcdf4']  \n",
    "    # Start of the date range for which to extract data: MM-DD-YYYY\n",
    "    startDate = next_image              \n",
    "    # End of the date range for which to extract data: MM-DD-YYYY. Set as d1 to use date of last image available in folder\n",
    "    endDate = d1                 \n",
    "    # Specify True for a recurring date range\n",
    "    recurring = False                   \n",
    "    # if recurring = True, set yearRange, change start/end date to MM-DD\n",
    "    #yearRange = [2000,2016]            \n",
    "\n",
    "    task = {\n",
    "        'task_type': task_type[1],\n",
    "        'task_name': task_name,\n",
    "        'params': {\n",
    "             'dates': [\n",
    "             {\n",
    "                 'startDate': startDate,\n",
    "                 'endDate': endDate\n",
    "             }],\n",
    "             'layers': prodLayer,\n",
    "             'output': {\n",
    "                     'format': {\n",
    "                             'type': outFormat[0]}, \n",
    "                             'projection': proj},\n",
    "             'geo': aoi_gc,\n",
    "        }\n",
    "    }\n",
    "    \n",
    "    # Post json to the API task service, return response as json\n",
    "    global task_response\n",
    "    task_response = r.post('{}task'.format(api), json=task, headers=head).json()  \n",
    "\n",
    "    # enable button to check status\n",
    "    checkStatusButton.disabled = False    \n",
    "\n",
    "# Checks status on submitted task\n",
    "def on_checkStatusButton_clicked(b):\n",
    "    task_id = task_response['task_id']\n",
    "    print('Status: %s' % r.get('{}task/{}'.format(api, task_id), headers=head).json()['status'])\n",
    "\n",
    "# Displays all available data queries\n",
    "def on_productsButton_clicked(b, user, password):\n",
    "\n",
    "    # Insert API URL, call login service, provide credentials & return json\n",
    "    token_response = r.post('{}login'.format(api), auth=(user.value, password.value)).json() \n",
    "    # Remove user and password information        \n",
    "    del user, password         \n",
    "    # Save login token to a variable\n",
    "    token = token_response['token']    \n",
    "    # Create a header to store token information, needed to submit a request\n",
    "    head = {'Authorization': 'Bearer {}'.format(token)}  \n",
    "\n",
    "    # Query task service, setting params and header \n",
    "    tasks_response = r.get('{}task'.format(api), headers=head).json() \n",
    "    # print table with download buttons\n",
    "    grid = ipywidgets.GridspecLayout(len(tasks_response)+1, 4)\n",
    "    grid[0,0] = ipywidgets.HTML('<h4>Name</h4>')\n",
    "    grid[0,1] = ipywidgets.HTML('<h4>Requested on</h4>')\n",
    "    grid[0,2] = ipywidgets.HTML('<h4>Status</h4>')\n",
    "    for i in range(len(tasks_response)):\n",
    "        grid[i+1,0] = ipywidgets.Label(tasks_response[i]['task_name'])\n",
    "        grid[i+1,1] = ipywidgets.Label(tasks_response[i]['created'])\n",
    "        grid[i+1,2] = ipywidgets.Label(tasks_response[i]['status'])\n",
    "        grid[i+1,3] = ipywidgets.Button(description = 'Download')\n",
    "        grid[i+1,3].on_click(functools.partial(on_downloadButton_clicked, tasks_response=tasks_response[i]))\n",
    "    display(grid)\n",
    "    \n",
    "# Downloads data\n",
    "def on_downloadButton_clicked(b, tasks_response):\n",
    "        \n",
    "    # Checks status of requested dataset\n",
    "    if tasks_response['status'] == 'done':    \n",
    "\n",
    "        task_id = tasks_response['task_id']\n",
    "        # Call API and return bundle contents for the task_id as json\n",
    "        bundle = r.get('{}bundle/{}'.format(api, task_id)).json()  \n",
    "        # Create empty dictionary\n",
    "        files = {}            \n",
    "        # Fill dictionary with file_id as keys and file_name as values\n",
    "        for f in bundle['files']: \n",
    "            files[f['file_id']] = f['file_name']                        \n",
    "        # Instantiates progress bar\n",
    "        progressBar = ipywidgets.IntProgress(description = 'Downloading %d images:' % len(files),\n",
    "                                             min         = 0,\n",
    "                                             max         = len(files),\n",
    "                                             layout      = ipywidgets.Layout(width='60%'),\n",
    "                                             style       = {'description_width': '25ex'})    \n",
    "        # display the bar\n",
    "        display(progressBar)                                           \n",
    "            \n",
    "        for f in files:\n",
    "            # Get a stream to the bundle file\n",
    "            dl = r.get('{}bundle/{}/{}'.format(api, task_id, f), stream=True)    \n",
    "            # Parse the name from Content-Disposition header \n",
    "            filename = os.path.basename(cgi.parse_header(dl.headers['Content-Disposition'])[1]['filename'])  \n",
    "            if fnmatch.fnmatch(filename, '*EVI*'):\n",
    "                filepath = os.path.join(path_evi, filename)\n",
    "            elif fnmatch.fnmatch(filename,'*pixel_reliability*'):\n",
    "                filepath = os.path.join(path_pr, filename)\n",
    "            elif fnmatch.fnmatch(filename,'*Quality*'):\n",
    "                filepath = os.path.join(path_qual, filename)\n",
    "            else:\n",
    "                filepath = os.path.join(input_path, filename)\n",
    "            # Write file to dest dir\n",
    "            with open(filepath, 'wb') as f:                                                                  \n",
    "                for data in dl.iter_content(chunk_size=8192): \n",
    "                    f.write(data) \n",
    "            progressBar.value += 1\n",
    "        print('Downloaded files can be found at: {}'.format(input_path))\n",
    "\n",
    "    else:\n",
    "        print('Data query not complete. Please check status.')\n",
    "\n",
    "\n",
    "\n",
    "####################################################\n",
    "####################### CODE #######################\n",
    "####################################################\n",
    "        \n",
    "# Checks whether data shall be downloaded\n",
    "if download.value == 'Yes':\n",
    "    # Creates interactive map for country/area selection\n",
    "    m = ipyleaflet.Map(zoom = 2, basemap = ipyleaflet.basemaps.CartoDB.Positron)\n",
    "    countries = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))\n",
    "    geo_data = ipyleaflet.GeoData(geo_dataframe = countries,\n",
    "                                  style={'opacity' : 0, 'fillOpacity' : 0},\n",
    "                                  hover_style={'fillColor': 'red' , 'fillOpacity': 0.2},\n",
    "                                  name = 'Countries')\n",
    "    m.add_layer(geo_data)\n",
    "    tile_ID = geo_data.on_click(click_handler)\n",
    "    html = ipywidgets.HTML('')\n",
    "    control = ipyleaflet.WidgetControl(widget=html, position='topright')\n",
    "    m.add_control(control)\n",
    "    display(m)\n",
    "\n",
    "    # Give the project a name. This will be used throughout the script.\n",
    "    # User-defined name of the task\n",
    "    task_name = study_area \n",
    "\n",
    "    # Change to working directory\n",
    "    os.chdir(input_path)                     \n",
    "    # Set the AρρEEARS API to a variable\n",
    "    api = 'https://lpdaacsvc.cr.usgs.gov/appeears/api/'  \n",
    "    # Input NASA Earthdata Login Username\n",
    "    user = ipywidgets.Text(description = 'Username')      \n",
    "    # Input NASA Earthdata Login Password\n",
    "    password = ipywidgets.Password(description = 'Password')  \n",
    "    start_date_picker = ipywidgets.DatePicker(description='Start Date:', value=datetime.datetime(2000,1,1).date())\n",
    "    requestButton = ipywidgets.Button(description = 'Request')\n",
    "    requestButton.on_click(functools.partial(on_requestButton_clicked, user=user, password=password))\n",
    "    productsButton = ipywidgets.Button(description = 'See Products')\n",
    "    productsButton.on_click(functools.partial(on_productsButton_clicked, user=user, password=password))\n",
    "    checkStatusButton = ipywidgets.Button(description = 'Check Status', disabled=True)\n",
    "    checkStatusButton.on_click(on_checkStatusButton_clicked)\n",
    "    display(ipywidgets.HBox([user, password, start_date_picker]))    \n",
    "    display(ipywidgets.HBox([requestButton, checkStatusButton, productsButton]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Processing\n",
    "\n",
    "<img src=\"https://github.com/vhertel/drought-monitoring/blob/master/resources/charts/chart2.png?raw=1\" width=\"1000\"/>\n",
    "\n",
    "### Standardized Vegetation Index (SVI)\n",
    "\n",
    "Originally building on the NDVI anomaly concept, the Standardized Vegetation Index (SVI) developed by Peters et al. (2002) describes the probability of variation from the normal NDVI over multiple years of data (e.g. 12 years) on a weekly time step. The SVI is a z-score deviation from the mean in units of the standard deviation, calculated from the NDVI or EVI values for each pixel location of a composite period for each year during a given reference period. The equation below shows the general calculation of the SVI\n",
    "\n",
    "$$z_{ijk} = \\frac{\\text{VI}_{ijk} - \\mu_{ij}}{\\sigma_{ij}}$$\n",
    "\n",
    "where $z_{ijk}$ is the z-value for the pixel $i$ during week $j$ for year $k$, $\\text{VI}_{ij}$ is the weekly $\\text{VI}$ value for pixel $i$ during week $j$ for year $k$ whereby both the NDVI or EVI can be utilized as $\\text{VI}$, $\\mu_{ij}$ is the mean for pixel $i$ during week $j$ over $n$ years, and $\\sigma_{ij}$ is the standard deviation of pixel $i$ during week $j$ over $n$ years. This recommended practice calculates the SVI based on the Enhanced Vegetation Index (EVI) which has some advantages compared to the NDVI such as an improved sensitivity over dense vegetation conditions and is less affected by aerosol influences. More detailed information on the Standardized Vegetation Index can be found [here](https://www.un-spider.org/advisory-support/recommended-practices/recommended-practice-agricultural-drought-monitoring-svi/in-detail/standardized-vegetation-index).\n",
    "\n",
    "### Vegetation Condition Index (VCI)\n",
    "\n",
    "Kogan proposed a Vegetation Condition Index (VCI) based on the relative Normalized Difference Vegetation Index (NDVI) change with respect to minimum historical NDVI value. The VCI therefore compares the current Vegetation Index (VI) such as NDVI or Enhanced Vegetation Index (EVI) to the values observed in the same period in previous years within a specific pixel. The VCI is calculated as shown below,\n",
    "\n",
    "$$ \\text{VCI}_{ijk} = \\frac{\\text{VI}_{ijk} - \\text{VI}_{i,min}}{\\text{VI}_{i,max} - \\text{VI}_{i,min}} \\cdot 100$$\n",
    "\n",
    "where $\\text{VCI}_{ijk}$ is the $\\text{VCI}$ value for the pixel $i$ during week/month/DOY $j$ for year $k$, $\\text{VI}_{ijk}$ is the weekly/monthly/DOYs $\\text{VI}$ value for pixel $i$ in week/month/DOY $j$ for year $k$ whereby both the NDVI or EVI can utilized as $\\text{VI}$, $\\text{VI}_{i,min}$ and $\\text{VI}_{i,max}$ are the multi-year minimum and maximum $\\text{VI}$, respectively, for pixel $i$.\n",
    "\n",
    "The resulting percentage of the observed value is situated between the extreme values (minimum and maximum) in the previous years. Lower and higher values, therefore, indicate bad and good vegetation state conditions, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Click to run\n",
    "\n",
    "####################################################\n",
    "####################### CODE #######################\n",
    "####################################################\n",
    "\n",
    "# The following regular expression will match 'doy' followed by 7 digits,\n",
    "# followed by 0 or more characters and ending in '.tif'\n",
    "# It is used for building the initial file list of tif images that have a valid\n",
    "# DOY structure\n",
    "re_doy = re.compile(r\".*doy(\\d{7}).*\\.tif$\")\n",
    "\n",
    "# Build initial file list of tuples containing the filename, year and day for\n",
    "# each file\n",
    "# e.g. [(<str> '/path/to/my_evi_doy2000193.tif', <int> 2000, <str> '193')]\n",
    "# Note: the data types are stated in '<>'\n",
    "evi_files = []\n",
    "\n",
    "# Create a list of files, which include the defined DOY in their filename\n",
    "# for the EVI data\n",
    "for _, _, files in os.walk(path_evi):\n",
    "    for file in files:\n",
    "        doy = get_doy(re_doy, file)\n",
    "        if doy is None:\n",
    "            continue\n",
    "        else:\n",
    "            evi_files.append((os.path.join(path_evi, file), doy[0], doy[1]))\n",
    "\n",
    "# Create a list of files, which include the defined DOY in their filename\n",
    "# for the pixel reliability data\n",
    "pr_files = []\n",
    "\n",
    "for _, _, files in os.walk(path_pr):\n",
    "    for file in files:\n",
    "        doy = get_doy(re_doy, file)\n",
    "        if doy is None:\n",
    "            continue\n",
    "        else:\n",
    "            pr_files.append((os.path.join(path_pr, file), doy[0], doy[1]))\n",
    "\n",
    "# Read an example file and define the shape of the data arrays\n",
    "# Get the first file of the file list as example file\n",
    "example_file = gdal.Open(evi_files[0][0])\n",
    "\n",
    "# Store necessary data from the example file\n",
    "geo_transform = example_file.GetGeoTransform()\n",
    "projection = example_file.GetProjection()\n",
    "x_size = example_file.RasterXSize\n",
    "y_size = example_file.RasterYSize\n",
    "\n",
    "del example_file  # Save RAM\n",
    "\n",
    "# Preparing map reshaping\n",
    "lon_start = geo_transform[0]\n",
    "lon_stop = (x_size * geo_transform[1]) + geo_transform[0]\n",
    "lon_step = geo_transform[1]\n",
    "lat_start = geo_transform[3]\n",
    "lat_stop = (y_size * geo_transform[5]) + geo_transform[3]\n",
    "lat_step = geo_transform[5]\n",
    "\n",
    "extent = (lon_start, lon_stop, lat_stop, lat_start)\n",
    "\n",
    "# Build a dictionary where keys are the available DOYs and values are sorted\n",
    "# dictionaries for evi files and pr files\n",
    "doy_dict = check_prepare_files(evi_files, pr_files)\n",
    "\n",
    "# The script creates some warnings due to the NA in the data.\n",
    "# They are ignored by executing this cell.\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# If tqdm is available, use it.\n",
    "# Note: both update_iter_desc and days_iterator are set in either case since it\n",
    "# allows to avoid having many if clauses throughout the program and reduces\n",
    "# code duplication.\n",
    "if tqdm is not None:\n",
    "    def update_iter_desc(days_iterator, desc):\n",
    "        days_iterator.set_description(desc)\n",
    "    days_iterator = tqdm(sorted(doy_dict), desc=\"Processing\",\n",
    "                         file=sys.stdout)\n",
    "else:\n",
    "    def update_iter_desc(days_iterator, desc):\n",
    "            print(desc)\n",
    "    days_iterator = sorted(doy_dict)\n",
    "\n",
    "# Begin iterating through days of the year for which we have data.\n",
    "for day in days_iterator:\n",
    "    # Provide a progress update\n",
    "    update_iter_desc(days_iterator, \"Processing DOY {}\".format(day))\n",
    "\n",
    "    # Get the necessary files lists from doy_dict\n",
    "    evi_files_day = doy_dict[day]['evi']\n",
    "    pr_files_day = doy_dict[day]['pr']\n",
    "\n",
    "    # Number of years for which we have data\n",
    "    num_years = len(evi_files_day)\n",
    "\n",
    "    # Create a zero-filled 3D numpy array based on the example file\n",
    "    # dimensions\n",
    "    array_size = (y_size, x_size, num_years)\n",
    "    try:\n",
    "        # Adjust the size of the cloud mask array in place in order to save\n",
    "        # RAM and speed up processing\n",
    "        cloud_mask.resize(array_size)\n",
    "    except NameError:\n",
    "        # If this is the first iteration, cloud_mask will be undefined and\n",
    "        # trying to resize it will throw a NameError. This error is caught\n",
    "        # and cloud_mask is instantiated\n",
    "        cloud_mask = np.zeros(array_size, dtype=np.float64)\n",
    "\n",
    "    # Reading the reliability data\n",
    "    cloud_mask = load_cloud_mask(pr_files_day, cloud_mask)\n",
    "\n",
    "    # Preparing the reliability data\n",
    "    cloud_mask = prepare_cloud_mask(cloud_mask)\n",
    "\n",
    "    # Reading the EVI  data\n",
    "    evi = load_evi(evi_files_day, cloud_mask)\n",
    "\n",
    "    # Preparing the EVI  data\n",
    "    evi = prepare_evi(evi)\n",
    "    \n",
    "    if index.value == 'SVI':\n",
    "        update_iter_desc(days_iterator, \"Calculating SVI DOY {}\".format(day))\n",
    "        # Calculate SVI\n",
    "        svi = calculate_svi(evi)\n",
    "\n",
    "        update_iter_desc(days_iterator, \"Generating images DOY {}\".format(day))\n",
    "        # Generating png images and geotiffs\n",
    "        generate_SVI_images(evi_files_day, svi, extent, geo_transform, projection)\n",
    "\n",
    "        # Remove references to avoid array resize errors in future loops\n",
    "        del(evi, svi)\n",
    "    elif index.value == 'VCI':\n",
    "        update_iter_desc(days_iterator, \"Calculating VCI DOY {}\".format(day))\n",
    "        # Calculate VCI\n",
    "        vci = calculate_vci(evi)\n",
    "\n",
    "        update_iter_desc(days_iterator, \"Generating images DOY {}\".format(day))\n",
    "        # Generating png images and geotiffs\n",
    "        generate_VCI_images(evi_files_day, vci, extent, geo_transform, projection)\n",
    "\n",
    "        # Remove references to avoid array resize errors in future loops\n",
    "        del(evi, vci)\n",
    "    \n",
    "# create gif\n",
    "file_list = os.listdir(path_png)\n",
    "s = sorted(file_list, key = lambda x: (x.split('_')[-1],x.split('.')[0]))\n",
    "images = []\n",
    "for file_name in s:\n",
    "    if file_name.endswith('.png'):\n",
    "        file_path = os.path.join(path_png, file_name)\n",
    "        images.append(imageio.imread(file_path))\n",
    "\n",
    "gif_path = os.path.join(output_path, study_area+'_svi_animation.gif')\n",
    "imageio.mimsave(gif_path, images)\n",
    "\n",
    "# create grid image\n",
    "first_year = int(s[0].split('.')[0].split('_')[-1])\n",
    "last_year = int(s[-1].split('.')[0].split('_')[-1])\n",
    "# correlation between DOY and grid row\n",
    "row_dict = {\n",
    "                '001': 0,\n",
    "                '017': 1,\n",
    "                '033': 2,\n",
    "                '049': 3,\n",
    "                '065': 4,\n",
    "                '081': 5,\n",
    "                '097': 6,\n",
    "                '113': 7,\n",
    "                '129': 8,\n",
    "                '145': 9,\n",
    "                '161': 10,\n",
    "                '177': 11,\n",
    "                '193': 12,\n",
    "                '209': 13,\n",
    "                '225': 14,\n",
    "                '241': 15,\n",
    "                '257': 16,\n",
    "                '273': 17,\n",
    "                '289': 18,\n",
    "                '305': 19,\n",
    "                '321': 20,\n",
    "                '337': 21,\n",
    "                '353': 22,\n",
    "}\n",
    "# correlation between year and grid column\n",
    "col_dict = {}\n",
    "for year in range(first_year, last_year+1):\n",
    "    col_dict[str(year)] = year - first_year\n",
    "\n",
    "# number of rows and columns\n",
    "nrows = len(row_dict)\n",
    "ncols = last_year-first_year+1\n",
    "\n",
    "# figure\n",
    "fig = plt.figure(figsize=(ncols, nrows))\n",
    "spec = gridspec.GridSpec(nrows=nrows, ncols=ncols, figure=fig, wspace=0.05, hspace=0.05)\n",
    "# inserts png images\n",
    "i=0\n",
    "for file_name in s:\n",
    "    row = row_dict[file_name.split('_')[2]]\n",
    "    col = col_dict[file_name.split('.')[0].split('_')[-1]]\n",
    "    ax1 = fig.add_subplot(spec[row, col], frameon=False)\n",
    "    ax1.set(aspect = 'equal')\n",
    "    ax1.set_xticks([])\n",
    "    ax1.set_yticks([])\n",
    "    ax1.imshow(images[i])\n",
    "    i = i + 1\n",
    "# inserts DOYs for first column\n",
    "j = 0\n",
    "for key, value in row_dict.items():\n",
    "    ax2 = fig.add_subplot(spec[j, 0], frameon=False)\n",
    "    ax2.set_ylabel('DOY %s' % key, rotation=90)\n",
    "    ax2.set_xticks([])\n",
    "    ax2.set_yticks([])\n",
    "    j = j + 1\n",
    "j = 0\n",
    "# inserts years for first row\n",
    "for key, value in col_dict.items():\n",
    "    ax3 = fig.add_subplot(spec[0, j], frameon=False)\n",
    "    ax3.set_title(key)\n",
    "    ax3.set_xticks([])\n",
    "    ax3.set_yticks([])\n",
    "    j = j + 1\n",
    "# saves figure\n",
    "grid_path = os.path.join(output_path, study_area+'_grid_figure.png')\n",
    "plt.savefig(grid_path, dpi=400, bbox_inches='tight', edgecolor='none', facecolor='w')"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}