2a_Zonal_stat.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div>\n",
    "<table style=\"width: 100%\">\n",
    "\t<tr>\n",
    "\t\t<td>\n",
    "\t\t<table style=\"width: 100%\">\n",
    "\t\t\t<tr>\n",
    "                <td ><center><font size=\"5\"><b>Module 49</b></font><center>\n",
    "                <center><font size=\"6\">Digital Innovations for Water Challenges</font><center></td>\n",
    "\t\t\t</tr>\n",
    "\t\t\t<tr>\n",
    "                <td><center><font size=\"14\">Notebook 2a</font><center></td>\n",
    "\t\t\t</tr>\n",
    "\t\t\t<tr>\n",
    "                <td><center><font size=\"6\"><b>Zonal Statistics</b></font><center></td>\n",
    "\t\t\t</tr>\n",
    "\t\t</table>\n",
    "\t\t</td>\n",
    "\t\t<td><center><img src='images\\ihe-delft-institute_unesco_fc-lr.jpg'></img></td>\n",
    "\t</tr>\n",
    "</table>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Table of contents\n",
    "1. [Learning objectives](#learningobs)\n",
    "2. [Introduction to Zonal Statistics](#introduction)\n",
    "3. [An example](#Example)\n",
    "\n",
    "# 1. Learning objectives<a name=\"learningobs\"></a>\n",
    "\n",
    "- Understand what Zonal Statistics is.\n",
    "- Compute Zonal Statistics from raster file using python. \n",
    "\n",
    "# 2. Introduction<a name=\"introduction\"></a>\n",
    "\n",
    "\n",
    "In water resources management we often work with spatial data (raster data). Many of the spatial data manipulation are easier with python than with GIS software. Think of repeatability and speed!<br>\n",
    "An example of a commonly used type of spatial data - the land cover map:\n",
    "<figure>\n",
    "<img src=\"./2a_Data/Upper Litani_LLC_2017_thumbnail.png\" style=\"width:75%\">\n",
    "</figure>\n",
    "\n",
    "\n",
    "#### Vector vs Raster\n",
    "<figure>\n",
    "<img src=\"./2a_Data/zstat2_thumbnail.png\" style=\"width:40%\">\n",
    "<figcaption align = \"center\">Vectors (first line) vs. rasters (second line) </figcaption>\n",
    "</figure>\n",
    "\n",
    "Raster data has its data represented in a grid matrix or cells organized into rows and columns. On the other hand, information in vector data utilizes vertices or sequential points. Polygons, lines, and points are used to represent data in vector formats.\n",
    "\n",
    "The type of data represented in raster data is continuous, while vector data is discrete. For instance, details concerning temperature, soil PH, flow, elevation, distance, and air pressure are examples of raster data. Examples of vector data include rivers, linear features, and administrative borders. Note that some raster data can be represented as vectors by simplifying it to isolines (commonly used for e.g. elevation or temperature).<br> The type of data to be represented determines the type of spatial data to be used. Some of spatial data formats are: Geotiff, Shapefile, Geojson, KML, Geopackage, Spatialite. You can read <a href=\"https://www.e-education.psu.edu/geog585/node/691\">here</a> and <a href=\"https://gisgeography.com/gis-formats/\">here</a> about spatial data formats.\n",
    "\n",
    "### Zonal statistics\n",
    "A zone is all the cells in a raster that have the same value, regardless of whether or not they are contiguous. Zones can also be in vector format, they will then be represented by pologons sharing one attribute.\n",
    "<figure>\n",
    "<img src=\"./2a_Data/zstat1_thumbnail.png\" style=\"width:40%\">\n",
    "<figcaption align = \"center\"> Four zones in this raster. </figcaption>\n",
    "</figure>\n",
    "\n",
    "We can calculate different statistics for each zone: Total, Mean, Median, Standard deviation, Min and max, etc.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 - An example <a name=\"Example\"></a>\n",
    "### 3.1 - Zonal stat from two raster files\n",
    "We will compute evapotranspiration per land cover class from two raster files. We need two raster files. One for the land cover class which specifies for each grid a value representing its land cover type. The other is evapotranspiration raster which specifies the monthly evapotranspiration value in mm for each grid.\n",
    "\n",
    "For this example, we use data for a small river basin called Upper Litani in Lebanon. The data are obtained from the  <a href=\"https://wapor.apps.fao.org/home/WAPOR_2/1\">WaPOR</a> database. In the figure below, the first layer represents the land cover class raster, the second layer resents the evapotranspiration raster and the bar chart represents the statistics we would like to calculate.\n",
    "\n",
    "<figure>\n",
    "<img src=\"./2a_Data/zonal_stat2_thumbnail.png\" style=\"width:30%\">\n",
    "<figcaption align = \"center\"> Calculate Zonal statistics. <a href=\"https://bikeshbade.com.np/tutorials/Detail/?title=Calculate%20Zonal%20Statistics%20and%20export%20as%20CSV&code=33\">Source</a></figcaption>\n",
    "\n",
    "    \n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read the raster files and the land cover class description\n",
    "The raster can be found in the data folder and the land cover class description are obtained from WaPOR database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import rioxarray as rio\n",
    "import os\n",
    "\n",
    "# path to the LCC data \n",
    "lcc_file = r'./2a_Data/lcc.tiff'\n",
    "# open lcc file \n",
    "LCC = rio.open_rasterio(lcc_file)\n",
    "LCC = LCC.squeeze(dim = 'band', drop = True)\n",
    "LCC = LCC.rio.write_crs(\"epsg:4326\", inplace=True)\n",
    "LCC = LCC.rename({'x': 'longitude','y': 'latitude'})\n",
    "\n",
    "# path to the evapotranspiration data \n",
    "et_file = r'./2a_Data/et_2017.tif'\n",
    "et = rio.open_rasterio(et_file)\n",
    "et = et.squeeze(dim = 'band', drop = True)\n",
    "et = et.rio.write_crs(\"epsg:4326\", inplace=True)\n",
    "et = et.rename({'x': 'longitude','y': 'latitude'})\n",
    "\n",
    "## land cover classes in WaPOR Landuse data\n",
    "lcc_dict = {\n",
    "                0: 'No data',\n",
    "                20:\t 'Shrubland',\n",
    "                30:\t 'Grassland', \n",
    "                41:\t 'Cropland, rainfed',\n",
    "                42:\t 'Cropland, irrigated', \n",
    "                43:\t 'Cropland, fallow', \n",
    "                50:\t 'Built-up', \n",
    "                60:\t 'Bare/sparse vegetation', \n",
    "                70:\t 'Permanent  snow/ice', \n",
    "                80:\t 'Water bodies', \n",
    "                81:\t 'Temporary water bodies', \n",
    "                90:\t 'Shrub, flooded', \n",
    "                111:  'Tree: closed, evergreen needle-leaved', \n",
    "                112: 'Tree: closed, evergreen broadleaved', \n",
    "                114: 'Tree: closed, deciduous broadleaved', \n",
    "                115: 'Tree: closed, mixed type', \n",
    "                116: 'Tree: closed, unknown type', \n",
    "                121: 'Tree: open, evergreen needle-leaved', \n",
    "                122: 'Tree: open, evergreen broadleaved', \n",
    "                123: 'Tree: open, deciduous needle-leaved', \n",
    "                124: 'Tree: open, deciduous broadleaved', \n",
    "                125: 'Tree: open, mixed type',\n",
    "                126: 'Tree: open, unknown type', \n",
    "                200: 'Sea wate'\n",
    "             }\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# All the land cover classes (lcc) may not be present in our raster file. so we need to filter the lcc which are in \n",
    "# our data and prepare a new lcc dict\n",
    "lcc_unq=np.unique(LCC.values[~np.isnan(LCC.values)])\n",
    "new_lcc_dict = {k: lcc_dict[k] for k in lcc_unq if k in lcc_dict}\n",
    "\n",
    "new_lcc_dict\n",
    "# change dict values to list\n",
    "lcc_values = list(new_lcc_dict.values())\n",
    "# change list ket to list\n",
    "lcc_keys = list(new_lcc_dict.keys())\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can loop through the lcc and calculate mean values of the evapotranspiration for each lcc\n",
    "et_lst=[] \n",
    "for code in new_lcc_dict.keys():\n",
    "    code=float(code)\n",
    "    \n",
    "    #mask the lcc type based on the code\n",
    "    mask = xr.where(LCC==code, 1,0)\n",
    "\n",
    "    #et \n",
    "    ds_lcc = et.where(mask, drop=True).mean()\n",
    "    code_et = ds_lcc.values.tolist()\n",
    "    et_lst.append(code_et)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Once we get the et per lcc, we can create a dataframe from the data and save it to a file and plot it for visulaization\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# create a dataframe from ET and LCC data\n",
    "lcc_keys = list(new_lcc_dict)\n",
    "data={'LU Code':lcc_keys,'Land Cover Class Description':lcc_values,\n",
    "        'Evapotranspiration (mm/year)':et_lst}\n",
    "df=pd.DataFrame(data)\n",
    "df = df.sort_values('Evapotranspiration (mm/year)',ascending=False)\n",
    "\n",
    "# write the dataframe to file\n",
    "df.to_csv(os.path.join(r'./2a_Data/ETperLCC_Kherlen.csv'),sep=',')\n",
    "\n",
    "# plot the data\n",
    "df.index = df['Land Cover Class Description']\n",
    "ax = df['Evapotranspiration (mm/year)'].plot.bar(rot=90, width=0.65, figsize=(10,5))\n",
    "ax.grid(visible=True, which='Major', color='0.4',linestyle='--', linewidth=0.25, zorder = 0) #Grid line\n",
    "                 \n",
    "ax.set_xlabel(\"Land Cover Class Description\")\n",
    "ax.set_ylabel('ET [mm/year]')\n",
    "plt.title('ET per land cover class for Upper Litani River Basin - 2009')\n",
    "plt.tight_layout()\n",
    "\n",
    "# save the plot to file \n",
    "plt.savefig(os.path.join(r'./2a_Data/ET_per_LCC_Kherlen_RB.png'), bbox_inches='tight', dpi=600)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compute P-ET per land cover class\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "\n",
    "**EXERCISE 2**:\n",
    "\n",
    "* Compute mean P-ET per land cover class for Upper Litani River Basin.\n",
    "    \n",
    "\n",
    "<details>\n",
    "  <summary>Hints</summary>\n",
    "\n",
    "* Follow the same step to read the precipitation data, compute precipitation (P) minus evapotranspiration (ET) and follow similar steps to compute mean P-ET.\n",
    "\n",
    "</details>\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}