lecture2_droughts.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "23047bb9-a4e6-45ce-b586-ecfd15327c37",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import ipywidgets as widgets\n",
    "import cartopy\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.io.img_tiles import OSM\n",
    "import matplotlib.pyplot as plt\n",
    "from utils import sampling\n",
    "from scipy.stats import theilslopes\n",
    "from ismn.interface import ISMN_Interface\n",
    "from pytesmo.df_metrics import ubrmsd\n",
    "from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n",
    "import utils\n",
    "from matplotlib.gridspec import GridSpec\n",
    "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
    "import datetime\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0701ab09-46d3-4a22-b64f-7463921b99c7",
   "metadata": {},
   "source": [
    "# Data sets for this lecture\n",
    "\n",
    "## Summary\n",
    "\n",
    "In this lecture we will compute the following soil moisture anomaly based drought indices using monthly soil moisture data from C3S SM (v202012, downloaded from DOI: 10.24381/cds.d7782f18, and v202212):\n",
    "\n",
    "* Soil Moisture anomalies\n",
    "* Anomaly Quantiles\n",
    "* Z-Scores\n",
    "\n",
    "We will also use pre-processed Soil Moisture Anomalies Stantardised Index (SMASI) data.\n",
    "\n",
    "In addition we will use ECMWF ERA-5 monthly precipitation data to compute:\n",
    "\n",
    "* Standardised Precipitation Index (SPI)\n",
    "\n",
    "Finally we will use soil moisture and Vegetation Optical Depth (VOD) time series to assess temporal lags between them."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a44d97a1-2917-44c2-9b4b-7e64623867e1",
   "metadata": {},
   "source": [
    "First we load our soil moisture stack. In this case use the same stack as the one from the previous lecture on satellite soil moisture retrieval. The data is produced as part of the Copernicus Climate Change Service (C3S), which combines retrievals from multiple passive sensors. The used subset contains monthly mean values for Europe from 1978 to September 2022.\n",
    "More information on this data can be found in the [C3S Data Store](https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-soil-moisture?tab=overview)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ffcfcfe1-15a9-4246-884c-301ec0c758a2",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
       " *\n",
       " */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
       "  --xr-background-color: var(--jp-layout-color0, white);\n",
       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
       "}\n",
       "\n",
       "html[theme=dark],\n",
       "body[data-theme=dark],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: rgba(255, 255, 255, 1);\n",
       "  --xr-font-color2: rgba(255, 255, 255, 0.54);\n",
       "  --xr-font-color3: rgba(255, 255, 255, 0.38);\n",
       "  --xr-border-color: #1F1F1F;\n",
       "  --xr-disabled-color: #515151;\n",
       "  --xr-background-color: #111111;\n",
       "  --xr-background-color-row-even: #111111;\n",
       "  --xr-background-color-row-odd: #313131;\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block !important;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-array-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: '►';\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: '▼';\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "  padding-bottom: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: '(';\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: ')';\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: ',';\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-index-preview {\n",
       "  grid-column: 2 / 5;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data,\n",
       ".xr-index-data {\n",
       "  display: none;\n",
       "  background-color: var(--xr-background-color) !important;\n",
       "  padding-bottom: 5px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data,\n",
       ".xr-index-data-in:checked ~ .xr-index-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-index-name div,\n",
       ".xr-index-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data,\n",
       ".xr-index-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2,\n",
       ".xr-no-icon {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt;\n",
       "Dimensions:  (lat: 160, lon: 240, time: 530)\n",
       "Coordinates:\n",
       "  * lat      (lat) float32 69.88 69.62 69.38 69.12 ... 30.88 30.62 30.38 30.12\n",
       "  * lon      (lon) float32 -9.875 -9.625 -9.375 -9.125 ... 49.38 49.62 49.88\n",
       "  * time     (time) datetime64[ns] 1978-11-01 1978-12-01 ... 2022-12-01\n",
       "Data variables:\n",
       "    sm       (time, lat, lon) float32 ...\n",
       "    nobs     (time, lat, lon) float32 ...\n",
       "    sm_anom  (time, lat, lon) float32 ...\n",
       "Attributes: (12/44)\n",
       "    title:                        C3S Surface Surface Soil Moisture merged PA...\n",
       "    institution:                  TU Wien (AUT); VanderSat B.V. (NL)\n",
       "    contact:                      cci_sm_contact@eodc.eu\n",
       "    source:                       LPRMv7/SMMR/Nimbus 7 L3 Surface Soil Moistu...\n",
       "    platform:                     Nimbus 7, DMSP, TRMM, AQUA, Coriolis, GCOM-...\n",
       "    processing_level:             Quality-controlled, super-collocated (L3S) ...\n",
       "    ...                           ...\n",
       "    id:                           C3S-SOILMOISTURE-L3S-SSMV-PASSIVE-MONTHLY-1...\n",
       "    history:                      2023-06-02T15:23:23.165973 mean calculated\n",
       "    date_created:                 2023-06-02T15:23:23Z\n",
       "    time_coverage_duration:       P1M\n",
       "    time_coverage_start:          1978-10-31T12:00:00Z\n",
       "    time_coverage_end:            1978-11-30T12:00:00Z</pre><div class='xr-wrap' style='display:none'><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-ea1c529c-27ce-44a8-8ce4-79f271687f78' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-ea1c529c-27ce-44a8-8ce4-79f271687f78' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span class='xr-has-index'>lat</span>: 160</li><li><span class='xr-has-index'>lon</span>: 240</li><li><span class='xr-has-index'>time</span>: 530</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-8a3034e1-6840-4aee-87fa-719d19b1955c' class='xr-section-summary-in' type='checkbox'  checked><label for='section-8a3034e1-6840-4aee-87fa-719d19b1955c' class='xr-section-summary' >Coordinates: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>lat</span></div><div class='xr-var-dims'>(lat)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>69.88 69.62 69.38 ... 30.38 30.12</div><input id='attrs-add22e25-9a6f-4557-b4c6-258a8dd0ec31' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-add22e25-9a6f-4557-b4c6-258a8dd0ec31' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-64dd26bd-a086-4968-99b5-9f1c7b8e2996' class='xr-var-data-in' type='checkbox'><label for='data-64dd26bd-a086-4968-99b5-9f1c7b8e2996' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>standard_name :</span></dt><dd>latitude</dd><dt><span>units :</span></dt><dd>degrees_north</dd><dt><span>valid_range :</span></dt><dd>[-90.  90.]</dd><dt><span>_CoordinateAxisType :</span></dt><dd>Lat</dd></dl></div><div class='xr-var-data'><pre>array([69.875, 69.625, 69.375, 69.125, 68.875, 68.625, 68.375, 68.125, 67.875,\n",
       "       67.625, 67.375, 67.125, 66.875, 66.625, 66.375, 66.125, 65.875, 65.625,\n",
       "       65.375, 65.125, 64.875, 64.625, 64.375, 64.125, 63.875, 63.625, 63.375,\n",
       "       63.125, 62.875, 62.625, 62.375, 62.125, 61.875, 61.625, 61.375, 61.125,\n",
       "       60.875, 60.625, 60.375, 60.125, 59.875, 59.625, 59.375, 59.125, 58.875,\n",
       "       58.625, 58.375, 58.125, 57.875, 57.625, 57.375, 57.125, 56.875, 56.625,\n",
       "       56.375, 56.125, 55.875, 55.625, 55.375, 55.125, 54.875, 54.625, 54.375,\n",
       "       54.125, 53.875, 53.625, 53.375, 53.125, 52.875, 52.625, 52.375, 52.125,\n",
       "       51.875, 51.625, 51.375, 51.125, 50.875, 50.625, 50.375, 50.125, 49.875,\n",
       "       49.625, 49.375, 49.125, 48.875, 48.625, 48.375, 48.125, 47.875, 47.625,\n",
       "       47.375, 47.125, 46.875, 46.625, 46.375, 46.125, 45.875, 45.625, 45.375,\n",
       "       45.125, 44.875, 44.625, 44.375, 44.125, 43.875, 43.625, 43.375, 43.125,\n",
       "       42.875, 42.625, 42.375, 42.125, 41.875, 41.625, 41.375, 41.125, 40.875,\n",
       "       40.625, 40.375, 40.125, 39.875, 39.625, 39.375, 39.125, 38.875, 38.625,\n",
       "       38.375, 38.125, 37.875, 37.625, 37.375, 37.125, 36.875, 36.625, 36.375,\n",
       "       36.125, 35.875, 35.625, 35.375, 35.125, 34.875, 34.625, 34.375, 34.125,\n",
       "       33.875, 33.625, 33.375, 33.125, 32.875, 32.625, 32.375, 32.125, 31.875,\n",
       "       31.625, 31.375, 31.125, 30.875, 30.625, 30.375, 30.125], dtype=float32)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>lon</span></div><div class='xr-var-dims'>(lon)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>-9.875 -9.625 ... 49.62 49.88</div><input id='attrs-74148a17-c0da-47f2-8139-34b016f62ff2' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-74148a17-c0da-47f2-8139-34b016f62ff2' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-52687827-5d2d-4561-bce0-d7a8086522a6' class='xr-var-data-in' type='checkbox'><label for='data-52687827-5d2d-4561-bce0-d7a8086522a6' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>standard_name :</span></dt><dd>longitude</dd><dt><span>units :</span></dt><dd>degrees_east</dd><dt><span>valid_range :</span></dt><dd>[-180.  180.]</dd><dt><span>_CoordinateAxisType :</span></dt><dd>Lon</dd></dl></div><div class='xr-var-data'><pre>array([-9.875, -9.625, -9.375, ..., 49.375, 49.625, 49.875], dtype=float32)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>time</span></div><div class='xr-var-dims'>(time)</div><div class='xr-var-dtype'>datetime64[ns]</div><div class='xr-var-preview xr-preview'>1978-11-01 ... 2022-12-01</div><input id='attrs-e0f74746-9ac7-4aa5-8201-7d069fa501cf' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-e0f74746-9ac7-4aa5-8201-7d069fa501cf' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-8d022c2a-4c05-4d3b-8042-8a9addb7e214' class='xr-var-data-in' type='checkbox'><label for='data-8d022c2a-4c05-4d3b-8042-8a9addb7e214' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>standard_name :</span></dt><dd>time</dd><dt><span>_CoordinateAxisType :</span></dt><dd>Time</dd></dl></div><div class='xr-var-data'><pre>array([&#x27;1978-11-01T00:00:00.000000000&#x27;, &#x27;1978-12-01T00:00:00.000000000&#x27;,\n",
       "       &#x27;1979-01-01T00:00:00.000000000&#x27;, ..., &#x27;2022-10-01T00:00:00.000000000&#x27;,\n",
       "       &#x27;2022-11-01T00:00:00.000000000&#x27;, &#x27;2022-12-01T00:00:00.000000000&#x27;],\n",
       "      dtype=&#x27;datetime64[ns]&#x27;)</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-71c7dbce-f409-4f5f-bafc-54a80ae5e3f4' class='xr-section-summary-in' type='checkbox'  checked><label for='section-71c7dbce-f409-4f5f-bafc-54a80ae5e3f4' class='xr-section-summary' >Data variables: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>sm</span></div><div class='xr-var-dims'>(time, lat, lon)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-89aac685-5358-4b52-9e54-cb4350574f07' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-89aac685-5358-4b52-9e54-cb4350574f07' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-c747df51-60d7-4594-85fa-7e83c27084ed' class='xr-var-data-in' type='checkbox'><label for='data-c747df51-60d7-4594-85fa-7e83c27084ed' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>dtype :</span></dt><dd>float32</dd><dt><span>units :</span></dt><dd>m3 m-3</dd><dt><span>valid_range :</span></dt><dd>[0. 1.]</dd><dt><span>long_name :</span></dt><dd>Volumetric Soil Moisture</dd><dt><span>ancillary_variables :</span></dt><dd>sm_uncertainty flag t0</dd><dt><span>_CoordinateAxes :</span></dt><dd>time lat lon</dd></dl></div><div class='xr-var-data'><pre>[20352000 values with dtype=float32]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>nobs</span></div><div class='xr-var-dims'>(time, lat, lon)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-dbf9e38e-b311-40bd-98d0-6bcdfdfa4141' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-dbf9e38e-b311-40bd-98d0-6bcdfdfa4141' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-84a21c34-e72e-4996-90b7-d2693f126d0d' class='xr-var-data-in' type='checkbox'><label for='data-84a21c34-e72e-4996-90b7-d2693f126d0d' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>Number of valid observation</dd><dt><span>CoordinateAxes :</span></dt><dd>time lat lon</dd></dl></div><div class='xr-var-data'><pre>[20352000 values with dtype=float32]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>sm_anom</span></div><div class='xr-var-dims'>(time, lat, lon)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-296d8bd5-e830-4ca3-9554-ecc9a4c46278' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-296d8bd5-e830-4ca3-9554-ecc9a4c46278' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-7ccde2e8-d7d6-433d-914a-e7d2eaba5830' class='xr-var-data-in' type='checkbox'><label for='data-7ccde2e8-d7d6-433d-914a-e7d2eaba5830' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>[20352000 values with dtype=float32]</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-d8112079-bc39-4f01-a31b-f84913a92a8d' class='xr-section-summary-in' type='checkbox'  ><label for='section-d8112079-bc39-4f01-a31b-f84913a92a8d' class='xr-section-summary' >Indexes: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-index-name'><div>lat</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-b11e1cee-82aa-4c2f-bbfb-8a0ea87b04ca' class='xr-index-data-in' type='checkbox'/><label for='index-b11e1cee-82aa-4c2f-bbfb-8a0ea87b04ca' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(Float64Index([69.875, 69.625, 69.375, 69.125, 68.875, 68.625, 68.375, 68.125,\n",
       "              67.875, 67.625,\n",
       "              ...\n",
       "              32.375, 32.125, 31.875, 31.625, 31.375, 31.125, 30.875, 30.625,\n",
       "              30.375, 30.125],\n",
       "             dtype=&#x27;float64&#x27;, name=&#x27;lat&#x27;, length=160))</pre></div></li><li class='xr-var-item'><div class='xr-index-name'><div>lon</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-9bb1a9c7-ace4-4cdb-bbc5-7c81b4706979' class='xr-index-data-in' type='checkbox'/><label for='index-9bb1a9c7-ace4-4cdb-bbc5-7c81b4706979' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(Float64Index([-9.875, -9.625, -9.375, -9.125, -8.875, -8.625, -8.375, -8.125,\n",
       "              -7.875, -7.625,\n",
       "              ...\n",
       "              47.625, 47.875, 48.125, 48.375, 48.625, 48.875, 49.125, 49.375,\n",
       "              49.625, 49.875],\n",
       "             dtype=&#x27;float64&#x27;, name=&#x27;lon&#x27;, length=240))</pre></div></li><li class='xr-var-item'><div class='xr-index-name'><div>time</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-c42d591d-b2e8-4f48-947a-03caf9aeb72c' class='xr-index-data-in' type='checkbox'/><label for='index-c42d591d-b2e8-4f48-947a-03caf9aeb72c' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(DatetimeIndex([&#x27;1978-11-01&#x27;, &#x27;1978-12-01&#x27;, &#x27;1979-01-01&#x27;, &#x27;1979-02-01&#x27;,\n",
       "               &#x27;1979-03-01&#x27;, &#x27;1979-04-01&#x27;, &#x27;1979-05-01&#x27;, &#x27;1979-06-01&#x27;,\n",
       "               &#x27;1979-07-01&#x27;, &#x27;1979-08-01&#x27;,\n",
       "               ...\n",
       "               &#x27;2022-03-01&#x27;, &#x27;2022-04-01&#x27;, &#x27;2022-05-01&#x27;, &#x27;2022-06-01&#x27;,\n",
       "               &#x27;2022-07-01&#x27;, &#x27;2022-08-01&#x27;, &#x27;2022-09-01&#x27;, &#x27;2022-10-01&#x27;,\n",
       "               &#x27;2022-11-01&#x27;, &#x27;2022-12-01&#x27;],\n",
       "              dtype=&#x27;datetime64[ns]&#x27;, name=&#x27;time&#x27;, length=530, freq=None))</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-c53177ef-c6af-44e8-abe1-df4525e4c1e6' class='xr-section-summary-in' type='checkbox'  ><label for='section-c53177ef-c6af-44e8-abe1-df4525e4c1e6' class='xr-section-summary' >Attributes: <span>(44)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'><dt><span>title :</span></dt><dd>C3S Surface Surface Soil Moisture merged PASSIVE Product</dd><dt><span>institution :</span></dt><dd>TU Wien (AUT); VanderSat B.V. (NL)</dd><dt><span>contact :</span></dt><dd>cci_sm_contact@eodc.eu</dd><dt><span>source :</span></dt><dd>LPRMv7/SMMR/Nimbus 7 L3 Surface Soil Moisture, Ancillary Params, and quality flags; LPRMv7/SSMI/F08, F11, F13 DMSP L3 Surface Soil Moisture, Ancillary Params, and quality flags; LPRMv7/TMI/TRMM L2 Surface Soil Moisture, Ancillary Params, and QC; LPRMv7/AMSR-E/Aqua L2B Surface Soil Moisture, Ancillary Params, and QC; LPRMv7/WINDSAT/CORIOLIS L2 Surface Soil Moisture, Ancillary Params, and QC; LPRMv7/AMSR2/GCOM-W1 L3 Surface Soil Moisture, Ancillary Params; LPRMv7/SMOS/MIRAS L3 Surface Soil Moisture, CATDS Level 3 Brightness Temperatures (L3TB) version 300 RE03 &amp; RE04; LPRMv7/SMAP_radiometer/SMAP L2 Surface Soil Moisture, Ancillary Params, and QC; LPRMv7/GMI/GPM L3 Surface Soil Moisture, Ancillary Params, and quality flags; LPRMv7/VIRR/FengYun-3B L3 Surface Soil Moisture, Ancillary Params, and quality flags; LPRMv7/VIRR/FengYun-3C L3 Surface Soil Moisture, Ancillary Params, and quality flags; LPRMv7/VIRR/FengYun-3D L3 Surface Soil Moisture, Ancillary Params, and quality flags;</dd><dt><span>platform :</span></dt><dd>Nimbus 7, DMSP, TRMM, AQUA, Coriolis, GCOM-W1, MIRAS, SMAP, GPM, FengYun-3B, FengYun-3C, FengYun-3D</dd><dt><span>processing_level :</span></dt><dd>Quality-controlled, super-collocated (L3S) Satellite Soil Moisture (SM) data from multiple sensors</dd><dt><span>sensor :</span></dt><dd>SMMR, SSM/I, TMI, AMSR-E, WindSat, AMSR2, SMOS, SMAP_radiometer, GMI, VIRR-3B, VIRR-3C, VIRR-3D</dd><dt><span>references :</span></dt><dd>https://doi.org/10.24381/cds.d7782f18; Dorigo, W.A., Wagner, W., Albergel, C., Albrecht, F.,  Balsamo, G., Brocca, L., Chung, D., Ertl, M., Forkel, M., Gruber, A., Haas, E., Hamer, D. P. Hirschi, M., Ikonen, J., De Jeu, R. Kidd, R. Lahoz, W., Liu, Y.Y., Miralles, D., Lecomte, P. (2017) ESA CCI Soil Moisture for improved Earth system understanding: State-of-the art and future directions. In Remote Sensing of Environment, 2017, ISSN 0034-4257, https://doi.org/10.1016/j.rse.2017.07.001; Gruber, A., Scanlon, T., van der Schalie, R., Wagner, W., Dorigo, W. (2019) Evolution of the ESA CCI Soil Moisture Climate Data Records and their underlying merging methodology. Earth System Science Data 11, 717-739, https://doi.org/10.5194/essd-11-717-2019; Gruber, A., Dorigo, W. A., Crow, W., Wagner W. (2017). Triple Collocation-Based Merging of Satellite Soil Moisture Retrievals. IEEE Transactions on Geoscience and Remote Sensing. PP. 1-13. https://doi.org/10.1109/TGRS.2017.2734070</dd><dt><span>product_version :</span></dt><dd>v202212.0.0</dd><dt><span>tracking_id :</span></dt><dd>11dfa6eb-c9c7-4163-9597-09930590dd05</dd><dt><span>Conventions :</span></dt><dd>CF-1.9</dd><dt><span>standard_name_vocabulary :</span></dt><dd>CF Standard Name Table v77</dd><dt><span>summary :</span></dt><dd>This dataset was produced with funding of the Copernicus Climate Change Service implemented by ECMWF through C3S 312a Lot 7 &amp; C3S2 312a Lot 4 Soil Moisture</dd><dt><span>keywords :</span></dt><dd>soil moisture, soil water content, microwave remote sensing, active, passive, combined, climate data record</dd><dt><span>naming_authority :</span></dt><dd>TU Wien, EODC</dd><dt><span>keywords_vocabulary :</span></dt><dd>NASA Global Change Master Directory (GCMD) Science Keywords</dd><dt><span>cdm_data_type :</span></dt><dd>Grid</dd><dt><span>comment :</span></dt><dd>This dataset was produced with funding of the Copernicus Climate Change Service implemented by ECMWF through C3S 312a Lot 7 &amp; C3S2 312a Lot 4 Soil Moisture</dd><dt><span>creator_name :</span></dt><dd>Earth Observation Data Centre (EODC)</dd><dt><span>creator_url :</span></dt><dd>https://eodc.eu/</dd><dt><span>creator_email :</span></dt><dd>C3S_SM_Science@eodc.eu</dd><dt><span>project :</span></dt><dd>C3S2 312a Lot 4 Soil Moisture</dd><dt><span>license :</span></dt><dd>Copernicus Data License</dd><dt><span>time_coverage_start_product :</span></dt><dd>19781101T000000Z</dd><dt><span>time_coverage_end_product :</span></dt><dd>20221231T235959Z</dd><dt><span>key_variables :</span></dt><dd>sm</dd><dt><span>time_coverage_resolution :</span></dt><dd>P1D</dd><dt><span>geospatial_lat_min :</span></dt><dd>-90.0</dd><dt><span>geospatial_lat_max :</span></dt><dd>90.0</dd><dt><span>geospatial_lon_min :</span></dt><dd>-180.0</dd><dt><span>geospatial_lon_max :</span></dt><dd>180.0</dd><dt><span>geospatial_vertical_min :</span></dt><dd>0.0</dd><dt><span>geospatial_vertical_max :</span></dt><dd>0.0</dd><dt><span>geospatial_lat_units :</span></dt><dd>degrees_north</dd><dt><span>geospatial_lon_units :</span></dt><dd>degrees_east</dd><dt><span>geospatial_lat_resolution :</span></dt><dd>0.25 degree</dd><dt><span>geospatial_lon_resolution :</span></dt><dd>0.25 degree</dd><dt><span>spatial_resolution :</span></dt><dd>25km</dd><dt><span>id :</span></dt><dd>C3S-SOILMOISTURE-L3S-SSMV-PASSIVE-MONTHLY-19781101000000-TCDR-v202212.0.0.nc</dd><dt><span>history :</span></dt><dd>2023-06-02T15:23:23.165973 mean calculated</dd><dt><span>date_created :</span></dt><dd>2023-06-02T15:23:23Z</dd><dt><span>time_coverage_duration :</span></dt><dd>P1M</dd><dt><span>time_coverage_start :</span></dt><dd>1978-10-31T12:00:00Z</dd><dt><span>time_coverage_end :</span></dt><dd>1978-11-30T12:00:00Z</dd></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (lat: 160, lon: 240, time: 530)\n",
       "Coordinates:\n",
       "  * lat      (lat) float32 69.88 69.62 69.38 69.12 ... 30.88 30.62 30.38 30.12\n",
       "  * lon      (lon) float32 -9.875 -9.625 -9.375 -9.125 ... 49.38 49.62 49.88\n",
       "  * time     (time) datetime64[ns] 1978-11-01 1978-12-01 ... 2022-12-01\n",
       "Data variables:\n",
       "    sm       (time, lat, lon) float32 ...\n",
       "    nobs     (time, lat, lon) float32 ...\n",
       "    sm_anom  (time, lat, lon) float32 ...\n",
       "Attributes: (12/44)\n",
       "    title:                        C3S Surface Surface Soil Moisture merged PA...\n",
       "    institution:                  TU Wien (AUT); VanderSat B.V. (NL)\n",
       "    contact:                      cci_sm_contact@eodc.eu\n",
       "    source:                       LPRMv7/SMMR/Nimbus 7 L3 Surface Soil Moistu...\n",
       "    platform:                     Nimbus 7, DMSP, TRMM, AQUA, Coriolis, GCOM-...\n",
       "    processing_level:             Quality-controlled, super-collocated (L3S) ...\n",
       "    ...                           ...\n",
       "    id:                           C3S-SOILMOISTURE-L3S-SSMV-PASSIVE-MONTHLY-1...\n",
       "    history:                      2023-06-02T15:23:23.165973 mean calculated\n",
       "    date_created:                 2023-06-02T15:23:23Z\n",
       "    time_coverage_duration:       P1M\n",
       "    time_coverage_start:          1978-10-31T12:00:00Z\n",
       "    time_coverage_end:            1978-11-30T12:00:00Z"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "c3s_europe = xr.open_dataset('./LTC_DATA/STACK_C3S-SOILMOISTURE_EUROPE_v202212_PASSIVE_MONTHLY.nc')\n",
    "display(c3s_europe)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c0f856e-0b54-4d65-98a6-decb4c60489c",
   "metadata": {},
   "source": [
    "In addition we load precipitation data from ECMWF's ERA5 reanalysis (called `total precipitation` or `tp`). We have downloaded monthly averages for Europe in a similar time range as the soil moisture data. We sum up the (average) rainfall amount for each month. More background information on how to correctly sum-up precipitation data is available in the [support portal](https://confluence.ecmwf.int/pages/viewpage.action?pageId=197702790). \n",
    "ERA5 monthly means are available through the [Copernicus Climate Data Store](https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels-monthly-means?tab=overview)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0da2468e-701d-4bb9-be6b-160c0a114523",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
       " *\n",
       " */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
       "  --xr-background-color: var(--jp-layout-color0, white);\n",
       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
       "}\n",
       "\n",
       "html[theme=dark],\n",
       "body[data-theme=dark],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: rgba(255, 255, 255, 1);\n",
       "  --xr-font-color2: rgba(255, 255, 255, 0.54);\n",
       "  --xr-font-color3: rgba(255, 255, 255, 0.38);\n",
       "  --xr-border-color: #1F1F1F;\n",
       "  --xr-disabled-color: #515151;\n",
       "  --xr-background-color: #111111;\n",
       "  --xr-background-color-row-even: #111111;\n",
       "  --xr-background-color-row-odd: #313131;\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block !important;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-array-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: '►';\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: '▼';\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "  padding-bottom: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: '(';\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: ')';\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: ',';\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-index-preview {\n",
       "  grid-column: 2 / 5;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data,\n",
       ".xr-index-data {\n",
       "  display: none;\n",
       "  background-color: var(--xr-background-color) !important;\n",
       "  padding-bottom: 5px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data,\n",
       ".xr-index-data-in:checked ~ .xr-index-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-index-name div,\n",
       ".xr-index-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data,\n",
       ".xr-index-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2,\n",
       ".xr-no-icon {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt;\n",
       "Dimensions:  (lon: 241, lat: 161, time: 548)\n",
       "Coordinates:\n",
       "  * lon      (lon) float32 -10.0 -9.75 -9.5 -9.25 -9.0 ... 49.25 49.5 49.75 50.0\n",
       "  * lat      (lat) float32 70.0 69.75 69.5 69.25 69.0 ... 30.75 30.5 30.25 30.0\n",
       "  * time     (time) datetime64[ns] 1978-01-01 1978-02-01 ... 2023-08-01\n",
       "Data variables:\n",
       "    tp       (time, lat, lon) float64 64.75 65.48 65.9 67.05 ... 0.0 0.0 0.0 0.0</pre><div class='xr-wrap' style='display:none'><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-0bcf3d97-b645-4cf7-a247-fba6adc12d30' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-0bcf3d97-b645-4cf7-a247-fba6adc12d30' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span class='xr-has-index'>lon</span>: 241</li><li><span class='xr-has-index'>lat</span>: 161</li><li><span class='xr-has-index'>time</span>: 548</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-28d396cd-8f8b-412b-b1c9-765937554818' class='xr-section-summary-in' type='checkbox'  checked><label for='section-28d396cd-8f8b-412b-b1c9-765937554818' class='xr-section-summary' >Coordinates: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>lon</span></div><div class='xr-var-dims'>(lon)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>-10.0 -9.75 -9.5 ... 49.75 50.0</div><input id='attrs-e50895f0-0e42-43f9-8b61-e36899ca49c7' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-e50895f0-0e42-43f9-8b61-e36899ca49c7' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-6ddc43b9-bab7-4431-9eb0-894470c48f77' class='xr-var-data-in' type='checkbox'><label for='data-6ddc43b9-bab7-4431-9eb0-894470c48f77' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>degrees_east</dd><dt><span>long_name :</span></dt><dd>longitude</dd></dl></div><div class='xr-var-data'><pre>array([-10.  ,  -9.75,  -9.5 , ...,  49.5 ,  49.75,  50.  ], dtype=float32)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>lat</span></div><div class='xr-var-dims'>(lat)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>70.0 69.75 69.5 ... 30.5 30.25 30.0</div><input id='attrs-328e5031-0049-4fae-9cf5-09f1ffed7940' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-328e5031-0049-4fae-9cf5-09f1ffed7940' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-f826cbac-9ca7-4dba-affb-e5d6829c9591' class='xr-var-data-in' type='checkbox'><label for='data-f826cbac-9ca7-4dba-affb-e5d6829c9591' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>degrees_north</dd><dt><span>long_name :</span></dt><dd>latitude</dd></dl></div><div class='xr-var-data'><pre>array([70.  , 69.75, 69.5 , 69.25, 69.  , 68.75, 68.5 , 68.25, 68.  , 67.75,\n",
       "       67.5 , 67.25, 67.  , 66.75, 66.5 , 66.25, 66.  , 65.75, 65.5 , 65.25,\n",
       "       65.  , 64.75, 64.5 , 64.25, 64.  , 63.75, 63.5 , 63.25, 63.  , 62.75,\n",
       "       62.5 , 62.25, 62.  , 61.75, 61.5 , 61.25, 61.  , 60.75, 60.5 , 60.25,\n",
       "       60.  , 59.75, 59.5 , 59.25, 59.  , 58.75, 58.5 , 58.25, 58.  , 57.75,\n",
       "       57.5 , 57.25, 57.  , 56.75, 56.5 , 56.25, 56.  , 55.75, 55.5 , 55.25,\n",
       "       55.  , 54.75, 54.5 , 54.25, 54.  , 53.75, 53.5 , 53.25, 53.  , 52.75,\n",
       "       52.5 , 52.25, 52.  , 51.75, 51.5 , 51.25, 51.  , 50.75, 50.5 , 50.25,\n",
       "       50.  , 49.75, 49.5 , 49.25, 49.  , 48.75, 48.5 , 48.25, 48.  , 47.75,\n",
       "       47.5 , 47.25, 47.  , 46.75, 46.5 , 46.25, 46.  , 45.75, 45.5 , 45.25,\n",
       "       45.  , 44.75, 44.5 , 44.25, 44.  , 43.75, 43.5 , 43.25, 43.  , 42.75,\n",
       "       42.5 , 42.25, 42.  , 41.75, 41.5 , 41.25, 41.  , 40.75, 40.5 , 40.25,\n",
       "       40.  , 39.75, 39.5 , 39.25, 39.  , 38.75, 38.5 , 38.25, 38.  , 37.75,\n",
       "       37.5 , 37.25, 37.  , 36.75, 36.5 , 36.25, 36.  , 35.75, 35.5 , 35.25,\n",
       "       35.  , 34.75, 34.5 , 34.25, 34.  , 33.75, 33.5 , 33.25, 33.  , 32.75,\n",
       "       32.5 , 32.25, 32.  , 31.75, 31.5 , 31.25, 31.  , 30.75, 30.5 , 30.25,\n",
       "       30.  ], dtype=float32)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>time</span></div><div class='xr-var-dims'>(time)</div><div class='xr-var-dtype'>datetime64[ns]</div><div class='xr-var-preview xr-preview'>1978-01-01 ... 2023-08-01</div><input id='attrs-6f3f4318-66e0-4e2a-89ad-bee09cd564ed' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-6f3f4318-66e0-4e2a-89ad-bee09cd564ed' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-88060384-c4fe-43ca-b767-555dd66b071a' class='xr-var-data-in' type='checkbox'><label for='data-88060384-c4fe-43ca-b767-555dd66b071a' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>time</dd></dl></div><div class='xr-var-data'><pre>array([&#x27;1978-01-01T00:00:00.000000000&#x27;, &#x27;1978-02-01T00:00:00.000000000&#x27;,\n",
       "       &#x27;1978-03-01T00:00:00.000000000&#x27;, ..., &#x27;2023-06-01T00:00:00.000000000&#x27;,\n",
       "       &#x27;2023-07-01T00:00:00.000000000&#x27;, &#x27;2023-08-01T00:00:00.000000000&#x27;],\n",
       "      dtype=&#x27;datetime64[ns]&#x27;)</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-ec4b7280-fa20-4b84-8c24-a9142af8659b' class='xr-section-summary-in' type='checkbox'  checked><label for='section-ec4b7280-fa20-4b84-8c24-a9142af8659b' class='xr-section-summary' >Data variables: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>tp</span></div><div class='xr-var-dims'>(time, lat, lon)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>64.75 65.48 65.9 ... 0.0 0.0 0.0</div><input id='attrs-862a0cc5-05e2-4bf9-b280-1738b462de27' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-862a0cc5-05e2-4bf9-b280-1738b462de27' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-d516e668-9ca8-494c-a84b-a903904820e2' class='xr-var-data-in' type='checkbox'><label for='data-d516e668-9ca8-494c-a84b-a903904820e2' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([[[ 64.74664927,  65.48381782,  65.89845896, ...,  70.72067165,\n",
       "          68.708848  ,  66.72774863],\n",
       "        [ 61.96698081,  62.70414567,  63.11878681, ...,  58.18906868,\n",
       "          55.37867594,  53.36685967],\n",
       "        [ 58.12763119,  57.43654537,  59.24872029, ...,  53.75078762,\n",
       "          50.78680301,  47.89961433],\n",
       "        ...,\n",
       "        [ 82.00833273,  68.5245688 ,  60.60017133, ...,  52.82933986,\n",
       "          57.55942035,  69.86065769],\n",
       "        [ 86.47733808,  82.42298126,  59.54051971, ...,  44.32137311,\n",
       "          48.3910588 ,  42.95456362],\n",
       "        [ 80.05795765,  61.87482643,  74.52930069, ...,  48.83642054,\n",
       "          46.76319635,  35.69053137]],\n",
       "\n",
       "       [[ 79.67599583,  80.59151077,  80.21698952, ...,  72.62943935,\n",
       "          68.52358437,  65.4164505 ],\n",
       "        [ 84.00381851,  84.91930676,  84.60028076, ...,  51.32332754,\n",
       "          51.43431139,  51.15687847],\n",
       "        [ 87.31902122,  87.95709991,  87.52709007, ...,  43.58322048,\n",
       "          43.68030548,  43.79128933],\n",
       "...\n",
       "        [  0.35320781,   0.29177032,   0.23033284, ...,   0.        ,\n",
       "           0.        ,   0.        ],\n",
       "        [  0.53752027,   0.8292906 ,   0.35320781, ...,   0.        ,\n",
       "           0.        ,   0.        ],\n",
       "        [  0.53752027,   2.54930928,   1.05968118, ...,   0.        ,\n",
       "           0.        ,   0.        ]],\n",
       "\n",
       "       [[ 55.82404137,  57.46726227,  56.88369298, ...,  47.42356157,\n",
       "          47.42356157,  47.65392351],\n",
       "        [ 55.99296582,  57.65154517,  57.05261743, ...,  53.21326411,\n",
       "          51.90789211,  52.92149425],\n",
       "        [ 54.16543245,  55.16367614,  55.16367614, ...,  58.12763119,\n",
       "          56.59189355,  55.04080117],\n",
       "        ...,\n",
       "        [  2.65682489,   2.24217959,   1.82753429, ...,   1.72001868,\n",
       "           1.95040926,   2.0732265 ],\n",
       "        [  2.0732265 ,   3.07147019,   3.96219827, ...,   0.29177032,\n",
       "           0.23033284,   0.29177032],\n",
       "        [  1.47432648,   3.48611549,   2.96395458, ...,   0.        ,\n",
       "           0.        ,   0.        ]]])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-5412e6cd-469e-4383-8245-899dea78bddd' class='xr-section-summary-in' type='checkbox'  ><label for='section-5412e6cd-469e-4383-8245-899dea78bddd' class='xr-section-summary' >Indexes: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-index-name'><div>lon</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-9e89f245-28e5-4b92-b9d2-16bf122a4d4c' class='xr-index-data-in' type='checkbox'/><label for='index-9e89f245-28e5-4b92-b9d2-16bf122a4d4c' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(Float64Index([-10.0, -9.75,  -9.5, -9.25,  -9.0, -8.75,  -8.5, -8.25,  -8.0,\n",
       "              -7.75,\n",
       "              ...\n",
       "              47.75,  48.0, 48.25,  48.5, 48.75,  49.0, 49.25,  49.5, 49.75,\n",
       "               50.0],\n",
       "             dtype=&#x27;float64&#x27;, name=&#x27;lon&#x27;, length=241))</pre></div></li><li class='xr-var-item'><div class='xr-index-name'><div>lat</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-a7e6e19e-9a58-43ea-877c-2478b9f4b6f0' class='xr-index-data-in' type='checkbox'/><label for='index-a7e6e19e-9a58-43ea-877c-2478b9f4b6f0' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(Float64Index([ 70.0, 69.75,  69.5, 69.25,  69.0, 68.75,  68.5, 68.25,  68.0,\n",
       "              67.75,\n",
       "              ...\n",
       "              32.25,  32.0, 31.75,  31.5, 31.25,  31.0, 30.75,  30.5, 30.25,\n",
       "               30.0],\n",
       "             dtype=&#x27;float64&#x27;, name=&#x27;lat&#x27;, length=161))</pre></div></li><li class='xr-var-item'><div class='xr-index-name'><div>time</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-fa3b2b5b-256f-4321-85e8-4d461977b0cc' class='xr-index-data-in' type='checkbox'/><label for='index-fa3b2b5b-256f-4321-85e8-4d461977b0cc' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(DatetimeIndex([&#x27;1978-01-01&#x27;, &#x27;1978-02-01&#x27;, &#x27;1978-03-01&#x27;, &#x27;1978-04-01&#x27;,\n",
       "               &#x27;1978-05-01&#x27;, &#x27;1978-06-01&#x27;, &#x27;1978-07-01&#x27;, &#x27;1978-08-01&#x27;,\n",
       "               &#x27;1978-09-01&#x27;, &#x27;1978-10-01&#x27;,\n",
       "               ...\n",
       "               &#x27;2022-11-01&#x27;, &#x27;2022-12-01&#x27;, &#x27;2023-01-01&#x27;, &#x27;2023-02-01&#x27;,\n",
       "               &#x27;2023-03-01&#x27;, &#x27;2023-04-01&#x27;, &#x27;2023-05-01&#x27;, &#x27;2023-06-01&#x27;,\n",
       "               &#x27;2023-07-01&#x27;, &#x27;2023-08-01&#x27;],\n",
       "              dtype=&#x27;datetime64[ns]&#x27;, name=&#x27;time&#x27;, length=548, freq=None))</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-6b03a5f3-fedc-49cb-a145-ba185b924a24' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-6b03a5f3-fedc-49cb-a145-ba185b924a24' class='xr-section-summary'  title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (lon: 241, lat: 161, time: 548)\n",
       "Coordinates:\n",
       "  * lon      (lon) float32 -10.0 -9.75 -9.5 -9.25 -9.0 ... 49.25 49.5 49.75 50.0\n",
       "  * lat      (lat) float32 70.0 69.75 69.5 69.25 69.0 ... 30.75 30.5 30.25 30.0\n",
       "  * time     (time) datetime64[ns] 1978-01-01 1978-02-01 ... 2023-08-01\n",
       "Data variables:\n",
       "    tp       (time, lat, lon) float64 64.75 65.48 65.9 67.05 ... 0.0 0.0 0.0 0.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "era5_europe = xr.open_dataset('./LTC_DATA/era5_tp.nc').sel(time=slice('1978-01-01',None))\n",
    "\n",
    "days_in_month = pd.DatetimeIndex(era5_europe.time).days_in_month.values\n",
    "era5_europe['tp'].values = era5_europe['tp'].values * 1000 * days_in_month.reshape(len(days_in_month), 1, 1)\n",
    "era5_europe = era5_europe.rename(dict(longitude=\"lon\", latitude=\"lat\"))\n",
    "display(era5_europe)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1d2fa28-e3c6-4435-b676-7fee22fa4b90",
   "metadata": {},
   "source": [
    "# Loading data for the following examples\n",
    "Droughts are extreme events in time. We therefore use time series as inputs to compute the indices. The following example allows you to extract the soil moisture and precipitation data for a chosen point in Europe. The left plot shows the location of the extracted time series.\n",
    "We extract a time series for both data sets from our stacks and combine them in a pandas DataFrame for further processing. The figure on the right shows the extracted values for the chosen location and time period.\n",
    "\n",
    "**Try:**\n",
    "\n",
    "* **Changing the coordinates and see what a time series for a different location than the default one looks like. Find a suitable time series for the rest of the lecture or keep the default one in the area West of Wroclaw with a mixture of cultivated land and forest patches over a gentle topography.**\n",
    "\n",
    "You can always reset the interactive controls by re-running the cell.\n",
    "\n",
    "<font color='red'>**Note: The data you finally choose here will be used in the rest of the notebook**</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "03b15f81-2126-4cfa-affe-455c9e3ccd16",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "616a897d5e1346209a4c487a2aa0fd09",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntRangeSlider(value=(1991, 2022), description='Time series range (year from, year to)',…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "extracted_ts = dict()\n",
    "\n",
    "slider=widgets.IntRangeSlider(min=1979, max=2022, value=[1991, 2022], step=1, description='Time series range (year from, year to)', \n",
    "                              continuous_update=True, style={'description_width': 'initial'}, layout=widgets.Layout(width='30%'))\n",
    "@widgets.interact(period=slider, Longitude=\"16.5\", Latitude=\"51.0\")\n",
    "def extract_ts(period, Longitude, Latitude):\n",
    "    lon, lat = float(Longitude), float(Latitude)\n",
    "    \n",
    "    # Extract a pandas time series at location\n",
    "    ts_sm = c3s_europe.sel(lon=lon, lat=lat, method='nearest') \\\n",
    "                   .to_pandas().rename(columns={'sm': 'soil_moisture'}) \\\n",
    "                   .loc[f\"{period[0]}-01-01\":f\"{period[1]}-12-31\", ['soil_moisture']]\n",
    "    \n",
    "    ts_tp = era5_europe.sel(lon=lon, lat=lat, method='nearest') \\\n",
    "               .to_pandas().rename(columns={'tp': 'precipitation'}) \\\n",
    "               .loc[f\"{period[0]}-01-01\":f\"{period[1]}-12-31\", ['precipitation']]\n",
    "    \n",
    "    ts = pd.concat([ts_sm, ts_tp], axis=1).resample('MS').mean()\n",
    "\n",
    "    # Set up figures\n",
    "    fig = plt.figure(figsize=(15,4), constrained_layout=True)\n",
    "    gs = fig.add_gridspec(1, 3)\n",
    "    map_ax = fig.add_subplot(gs[0, 0], projection=ccrs.PlateCarree())\n",
    "    ts_ax = fig.add_subplot(gs[0, 1:])\n",
    "    bar_ax = ts_ax.twinx()\n",
    "\n",
    "    # Plot overview map\n",
    "    map_ax.set_extent([-14, 46, 29, 71])\n",
    "    map_ax.add_feature(cartopy.feature.LAND, zorder=0)\n",
    "    map_ax.add_feature(cartopy.feature.BORDERS)\n",
    "    map_ax.coastlines()\n",
    "    map_ax.plot([lon], [lat], 'ro', transform=ccrs.PlateCarree())\n",
    "    map_ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)\n",
    "    \n",
    "    # Plot time series\n",
    "    ts['soil_moisture'].plot(style='-o', ax=ts_ax, label='soil moisture')\n",
    "    ts['precipitation'].plot.area(ax=bar_ax, alpha=0.25, label='precipitation', linewidth=0)\n",
    "    \n",
    "    ts_ax.set_xlabel('Time [Year]')\n",
    "    ts_ax.set_ylabel(f'Soil moisture $[m³/m³]$')\n",
    "    bar_ax.set_ylabel(f'Precipitation $[mm]$')\n",
    "    ts_ax.set_title('Satellite Soil Moisture and ERA5 Precipitation')\n",
    "    \n",
    "    ts_ax.legend(loc='upper left')\n",
    "    bar_ax.legend(loc='upper right')\n",
    "    \n",
    "    # Make data available for the next example\n",
    "    global extracted_ts\n",
    "    extracted_ts = dict(data=ts.copy(), lon=lon, lat=lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "815bc02a-29fe-4421-9a72-7633784d9e0d",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Compute anomalies for the chosen location\n",
    "Here we use the time series that was extracted in the previous example.\n",
    "\n",
    "We now compute the average soil moisture conditions---**climatology** ($\\overline{SM}$)---for each month using data from all years within the selected reference period. \n",
    "The deviation of the absolute values ($SM$) from this climatological reference are the soil moisture **anomalies**. Values close to 0 therefore indicate normal conditions, while negative anomalies indicate drier than usual conditions (droughts). In case of strong precipitation events, positive anomalies are found.\n",
    "\n",
    "$\\Large SMA = SM_{k,i} - \\overline{SM_i}$\n",
    "\n",
    "**Try:**\n",
    "    \n",
    "* **Changing the baseline / reference period to see if / how it affects the climatology computation.**\n",
    "* **Changing the amount of smoothing applied to the climatology.**\n",
    "\n",
    "*Note: You can always reset the interactive controls by re-running the cell.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "8339eec9-932b-4abc-85c9-82700c61a184",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "46bd905c484f4d47842e2f1f858da46e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntRangeSlider(value=(1991, 2022), description='Reference period (year from, year to):',…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "anomaly_ts = None\n",
    "\n",
    "@widgets.interact(baseline=widgets.IntRangeSlider(min=1991, max=2022, value=[1991, 2022], step=1, style={'description_width': 'initial'}, \n",
    "                                                  description='Reference period (year from, year to):', layout=widgets.Layout(width='30%')), \n",
    "                  smoothing=widgets.Dropdown(options=['strong', 'medium', 'none'], value='medium', style={'description_width': 'initial'}, description='Smoothing of climatology'))\n",
    "def plot_components(baseline, smoothing):\n",
    "    \"\"\"\n",
    "    Compute climatology and anomalies for the loaded soil moisture time series\n",
    "    \"\"\"\n",
    "    loc = f\"Lon: {extracted_ts['lon']}°, Lat: {extracted_ts['lat']}°\" \n",
    "    \n",
    "    fig, axs = plt.subplots(3, 1, figsize=(10, 7))\n",
    "    ts = extracted_ts['data'].copy().loc['1991-01-01':,:]\n",
    "    \n",
    "    clim_data = ts['soil_moisture'].loc[f'{baseline[0]}-01-01':f'{baseline[1]}-12-31']\n",
    "    clim_std = pd.Series(clim_data.groupby(clim_data.index.month).std(), name='climatology_std')      \n",
    "    clim_mean = pd.Series(clim_data.groupby(clim_data.index.month).mean(), name='climatology')\n",
    "    \n",
    "    if smoothing.lower() == 'medium':\n",
    "        window = 3\n",
    "    elif smoothing.lower() == 'strong':\n",
    "        window = 5\n",
    "    else:\n",
    "        window = 0\n",
    "    \n",
    "    if window != 0:\n",
    "        clim_mean = clim_mean.rolling(window, min_periods=1).mean()\n",
    "        \n",
    "    ts = ts.join(on=ts.index.month, other=clim_mean)\n",
    "    ts = ts.join(on=ts.index.month, other=clim_std)\n",
    "    \n",
    "    ts['anomaly'] = ts['soil_moisture'] - ts['climatology']\n",
    "    \n",
    "    ts['soil_moisture'].plot(ax=axs[0], title=f'Soil Moisture (absolute) at {loc}', ylabel=f'SM $[m³/m³]$', xlabel='Time [Year]')\n",
    "    for i, g in clim_data.groupby(clim_data.index.year):\n",
    "        axs[1].plot(range(1,13), g.values, alpha=0.2)\n",
    "        \n",
    "    clim_mean.plot(ax=axs[1], color='blue', title=f'Climatology at {loc}', ylabel='SM $[m³/m³]$', xlabel='Time [Month]', label='mean')\n",
    "    clim_std.plot(ax=axs[1], label='$\\sigma$')\n",
    "    axs[1].legend()\n",
    "    \n",
    "    global anomaly_ts\n",
    "    anomaly_ts = ts.copy()\n",
    "    \n",
    "    axs[2].axhline(0, color='k')\n",
    "    axs[2].fill_between(ts['anomaly'].index,ts['anomaly'].values,where=ts['anomaly'].values>=0, color='blue')\n",
    "    axs[2].fill_between(ts['anomaly'].index,ts['anomaly'].values,where=ts['anomaly'].values<0, color='red')\n",
    "    axs[2].set_ylabel('SM $[m³/m³]$')\n",
    "    axs[2].set_xlabel('Time [Year]')\n",
    "    axs[2].set_title(f\"SM Anomalies at {loc}\")\n",
    "    plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37bb59bc",
   "metadata": {},
   "source": [
    "Spatial SM anomalies give a good representation of drought events and their evolution. Below monthly SM anomalies can be plotted for a particular year and country in Europe, from the C3S v202212 data set.\n",
    "\n",
    "**Try:**\n",
    "\n",
    "* **Changing the year and study area (e.g., the 2022 summer drought in Italy) to plot soil moisture anomalies for.**\n",
    "* **Finding patterns in the data, explaining outliers or data gaps.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8650c263",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "08a9c8a23ea5491a81e841b54014b00a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(SelectionSlider(continuous_update=False, description='Select a year:', index=12, options…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mask = xr.open_dataset('./grid_eu.nc')['country'].sortby('lat', ascending=False)\n",
    "lut = dict(zip(mask.attrs['flag_meanings'], mask.attrs['flag_values']))\n",
    "\n",
    "def load_study_area(country, ds, var='sm', time_from='2010-01-01'):\n",
    "    subset = ds.sel(time=slice(time_from, None))[var]\n",
    "    sel = mask.where(mask==lut[country], drop=True)\n",
    "    subset = subset.sel(lat=sel.lat, lon=sel.lon, tolerance=.25, method=\"nearest\")\n",
    "    subset.values[:, ~np.isfinite(sel)] = np.nan\n",
    "    return subset\n",
    "\n",
    "country_data = dict()\n",
    "\n",
    "dropdown = widgets.Dropdown(options=utils.europe_countries, value='Poland',description='Country:')\n",
    "slider = widgets.SelectionSlider(options=np.arange(2010,2024), value=2022, description='Select a year:', continuous_update=False)\n",
    "@widgets.interact(year=slider, country=dropdown)\n",
    "def plot_months(year=2022, country='Poland'):\n",
    "    global country_data\n",
    "    \n",
    "    if not (country in country_data.keys()):\n",
    "        country_data = {country: (load_study_area(country, c3s_europe, var='sm_anom'))}\n",
    "\n",
    "    dat = country_data[country].sel(time=slice(f'{year}-01-01', f'{year}-12-31'))\n",
    "    p = dat.plot(transform=ccrs.PlateCarree(),\n",
    "                 col='time', col_wrap=4, vmin=-0.15, vmax=0.15, cmap=plt.get_cmap('RdBu'),\n",
    "                 subplot_kws={'projection': ccrs.PlateCarree()})\n",
    "    for ax in p.axs.flat:\n",
    "        ax.axes.add_feature(cartopy.feature.BORDERS, linewidth=0.5, zorder=2)\n",
    "        ax.axes.add_feature(cartopy.feature.RIVERS, zorder=2)        \n",
    "        ax.axes.coastlines()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97ed3492-c3ca-4756-8d1c-8b4a28092768",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Mapped Quantiles\n",
    "\n",
    "A simple way to classify potential droughts is by applying thresholds based on quantiles of anomalies in a given time series. Extremely low anomaly values therefore correspond to droughts, while anomalies around and above zero do not. \n",
    "\n",
    "Long time series make the analysis more robust, as the sampling becomes denser. The default values are taken from literature.\n",
    "\n",
    "**Try:**\n",
    "\n",
    "* **Changing the theshold values for the quantile-based discretization function and find a classification for the time series you have chosen above.**\n",
    "\n",
    "Which areas have been most affected by the 2022 summer drought?\n",
    "\n",
    "*Note: You can always reset the interactive controls and change the country by re-running the cell.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "56cfa419-ba67-4483-aab0-67c1a8da5890",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "329d4bd69bda4c409e13ca3bcafd05ba",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.35, description=' [quantile]', layout=Layout(width='30%'), max=0.5, …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = anomaly_ts.copy()\n",
    "\n",
    "drought_levels = {0: 'SM deficit', 1: 'Mild Drought', 2: 'Moderate Drought', \n",
    "                  3: 'Significant Drought', 4: 'Severe Drought', 5: 'Extreme Drought', }\n",
    "\n",
    "@widgets.interact(\n",
    "mild=widgets.FloatSlider(value=0.35, min=0, max=0.5, step=0.05, description=' [quantile]', style={'description_width': 'initial'},layout=widgets.Layout(width='30%')),\n",
    "moderate=widgets.FloatSlider(value=0.3, min=0, max=0.5, step=0.05, description=' [quantile]', style={'description_width': 'initial'},layout=widgets.Layout(width='30%')),\n",
    "significant=widgets.FloatSlider(value=0.25, min=0, max=0.5, step=0.05, description=' [quantile]', style={'description_width': 'initial'},layout=widgets.Layout(width='30%')),\n",
    "severe=widgets.FloatSlider(value=0.15, min=0, max=0.5, step=0.05, description=' [quantile]', style={'description_width': 'initial'},layout=widgets.Layout(width='30%')),\n",
    "extreme=widgets.FloatSlider(value=0.1, min=0, max=0.5, step=0.05, description=' [quantile]', style={'description_width': 'initial'}, layout=widgets.Layout(width='30%')),\n",
    "month=widgets.SelectionSlider(options=np.arange(1, 13), value=6, description='Select a 2022 month to plot:', continuous_update=False, style={'description_width': 'initial'}),\n",
    "continuous_update=False)\n",
    "def plot(mild=0.2, moderate=0.1, significant=0.05, severe=0.02, extreme=0.01, month=6):\n",
    "    bins = df['anomaly'].quantile([0, extreme, severe, significant, moderate, mild, .5]).values\n",
    "\n",
    "    fig = plt.figure(figsize=(12, 12))\n",
    "    gs = GridSpec(2, 1, height_ratios=(1,4), hspace=.3)\n",
    "    \n",
    "    ax = fig.add_subplot(gs[0, 0],)\n",
    "    ax.axhline(0, color='k', linestyle='--')\n",
    "    ax.axhspan(bins[4], bins[5], color='#fee5d9', alpha=0.5, label='mild')\n",
    "    ax.axhspan(bins[3], bins[4], color='#fcbba1', alpha=0.5, label='moderate')\n",
    "    ax.axhspan(bins[2], bins[3], color='#fc9272', alpha=0.5, label='significant')\n",
    "    ax.axhspan(bins[1], bins[2], color='#fb6a4a', alpha=0.5, label='severe')\n",
    "    ax.axhspan(bins[0], bins[1], color='#de2d26', alpha=0.5, label='extreme')\n",
    "    \n",
    "    p = df['anomaly'].plot(marker='o', color='black', ax=ax, ylabel='SM Anomaly $[m³/m³]$', \n",
    "                           title=f\"Soil Moisture Anomaly drought classification at Lon: {extracted_ts['lon']}°, Lat: {extracted_ts['lat']}°\")\n",
    "    plt.legend()\n",
    "    \n",
    "    ax_map = fig.add_subplot(gs[1, 0], projection=ccrs.PlateCarree())\n",
    "    ax_map.axes.add_feature(cartopy.feature.BORDERS, linewidth=0.5, zorder=2)\n",
    "    ax_map.axes.add_feature(cartopy.feature.RIVERS, zorder=2) \n",
    "    ax_map.axes.coastlines()\n",
    "    \n",
    "    gl = ax_map.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)\n",
    "\n",
    "    gl.top_labels = False\n",
    "    gl.right_labels = False\n",
    "    gl.xformatter = LONGITUDE_FORMATTER\n",
    "    gl.yformatter = LATITUDE_FORMATTER\n",
    "    \n",
    "    selected_country = list(country_data.keys())[0]\n",
    "    month_data = country_data[selected_country].sel(time=datetime.datetime(2022, month, 1))\n",
    "    mapplot = month_data.plot(transform=ccrs.PlateCarree(), levels=np.round(bins, 2), cmap=plt.get_cmap('hot'), \n",
    "        subplot_kws={'projection': ccrs.PlateCarree()}, add_colorbar=False)\n",
    "    \n",
    "    binned = np.digitize(month_data, bins)\n",
    "    print(\"Percentage area of drough levels [%]:\"\n",
    "          \"\\n=====================================\")\n",
    "    for b in range(len(drought_levels)):\n",
    "        percent_area = np.sum(binned==(b + 1))/binned.size * 100\n",
    "        print(f\"{drought_levels[b]} {percent_area:.1f}\")\n",
    "    \n",
    "    cax = fig.add_axes([ax_map.get_position().x0, ax_map.get_position().y0-.05, ax_map.get_position().width, 0.02])\n",
    "    cbar = fig.colorbar(mapplot, ax=ax, cax=cax, orientation=\"horizontal\")\n",
    "    cbar.set_label('Selected drought levels (in SM anomalies, m³/m³)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "594a0d04-122d-4869-9574-993afb464005",
   "metadata": {},
   "source": [
    "## Z-Scores\n",
    "Z-scores are a way of standardising values from different normal distributions. Z-scores show the number of standard deviations from the mean of the sample.\n",
    "\n",
    "$\\Large Z = \\frac{x - \\mu}{\\sigma}$\n",
    "\n",
    "In the case of soil moisture anomalies, the distributions of values within a year can vary, e.g., due to distinct seasonal hydrological regimes, or due to differences in data coverage (i.e. a different distribution for values in summer than for those in winter). This can affect the computed climatology and therefore the derived anomalies. By taking into account the standard deviation of samples used to compute the climatologies, we can normalize the anomalies and reduce intra-annual biases.\n",
    "\n",
    "$\\Large Z_{k,i} = \\frac{SM_{k,i} - \\overline{SM_i}}{\\sigma_{i}}$\n",
    "\n",
    "The default values in the example are taken from literature.\n",
    "\n",
    "**Try:**\n",
    "\n",
    "* **Compare the Z-score time series to the anomaly time series above.**\n",
    "* **Change the threshold values for the quantile-based discretization function find a meaningful classification for the time series you have chosen above.**\n",
    "* **Extract a time series in areas where data coverage varies strongly within a year and compare the output.**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "43b4a52f-8aff-4c70-85ff-c46ea98d2517",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a63acf7fdf974131abe822ae550e0700",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.0, description='Moderate Drought $[\\\\sigma]$', max=10.0, step=0.25, …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = anomaly_ts.copy()\n",
    "\n",
    "@widgets.interact(\n",
    "moderate=widgets.FloatSlider(value=1, min=0, max=10, step=0.25, description='Moderate Drought $[\\sigma]$', style={'description_width': 'initial'},),\n",
    "extreme=widgets.FloatSlider(value=2, min=0, max=10, step=0.25, description='Extreme Drought $[\\sigma]$', style={'description_width': 'initial'},),\n",
    "continuous_update=False)\n",
    "def plot(moderate, extreme):\n",
    "    df.loc[:, 'zscore'] = (df['soil_moisture'] - df['climatology']) / df['climatology_std']\n",
    "    bins = [df['zscore'].min(), -extreme, -moderate, 0]\n",
    "\n",
    "    fig, ax = plt.subplots(1,1, figsize=(15, 5))\n",
    "    ax_anom = ax.twinx()\n",
    "    ax.axhline(0, color='k', linestyle='--')\n",
    "    ax.axhspan(bins[1], bins[2], color='#fc9272', alpha=0.5, label='moderate')\n",
    "    ax.axhspan(bins[0], bins[1], color='#de2d26', alpha=0.5, label='extreme')\n",
    "    df['zscore'].plot(marker='o', color='black', ax=ax, ylabel='z score [-]', \n",
    "                      title=f\"Soil Moisture Z-Score drought classification (Lon: {extracted_ts['lon']}°, Lat: {extracted_ts['lat']}°)\")\n",
    "    ax_anom.plot(df['anomaly'].index, df['anomaly'], color='blue', linewidth=1, label=\"anomaly\")\n",
    "    ax_anom.set_ylabel('SM Anomaly [m³/m³]')\n",
    "    \n",
    "    ax.legend(loc=\"upper left\")\n",
    "    ax_anom.legend(loc=\"upper right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c57a085-1fc6-41cf-9b45-f28046710307",
   "metadata": {},
   "source": [
    "## SMASI\n",
    "\n",
    "One index that was shown in the lecture is the Standardised Vegetation Optical Depth Index [(Moesinger et al., 2022)](https://doi.org/10.5194/bg-19-5107-2022).\n",
    "\n",
    "Similarly, a Soil Moisture Anomalies Standardized Index (SMASI) can be obtained with a similar method. Compared to the C3S SM-derived indices shown before, SMASI is more robust to noise variations between merging periods (i.e., variations of the number of sensors in the merging). \n",
    "\n",
    "The main innovation is represented by the fact that the empirical multivariate distribution of the indices is scaled to a uniform distribution, which renders the time series more consistent in time. The principle is shown in Figure 2 of Moesinger et al. (2022):\n",
    "\n",
    "<img src=\"img/svodi_bivariate_distribution.png\" width=\"90%\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a18ee82e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
       " *\n",
       " */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
       "  --xr-background-color: var(--jp-layout-color0, white);\n",
       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
       "}\n",
       "\n",
       "html[theme=dark],\n",
       "body[data-theme=dark],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: rgba(255, 255, 255, 1);\n",
       "  --xr-font-color2: rgba(255, 255, 255, 0.54);\n",
       "  --xr-font-color3: rgba(255, 255, 255, 0.38);\n",
       "  --xr-border-color: #1F1F1F;\n",
       "  --xr-disabled-color: #515151;\n",
       "  --xr-background-color: #111111;\n",
       "  --xr-background-color-row-even: #111111;\n",
       "  --xr-background-color-row-odd: #313131;\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block !important;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-array-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: '►';\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: '▼';\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "  padding-bottom: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: '(';\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: ')';\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: ',';\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-index-preview {\n",
       "  grid-column: 2 / 5;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data,\n",
       ".xr-index-data {\n",
       "  display: none;\n",
       "  background-color: var(--xr-background-color) !important;\n",
       "  padding-bottom: 5px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data,\n",
       ".xr-index-data-in:checked ~ .xr-index-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-index-name div,\n",
       ".xr-index-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data,\n",
       ".xr-index-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2,\n",
       ".xr-no-icon {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt;\n",
       "Dimensions:         (lon: 240, time: 365, lat: 160)\n",
       "Coordinates:\n",
       "  * lon             (lon) float64 -9.875 -9.625 -9.375 ... 49.38 49.62 49.88\n",
       "  * time            (time) datetime64[ns] 2022-01-01 2022-01-02 ... 2022-12-31\n",
       "  * lat             (lat) float64 69.88 69.62 69.38 69.12 ... 30.62 30.38 30.12\n",
       "    reference_time  (time) datetime64[ns] ...\n",
       "Data variables:\n",
       "    SMASI           (time, lat, lon) float64 ...</pre><div class='xr-wrap' style='display:none'><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-b1b7612c-6e2b-4006-a13e-35fd41d84fea' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-b1b7612c-6e2b-4006-a13e-35fd41d84fea' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span class='xr-has-index'>lon</span>: 240</li><li><span class='xr-has-index'>time</span>: 365</li><li><span class='xr-has-index'>lat</span>: 160</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-446992a3-8d3b-4d2d-83b2-2608a70f3262' class='xr-section-summary-in' type='checkbox'  checked><label for='section-446992a3-8d3b-4d2d-83b2-2608a70f3262' class='xr-section-summary' >Coordinates: <span>(4)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>lon</span></div><div class='xr-var-dims'>(lon)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>-9.875 -9.625 ... 49.62 49.88</div><input id='attrs-5993b989-4412-4a36-87a8-95de5663f3a3' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-5993b989-4412-4a36-87a8-95de5663f3a3' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-ceef539b-294a-42de-bec0-ff2b02e4c441' class='xr-var-data-in' type='checkbox'><label for='data-ceef539b-294a-42de-bec0-ff2b02e4c441' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([-9.875, -9.625, -9.375, ..., 49.375, 49.625, 49.875])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>time</span></div><div class='xr-var-dims'>(time)</div><div class='xr-var-dtype'>datetime64[ns]</div><div class='xr-var-preview xr-preview'>2022-01-01 ... 2022-12-31</div><input id='attrs-9b2df63b-5b35-43e1-a049-352266015838' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-9b2df63b-5b35-43e1-a049-352266015838' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-d37e75e4-6b18-4379-9a22-ac9e429b0234' class='xr-var-data-in' type='checkbox'><label for='data-d37e75e4-6b18-4379-9a22-ac9e429b0234' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([&#x27;2022-01-01T00:00:00.000000000&#x27;, &#x27;2022-01-02T00:00:00.000000000&#x27;,\n",
       "       &#x27;2022-01-03T00:00:00.000000000&#x27;, ..., &#x27;2022-12-29T00:00:00.000000000&#x27;,\n",
       "       &#x27;2022-12-30T00:00:00.000000000&#x27;, &#x27;2022-12-31T00:00:00.000000000&#x27;],\n",
       "      dtype=&#x27;datetime64[ns]&#x27;)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>lat</span></div><div class='xr-var-dims'>(lat)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>69.88 69.62 69.38 ... 30.38 30.12</div><input id='attrs-fc2c2069-969a-479b-8e51-1911912b2b69' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-fc2c2069-969a-479b-8e51-1911912b2b69' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-49389bcd-f4f7-44af-bc37-d3da9250fd5f' class='xr-var-data-in' type='checkbox'><label for='data-49389bcd-f4f7-44af-bc37-d3da9250fd5f' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([69.875, 69.625, 69.375, 69.125, 68.875, 68.625, 68.375, 68.125, 67.875,\n",
       "       67.625, 67.375, 67.125, 66.875, 66.625, 66.375, 66.125, 65.875, 65.625,\n",
       "       65.375, 65.125, 64.875, 64.625, 64.375, 64.125, 63.875, 63.625, 63.375,\n",
       "       63.125, 62.875, 62.625, 62.375, 62.125, 61.875, 61.625, 61.375, 61.125,\n",
       "       60.875, 60.625, 60.375, 60.125, 59.875, 59.625, 59.375, 59.125, 58.875,\n",
       "       58.625, 58.375, 58.125, 57.875, 57.625, 57.375, 57.125, 56.875, 56.625,\n",
       "       56.375, 56.125, 55.875, 55.625, 55.375, 55.125, 54.875, 54.625, 54.375,\n",
       "       54.125, 53.875, 53.625, 53.375, 53.125, 52.875, 52.625, 52.375, 52.125,\n",
       "       51.875, 51.625, 51.375, 51.125, 50.875, 50.625, 50.375, 50.125, 49.875,\n",
       "       49.625, 49.375, 49.125, 48.875, 48.625, 48.375, 48.125, 47.875, 47.625,\n",
       "       47.375, 47.125, 46.875, 46.625, 46.375, 46.125, 45.875, 45.625, 45.375,\n",
       "       45.125, 44.875, 44.625, 44.375, 44.125, 43.875, 43.625, 43.375, 43.125,\n",
       "       42.875, 42.625, 42.375, 42.125, 41.875, 41.625, 41.375, 41.125, 40.875,\n",
       "       40.625, 40.375, 40.125, 39.875, 39.625, 39.375, 39.125, 38.875, 38.625,\n",
       "       38.375, 38.125, 37.875, 37.625, 37.375, 37.125, 36.875, 36.625, 36.375,\n",
       "       36.125, 35.875, 35.625, 35.375, 35.125, 34.875, 34.625, 34.375, 34.125,\n",
       "       33.875, 33.625, 33.375, 33.125, 32.875, 32.625, 32.375, 32.125, 31.875,\n",
       "       31.625, 31.375, 31.125, 30.875, 30.625, 30.375, 30.125])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>reference_time</span></div><div class='xr-var-dims'>(time)</div><div class='xr-var-dtype'>datetime64[ns]</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-0b41c398-aee3-4ad1-8785-d1609370127b' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-0b41c398-aee3-4ad1-8785-d1609370127b' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-6877c9ae-13a5-4710-8821-ba17019e64b7' class='xr-var-data-in' type='checkbox'><label for='data-6877c9ae-13a5-4710-8821-ba17019e64b7' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>[365 values with dtype=datetime64[ns]]</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-8121ddc4-4eee-4454-8592-47f5585814fb' class='xr-section-summary-in' type='checkbox'  checked><label for='section-8121ddc4-4eee-4454-8592-47f5585814fb' class='xr-section-summary' >Data variables: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>SMASI</span></div><div class='xr-var-dims'>(time, lat, lon)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-ae6567ac-13d1-4e95-ad94-b20b3f23a055' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-ae6567ac-13d1-4e95-ad94-b20b3f23a055' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-e1cc21cf-2fe3-4602-b9e8-422c015eb097' class='xr-var-data-in' type='checkbox'><label for='data-e1cc21cf-2fe3-4602-b9e8-422c015eb097' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>Soil Moisture Anomaly Standardised Index</dd></dl></div><div class='xr-var-data'><pre>[14016000 values with dtype=float64]</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-582b950f-658e-4bff-80c1-74a06477569c' class='xr-section-summary-in' type='checkbox'  ><label for='section-582b950f-658e-4bff-80c1-74a06477569c' class='xr-section-summary' >Indexes: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-index-name'><div>lon</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-2e1226ae-bf10-49e6-8e02-6801d0193e2c' class='xr-index-data-in' type='checkbox'/><label for='index-2e1226ae-bf10-49e6-8e02-6801d0193e2c' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(Float64Index([-9.875, -9.625, -9.375, -9.125, -8.875, -8.625, -8.375, -8.125,\n",
       "              -7.875, -7.625,\n",
       "              ...\n",
       "              47.625, 47.875, 48.125, 48.375, 48.625, 48.875, 49.125, 49.375,\n",
       "              49.625, 49.875],\n",
       "             dtype=&#x27;float64&#x27;, name=&#x27;lon&#x27;, length=240))</pre></div></li><li class='xr-var-item'><div class='xr-index-name'><div>time</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-f7db1e2e-fa3f-48a9-b993-e46f21f07cfb' class='xr-index-data-in' type='checkbox'/><label for='index-f7db1e2e-fa3f-48a9-b993-e46f21f07cfb' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(DatetimeIndex([&#x27;2022-01-01&#x27;, &#x27;2022-01-02&#x27;, &#x27;2022-01-03&#x27;, &#x27;2022-01-04&#x27;,\n",
       "               &#x27;2022-01-05&#x27;, &#x27;2022-01-06&#x27;, &#x27;2022-01-07&#x27;, &#x27;2022-01-08&#x27;,\n",
       "               &#x27;2022-01-09&#x27;, &#x27;2022-01-10&#x27;,\n",
       "               ...\n",
       "               &#x27;2022-12-22&#x27;, &#x27;2022-12-23&#x27;, &#x27;2022-12-24&#x27;, &#x27;2022-12-25&#x27;,\n",
       "               &#x27;2022-12-26&#x27;, &#x27;2022-12-27&#x27;, &#x27;2022-12-28&#x27;, &#x27;2022-12-29&#x27;,\n",
       "               &#x27;2022-12-30&#x27;, &#x27;2022-12-31&#x27;],\n",
       "              dtype=&#x27;datetime64[ns]&#x27;, name=&#x27;time&#x27;, length=365, freq=None))</pre></div></li><li class='xr-var-item'><div class='xr-index-name'><div>lat</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-525247f3-422d-4b61-b59d-607bfe025752' class='xr-index-data-in' type='checkbox'/><label for='index-525247f3-422d-4b61-b59d-607bfe025752' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(Float64Index([69.875, 69.625, 69.375, 69.125, 68.875, 68.625, 68.375, 68.125,\n",
       "              67.875, 67.625,\n",
       "              ...\n",
       "              32.375, 32.125, 31.875, 31.625, 31.375, 31.125, 30.875, 30.625,\n",
       "              30.375, 30.125],\n",
       "             dtype=&#x27;float64&#x27;, name=&#x27;lat&#x27;, length=160))</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-8a9279ff-fa3f-47e7-8157-babb9bdd607b' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-8a9279ff-fa3f-47e7-8157-babb9bdd607b' class='xr-section-summary'  title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:         (lon: 240, time: 365, lat: 160)\n",
       "Coordinates:\n",
       "  * lon             (lon) float64 -9.875 -9.625 -9.375 ... 49.38 49.62 49.88\n",
       "  * time            (time) datetime64[ns] 2022-01-01 2022-01-02 ... 2022-12-31\n",
       "  * lat             (lat) float64 69.88 69.62 69.38 69.12 ... 30.62 30.38 30.12\n",
       "    reference_time  (time) datetime64[ns] ...\n",
       "Data variables:\n",
       "    SMASI           (time, lat, lon) float64 ..."
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "smasi_europe = xr.open_dataset('./LTC_DATA/STACK_C3S-SOILMOISTURE_EUROPE_SMASI_5x_DAILY.nc')\n",
    "display(smasi_europe)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1816afe3",
   "metadata": {},
   "source": [
    "The SMASI data set used here contains a harmonised index from 5 satellites, covering the 2007-2022 period. The main differences (synthetically) are that:\n",
    "\n",
    "* The standardization is performed on a day-of-year basis, with a rolling 90-days window.\n",
    "* The multivariate distribution is scaled, as mentioned above\n",
    "\n",
    "The plot below is showing the percentag area of Poland with a under(over)-saturated SM, corresponding to SMASI values below(over) the selected threshold. For refernce, the median (smoothed) SMASI over the area is given.\n",
    "\n",
    "**Try:**\n",
    "\n",
    "* **Changing the magnitude of the index threshold to see the affected area in Poland**\n",
    "* **Changing the level of smoothing to see if intra/inter-annual patterns emerge**\n",
    "\n",
    "Which years show the most extreme events? Can you see temporal patterns and are they as you would expect them? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "98fc77fa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4d5ed4d4b66745afbedfaa4fb954ca86",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.5, continuous_update=False, description='Index threshold [-]:', layo…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "slider = widgets.FloatSlider(min=.5, max=2.5, value=1.5, step=0.5, style={'description_width': 'initial'}, description='Index threshold [-]:', layout=widgets.Layout(width='30%'), continuous_update=False)\n",
    "slider2 = widgets.IntSlider(min=1, max=36, value=4, step=2, style={'description_width': 'initial'}, description='Median smoothing (# months):', layout=widgets.Layout(width='30%'), continuous_update=False)\n",
    "@widgets.interact(thresholds=slider, smoothing=slider2)\n",
    "def plot(thresholds, smoothing):\n",
    "    df = pd.read_csv(f'./LTC_DATA/mean_ts/smasi_poland.csv', index_col=0, parse_dates=True)\n",
    "    df = df.rolling(smoothing * 30, min_periods=0).median()\n",
    "    fig, ax = plt.subplots(1,1, figsize=(15, 7))\n",
    "    \n",
    "    ax.fill_between(df.index, -df[f\"SMASI_below_{.5:.1f}\"], df[f\"SMASI_above_{.5:.1f}\"], color=\"#F9AFBC\", label=f\"|SMASI| > {.5}\")\n",
    "    ax.fill_between(df.index, -df[f\"SMASI_below_{thresholds:.1f}\"], df[f\"SMASI_above_{thresholds:.1f}\"], color=\"#6495ED\", label=f\"|SMASI| > {thresholds}\")\n",
    "    ax.set_ylabel(\"Percentage area (%)\")\n",
    "    ax.set_ylim([-100, 100])\n",
    "    ax.set_yticks(np.linspace(-100, 100, 11))\n",
    "    ax.set_yticklabels(np.abs(np.linspace(-100, 100, 11)).astype(int))\n",
    "    \n",
    "    ax_mean = ax.twinx()\n",
    "    df[\"SMASI\"].plot(ax=ax_mean, c=\"blue\", linestyle=\"--\", linewidth=1.5, label=\"SMASI smoothed median\")\n",
    "    ax_mean.set_ylim([-2, 2])\n",
    "    \n",
    "    ax.set_xlim([datetime.datetime(2010,1,1), datetime.datetime(2023,1,1)])\n",
    "    \n",
    "    ax.set_title(\"Under/Over-saturated area in Poland\")\n",
    "    ax.hlines(0, datetime.datetime(2010,1,1), datetime.datetime(2023,1,1), color=\"r\", lw=1)\n",
    "    ax.legend(loc=\"lower right\")\n",
    "    ax_mean.legend(loc=\"upper left\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4915bb24",
   "metadata": {},
   "source": [
    "We can now compare SMASI to the SM Anomalies over a selected country area to see what is their agreement. The values are averaged in the Latitude dimension to ease the visualization.\n",
    "\n",
    "Bear in mind that SMASI is provided in daily resolution, whereas the C3S SM data used here is at monthly intervals.\n",
    "\n",
    "**Try:**\n",
    "\n",
    "* **Changing the area below to understand whether SMASI and regular SM anomalies agree with each other.**\n",
    "\n",
    "Where do the differences between the two originate?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "8a23467d",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "143d01a8ba564f85a74c24f1547eb465",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='Country:', index=27, options=('Albania', 'Austria', 'Belarus', 'Be…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dropdown = widgets.Dropdown(options=utils.europe_countries, value='Poland',description='Country:')\n",
    "slider = widgets.SelectionSlider(options=np.arange(2010,2024), value=2022, description='Select a year:', continuous_update=False)\n",
    "@widgets.interact(year=slider, country=dropdown)\n",
    "def plot_hovmoeller(country='Poland'):\n",
    "    global country_data\n",
    "    \n",
    "    if not (country in country_data.keys()):\n",
    "        country_data = {country: (load_study_area(country, c3s_europe, var='sm_anom'))}\n",
    "    \n",
    "    dat = country_data[country]\n",
    "    dat_2022 = dat.sel(time=slice(datetime.datetime(2021,12,31), datetime.datetime(2023,1,1)))\n",
    "    smasi_dat = load_study_area(country, smasi_europe, var='SMASI')\n",
    "    dat_2022 = dat_2022.reindex_like(smasi_dat, method=\"ffill\")  # match the time of the SMASI and SM Anom data sets\n",
    "    \n",
    "    fig = plt.figure(figsize=(12, 10))\n",
    "    gs = GridSpec(3, 1, height_ratios=(1,3,3), hspace=.2)\n",
    "    \n",
    "    ax = fig.add_subplot(gs[0, 0],)\n",
    "    ax.set_xticks([])\n",
    "    ax_hov_smasi = fig.add_subplot(gs[1, 0],)\n",
    "    ax_hov_smasi.set_xticks([])    \n",
    "    ax_hov_anom = fig.add_subplot(gs[2, 0],)\n",
    "    \n",
    "    mean = dat_2022.mean((\"lon\", \"lat\"))\n",
    "    std = dat_2022.std((\"lon\", \"lat\"))\n",
    "    \n",
    "    mean_smasi = smasi_dat.mean((\"lon\", \"lat\"))\n",
    "    std_smasi = smasi_dat.std((\"lon\", \"lat\"))\n",
    "    \n",
    "    ax.plot(mean.time, mean, label=\"SM Anomalies\")\n",
    "    ax.fill_between(mean.time, mean - std, mean + std, alpha=.4)\n",
    "    \n",
    "    ax_anom = ax.twinx()\n",
    "    ax_anom.plot(mean_smasi.time, mean_smasi, label=\"SMASI\", color=\"orange\")\n",
    "    ax_anom.fill_between(mean_smasi.time, mean_smasi - std_smasi, mean_smasi + std_smasi, alpha=.4, color=\"orange\")\n",
    "    ax.legend(loc=\"upper left\")\n",
    "    ax_anom.legend(loc=\"upper right\")\n",
    "    \n",
    "    ax.set_ylabel(\"SM Anomalies [m³/m³]\")\n",
    "    ax_anom.set_ylabel('SMASI [-]')\n",
    "    \n",
    "    smasi_plot = smasi_dat.mean(\"lon\").transpose(\"lat\", \"time\").plot(\n",
    "        ax=ax_hov_smasi, add_colorbar=False, cmap=\"RdBu\")\n",
    "    ax_hov_smasi.set_xlabel('')\n",
    "    ax_hov_smasi.set_ylabel('Latitude (deg.)')\n",
    "    \n",
    "    anom_plot = dat_2022.mean(\"lon\").transpose(\"lat\", \"time\").plot(\n",
    "        ax=ax_hov_anom, add_colorbar=False, cmap=\"RdBu\")\n",
    "    ax_hov_anom.set_ylabel('Latitude (deg.)')\n",
    "    \n",
    "    ax_hov_anom.set_title(\"SM Anomalies 2022 evolution\")\n",
    "    ax_hov_smasi.set_title(\"SMASI 2022 evolution\")\n",
    "    \n",
    "    cax = fig.add_axes([ax_hov_anom.get_position().x0, ax_hov_anom.get_position().y0-.1, ax_hov_anom.get_position().width/3, 0.02])\n",
    "    cbar = fig.colorbar(smasi_plot, ax=ax, cax=cax, orientation=\"horizontal\")\n",
    "    cbar.set_label('SMASI [-]')\n",
    "    \n",
    "    cax_anom = fig.add_axes([ax_hov_anom.get_position().x0 + ax_hov_anom.get_position().width/3*2, ax_hov_anom.get_position().y0-.1, ax_hov_anom.get_position().width/3, 0.02])\n",
    "    cbar_anom = fig.colorbar(anom_plot, ax=ax, cax=cax_anom, orientation=\"horizontal\")\n",
    "    cbar_anom.set_label('SM Anomalies [m³/m³]')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9073ec27-acf5-4922-aa7b-e441b246b302",
   "metadata": {},
   "source": [
    "## Standardized Precipitation Index\n",
    "\n",
    "Here we compute SPI values for the previously extracted precipitation time series (see first cell in this chapter).\n",
    "SPI can be computed for different integration time lengths.\n",
    "\n",
    "SPI values can be aggregated to get an estimate for the duration, severity and intensity of a drought.\n",
    "\n",
    "1) **Duration**: Number of consecutive time stamps where SPI < 0\n",
    "2) **Severity**: $\\sum{SPI}$ over the duration of a drought\n",
    "3) **Intensity**: $ \\frac{Severity}{Duration}$\n",
    "\n",
    "**Try:**\n",
    "* **Changing the integration time length and see how it affects the SPI time series**\n",
    "* **Choosing a different location in the data extraction cell above, re-run this cell and see how SPI changes**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "382f6cb7-8219-436b-b75c-d8f1fd79b77c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "690b5a123e3e4ebba75cacd13e7f6336",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='Integration Time [Months]:', index=4, options=(0, 1, 2, 3, 12, 48)…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import gamma, norm\n",
    "\n",
    "ts = extracted_ts['data'].copy().to_xarray()\n",
    "\n",
    "def calc_zscore(ts):\n",
    "    fit_alpha, fit_loc, fit_beta = gamma.fit(ts)\n",
    "    pvals = gamma(fit_alpha, loc=fit_loc, scale=fit_beta).cdf(ts)\n",
    "    zscore = norm.ppf(pvals)\n",
    "    return zscore\n",
    "\n",
    "def calc_spi_integrationtime(ts, integrationtime=None):\n",
    "    if (integrationtime is not None) and (integrationtime!=0):\n",
    "        ts = ts.rolling(time=integrationtime).mean().dropna(\"time\")\n",
    "    spi = ts.groupby('time.month').map(calc_zscore)\n",
    "    return spi\n",
    "\n",
    "def index_metrics(index):\n",
    "    duration = [0]\n",
    "    severity = [0]\n",
    "\n",
    "    for i in index:\n",
    "        if i < 0:\n",
    "            duration.append(duration[-1] + 1)\n",
    "            severity.append(severity[-1] + i)\n",
    "        else:\n",
    "            duration.append(0)\n",
    "            severity.append(0)\n",
    "    \n",
    "    duration = np.array(duration[1:])\n",
    "    severity = np.array(severity[1:])\n",
    "    \n",
    "    return duration, severity\n",
    "\n",
    "\n",
    "@widgets.interact(integration_time=widgets.Dropdown(options=[0, 1, 2, 3, 12, 48], value=12, description='Integration Time [Months]:', style={'description_width': 'initial'}))\n",
    "def plot_spi(integration_time):\n",
    "    spi = pd.DataFrame(data={'index': calc_spi_integrationtime(ts['precipitation'], integration_time).to_pandas()})\n",
    "    \n",
    "    duration, severity = index_metrics(spi['index'])\n",
    "    \n",
    "    spi.loc[:, 'duration'] = duration\n",
    "    spi.loc[:, 'severity'] = severity\n",
    "    spi.loc[:, 'intensity'] = spi['severity'] / spi['duration']\n",
    "    \n",
    "    fig, axs = plt.subplots(2, 2, figsize=(12, 5))\n",
    "    \n",
    "    axs[0, 0].fill_between(spi['index'].index, spi['index'].values, where=spi['index'].values>=0, color='blue')\n",
    "    axs[0, 0].fill_between(spi['index'].index, spi['index'].values, where=spi['index'].values<0, color='red')\n",
    "    axs[0, 0].set_ylabel('SPI[-]')\n",
    "    axs[0, 0].set_title(f\"SPI-{integration_time} at Lon: {extracted_ts['lon']}°, Lat: {extracted_ts['lat']}°\")\n",
    "    \n",
    "    spi['intensity'].fillna(0).plot(title=f'Intensity', ax=axs[1, 0], ylabel='Severity / Duration')\n",
    "    spi['duration'].plot(title=f'Duration', ax=axs[0, 1], ylabel='months')\n",
    "    spi['severity'].plot(title=f'Severity', ax=axs[1, 1], ylabel=\"$\\sum{SPI}$\")\n",
    "    plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "865c3ae2",
   "metadata": {},
   "source": [
    "# SMASI and SPI comparison over Poland\n",
    "\n",
    "The various drought indicators (Duration, Intensity, Severity) can be computed using a standardized SM anomaly index like SMASI, the same way it is done for precipitation. The relationship of the two can be explored by considering the averaged monthly SMASI and SPI values over Poland for the 2011-2022 period.\n",
    "\n",
    "**Try:**\n",
    "* **Looking at the various drought indicators for the two indices**\n",
    "\n",
    "Are the two indices comparable? What can be observed on the precipitation/soil moisture deficit relationship?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6682c9d6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2ee42f00dbca4eae809347d0781dafb7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='Drought indicator:', options=('Duration', 'Intensity', 'Severity')…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@widgets.interact(indicator=widgets.Dropdown(options=[\"Duration\", \"Intensity\", \"Severity\"], value='Duration',description='Drought indicator:'))\n",
    "def plot_spi_smasi_comparison(indicator):\n",
    "    tp_poland = pd.read_csv(f'./LTC_DATA/mean_ts/tp_poland.csv', index_col=0, parse_dates=True)[\"tp\"].to_xarray()\n",
    "    spi_df = pd.DataFrame(data={'index': calc_spi_integrationtime(tp_poland, 1).to_pandas()})  \n",
    "    \n",
    "    smasi_df = pd.read_csv(f'./LTC_DATA/mean_ts/smasi_poland.csv', index_col=0, parse_dates=True).resample(\"1M\").mean()\n",
    "    \n",
    "    duration, severity = index_metrics(spi_df[\"index\"])\n",
    "    \n",
    "    spi_df.loc[:, 'duration'] = duration\n",
    "    spi_df.loc[:, 'severity'] = severity\n",
    "    spi_df.loc[:, 'intensity'] = spi_df['severity'] / spi_df['duration']\n",
    "    \n",
    "    duration, severity = index_metrics(smasi_df[\"SMASI\"])\n",
    "    \n",
    "    smasi_df.loc[:, 'duration'] = duration\n",
    "    smasi_df.loc[:, 'severity'] = severity\n",
    "    smasi_df.loc[:, 'intensity'] = smasi_df['severity'] / smasi_df['duration']\n",
    "    \n",
    "    units_lut = dict(duration=\"Months\", intensity=\"$\\sum{SPI}$\", severity=\"Severity / Duration\")\n",
    "    \n",
    "    fig, axs = plt.subplots(1, 1, figsize=(12, 4))    \n",
    "    spi_df[indicator.lower()].fillna(0).plot(title=indicator, ax=axs, ylabel=units_lut[indicator.lower()], label=\"SPI\")\n",
    "    smasi_df[indicator.lower()].fillna(0).plot(title=indicator, ax=axs, label=\"SMASI\")\n",
    "    \n",
    "    axs.set_xlim([datetime.datetime(2010,1,1), datetime.datetime(2023,1,1)])\n",
    "    \n",
    "    plt.legend()\n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "def578cc-0d87-4e97-a5e6-f002901b5be5",
   "metadata": {},
   "source": [
    "# Soil Moisture / VOD Lag analysis\n",
    "\n",
    "For this final example use pre-processed, daily, average VOD and soil moisture time series for larger study areas (bounding boxes around Poland, France, Italy, Spain and Europe). We then compute cross-correlation values between the soil moisture and VOD time series with a number of lags applied to the soil moisture data. \n",
    "\n",
    "<img src=\"img/study_areas.png\" width=\"50%\"/>\n",
    "\n",
    "**Try:**\n",
    "* **Changing the study area and see if you find differences in the computed cross correlations**.\n",
    "* **Changing the range of lags that is applied to the soil moisture data (positive values refer to shifts forward in time)**\n",
    "* **What causes these patterns and when do they repeat?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "08936894-1ece-4262-9417-a5a140615a19",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "764c97b6bff14dcfbc965d13fb2a3500",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='Study area:', index=4, options=('Europe', 'France', 'Spain', 'Ital…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def crosscorr(vod, sm, lag=0):\n",
    "    return vod.corr(sm.shift(lag))\n",
    "\n",
    "slider = widgets.IntRangeSlider(min=-366, max=366, value=[0, 200], step=1, style={'description_width': 'initial'}, description='Lags [days] (min, max):', layout=widgets.Layout(width='30%'), continuous_update=False)\n",
    "dropdown = widgets.Dropdown(options=['Europe', 'France', 'Spain', 'Italy', 'Poland'], value='Poland', description='Study area:', style={'description_width': 'initial'})\n",
    "@widgets.interact(study_area=dropdown,lags=slider)\n",
    "def plot(study_area, lags, Anomalies: bool = False):\n",
    "    vod_var, sm_var = 'anomaly' if Anomalies else 'vod_K', 'sm_anomaly' if Anomalies else 'sm'\n",
    "    df_vod = pd.read_csv(f'./LTC_DATA/mean_ts/vod_{study_area.lower()}.csv', index_col=0, parse_dates=True)[vod_var]\n",
    "    df_sm = pd.read_csv(f'./LTC_DATA/mean_ts/c3s_{study_area.lower()}.csv', index_col=0, parse_dates=True)[sm_var]\n",
    "    df = pd.concat([df_vod.loc['1988-01-01':'2020-12-31'], df_sm['1988-01-01':'2020-12-31']], axis=1)\n",
    "    x = np.arange(lags[0], lags[1], 1)\n",
    "    plt.plot(x, np.array([crosscorr(df[vod_var], df[sm_var], lag=lag) for lag in x]))\n",
    "    plt.title(f\"SM / VOD {'anomaly' if Anomalies else ''} cross-correlation ({study_area})\")\n",
    "    plt.xlabel('Lag [days]')\n",
    "    plt.ylabel('R[-]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c68a7619",
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}