Use_Case.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9acfbb57-185f-49d9-8ffb-d3e0db5b8cc0",
   "metadata": {},
   "source": [
    "![](./images/Title_ex_6.6.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce8c16c3-4e2f-46e2-a704-c550b831d02a",
   "metadata": {},
   "source": [
    "# 1. Introduction\n",
    "\n",
    "In the previous exercises we explored different online data portals and analysis techniques to investigate vegetation and droughts. In this notebook we will put everything together in one use case and look at how vegetation production in a certain agricultural region is impacted by drought and how many people are likely to be directly impacted by that drought. We will be using the **Vegetation Health Index** from the [ASIS data portal](https://www.fao.org/giews/earthobservation/asis/index_1.jsp?lang=en) to look into a specific drought and its impact on vegetation.\n",
    "For more details on the methodology used to derive this index, have a look at [this paper](https://www.tandfonline.com/doi/full/10.1080/01431161.2015.1126378?scroll=top&needAccess=true).\n",
    "\n",
    "The **Vegetation Health Index (VHI)** illustrates the severity of drought based on the vegetation health and the influence of temperature on plant conditions. The VHI is a composite index, combining the Vegetation Condition Index (VCI) and the Temperature Condition Index (TCI). VCI evaluates the current vegetation health in comparison to historical trends by relating the current dekadal NDVI to its long-term minimimum and maximum value. The TCI is calculated using a similar equation to the VCI, but relates the current temperature to the long-term maximum and minimum , as it is assumed that higher temperatures tend to cause a deterioration in vegetation conditions. A decrease in the VHI would, for example, indicate relatively poor vegetation conditions and warmer temperatures, signifying stressed vegetation conditions, and over a longer period would be indicative of drought.\n",
    "\n",
    "Within ASIS, VHI is available as a dekadal, monthly and annual metric. An interesting derived product is the **Mean Vegetation Health Index (MVHI)**, available per dekad and per season. For each pixel, this MVHI represents the temporally averaged VHI starting from the start of season up to the requested dekad and provides a good indication on drought conditions during a growing season of interest.\n",
    "\n",
    "In the scope of this use case, we will only focus on one (the first) growing season. During a typical growing season, a crop is most sensitive to a drought event during the first half of the season, when the crop is fully developing. Therefore, we will zoom in on the first half of the growing season and derive the MVHI for each pixel at the moment of the peak of the season (this moment could be different for different pixels!). For each pixel, we will be needing information on the timing of the peak of season, which we will get through the WaPOR phenology products.\n",
    "\n",
    "Based on the MVHI value, the growing season can be classified into different drought severity classes:\n",
    "\n",
    "- Class 1: MVHI < 0.25    : Extreme\n",
    "- Class 2: MVHI 0.25-0.35 : Severe\n",
    "- Class 3: MVHI 0.35-0.38 : Moderate \n",
    "- Class 4: MVHI 0.38-0.42 : Mild\n",
    "- Class 5: MVHI 0.42-1    : None\n",
    "- Class 6: MVHI > 1       : Insufficient data\n",
    "\n",
    "Once we derived these drought severity classes for our area of interest, we will be checking the impact of the drought on vegetation productivity in the most affected area by inspecting the biomass productivity metrics as provided through the WaPOR portal.\n",
    "\n",
    "Finally, we will also be looking at the number of people affected by drought. For this, we can make use of the WPOP database, providing numbers of population density throughout the world and also available through the FAO API.\n",
    "\n",
    "\n",
    "\n",
    "**Throughout the notebook, we ask a series of QUESTIONS. Reflect upon these, collect your answers in one WORD document and submit to the teaching staff. This part is mandatory in order to pass the course.**\n",
    "\n",
    "\n",
    "\n",
    "Notebook outline: \n",
    "- [1. Introduction](#1.-Introduction)\n",
    "- [2. Libraries](#2.-Import-the-necessary-libraries)\n",
    "- [3. Choose a drought event and location](#3.-Choose-a-drought-event-and-location)\n",
    "- [4. Extract land cover information](#4.-Land-cover-mask)\n",
    "- [5. Extract phenology](#5.-Extract-phenology)\n",
    "- [6. Get mean VHI](#6.-Get-mean-VHI)\n",
    "- [7. Check precipitation (and evapotranspiration) patterns](#7.-Check-precipitation-(and-evapotranspiration)-patterns)\n",
    "- [8. Impact of drought on vegetation](#8.-Impact-of-drought-on-vegetation)\n",
    "- [9. How many people are affected?](#9.-How-many-people-are-affected?)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "020d081e-6526-42f7-a507-a478066a92f3",
   "metadata": {},
   "source": [
    "# 2. Import the necessary libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31cc0308-b5ee-4304-a280-17aad3133b4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "import json\n",
    "import datetime\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import geopandas as gpd\n",
    "import os \n",
    "import json\n",
    "from pathlib import Path\n",
    "import rioxarray\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "from osgeo import gdal\n",
    "import glob\n",
    "from shapely.geometry import Polygon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d6f56c6-25e1-4c9a-b29c-f8c7905aece6",
   "metadata": {},
   "source": [
    "# 3. Choose a drought event and location\n",
    "\n",
    "First we need to choose a drought event of which we want to investigate the impact from. NASA's GRACE sensor provides coarse resolution information on Surface Soil Moisture, Root Zone Soil Moisture and Groundwater levels. This information can be accessed [HERE](https://nasagrace.unl.edu/Archive.aspx) to get a first impression on when and where dry conditions are occuring.\n",
    "\n",
    "Based on this information, we have already identified a potentially interesting drought event in the 2015 growing season of the Gadaref region in Sudan. Compared to 2015, 2018 has been a relatively wet year in the same region. Throughout this exercise we will compare both periods in this region and check whether we can see any impact of the drought on vegetation productivity.\n",
    "\n",
    "To delineate the region of interest to which we are going to clip our data we can make use of the FAO Global Administrative Unit Layer (GAUL) database which provides shapefiles for each departement within all countries."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8d2e7d2-11d6-465b-9018-0e7f2d23bc83",
   "metadata": {},
   "source": [
    "The GAUL dataset can be found in the folder ./data/GAUL. Different levels of detail exist: level 0 represents the country borders, level 1 represents the departement borders and level 2 (not provided here) goes up to province level. Feel free to inspect these shapefiles in QGIS to get familiar with the data.\n",
    "\n",
    "From the level 1 dataset, we will first isolate our country of interest (Sudan)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6d875ea-fcfc-49b8-9d3b-7d2a22b8df62",
   "metadata": {},
   "outputs": [],
   "source": [
    "country = 'Sudan'\n",
    "gaul1_file = './data/GAUL/GAUL1_thinned.shp'\n",
    "\n",
    "gaul1 = gpd.read_file(gaul1_file)\n",
    "\n",
    "sudan = gaul1.loc[gaul1['ADM0_NAME'] == country]\n",
    "\n",
    "outdir = Path('./results')\n",
    "outdir.mkdir(exist_ok=True)\n",
    "\n",
    "outfile = str(outdir / 'Sudan.gpkg')\n",
    "sudan.to_file(outfile, driver='GPKG')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16c86c8c-6193-401a-a771-222550cb3287",
   "metadata": {},
   "source": [
    "Plot the country of interest and look at which departments are available:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eddb467c-6a31-4b25-93a4-7932c3a52c7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "sudan.plot()\n",
    "sudan"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4883808a-f222-49ae-a0c6-84b76cfd3362",
   "metadata": {},
   "source": [
    "Now create a new shapefile with just the departement of Gadaref. This will be our Region of Interest (ROI) for the remainder of this exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1174cf55-1afe-419d-8751-cbe48d835f4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "department = 'Gadaref'\n",
    "\n",
    "# complete the following line of code\n",
    "roi = sudan.loc[sudan[xxx] == xxx]\n",
    "\n",
    "# save result\n",
    "roifile = str(outdir / 'ROI.gpkg')\n",
    "roi.to_file(outfile, driver='GPKG')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "598db987-b4c4-40c3-b89e-8472008effee",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 4. Land cover mask\n",
    "\n",
    "In order to focus our analysis just on cropland and grassland, we need land cover information for our region of interest."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e934d2d7-3ad5-4f72-a2b6-d3f6d8d3aedf",
   "metadata": {},
   "source": [
    "## 4.1 Download data from WaPOR\n",
    "\n",
    "We will use the Copernicus Global Land Cover product available through the WaPOR data portal.\n",
    "\n",
    "Download the **100m version** of this product for our region of interest for both the wet (2018) and dry (2015) years.\n",
    "\n",
    "Make use of the **request_tif_files** function from earlier exercises, but this time use the ROI shapefile as input for the function.\n",
    "\n",
    "TIP: use \"ADM1_NAME\" as value for the \"labelcol\" argument of the function.\n",
    "\n",
    "Upload the two products to a new folder ./results/LC and name them as follows:\n",
    "- L2_LCC_dry_clipped.tif\n",
    "- L2_LCC_wet_clipped.tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fd32f62",
   "metadata": {},
   "outputs": [],
   "source": [
    "def request_tif_files(workspace, product_code, crs,\n",
    "                      outputfilename, dimensions,\n",
    "                      shapefile=None, labelcol=None,\n",
    "                      bbox=None, APIToken=None):\n",
    "    '''\n",
    "    Function to extract .tif files for one or multiple polygons from the FAO portal.\n",
    "    \n",
    "    Parameters:\n",
    "    - workspace: workspace in the FAO catalog where the product of interest is located.\n",
    "    - product_code: specific product you wish to download.\n",
    "   \n",
    "    - crs: EPSG code of projection system (int)\n",
    "    - dimensions: list of dimensions specifying the characteristics of the exact product\n",
    "      you wish to download.\n",
    "    \n",
    "    For spatial information, you can either provide a shapefile or a series of bounding boxes:\n",
    "    - shapefile: path to shapefile containing one or multiple polygons for which data\n",
    "      needs to be extracted.\n",
    "    - labelcol: this parameter is only required if \"shapefile\" is provided and represents\n",
    "      the name of the column in the shapefile which serves as a unique identifier for\n",
    "      the polygons.\n",
    "    - bbox: dictionary with keys representing label of the bbox\n",
    "        and values the actual bounding box in the CRS defined by the next parameter and\n",
    "        in the format [xmin, ymin, xmax, ymax].\n",
    "        \n",
    "    The function will return a list of download links.\n",
    "    '''\n",
    "    \n",
    "    # define necessary paths to the service\n",
    "    path_query=r'https://io.apps.fao.org/gismgr/api/v1/query/'\n",
    "    path_sign_in=r'https://io.apps.fao.org/gismgr/api/v1/iam/sign-in/'\n",
    "\n",
    "    if APIToken is None:\n",
    "      # request API token\n",
    "      APIToken = input('Your API token: ')\n",
    "    # get access token\n",
    "    resp_signin = requests.post(path_sign_in,\n",
    "                              headers={'X-GISMGR-API-KEY': APIToken})\n",
    "    resp_signin = resp_signin.json()\n",
    "    AccessToken = resp_signin['response']['accessToken']\n",
    "    \n",
    "    # get datacube measure\n",
    "    cube_url = f'https://io.apps.fao.org/gismgr/api/v1/catalog/workspaces/{workspace}/cubes/{product_code}/measures'\n",
    "    resp = requests.get(cube_url).json()\n",
    "    measure = resp['response']['items'][0]['code']\n",
    "    print('MEASURE: ', measure)\n",
    "    \n",
    "    # prepare spatial information, either from bbox or shapefile...\n",
    "    crs = f'EPSG:{crs}'\n",
    "    shapes = {}\n",
    "    if bbox is not None:\n",
    "        npoly = len(bbox.keys())\n",
    "        for name, box in bbox.items():\n",
    "            xmin, ymin, xmax, ymax = box[0], box[1], box[2], box[3]\n",
    "            shape = [\n",
    "                    [xmin,ymin],\n",
    "                    [xmin,ymax],\n",
    "                    [xmax,ymax],\n",
    "                    [xmax,ymin],\n",
    "                    [xmin,ymin]\n",
    "                    ]\n",
    "            shapes[name] = [shape]\n",
    "            \n",
    "    elif shapefile is not None:\n",
    "      polys = gpd.read_file(shapefile)\n",
    "      npoly = len(polys)\n",
    "      names = polys[labelcol].values\n",
    "      geometry = json.loads(polys.geometry.to_json())['features']\n",
    "      for k, name in enumerate(names):\n",
    "        shapes[name] = geometry[k]['geometry']['coordinates']\n",
    "\n",
    "    else:\n",
    "        raise ValueError('No spatial information provided!')\n",
    "    \n",
    "    print('Sending query for each polygon...')\n",
    "    output = []\n",
    "    i = 0    \n",
    "    for label, geom in shapes.items():\n",
    "      \n",
    "      print(f'Processing polygon {i+1} / {npoly}...')\n",
    "      \n",
    "      # build query\n",
    "      query = {\n",
    "      \"type\": \"CropRaster\",\n",
    "      \"params\": {\n",
    "        \"properties\": {\n",
    "          \"outputFileName\": outputfilename,\n",
    "          \"cutline\": True,\n",
    "          \"tiled\": True,\n",
    "          \"compressed\": True,\n",
    "          \"overviews\": True\n",
    "        },\n",
    "        \"cube\": {\n",
    "          \"code\": product_code,\n",
    "          \"workspaceCode\": workspace,\n",
    "          \"language\": \"en\"\n",
    "        },\n",
    "        \"dimensions\": dimensions,\n",
    "        \"measures\": [\n",
    "          measure\n",
    "        ],\n",
    "        \"shape\": {\n",
    "          \"type\": \"Polygon\",\n",
    "          \"properties\": {\n",
    "              \"name\": crs},\n",
    "          \"coordinates\": geom}}}\n",
    "      \n",
    "      # Post query\n",
    "      resp_query = requests.post(path_query,\n",
    "                                 headers={'Authorization':'Bearer {0}'.format(AccessToken)},\n",
    "                                 json=query)\n",
    "      resp_query = resp_query.json()\n",
    "      job_url = resp_query['response']['links'][0]['href']\n",
    "      \n",
    "      # Get results\n",
    "      j = 0\n",
    "      print('RUNNING', end=\" \")\n",
    "      while j == 0:        \n",
    "          resp = requests.get(job_url)\n",
    "          resp = resp.json()\n",
    "          if resp['response']['status']=='RUNNING':\n",
    "              print('.', end =\" \")\n",
    "          if resp['response']['status']=='COMPLETED':\n",
    "              results = resp['response']['output']\n",
    "              out = results['downloadUrl']\n",
    "              j = 1\n",
    "          if resp['response']['status']=='COMPLETED WITH ERRORS':\n",
    "              print(resp['response']['log'])\n",
    "              out = None\n",
    "              j = 1\n",
    "      \n",
    "      output.append(out)\n",
    "\n",
    "    return output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f4670be",
   "metadata": {},
   "outputs": [],
   "source": [
    "< your code here >"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "bec48d86-1524-49da-b9ec-03f311fbe485",
   "metadata": {},
   "source": [
    "QUESTION: does the land cover information available through WaPOR already reflect the difference in climatological conditions between both years? If yes, how?\n",
    "\n",
    "Hint: Open the downloaded products in QGIS and use the LC_style_CGLOPS.qml style file to properly visualize the land cover data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cd524b4-01d5-4069-9a3a-3c4dd764619a",
   "metadata": {},
   "source": [
    "## 4.2 Create the mask\n",
    "\n",
    "We only want to focus our analysis on **grassland and cropland**. Create a mask based on both land cover products containing \"1\" for pixels of interest and NaN for other pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a9931ff-2677-4914-9ba1-0378c59d284f",
   "metadata": {},
   "outputs": [],
   "source": [
    "infiles = {'wet': './results/LC/L2_LCC_wet_clipped.tif',\n",
    "           'dry': './results/LC/L2_LCC_dry_clipped.tif'}\n",
    "\n",
    "for t, infile in infiles.items():\n",
    "    lcc = rioxarray.open_rasterio(infile)\n",
    "    \n",
    "    # fill in the right values at the location of the xxx\n",
    "    \n",
    "    lcc_masked = lcc.where(lcc.isin([xxx]))\n",
    "    lcc_binary = lcc_masked.where(np.isnan(lcc_masked), 1)\n",
    "    # export final mask to tiff file\n",
    "    outfile = f'./results/LC/L2_LCC_{t}_masked.tif'\n",
    "    lcc_binary.rio.to_raster(outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7560a68-8421-4e41-b088-b3c420622d0f",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 5. Extract phenology"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33697d54-163c-4f8c-89be-e9c9ca0bb3a2",
   "metadata": {},
   "source": [
    "Phenology metrics will be extracted from the WaPOR database. \n",
    "\n",
    "For this exercise we are interested in the maximum of season (MOS), as this will define for each pixel up till when we need to compute the Mean Vegetation Health Index.\n",
    "\n",
    "Download this indicator for our region of interest and both the wet and dry year. Do this for both season 1 and 2 (S1 and S2).\n",
    "\n",
    "Make use of the **request_tif_files** function from earlier exercises, but this time use the ROI shapefile as input for the function.\n",
    "\n",
    "In the code below, make sure to complete everything indicated by xxx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbf93d5d-fc7c-40ff-af6a-22ce1e97e49b",
   "metadata": {},
   "outputs": [],
   "source": [
    "product_code = xxx\n",
    "workspace = 'WAPOR_2'\n",
    "\n",
    "years = {'dry': \"[2015-01-01,2016-01-01)\",\n",
    "         'wet': \"[2018-01-01,2019-01-01)\"}\n",
    "stages = [xxx]\n",
    "seasons = [\"S1\", \"S2\"]\n",
    "\n",
    "labelcol = 'ADM1_NAME'\n",
    "\n",
    "crs = 4326\n",
    "\n",
    "for season in seasons:\n",
    "    for name, year in years.items():\n",
    "        for stage in stages:\n",
    "            outputfilename = f'L2_PHE_{name}_{stage}_{season}.tif'\n",
    "            dimensions = [{\"code\": \"SEASON\",\n",
    "                           \"values\": [season]},\n",
    "                          {\"code\": \"STAGE\",\n",
    "                           \"values\": [stage]},\n",
    "                          {\"code\": \"YEAR\",\n",
    "                           \"values\": [year]}]\n",
    "            output = request_tif_files(workspace, product_code,\n",
    "                                       crs, outputfilename,\n",
    "                                       dimensions, shapefile=roifile,\n",
    "                                       labelcol=labelcol)\n",
    "            print(output)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d418e68-63f8-4a16-947b-6c8603b29c4c",
   "metadata": {},
   "source": [
    "Upload the downloaded files into a designated folder *phenology* in your results directory (./results/phenology)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d41b6467",
   "metadata": {},
   "source": [
    "## 5.2 Check results visually and extract relevant information\n",
    "\n",
    "Now visualize MOS for the dry season, starting with S2.\n",
    "\n",
    "Note that MOS is expressed in dekad number, counting from the year before the year of interest, meaning:\n",
    "- dekad 1-36 is from the year before the year of interest \n",
    "- dekad 37-72 is from the year of interest \n",
    "- dekad 72-108 is from the year after the year of interest\n",
    "\n",
    "This also means that any number > 108 should be masked out!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70d0f12c-1a9d-4c10-9402-97766bc3ba16",
   "metadata": {},
   "outputs": [],
   "source": [
    "phe_dry_mos_img_s2 = rioxarray.open_rasterio(str(outdir / 'phenology' / 'L2_PHE_dry_MOS_S2.tif'))\n",
    "phe_dry_mos_img_s2 = phe_dry_mos_img_s2.where(phe_dry_mos_img_s2 <= 108)\n",
    "phe_dry_mos_img_s2.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b62bb14",
   "metadata": {},
   "source": [
    "As you can see in the result above, there are few pixels in the study area having a second growing season.\n",
    "In order to not make this exercise overly complex, we will ignore pixels with a second growing season.\n",
    "Mask these pixels in the S1 MOS product:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fbe4ceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "phe_dry_mos_img_s1 = rioxarray.open_rasterio(str(outdir / 'phenology' / 'L2_PHE_dry_MOS_S1.tif'))\n",
    "phe_dry_mos_img_s1 = phe_dry_mos_img_s1.where(phe_dry_mos_img_s1 <= 108)\n",
    "\n",
    "# mask phe_dry_mos_img_s1 using phe_dry_mos_img_s2\n",
    "phe_dry_mos_img_s1 = phe_dry_mos_img_s1.where(np.isnan(phe_dry_mos_img_s2))\n",
    "phe_dry_mos_img_s1.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "724a45a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check minimum and maximum values of MOS in the dry year:\n",
    "print(np.nanmin(phe_dry_mos_img_s1.values))\n",
    "print(np.nanmax(phe_dry_mos_img_s1.values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a37c3311",
   "metadata": {},
   "outputs": [],
   "source": [
    "# note that some pixels have a peak of season which is outside of our year of interest (MOS < 36 or MOS > 72)\n",
    "# for safety, we will exclude these pixels and also exclude pixels with MOS very early in the year from our analysis\n",
    "phe_dry_mos_img_s1 = phe_dry_mos_img_s1.where(xr.ufuncs.logical_and(phe_dry_mos_img_s1 > 45, phe_dry_mos_img_s1 <= 72))\n",
    "phe_dry_mos_img_s1.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dc3a48a",
   "metadata": {},
   "source": [
    "Now repeat these steps for the wet year as well..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1394a16e",
   "metadata": {},
   "outputs": [],
   "source": [
    "< your code here >"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecabb288",
   "metadata": {},
   "source": [
    "QUESTION: can you already see a difference between the wet and dry year in terms of timing of MOS, caused by the drought? Please describe!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "182850b6-1a3a-4e21-8ba0-147ae41b7994",
   "metadata": {},
   "source": [
    "## 5.3 Update LC mask based on phenology\n",
    "\n",
    "Extend the LC mask with the pixels masked out in the previous step of the analysis (i.e. pixels having two growing seasons and/or no clear seasonality in our year of interest)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db61f68e",
   "metadata": {},
   "outputs": [],
   "source": [
    "lcc_mask_dry = rioxarray.open_rasterio(str(outdir / 'LC' / 'L2_LCC_dry_masked.tif'))\n",
    "lcc_mask_dry = lcc_mask_dry.where(~np.isnan(phe_dry_mos_img_s1.values))\n",
    "lcc_mask_dry.plot()\n",
    "plt.show()\n",
    "outfile = str(outdir / 'LC' / 'L2_LCC_dry_masked_pheno.tif')\n",
    "lcc_mask_dry.rio.to_raster(outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "040cac11",
   "metadata": {},
   "source": [
    "Do the same for the wet year..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14c8389b",
   "metadata": {},
   "outputs": [],
   "source": [
    "< your code here >"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcc14211-0b74-4257-aeaf-314e13322040",
   "metadata": {},
   "source": [
    "## 5.3 Apply the updated mask to all phenology layers\n",
    "\n",
    "Now we apply the extended land cover mask to the MOS data for season 1 in order to make sure pixels which we are not interested in, will not be processed in the future...\n",
    "\n",
    "Meanwhile, we also make sure the phenology values range from 1 - 36 as we only focus on our year of interest...\n",
    "You can make use of the function below to adjust the MOS of Season 1 for both the wet and dry years."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9c55fdd-7dd9-4704-98c6-b2407398a2ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mask_phenology(input_path, mask):\n",
    "    '''\n",
    "    Little function to execute a recurrent operation faster\n",
    "    '''\n",
    "    \n",
    "    src = rioxarray.open_rasterio(input_path)\n",
    "    src = src.where(~np.isnan(mask.values))\n",
    "    src = src - 36\n",
    "    output_path = input_path[:-4] + '_masked.tif'\n",
    "    src.rio.to_raster(output_path)\n",
    "    \n",
    "    return src"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09b1eda9-502e-4000-8209-8fca66a16e7f",
   "metadata": {},
   "source": [
    "Now apply that function to both the dry and wet years:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0f36463-854c-41af-8884-57cd5a6abc9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "infile = str(outdir / 'phenology' / 'L2_PHE_dry_MOS_S1.tif')\n",
    "# apply the mask_phenology function here...\n",
    "phe_dry_mos = xxx\n",
    "\n",
    "infile = str(outdir / 'phenology' / 'L2_PHE_wet_MOS_S1.tif')\n",
    "# apply the mask_phenology function here...\n",
    "phe_wet_mos = xxx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0837030-bfe9-4f9c-87ee-3dc7fffa8df1",
   "metadata": {},
   "source": [
    "# 6. Get mean VHI\n",
    "\n",
    "We will now extract the Mean VHI for both years for our region of interest through the FAO map catalog..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c8a9a39-66e6-458d-af04-3933420fde20",
   "metadata": {},
   "source": [
    "## 6.1 Extract the mean VHI dekadal from the FAO map catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3700f795-7b5a-4ff0-bf33-5fe4765cc036",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "workspace = 'ASIS'\n",
    "overview = False\n",
    "paged = False\n",
    "cubes_url = f'https://io.apps.fao.org/gismgr/api/v1/catalog/workspaces/{workspace}/cubes?overview={overview}&paged={paged}'\n",
    "resp = requests.get(cubes_url).json()\n",
    "catalog = pd.DataFrame.from_dict(resp['response'])\n",
    "catalog"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac40a121-8b75-4e2e-9eb0-95663b762b38",
   "metadata": {},
   "source": [
    "We are interested in the Mean Vegetation Health Index on dekadal scale. The mean VHI averages the VHI from start of the season 1 to the dekad given in the query payload. In order to download the right data we need to know the range of MOS in both the wet and dry years..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70069417-1894-4dd3-a2e8-6a9a499450ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('MOS min dry: ', np.nanmin(phe_dry_mos))\n",
    "print('MOS max dry: ', np.nanmax(phe_dry_mos))\n",
    "print('MOS min wet: ', np.nanmin(phe_wet_mos))\n",
    "print('MOS max wet: ', np.nanmax(phe_wet_mos))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08dcc998-bd1a-4470-a1a8-c1f90c26505f",
   "metadata": {},
   "source": [
    "So basically we need to dowload data from dekad 10 onwards (month april) up till the end of the year...\n",
    "\n",
    "In order to link the dekad number with the files to download, we can prepare the following list, thereby making sure the dekad is presented in the format expected by the API:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c373fc38-c238-4200-9a26-a0d26de79151",
   "metadata": {},
   "outputs": [],
   "source": [
    "# dekad numbers range between 10 and 36\n",
    "dekad_nr = list(range(10, 37))\n",
    "years = [2015, 2018]\n",
    "dekad_df = []\n",
    "\n",
    "for year in years:\n",
    "    dekad_list = []\n",
    "    for month in range(4, 13):\n",
    "        if month < 9:\n",
    "            dekad_list.append(f'[{year}-0{month}-01,{year}-0{month}-11)')\n",
    "            dekad_list.append(f'[{year}-0{month}-11,{year}-0{month}-21)')\n",
    "            dekad_list.append(f'[{year}-0{month}-21,{year}-0{month+1}-01)')\n",
    "        elif month == 9:\n",
    "            dekad_list.append(f'[{year}-0{month}-01,{year}-0{month}-11)')\n",
    "            dekad_list.append(f'[{year}-0{month}-11,{year}-0{month}-21)')\n",
    "            dekad_list.append(f'[{year}-0{month}-21,{year}-{month+1}-01)')\n",
    "        elif month == 12:\n",
    "            dekad_list.append(f'[{year}-{month}-01,{year}-{month}-11)')\n",
    "            dekad_list.append(f'[{year}-{month}-11,{year}-{month}-21)')\n",
    "            dekad_list.append(f'[{year}-{month}-21,{year+1}-0{month-11}-01)')\n",
    "        else:\n",
    "            dekad_list.append(f'[{year}-{month}-01,{year}-{month}-11)')\n",
    "            dekad_list.append(f'[{year}-{month}-11,{year}-{month}-21)')\n",
    "            dekad_list.append(f'[{year}-{month}-21,{year}-{month+1}-01)')\n",
    "            \n",
    "    intermediate_dict = {'dekad_no': dekad_nr, 'dekad': dekad_list,\n",
    "                         'year': [year] * len(dekad_nr)}\n",
    "    dekad_df.append(pd.DataFrame(intermediate_dict))\n",
    "\n",
    "dekad_df = pd.concat(dekad_df, axis=0, ignore_index=True)\n",
    "dekad_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c6a91a2-0164-4ea2-8b93-3f3fe03a9d7a",
   "metadata": {},
   "source": [
    "The following cell demonstrates how the download of the MVHI files was done.\n",
    "As it takes a long time to download all, we have prepared the data for you in the ./data/mvhi folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b94809fa-b9a9-4bc1-9c2d-551ecd64e446",
   "metadata": {},
   "outputs": [],
   "source": [
    "product_code = 'MVHI_D'\n",
    "workspace = 'ASIS'\n",
    "\n",
    "season = \"S1\"\n",
    "land = \"C\"  \n",
    "\n",
    "crs = 4326\n",
    "\n",
    "api_token = '4c95eef5536df1786815b04ed92f8a9e7a59bb20976e7e832db2eeef8c4cafd53b0b74d222614556'\n",
    "\n",
    "for i, row in dekad_df.iterrows():\n",
    "    \n",
    "    year = row.year\n",
    "    dekad_no = row.dekad_no\n",
    "    dekad = row.dekad\n",
    "    \n",
    "    dimensions = [{\"code\": \"DEKAD\",\n",
    "                   \"values\": [dekad]},\n",
    "                  { \"code\": \"SEASON\",\n",
    "                   \"values\": [season]},\n",
    "                  {\"code\": \"LAND\",\n",
    "                   \"values\": [land]}]\n",
    "    \n",
    "    # indicate how you want to name your rasterfile\n",
    "    outputfilename = f'MVHI_D_{year}_{dekad_no}.tif' \n",
    "    \n",
    "    output = request_tif_files(workspace, product_code, crs,\n",
    "                      outputfilename, dimensions,\n",
    "                      shapefile=roifile, labelcol=labelcol,\n",
    "                      APIToken=api_token)\n",
    "    print(output)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d26b73a7-51ea-435f-a402-a72260b2578c",
   "metadata": {},
   "source": [
    "Before we can work properly with the MVHI maps, we need to warp them so their resolution matches with the LC mask."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d776f7d-b872-43bf-bb57-81ac1e00ae7f",
   "metadata": {},
   "source": [
    "## 6.2 Warp MVHI\n",
    "\n",
    "Again to optimize the workflow and automate the warping process a for loop is applied over all the files present in the inputfolder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c68f1283-5bf8-4aeb-afdf-a338c6f0fcfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define input_folder\n",
    "input_folder = Path(r'./data/mvhi' )\n",
    "# get list of tif files in the input_folder\n",
    "input_fhs = sorted(glob.glob(str(input_folder / '*.tif'))) \n",
    "\n",
    "# define the file of which we want the projection/resolution (destination info):\n",
    "dst_info = gdal.Info(gdal.Open(r'./results/phenology/L2_PHE_dry_MOS_S1_masked.tif'), format='json')\n",
    "\n",
    "# define output folder:\n",
    "output_folder = Path(r'./results/mvhi/warped')\n",
    "# create output_folder if it does not exist\n",
    "if not os.path.exists(output_folder): \n",
    "    os.makedirs(output_folder)\n",
    "    \n",
    "# loop over all files:\n",
    "for input_fh in input_fhs:\n",
    "    src_info = gdal.Info(gdal.Open(input_fh), format='json')\n",
    "    filename = os.path.basename(input_fh)\n",
    "    output_file = str(output_folder / f'{filename[:-4]}_warped.tif')\n",
    "    \n",
    "    output = gdal.Warp(output_file, input_fh, format='GTiff',\n",
    "              srcSRS=src_info['coordinateSystem']['wkt'],\n",
    "              dstSRS=dst_info['coordinateSystem']['wkt'],\n",
    "              width=dst_info['size'][0],\n",
    "              height=dst_info['size'][1],\n",
    "              outputBounds=(dst_info['cornerCoordinates']['lowerLeft'][0],\n",
    "                            dst_info['cornerCoordinates']['lowerLeft'][1],\n",
    "                            dst_info['cornerCoordinates']['upperRight'][0],\n",
    "                            dst_info['cornerCoordinates']['upperRight'][1]),\n",
    "              outputBoundsSRS=dst_info['coordinateSystem']['wkt'],\n",
    "              resampleAlg='near')\n",
    "    # (use the following line to make sure the output is written to disk correctly)\n",
    "    output = None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46cb612b-660d-4490-bbd2-c62a345c15bb",
   "metadata": {},
   "source": [
    "## 6.3 Compute MVHI for the region of interest\n",
    "\n",
    "Now we want for each pixel the Mean VHI over the first half of the growing season. We simply need for each pixel the value of the MVHI which corresponds with the dekad at which the season peaks. This can be achieved by looping over the different dekads and assigning the right values to the right pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9491138-e738-4c8e-8b81-3b6229f2c499",
   "metadata": {},
   "outputs": [],
   "source": [
    "def produce_mvhi(mos_ds, mvhi_dir, year, outfile):\n",
    "    \n",
    "    # define final result\n",
    "    output = xr.zeros_like(mos_ds)\n",
    "    \n",
    "    # check which unique MOS values occur, ignoring NaN\n",
    "    mos_val = list(np.unique(mos_ds.values[~np.isnan(mos_ds.values)]))\n",
    "    \n",
    "    for mos in mos_val:\n",
    "        mos = int(mos)\n",
    "        mvhi_mos = rioxarray.open_rasterio(str(mvhi_dir / f'MVHI_D_{year}_{mos}_warped.tif'))\n",
    "        # get MVHI values where phe_dry_mos = dekad\n",
    "        mask = np.where(mos_ds != mos, 0, mvhi_mos)\n",
    "        # add to final result\n",
    "        output = output + mask\n",
    "    \n",
    "    # make sure result contains NaN where MOS == np.NaN\n",
    "    output = output.where(~np.isnan(mos_ds))\n",
    "    # save to .tif file\n",
    "    output.rio.to_raster(outfile)\n",
    "    \n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "158cfc92-78c2-43dc-a369-6e112c54effa",
   "metadata": {},
   "source": [
    "Apply function to wet and dry year:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f9eeb1b-b6e8-4ca1-a3f1-146642cf23b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "outfile = str(output_folder / 'Final_mvhi_dry.tif')\n",
    "# apply the produce_mvhi function here...\n",
    "mvhi_dry = xxx\n",
    "\n",
    "outfile = str(output_folder / 'Final_mvhi_wet.tif')\n",
    "# apply the produce_mvhi function here...\n",
    "mvhi_wet = xxx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07e3b46b-9d7a-4b35-8ceb-3f9b93d1adb9",
   "metadata": {},
   "source": [
    "## 6.4 Classify MVHI values into drought severity classes\n",
    "\n",
    "Our map represents continuous values of the Mean VHI. To simplify and convert them to a drought intensity map, we can reclass the array. The classes are based on [this FAO website](https://www.fao.org/giews/earthobservation/asis/index_1.jsp?lang=en).\n",
    "\n",
    "- Class 1: MVHI < 0.25    : Extreme\n",
    "- Class 2: MVHI 0.25-0.35 : Severe\n",
    "- Class 3: MVHI 0.35-0.38 : Moderate \n",
    "- Class 4: MVHI 0.38-0.42 : Mild\n",
    "- Class 5: MVHI 0.42-1    : None\n",
    "- Class 6: MVHI > 1       : Insufficient data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4754b4c4-f8fc-4987-9e1c-b855695b0d9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def classify_drought(mvhi, outfile):\n",
    "    \n",
    "    # define output\n",
    "    output = xr.zeros_like(mvhi)\n",
    "    class1 = np.where((mvhi < 0.25) & (mvhi > 0), 1, output)\n",
    "    class2 = np.where((mvhi < 0.35) & (mvhi >= 0.25), 2, output)\n",
    "    class3 = np.where((mvhi < 0.38) & (mvhi >= 0.35), 3, output) \n",
    "    class4 = np.where((mvhi < 0.42) & (mvhi >= 0.38), 4, output) \n",
    "    class5 = np.where((mvhi <= 1) & (mvhi >= 0.42), 5, output) \n",
    "    class6 = np.where(mvhi > 1, 6, output) \n",
    "    \n",
    "    output.values = class1 + class2 + class3 + class4 + class5 + class6\n",
    "    output = output.where(output != 0)\n",
    "\n",
    "    output.rio.to_raster(outfile)\n",
    "    \n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e63de163-19e3-43d9-b188-7e59515beb39",
   "metadata": {},
   "source": [
    "Apply this function for both wet and dry years and plot the results...\n",
    "\n",
    "QUESTION: Based on the figures you produce, where in the Gadaref province do you expect to see the highest impact of the drought on the vegetation? Please describe. Add the figures of MVHI for both the wet and dry years in your report!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0270c17-da79-4560-871e-2cb3a407e5dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "drought_dir = outdir / 'drought_severity'\n",
    "drought_dir.mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "outfile = str(drought_dir / 'drought_wet.tif')\n",
    "# apply the classify_drought function here...\n",
    "drought_wet = xxx\n",
    "drought_wet.plot(cmap=\"RdBu\")\n",
    "plt.title(\"Drought Intenisty Classes Wet Year\")\n",
    "plt.ylabel(\"latitude\")\n",
    "plt.xlabel(\"longitude\")\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "outfile = str(drought_dir / 'drought_dry.tif')\n",
    "# apply the classify_drought function here...\n",
    "drought_dry = xxx\n",
    "drought_dry.plot(cmap=\"RdBu\")\n",
    "plt.title(\"Drought Intenisty Classes Dry Year\")\n",
    "plt.ylabel(\"latitude\")\n",
    "plt.xlabel(\"longitude\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc80cb82-b218-4120-a15f-1158e36c5a63",
   "metadata": {},
   "source": [
    "# 7. Check precipitation (and evapotranspiration) patterns\n",
    "\n",
    "QUESTION: Request precipitation and evapotranspiration data at 250 m resolution from the WaPOR data portal for a pixel clearly affected by the drought (based on the drought intensity maps you just produced). Compare the patterns between the dry and wet year. Make sure to add the figures to your report.\n",
    "\n",
    "Make use of the **get_pixel_timeseries** function adopted from an earlier exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3260718b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_pixel_timeseries(workspace, product_code,\n",
    "                         coordinates, crs,\n",
    "                         start_date, end_date):\n",
    "    '''\n",
    "    Function to extract one or multiple pixel time series from the FAO portal.\n",
    "    Parameters:\n",
    "    - workspace: workspace in the FAO catalog where the product of interest is located.\n",
    "    - product_code: specific product you wish to download.\n",
    "    - coordinates: dictionary with keys representing label of the point\n",
    "        and values the x,y coordinate pair in the CRS defined by the next parameter.\n",
    "        In case of Geographic CRS, the order should be [lat, lon].\n",
    "    - crs: EPSG code of projection system (int)\n",
    "    - start_date: start date of the time series.\n",
    "        Should be str in the following format: \"yyyy-mm-dd\".\n",
    "    - end_date: end date of the time series.\n",
    "        Should be str in the following format: \"yyyy-mm-dd\".\n",
    "        \n",
    "    The function will return a dataframe with time as index and each column\n",
    "    representing the query result for an individual point.\n",
    "    '''\n",
    "    \n",
    "    print('Making some preparations...')\n",
    "    \n",
    "    path_query = r'https://io.apps.fao.org/gismgr/api/v1/query/'\n",
    "    \n",
    "    crs = f'EPSG:{crs}'\n",
    "    \n",
    "    # get datacube measure\n",
    "    cube_url = f'https://io.apps.fao.org/gismgr/api/v1/catalog/workspaces/{workspace}/cubes/{product_code}/measures'\n",
    "    resp = requests.get(cube_url).json()\n",
    "    measure = resp['response']['items'][0]['code']\n",
    "    print('MEASURE: ', measure)\n",
    "    \n",
    "    # get datacube time dimension\n",
    "    cube_url = f'https://io.apps.fao.org/gismgr/api/v1/catalog/workspaces/{workspace}/cubes/{product_code}/dimensions'\n",
    "    resp = requests.get(cube_url).json()\n",
    "    items = pd.DataFrame.from_dict(resp['response']['items'])\n",
    "    dimension = items[items.type=='TIME']['code'].values[0]\n",
    "    print('DIMENSION: ', dimension)\n",
    "    \n",
    "    print('Sending query for each point...')\n",
    "    npoints = len(coordinates.keys())\n",
    "    i = 0    \n",
    "    for label, coo in coordinates.items():\n",
    "        \n",
    "        print(f'Processing point {i+1} / {npoints}...')\n",
    "        query = {\n",
    "        \"type\": \"PixelTimeSeries\",\n",
    "        \"params\": {\n",
    "            \"cube\": {\n",
    "            \"code\": product_code,\n",
    "            \"workspaceCode\": workspace,\n",
    "            \"language\": \"en\"\n",
    "            },\n",
    "            \"dimensions\": [\n",
    "            {\n",
    "                \"code\": dimension,\n",
    "                \"range\": f\"[{start_date},{end_date})\"\n",
    "            }\n",
    "            ],\n",
    "            \"measures\": [\n",
    "            measure\n",
    "            ],\n",
    "            \"point\": {\n",
    "            \"crs\": crs,            \n",
    "            \"x\":coo[1],\n",
    "            \"y\":coo[0]\n",
    "            }\n",
    "        }\n",
    "        }\n",
    "        \n",
    "        resp_query = requests.post(path_query, json=query)\n",
    "        resp = resp_query.json()\n",
    "        results = resp['response']\n",
    "        df = pd.DataFrame(results['items'], columns=['Time', label])\n",
    "        df.index = pd.to_datetime(df['Time'], format='%Y-%m-%d')\n",
    "        df = df.drop(columns=['Time'])\n",
    "        if i == 0:\n",
    "            df_fin = df\n",
    "        else:\n",
    "            df_fin = pd.merge(df_fin, df, left_index=True, right_index=True)  \n",
    "        \n",
    "        i += 1\n",
    "        print('Query successful!')\n",
    "        \n",
    "    return df_fin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da83e432",
   "metadata": {},
   "outputs": [],
   "source": [
    "# DRY YEAR \n",
    "\n",
    "# PRECIPITATION\n",
    "crs = 4326 \n",
    "coordinates = {'POI': xxx}\n",
    "\n",
    "product_code = xxx\n",
    "workspace = 'WAPOR_2'\n",
    "\n",
    "start_date = \"2015-01-01\"\n",
    "end_date = \"2015-12-31\"\n",
    "\n",
    "pcp_dry = get_pixel_timeseries(workspace, product_code,\n",
    "                         coordinates, crs,\n",
    "                         start_date, end_date)\n",
    "pcp_dry = pcp_dry.rename(columns={'POI': 'pcp'})\n",
    "\n",
    "\n",
    "# EVAPOTRANSPIRATION\n",
    "product_code = xxx\n",
    "\n",
    "aeti_dry = get_pixel_timeseries(workspace, product_code,\n",
    "                         coordinates, crs,\n",
    "                         start_date, end_date)\n",
    "aeti_dry = aeti_dry.rename(columns={'POI': 'aeti'})\n",
    "\n",
    "data_dry = pd.merge(pcp_dry, aeti_dry, right_index=True, left_index=True)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "data_dry.plot(ax=ax)\n",
    "ax.set_title('Dry year')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ebf52e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# WET YEAR \n",
    "\n",
    "< your code here >"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7214d68-9a35-47f3-900f-b65e88497734",
   "metadata": {},
   "source": [
    "# 8. Impact of drought on vegetation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2e08c96-4f71-4d1b-9492-98fc5f8ad142",
   "metadata": {},
   "source": [
    "## 8.1 Net Primary Productivity throughout the season"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21bbd8aa-2be1-4e7e-b3fa-a8dc3bd9671c",
   "metadata": {},
   "source": [
    "Retrieve the NPP timeseries at 250m resolution from the WaPOR data portal for the given location, for both the wet and dry years.\n",
    "\n",
    "Location: [13.4712, 35.3531] (Lat, Lon)\n",
    "\n",
    "Make use of the **get_pixel_timeseries** function from an earlier exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fadb90c-e573-4184-bf35-272cd6c6a950",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify parameters for API query\n",
    "point = {'POI': xxx} \n",
    "crs = 4326\n",
    "product_code = xxx\n",
    "workspace = 'WAPOR_2'\n",
    "start_date_dry = \"2015-01-01\"\n",
    "end_date_dry = \"2015-12-31\"\n",
    "start_date_wet = \"2018-01-01\"\n",
    "end_date_wet = \"2018-12-31\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c499e066-a03e-4fea-8191-6f1ac6ffa95d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dry year\n",
    "df_dry = get_pixel_timeseries(workspace, product_code,\n",
    "                              point, crs,\n",
    "                              start_date_dry, end_date_dry)\n",
    "# add dekad number\n",
    "df_dry['dekad'] = range(1, len(df_dry)+1)\n",
    "\n",
    "# Wet year\n",
    "df_wet = xxx\n",
    "\n",
    "# add dekad number\n",
    "df_wet['dekad'] = range(1, len(df_wet)+1)\n",
    "\n",
    "npp_df = pd.merge(df_dry, df_wet, left_on='dekad', right_on='dekad',\n",
    "                  suffixes=('_dry', '_wet'))\n",
    "npp_df.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c170a0db",
   "metadata": {},
   "source": [
    "Now plot both timeseries in one graph.\n",
    "\n",
    "QUESTION: Discuss how the drought had an impact on vegetation productivity for this specific location. Please include the figure in your answer to this question, as well as the product_code you used for the query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3ca3433-1a65-4add-a570-4ed960efce81",
   "metadata": {},
   "outputs": [],
   "source": [
    "figs, ax = plt.subplots(1, figsize=(12, 6))\n",
    "\n",
    "npp_df['test_dry'].plot(ax=ax,label='Dry', marker='o')\n",
    "npp_df['test_wet'].plot(ax=ax,label='Wet', marker='o')\n",
    "\n",
    "ax.set_ylabel('NPP [gC/m2/day]')\n",
    "ax.legend(loc='upper right', ncol=5)\n",
    "ax.set_title('NPP time series example for dry and wet year')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7457749a-14d2-44c1-b0bd-03aa40f140f7",
   "metadata": {},
   "source": [
    "## 8.2 Total Biomassa Production during the season\n",
    "\n",
    "Let's now check the impact on total biomass production across the entire FIRST season, again for the location specified above.\n",
    "Retrieve the Total Biomass Production product from the WaPOR data portal, this time at 100 m resolution\n",
    "\n",
    "Since we now need a seasonal product, we need to adjust the **get_pixel_timeseries** function from the earlier exercise a bit..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b389733",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_pixel_seasonal_value(workspace, product_code,\n",
    "                         coordinates, crs,\n",
    "                         dimensions):\n",
    "    '''\n",
    "    Function to extract one or multiple pixel time series from the FAO portal.\n",
    "    Parameters:\n",
    "    - workspace: workspace in the FAO catalog where the product of interest is located.\n",
    "    - product_code: specific product you wish to download.\n",
    "    - coordinates: dictionary with keys representing label of the point\n",
    "        and values the x,y coordinate pair in the CRS defined by the next parameter.\n",
    "        In case of Geographic CRS, the order should be [lat, lon].\n",
    "    - crs: EPSG code of projection system (int)\n",
    "    - dimensions: specifies for which year, season etc. the product needs to be downloaded\n",
    "        \n",
    "    The function will return a dataframe with time as index and each column\n",
    "    representing the query result for an individual point.\n",
    "    '''\n",
    "    \n",
    "    print('Making some preparations...')\n",
    "    \n",
    "    path_query = r'https://io.apps.fao.org/gismgr/api/v1/query/'\n",
    "    \n",
    "    crs = f'EPSG:{crs}'\n",
    "    \n",
    "    # get datacube measure\n",
    "    cube_url = f'https://io.apps.fao.org/gismgr/api/v1/catalog/workspaces/{workspace}/cubes/{product_code}/measures'\n",
    "    resp = requests.get(cube_url).json()\n",
    "    measure = resp['response']['items'][0]['code']\n",
    "    print('MEASURE: ', measure)\n",
    "    \n",
    "    # get datacube time dimension\n",
    "    cube_url = f'https://io.apps.fao.org/gismgr/api/v1/catalog/workspaces/{workspace}/cubes/{product_code}/dimensions'\n",
    "    resp = requests.get(cube_url).json()\n",
    "    items = pd.DataFrame.from_dict(resp['response']['items'])\n",
    "    dimension = items[items.type=='TIME']['code'].values[0]\n",
    "    print('DIMENSION: ', dimension)\n",
    "    \n",
    "    print('Sending query for each point...')\n",
    "    npoints = len(coordinates.keys())\n",
    "    i = 0    \n",
    "    for label, coo in coordinates.items():\n",
    "        \n",
    "        print(f'Processing point {i+1} / {npoints}...')\n",
    "        query = {\n",
    "        \"type\": \"PixelTimeSeries\",\n",
    "        \"params\": {\n",
    "            \"cube\": {\n",
    "            \"code\": product_code,\n",
    "            \"workspaceCode\": workspace,\n",
    "            \"language\": \"en\"\n",
    "            },\n",
    "            \"dimensions\": dimensions,\n",
    "            \"measures\": [\n",
    "            measure\n",
    "            ],\n",
    "            \"point\": {\n",
    "            \"crs\": crs,            \n",
    "            \"x\":coo[1],\n",
    "            \"y\":coo[0]\n",
    "            }\n",
    "        }\n",
    "        }\n",
    "        \n",
    "        resp_query = requests.post(path_query, json=query)\n",
    "        resp = resp_query.json()\n",
    "        results = resp['response']\n",
    "        df = pd.DataFrame(results['items'], columns=['Time', label])\n",
    "        df.index = pd.to_datetime(df['Time'], format='%Y-%m-%d')\n",
    "        df = df.drop(columns=['Time'])\n",
    "        if i == 0:\n",
    "            df_fin = df\n",
    "        else:\n",
    "            df_fin = pd.merge(df_fin, df, left_index=True, right_index=True)  \n",
    "        \n",
    "        i += 1\n",
    "        print('Query successful!')\n",
    "        \n",
    "    return df_fin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f11d7757-a35f-46eb-80c4-e1db9bbb760e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify parameters for API query\n",
    "point = {'POI': xxx} \n",
    "crs = 4326\n",
    "product_code = xxx\n",
    "workspace = 'WAPOR_2'\n",
    "year_dry = '[2015-01-01,2016-01-01)'\n",
    "year_wet = '[2018-01-01,2019-01-01)'\n",
    "season = xxx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f463f228",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dry year\n",
    "dimensions = [{\"code\": \"YEAR\",\n",
    "               \"values\": [year_dry]},\n",
    "              {\"code\": \"SEASON\",\n",
    "               \"values\": [season]}]\n",
    "tbp_dry = get_pixel_seasonal_value(workspace, product_code,\n",
    "                              point, crs,\n",
    "                              dimensions)\n",
    "\n",
    "# Wet year\n",
    "dimensions = [{\"code\": \"YEAR\",\n",
    "               \"values\": [year_wet]},\n",
    "              {\"code\": \"SEASON\",\n",
    "               \"values\": [season]}]\n",
    "tbp_wet = get_pixel_seasonal_value(workspace, product_code,\n",
    "                              point, crs,\n",
    "                              dimensions)\n",
    "\n",
    "TBP_df = pd.concat([tbp_dry, tbp_wet], ignore_index=False)\n",
    "TBP_df.index = [\"Dry\", \"Wet\"]\n",
    "TBP_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f07fce5",
   "metadata": {},
   "source": [
    "Now plot the values in one bar graph for quick comparison.\n",
    "\n",
    "QUESTION: Does this result confirm your earlier conclusion based on NPP? Does this make sense, given how total biomass production in WaPOR is derived? Please include the figure in your answer to this question, as well as the product_code you used in the query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e773957-4a69-4911-bcc3-17a46a8c2d90",
   "metadata": {},
   "outputs": [],
   "source": [
    "figs, ax = plt.subplots(1, figsize=(12, 6))\n",
    "TBP_df['POI'].plot.bar(ax=ax, label = 'TBP')\n",
    "ax.set_ylabel('TBP [kg/ha]')\n",
    "ax.legend(loc='upper right', ncol=5)\n",
    "ax.set_title('TBP druing season 1 in a dry and wet year')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5d30834-eae3-4e6f-a164-f82a7db6dc1a",
   "metadata": {},
   "source": [
    "# 9. How many people are affected by the drought? "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8183ca34-4b3b-46c4-9043-5beee28e7cd6",
   "metadata": {},
   "source": [
    "## 9.1 Get population count map\n",
    "\n",
    "From the FAO catalog, retrieve the population count for our region of interest and for the dry year only.\n",
    "\n",
    "Population count data is included in the \"WPOP\" dataset.\n",
    "\n",
    "Make use of the **request_tif_files** function from earlier exercise, but this time use the ROI shapefile as input for the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13bd7190-c8e0-4968-9380-dcef72e89873",
   "metadata": {},
   "outputs": [],
   "source": [
    "workspace = 'WPOP'\n",
    "overview = False\n",
    "paged = False\n",
    "cubes_url = f'https://io.apps.fao.org/gismgr/api/v1/catalog/workspaces/{workspace}/cubes?overview={overview}&paged={paged}'\n",
    "resp = requests.get(cubes_url).json()\n",
    "catalog = pd.DataFrame.from_dict(resp['response'])\n",
    "catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9915884-6a40-4de3-920e-02399e1550cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "product_code = xxx\n",
    "\n",
    "# Specify years\n",
    "years = {'dry': '[2015-01-01,2016-01-01)'}\n",
    "\n",
    "# specify output directory\n",
    "wpop_dir = outdir / 'wpop'\n",
    "wpop_dir.mkdir(exist_ok=True, parents=True)\n",
    "\n",
    "for label, year in years.items():\n",
    "    dimensions = [{\"code\": \"YEAR\",\n",
    "               \"values\": [year]}]\n",
    "    outputfilename = f'WPOP_count_{label}.tif'\n",
    "    output = request_tif_files(workspace, product_code,\n",
    "                               crs, outputfilename,\n",
    "                               dimensions, shapefile=roifile,\n",
    "                               labelcol=labelcol)\n",
    "    print(output)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "79a84658-f64e-4d67-835e-1328b39f4ac2",
   "metadata": {},
   "source": [
    "Download the files and upload into a wpop folder in your output directory."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e2e5d95-3b78-4986-9978-d4263ba71f88",
   "metadata": {},
   "source": [
    "## 9.2 Warp the population data to the resolution of our drought analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b68062a-3dbc-4a7a-884f-158016dd09dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the file of which we want the projection/resolution (destination info):\n",
    "dst_info = gdal.Info(gdal.Open(r'./results/mvhi/warped/Final_mvhi_dry.tif'), format='json')\n",
    "\n",
    "wpopfile = str(wpop_dir / 'WPOP_count_dry.tif')\n",
    "src_info = gdal.Info(gdal.Open(wpopfile), format='json')\n",
    "filename = os.path.basename(wpopfile)\n",
    "output_file = str(wpop_dir / (filename[:-4] + '_warped.tif'))\n",
    "    \n",
    "output = gdal.Warp(xxx, xxx, format='GTiff',\n",
    "          srcSRS=src_info['coordinateSystem']['wkt'],   \n",
    "          dstSRS=dst_info['coordinateSystem']['wkt'],  \n",
    "          srcNodata=src_info['bands'][0]['noDataValue'],\n",
    "          dstNodata=src_info['bands'][0]['noDataValue'],   \n",
    "          width=dst_info['size'][0],\n",
    "          height=dst_info['size'][1],\n",
    "          outputBounds=(dst_info['cornerCoordinates']['lowerLeft'][0],              \n",
    "                        dst_info['cornerCoordinates']['lowerLeft'][1],\n",
    "                        dst_info['cornerCoordinates']['upperRight'][0],\n",
    "                        dst_info['cornerCoordinates']['upperRight'][1]),\n",
    "          outputBoundsSRS=dst_info['coordinateSystem']['wkt'],\n",
    "          resampleAlg='near')\n",
    "output = None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d4eb7fd-9aa3-4688-84ab-059c7c702a99",
   "metadata": {},
   "source": [
    "WPOP_T data represents a people count per pixel. When warping and adjusting the pixel size, this needs to be accounted for. You can quickly calculate a conversion factor based on the sizes of the arrays, as they both represent the exact same area:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa5e436e-1f07-44f5-b573-eb545a15994e",
   "metadata": {},
   "outputs": [],
   "source": [
    "WPOP_arr = rioxarray.open_rasterio(wpopfile)\n",
    "WPOP_pixcount = len(WPOP_arr['y']) * len(WPOP_arr['x'])\n",
    "print(WPOP_pixcount)\n",
    "\n",
    "drought_arr = rioxarray.open_rasterio('./results/drought_severity/drought_dry.tif')\n",
    "drought_pixcount = xxx\n",
    "print(drought_pixcount)\n",
    "\n",
    "conversion = drought_pixcount / WPOP_pixcount\n",
    "print(conversion)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a0a1d6f-eb9e-491a-914c-8e5c53b83d3e",
   "metadata": {},
   "source": [
    "Apply the conversion factor to the data and get the people count where the drought was extreme or severe.\n",
    "\n",
    "Make sure to first filter the WPOP data to get rid of the nodata value!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2482fab1-4f98-408a-8add-60f022e883ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "WPOP_tot = rioxarray.open_rasterio(output_file)\n",
    "nodata = -9999\n",
    "WPOP_tot = WPOP_tot.where(WPOP_tot != nodata)\n",
    "WPOP_tot = WPOP_tot / conversion\n",
    "\n",
    "WPOP_count_arr = np.where((drought_arr.values != 1) & (drought_arr.values != 2), 0, WPOP_tot.values)\n",
    "\n",
    "total_people = int(np.nansum(WPOP_tot))\n",
    "people_affected = int(np.nansum(WPOP_count_arr))\n",
    "print(f'In total {people_affected} people were directly affected by this drought, in a region occupied by a total of {total_people} inhabitants.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a616c99",
   "metadata": {},
   "source": [
    "QUESTION: report the total number of people directly affected by the drought."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.12 ('eoafrica')",
   "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"
  },
  "vscode": {
   "interpreter": {
    "hash": "243fe273f189138c19d68eb82425e8e07c0bc17275fde91d33ad51583df3419f"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}