notebooks/Sentinel-2/sentinel2_leaf_area_index.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Leaf Area Index Retrieval\n",
    "\n",
    "Author: [Lukas Valentin Graf](https://github.com/lukasValentin/lukasValentin) (2023)\n",
    "\n",
    "## Learning Objectives\n",
    "\n",
    "In this notebook you will learn how to\n",
    "\n",
    "* retrieve the green leaf area index from Sentinel-2 imagery using PROSAIL\n",
    "* get some background information about the green leaf area index.\n",
    "\n",
    "## Tutorial Pre-Requisites\n",
    "\n",
    "This tutorial is advanced level. If you are unfamiliar with `EOdal`, better start [with the basic introduction](/notebooks/General/Loading_Images.ipynb) and then continue [here](/notebooks/General/EOdal_Mapper.ipynb).\n",
    "\n",
    "\n",
    "## Background\n",
    "The leaf area index (LAI) - defined as ratio of one-sided leaf area to unit ground area - is an important functional trait that is a proxy for the net assimilation rate in crops [Watson, 1958](https://www.jstor.org/stable/42907407). LAI can be split in two components: Green LAI (GLAI) and brown LAI (BLAI). GLAI and BLAI form the total LAI. GLAI is the LAI of the photosynthetically-active leaf area, while BLAI refers to non-photosynthetically leaf area, i.e., senescent or necrotic leaves.\n",
    "\n",
    "LAI correlates linearly with the fraction of incident solar radition intercepted by a canopy [Monsi and Saeki, 1953](https://doi.org/10.1093/aob/mci052) and can therefore be estimated from remotely sensed imagery using [canopy radiative transfer models (RTM)](https://edepot.wur.nl/210943).\n",
    "\n",
    "In this notebook, we will present a methodology using the four-stream RTM [PROSAIL](http://teledetection.ipgp.jussieu.fr/prosail/) to retrieve the GLAI for a winter wheat canopy during the stem elongation (BBCH 30) to heading phase (BBCH 59) from Sentinel-2 imagery. The methodology follows the work by [Graf et al. (2023)](https://doi.org/10.1016/j.rse.2023.113860)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Install required software\n",
    "\n",
    "First, we need to install [glai_processor software](https://github.com/terensis/GLAIProcessor_Optical_GPL3) that handles the actual GLAI retrieval from Sentinel-2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[33mDEPRECATION: Loading egg at /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages/crop_phenology-0.1-py3.11.egg is deprecated. pip 23.3 will enforce this behaviour change. A possible replacement is to use pip for package installation..\u001b[0m\u001b[33m\n",
      "\u001b[0mCollecting git+https://github.com/terensis/GLAIProcessor_Optical_GPL3\n",
      "  Cloning https://github.com/terensis/GLAIProcessor_Optical_GPL3 to /tmp/pip-req-build-esgeu3s_\n",
      "  Running command git clone --filter=blob:none --quiet https://github.com/terensis/GLAIProcessor_Optical_GPL3 /tmp/pip-req-build-esgeu3s_\n",
      "  Resolved https://github.com/terensis/GLAIProcessor_Optical_GPL3 to commit 34d4b0ef3d331422091332ff1c2a0c982676783c\n",
      "  Installing build dependencies ... \u001b[?25ldone\n",
      "\u001b[?25h  Getting requirements to build wheel ... \u001b[?25ldone\n",
      "\u001b[?25h  Preparing metadata (pyproject.toml) ... \u001b[?25ldone\n",
      "\u001b[?25hCollecting eodal@ git+https://github.com/EOA-team/eodal (from glai-processor==0.1)\n",
      "  Cloning https://github.com/EOA-team/eodal to /tmp/pip-install-982qcv5b/eodal_5d944b2a5f044b7ba59ef77ef5b67593\n",
      "  Running command git clone --filter=blob:none --quiet https://github.com/EOA-team/eodal /tmp/pip-install-982qcv5b/eodal_5d944b2a5f044b7ba59ef77ef5b67593\n",
      "  Resolved https://github.com/EOA-team/eodal to commit 3796df3cd64898d7516f1392151dcce70c96d0e9\n",
      "  Installing build dependencies ... \u001b[?25ldone\n",
      "\u001b[?25h  Getting requirements to build wheel ... \u001b[?25ldone\n",
      "\u001b[?25h  Preparing metadata (pyproject.toml) ... \u001b[?25ldone\n",
      "\u001b[?25hCollecting rtm-inv@ git+https://github.com/EOA-team/rtm_inv (from glai-processor==0.1)\n",
      "  Cloning https://github.com/EOA-team/rtm_inv to /tmp/pip-install-982qcv5b/rtm-inv_a9bbc8f6be9d4ea583b62d590eea588c\n",
      "  Running command git clone --filter=blob:none --quiet https://github.com/EOA-team/rtm_inv /tmp/pip-install-982qcv5b/rtm-inv_a9bbc8f6be9d4ea583b62d590eea588c\n",
      "  Resolved https://github.com/EOA-team/rtm_inv to commit e03d42619a3314d79d0ecdb22d3a69f072bee0fa\n",
      "  Installing build dependencies ... \u001b[?25ldone\n",
      "\u001b[?25h  Getting requirements to build wheel ... \u001b[?25ldone\n",
      "\u001b[?25h  Installing backend dependencies ... \u001b[?25ldone\n",
      "\u001b[?25h  Preparing metadata (pyproject.toml) ... \u001b[?25ldone\n",
      "\u001b[?25hRequirement already satisfied: rio-cogeo in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from glai-processor==0.1) (5.0.0)\n",
      "Requirement already satisfied: pandas in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.0.2)\n",
      "Requirement already satisfied: geopandas in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.13.2)\n",
      "Requirement already satisfied: sqlalchemy in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.4.50)\n",
      "Requirement already satisfied: geoalchemy2 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.13.3)\n",
      "Requirement already satisfied: numpy in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.24.4)\n",
      "Requirement already satisfied: pydantic-settings in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.0.3)\n",
      "Requirement already satisfied: pyproj in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (3.6.0)\n",
      "Requirement already satisfied: rasterio in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.3.7)\n",
      "Requirement already satisfied: matplotlib in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (3.7.1)\n",
      "Requirement already satisfied: lxml in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (4.9.2)\n",
      "Requirement already satisfied: joblib in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.2.0)\n",
      "Requirement already satisfied: psycopg2-binary in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.9.6)\n",
      "Requirement already satisfied: opencv-python-headless in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (4.7.0.72)\n",
      "Requirement already satisfied: requests in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.31.0)\n",
      "Requirement already satisfied: xarray in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2023.6.0)\n",
      "Requirement already satisfied: rtree in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.0.1)\n",
      "Requirement already satisfied: zarr in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.15.0)\n",
      "Requirement already satisfied: numba in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.57.1)\n",
      "Requirement already satisfied: pystac-client in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.7.2)\n",
      "Requirement already satisfied: sentinelsat in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.2.1)\n",
      "Requirement already satisfied: python-dotenv in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.0.0)\n",
      "Requirement already satisfied: planetary-computer>=1.0.0 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.0.0)\n",
      "Requirement already satisfied: alive-progress in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (3.1.4)\n",
      "Requirement already satisfied: rasterstats in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.19.0)\n",
      "Requirement already satisfied: pyyaml in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (6.0)\n",
      "Requirement already satisfied: myst-parser in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.0.0)\n",
      "Requirement already satisfied: pyotp in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.9.0)\n",
      "Requirement already satisfied: click>=7.0 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rio-cogeo->glai-processor==0.1) (8.1.3)\n",
      "Requirement already satisfied: morecantile<6.0,>=5.0 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rio-cogeo->glai-processor==0.1) (5.0.0)\n",
      "Requirement already satisfied: pydantic~=2.0 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rio-cogeo->glai-processor==0.1) (2.4.2)\n",
      "Requirement already satisfied: prosail in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rtm-inv@ git+https://github.com/EOA-team/rtm_inv->glai-processor==0.1) (2.0.5)\n",
      "Requirement already satisfied: spectral in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rtm-inv@ git+https://github.com/EOA-team/rtm_inv->glai-processor==0.1) (0.23.1)\n",
      "Requirement already satisfied: spart-python in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rtm-inv@ git+https://github.com/EOA-team/rtm_inv->glai-processor==0.1) (1.0.4)\n",
      "Requirement already satisfied: openpyxl in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rtm-inv@ git+https://github.com/EOA-team/rtm_inv->glai-processor==0.1) (3.1.2)\n",
      "Requirement already satisfied: lhsmdu in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rtm-inv@ git+https://github.com/EOA-team/rtm_inv->glai-processor==0.1) (1.1)\n",
      "Requirement already satisfied: attrs in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from morecantile<6.0,>=5.0->rio-cogeo->glai-processor==0.1) (23.1.0)\n",
      "Requirement already satisfied: pystac>=1.0.0 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from planetary-computer>=1.0.0->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.7.3)\n",
      "Requirement already satisfied: pytz>=2020.5 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from planetary-computer>=1.0.0->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2023.3)\n",
      "Requirement already satisfied: packaging in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from planetary-computer>=1.0.0->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (23.1)\n",
      "Requirement already satisfied: annotated-types>=0.4.0 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from pydantic~=2.0->rio-cogeo->glai-processor==0.1) (0.6.0)\n",
      "Requirement already satisfied: pydantic-core==2.10.1 in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from pydantic~=2.0->rio-cogeo->glai-processor==0.1) (2.10.1)\n",
      "Requirement already satisfied: typing-extensions>=4.6.1 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from pydantic~=2.0->rio-cogeo->glai-processor==0.1) (4.6.3)\n",
      "Requirement already satisfied: certifi in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from pyproj->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2023.5.7)\n",
      "Requirement already satisfied: python-dateutil>=2.8.2 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from pystac-client->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.8.2)\n",
      "Requirement already satisfied: affine in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rasterio->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.4.0)\n",
      "Requirement already satisfied: cligj>=0.5 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rasterio->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.7.2)\n",
      "Requirement already satisfied: snuggs>=1.4.1 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rasterio->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.4.7)\n",
      "Requirement already satisfied: click-plugins in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rasterio->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.1.1)\n",
      "Requirement already satisfied: setuptools in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from rasterio->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (65.5.1)\n",
      "Requirement already satisfied: charset-normalizer<4,>=2 in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from requests->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (3.1.0)\n",
      "Requirement already satisfied: idna<4,>=2.5 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from requests->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (3.4)\n",
      "Requirement already satisfied: urllib3<3,>=1.21.1 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from requests->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.0.3)\n",
      "Requirement already satisfied: about-time==4.2.1 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from alive-progress->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (4.2.1)\n",
      "Requirement already satisfied: grapheme==0.6.0 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from alive-progress->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.6.0)\n",
      "Requirement already satisfied: greenlet!=0.4.17 in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from sqlalchemy->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.0.2)\n",
      "Requirement already satisfied: fiona>=1.8.19 in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from geopandas->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.9.4.post1)\n",
      "Requirement already satisfied: shapely>=1.7.1 in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from geopandas->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.0.1)\n",
      "Requirement already satisfied: tzdata>=2022.1 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from pandas->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2023.3)\n",
      "Requirement already satisfied: scipy in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from lhsmdu->rtm-inv@ git+https://github.com/EOA-team/rtm_inv->glai-processor==0.1) (1.11.0)\n",
      "Requirement already satisfied: contourpy>=1.0.1 in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from matplotlib->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.1.0)\n",
      "Requirement already satisfied: cycler>=0.10 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from matplotlib->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.11.0)\n",
      "Requirement already satisfied: fonttools>=4.22.0 in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from matplotlib->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (4.40.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from matplotlib->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.4.4)\n",
      "Requirement already satisfied: pillow>=6.2.0 in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from matplotlib->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (9.5.0)\n",
      "Requirement already satisfied: pyparsing>=2.3.1 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from matplotlib->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (3.1.0)\n",
      "Requirement already satisfied: docutils<0.21,>=0.16 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.20.1)\n",
      "Requirement already satisfied: jinja2 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (3.1.2)\n",
      "Requirement already satisfied: markdown-it-py~=3.0 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (3.0.0)\n",
      "Requirement already satisfied: mdit-py-plugins~=0.4 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.4.0)\n",
      "Requirement already satisfied: sphinx<8,>=6 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (7.2.6)\n",
      "Requirement already satisfied: llvmlite<0.41,>=0.40.0dev0 in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from numba->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.40.1)\n",
      "Requirement already satisfied: et-xmlfile in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from openpyxl->rtm-inv@ git+https://github.com/EOA-team/rtm_inv->glai-processor==0.1) (1.1.0)\n",
      "Requirement already satisfied: pytest in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from prosail->rtm-inv@ git+https://github.com/EOA-team/rtm_inv->glai-processor==0.1) (7.4.0)\n",
      "Requirement already satisfied: simplejson in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from rasterstats->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (3.19.1)\n",
      "Requirement already satisfied: html2text in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sentinelsat->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2020.1.16)\n",
      "Requirement already satisfied: geojson>=2 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sentinelsat->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (3.0.1)\n",
      "Requirement already satisfied: tqdm>=4.58 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sentinelsat->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (4.65.0)\n",
      "Requirement already satisfied: geomet in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sentinelsat->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.0.0)\n",
      "Requirement already satisfied: asciitree in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from zarr->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.3.3)\n",
      "Requirement already satisfied: fasteners in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from zarr->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.18)\n",
      "Requirement already satisfied: numcodecs>=0.10.0 in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from zarr->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.11.0)\n",
      "Requirement already satisfied: six in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from fiona>=1.8.19->geopandas->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.16.0)\n",
      "Requirement already satisfied: mdurl~=0.1 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from markdown-it-py~=3.0->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.1.2)\n",
      "Requirement already satisfied: entrypoints in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from numcodecs>=0.10.0->zarr->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.4)\n",
      "Requirement already satisfied: jsonschema>=4.0.1 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from pystac>=1.0.0->planetary-computer>=1.0.0->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (4.17.3)\n",
      "Requirement already satisfied: sphinxcontrib-applehelp in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sphinx<8,>=6->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.0.7)\n",
      "Requirement already satisfied: sphinxcontrib-devhelp in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sphinx<8,>=6->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.0.5)\n",
      "Requirement already satisfied: sphinxcontrib-jsmath in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sphinx<8,>=6->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.0.1)\n",
      "Requirement already satisfied: sphinxcontrib-htmlhelp>=2.0.0 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sphinx<8,>=6->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.0.4)\n",
      "Requirement already satisfied: sphinxcontrib-serializinghtml>=1.1.9 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sphinx<8,>=6->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.1.9)\n",
      "Requirement already satisfied: sphinxcontrib-qthelp in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sphinx<8,>=6->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.0.6)\n",
      "Requirement already satisfied: Pygments>=2.14 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sphinx<8,>=6->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.15.1)\n",
      "Requirement already satisfied: snowballstemmer>=2.0 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sphinx<8,>=6->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.2.0)\n",
      "Requirement already satisfied: babel>=2.9 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sphinx<8,>=6->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.13.1)\n",
      "Requirement already satisfied: alabaster<0.8,>=0.7 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sphinx<8,>=6->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.7.13)\n",
      "Requirement already satisfied: imagesize>=1.3 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from sphinx<8,>=6->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (1.4.1)\n",
      "Requirement already satisfied: MarkupSafe>=2.0 in /mnt/ides/Lukas/venvs/GeoPython/lib64/python3.11/site-packages (from jinja2->myst-parser->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (2.1.3)\n",
      "Requirement already satisfied: iniconfig in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from pytest->prosail->rtm-inv@ git+https://github.com/EOA-team/rtm_inv->glai-processor==0.1) (2.0.0)\n",
      "Requirement already satisfied: pluggy<2.0,>=0.12 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from pytest->prosail->rtm-inv@ git+https://github.com/EOA-team/rtm_inv->glai-processor==0.1) (1.2.0)\n",
      "Requirement already satisfied: pyrsistent!=0.17.0,!=0.17.1,!=0.17.2,>=0.14.0 in /mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages (from jsonschema>=4.0.1->pystac>=1.0.0->planetary-computer>=1.0.0->eodal@ git+https://github.com/EOA-team/eodal->glai-processor==0.1) (0.19.3)\n",
      "\n",
      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.2.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m23.3.1\u001b[0m\n",
      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "!pip install git+https://github.com/terensis/GLAIProcessor_Optical_GPL3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How it works\n",
    "\n",
    "Three steps are required to get to the GLAI\n",
    "\n",
    "1. We have to fetch the Sentinel-2 surface reflectance data. We will use the 10 m bands B02 (blue), B03 (green), B04 (red), and B08 (near-infraread 1). Furthermore, we need to retrieve the angular information about illumination and viewing geometries as this required input by PROSAIL. Fortunately, this is already built-in into the `glai_processor`\n",
    "\n",
    "2. Once we have the image data, we run the RTM in so-called \"forward mode\" to simulate synthetic Sentinel-2 spectra using a range of leaf and canopy traits (among which we also have the GLAI). A full list of all traits required and the statistical distributions used, can be found in [Graf et al. (2023)](https://doi.org/10.1016/j.rse.2023.113860). The simulated spectra are stored in lookup-tables (essentially `pandas.DataFrame` objects). There is one lookup-table for each Sentinel-2 scene as the angles are slightly different for each scene.\n",
    "\n",
    "3. We use the real Sentinel-2 images and the simulated spectra to do the actual inversion. For the inversion, there is no analytical solution. Therefore, we have to compare \"real\" and simulated Sentinel-2 spectra stored in the lookup-tables per scene by means of a cost function (in our case simply the root-mean-squared error) to obtain estimates of the GLAI per pixel.\n",
    "\n",
    "The result is a GLAI map for each Sentinel-2 scene in 10 m spatial resolution.\n",
    "\n",
    "We define two functions first: The first function for fetching the Sentinel-2 data using `eodal`, and the second function for doing the actual inversion.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EOdal version:  0.2.3\n"
     ]
    }
   ],
   "source": [
    "# necessary imports\n",
    "import eodal\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "\n",
    "from datetime import datetime\n",
    "from eodal.config import get_settings\n",
    "from eodal.core.raster import RasterCollection, SceneProperties\n",
    "from eodal.core.scene import SceneCollection\n",
    "from eodal.core.sensors import Sentinel2\n",
    "from eodal.mapper.feature import Feature\n",
    "from eodal.mapper.filter import Filter\n",
    "from eodal.mapper.mapper import MapperConfigs\n",
    "from glai_processor.fetch_satellite_data import (\n",
    "    fetch_data,\n",
    "    preprocess_sentinel2_scene)\n",
    "from glai_processor.inversion import invert\n",
    "from glai_processor.utils import (\n",
    "    load_angles,\n",
    "    FPATH_SRF,\n",
    "    PLATFORMS,\n",
    "    SAMPLING_METHOD)\n",
    "from pathlib import Path\n",
    "from rtm_inv.core.lookup_table import generate_lut\n",
    "\n",
    "settings = get_settings()\n",
    "settings.USE_STAC = True\n",
    "\n",
    "print('EOdal version: ', eodal.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Fetch Sentinel-2 data\n",
    "\n",
    "We focus on a single winter wheat parcel (2022 growing seaon) in Western Switzerland. The parcel is the same as used [here](sentinel2_crop_growing_conditions.ipynb). We will retrieve data for May 2022, only, to keep the processing requirements low. However, the methodology can be extended to the whole growing season and other geographic regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_s2_data(\n",
    "    start_date: datetime,\n",
    "    end_date: datetime,\n",
    "    fpath_feature: Path,\n",
    "    output_dir: Path,\n",
    "    band_selection: list[str] = ['blue', 'green', 'red', 'nir_1']\n",
    ") -> None:\n",
    "    \"\"\"\n",
    "    Fetch the Sentinel-2 data (VIS + NIR) in 10 m resolution\n",
    "    and filter out all SCL classes but 4 and 5.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    start_date : datetime\n",
    "        Defines the start date.\n",
    "    harvest_date : datetime\n",
    "        Defines the end date.\n",
    "    fpath_feature : Path\n",
    "        Path to the feature file, i.e., the field parcel, area\n",
    "        of interest for which to fetch the data.\n",
    "    output_dir : Path\n",
    "        Path to the output directory.\n",
    "    band_selection : list[str], optional\n",
    "        List of Sentinel-2 bands to be fetched, by default\n",
    "        ['blue', 'green', 'red', 'nir_1']\n",
    "    \"\"\"\n",
    "    feature = Feature.from_geoseries(\n",
    "        gpd.read_file(fpath_feature).dissolve().geometry\n",
    "    )\n",
    "    metadata_filters: list[Filter] = [\n",
    "        Filter('cloudy_pixel_percentage', '<', 50),\n",
    "        Filter('processing_level', '==', 'Level-2A')\n",
    "    ]\n",
    "\n",
    "    # define the scene kwargs\n",
    "    scene_kwargs = {\n",
    "        'scene_constructor': Sentinel2.from_safe,\n",
    "        'scene_constructor_kwargs': {'band_selection': band_selection.copy(),\n",
    "                                     'apply_scaling': False,\n",
    "                                     'read_scl': True},\n",
    "        'scene_modifier': preprocess_sentinel2_scene,\n",
    "        'scene_modifier_kwargs': {'target_resolution': 10}\n",
    "    }\n",
    "    # define the mapper configs\n",
    "    mapper_configs = MapperConfigs(\n",
    "        collection='sentinel2-msi',\n",
    "        time_start=start_date,\n",
    "        time_end=end_date,\n",
    "        feature=feature,\n",
    "        metadata_filters=metadata_filters)\n",
    "\n",
    "    # fetch the data\n",
    "    fetch_data(\n",
    "        mapper_configs=mapper_configs,\n",
    "        scene_kwargs=scene_kwargs,\n",
    "        band_selection=band_selection,\n",
    "        output_dir=output_dir\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-11-03 15:11:51,805 eodal        INFO     Starting extraction of sentinel2 scenes\n",
      "2023-11-03 15:11:54,794 eodal        INFO     Finished extraction of sentinel2 scenes\n",
      "/mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages/eodal/core/raster.py:1777: RasterioDeprecationWarning: Source dataset should be opened in read-only mode. Use of datasets opened in modes other than 'r' will be disallowed in a future version.\n",
      "  cog_translate(\n",
      "/mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages/eodal/core/raster.py:1777: RasterioDeprecationWarning: Source dataset should be opened in read-only mode. Use of datasets opened in modes other than 'r' will be disallowed in a future version.\n",
      "  cog_translate(\n",
      "/mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages/eodal/core/raster.py:1777: RasterioDeprecationWarning: Source dataset should be opened in read-only mode. Use of datasets opened in modes other than 'r' will be disallowed in a future version.\n",
      "  cog_translate(\n"
     ]
    }
   ],
   "source": [
    "# set the paths to the data\n",
    "base_dir = Path(os.path.dirname(os.path.realpath(\"__file__\"))).parent.parent\n",
    "\n",
    "# path to the field parcel for which to fetch the data\n",
    "fpath_feature = base_dir.joinpath('data/sample_polygons/Parzelle35.shp')\n",
    "\n",
    "# path to the output directory. All fetched data will be stored here\n",
    "output_dir = base_dir.joinpath('results')\n",
    "output_dir.mkdir(exist_ok=True)\n",
    "\n",
    "# start and end date for the data fetch\n",
    "start_date = datetime(2022, 5, 1)\n",
    "end_date = datetime(2022, 5, 30)\n",
    "\n",
    "extract_s2_data(\n",
    "    start_date=start_date,\n",
    "    end_date=end_date,\n",
    "    fpath_feature=fpath_feature,\n",
    "    output_dir=output_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2 and 3 - Generate the lookup tables and perform the inversion\n",
    "\n",
    "This step can take a while as the RTM has to be run multiple times in forward mode (20000 times in the case per scene) to generate the synthetic spectra. The actual inversion is also performed. The result is a \"*_traits.tiff\" file in the [result folder](/results/) that has been generated in step 1 and filled with GeoTiff files containing the Sentinel-2 data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def retrieve_traits(\n",
    "    rtm_params: Path,\n",
    "    output_dir: Path,\n",
    "    traits: list[str] = ['lai', 'cab'],\n",
    "    band_selection_lut: list[str] = ['B02', 'B03', 'B04', 'B08'],\n",
    "    band_selection_srf: list[str] = ['blue', 'green', 'red', 'nir_1'],\n",
    "    n_solutions: int = 2000,\n",
    "    lut_size: int = 20000\n",
    ") -> None:\n",
    "    \"\"\"\n",
    "    Retrieve the traits from the Sentinel-2 data using the\n",
    "    lookup-table based inversion of the PROSAIL RTM.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    rtm_params : Path\n",
    "        Path to the RTM parameters\n",
    "    output_dir : Path\n",
    "        Path to the output directory\n",
    "    traits : list[str], optional\n",
    "        List of traits to be retrieved, by default ['lai', 'cab']\n",
    "    band_selection_lut : list[str], optional\n",
    "        List of Sentinel-2 bands to be used for the lookup-tables,\n",
    "        by default ['B02', 'B03', 'B04', 'B08']\n",
    "    band_selection_srf : list[str], optional\n",
    "        List of Sentinel-2 bands to be used for the spectral response\n",
    "        functions, by default ['blue', 'red', 'green', 'nir_1']\n",
    "    n_solutions : int, optional\n",
    "        Number of solutions to be retrieved, by default 2000.\n",
    "    lut_size : int, optional\n",
    "        Size of the lookup-tables, by default 20000\n",
    "    \"\"\"\n",
    "    # iterate over all angle files\n",
    "    for yaml_file in output_dir.glob('S2*_angles.yaml'):\n",
    "        # load the angles\n",
    "        angles = load_angles(yaml_file)\n",
    "\n",
    "        fname_lut = yaml_file.name.replace(\n",
    "            'angles.yaml', 'lut.pkl'\n",
    "        )\n",
    "\n",
    "        # GeoTiff file containing the spectral values\n",
    "        band_names_file = '-'.join(band_selection_srf) + '.tiff'\n",
    "        fname_srf = yaml_file.name.replace(\n",
    "            'angles.yaml', band_names_file)\n",
    "        fpath_srf = output_dir.joinpath(fname_srf)\n",
    "        # get the platform from the file name\n",
    "        platform = PLATFORMS[fpath_srf.name.split('_')[0]]\n",
    "\n",
    "        # generate the lookup-tables using the methodology from\n",
    "        # Graf et al. (2023, RSE, https://doi.org/10.1016/j.rse.2023.113860)\n",
    "        if not output_dir.joinpath(fname_lut).exists():\n",
    "            lut_srf = generate_lut(\n",
    "                sensor=platform,\n",
    "                lut_params=rtm_params,\n",
    "                solar_zenith_angle=angles['solar_zenith_angle'],\n",
    "                viewing_zenith_angle=angles['viewing_zenith_angle'],\n",
    "                solar_azimuth_angle=angles['solar_azimuth_angle'],\n",
    "                viewing_azimuth_angle=angles['viewing_azimuth_angle'],\n",
    "                lut_size=lut_size,\n",
    "                sampling_method=SAMPLING_METHOD,\n",
    "                fpath_srf=FPATH_SRF,\n",
    "                remove_invalid_green_peaks=True,\n",
    "                linearize_lai=False\n",
    "            )\n",
    "        else:\n",
    "            lut_srf = pd.read_pickle(output_dir.joinpath(fname_lut))\n",
    "\n",
    "        # save the lookup-tables\n",
    "        lut_srf.dropna(inplace=True)\n",
    "        fpath_lut = output_dir.joinpath(fname_lut)\n",
    "        lut_srf.to_pickle(fpath_lut)\n",
    "\n",
    "        # invert the RTM\n",
    "        invert(\n",
    "            fpath_lut=fpath_lut,\n",
    "            fpath_srf=fpath_srf,\n",
    "            output_dir=output_dir,\n",
    "            band_selection_lut=band_selection_lut,\n",
    "            band_selection_srf=band_selection_srf,\n",
    "            traits=traits,\n",
    "            n_solutions=n_solutions\n",
    "        )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the generation of the lookup-table and the subsequent inversion, we stick to defaults. However, this might not be the optimum settings and are mainly used for pure demonstration purposes.\n",
    "\n",
    "In the below you see the leaf and canopy parameters and how they are varied using Gaussian or Uniform distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Parameter</th>\n",
       "      <th>Min</th>\n",
       "      <th>Max</th>\n",
       "      <th>Mode</th>\n",
       "      <th>Std</th>\n",
       "      <th>Distribution</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>n</td>\n",
       "      <td>1.00</td>\n",
       "      <td>2.50</td>\n",
       "      <td>1.50</td>\n",
       "      <td>0.20</td>\n",
       "      <td>Gaussian</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>cab</td>\n",
       "      <td>0.00</td>\n",
       "      <td>80.00</td>\n",
       "      <td>50.00</td>\n",
       "      <td>40.00</td>\n",
       "      <td>Gaussian</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>car</td>\n",
       "      <td>0.00</td>\n",
       "      <td>15.00</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Gaussian</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>ant</td>\n",
       "      <td>0.00</td>\n",
       "      <td>5.00</td>\n",
       "      <td>2.00</td>\n",
       "      <td>0.80</td>\n",
       "      <td>Gaussian</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>cbrown</td>\n",
       "      <td>0.00</td>\n",
       "      <td>1.00</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Uniform</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>cw</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.07</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.02</td>\n",
       "      <td>Gaussian</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>cm</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.01</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Uniform</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>lai</td>\n",
       "      <td>0.00</td>\n",
       "      <td>8.00</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Uniform</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>lidfa</td>\n",
       "      <td>30.00</td>\n",
       "      <td>70.00</td>\n",
       "      <td>45.00</td>\n",
       "      <td>15.00</td>\n",
       "      <td>Uniform</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>lidfb</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Constant</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>hspot</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.50</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Uniform</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>rsoil</td>\n",
       "      <td>0.00</td>\n",
       "      <td>1.00</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Uniform</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>psoil</td>\n",
       "      <td>0.00</td>\n",
       "      <td>1.00</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Uniform</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>tts</td>\n",
       "      <td>25.90</td>\n",
       "      <td>25.90</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Constant</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>tto</td>\n",
       "      <td>5.90</td>\n",
       "      <td>5.90</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Constant</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>psi</td>\n",
       "      <td>131.10</td>\n",
       "      <td>131.10</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Constant</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>typelidf</td>\n",
       "      <td>2.00</td>\n",
       "      <td>2.00</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Constant</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Parameter     Min     Max   Mode    Std Distribution\n",
       "0          n    1.00    2.50   1.50   0.20     Gaussian\n",
       "1        cab    0.00   80.00  50.00  40.00     Gaussian\n",
       "2        car    0.00   15.00    NaN    NaN     Gaussian\n",
       "3        ant    0.00    5.00   2.00   0.80     Gaussian\n",
       "4     cbrown    0.00    1.00    NaN    NaN      Uniform\n",
       "5         cw    0.00    0.07   0.04   0.02     Gaussian\n",
       "6         cm    0.00    0.01    NaN    NaN      Uniform\n",
       "7        lai    0.00    8.00    NaN    NaN      Uniform\n",
       "8      lidfa   30.00   70.00  45.00  15.00      Uniform\n",
       "9      lidfb    0.00    0.00    NaN    NaN     Constant\n",
       "10     hspot    0.01    0.50    NaN    NaN      Uniform\n",
       "11     rsoil    0.00    1.00    NaN    NaN      Uniform\n",
       "12     psoil    0.00    1.00    NaN    NaN      Uniform\n",
       "13       tts   25.90   25.90    NaN    NaN     Constant\n",
       "14       tto    5.90    5.90    NaN    NaN     Constant\n",
       "15       psi  131.10  131.10    NaN    NaN     Constant\n",
       "16  typelidf    2.00    2.00    NaN    NaN     Constant"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# we need to the RTM which leaf and canopy parameter combinations to use.\n",
    "# Here, we use the settings suggested by Danner et al. \n",
    "rtm_params = pd.read_csv('https://raw.githubusercontent.com/EOA-team/sentinel2_crop_traits/main/src/lut_params/prosail_danner-etal_all_phases.csv')\n",
    "rtm_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages/eodal/core/raster.py:1777: RasterioDeprecationWarning: Source dataset should be opened in read-only mode. Use of datasets opened in modes other than 'r' will be disallowed in a future version.\n",
      "  cog_translate(\n",
      "/mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages/openpyxl/worksheet/_reader.py:329: UserWarning: Unknown extension is not supported and will be removed\n",
      "  warn(msg)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulated spectrum 500/20000\n",
      "Simulated spectrum 1000/20000\n",
      "Simulated spectrum 1500/20000\n",
      "Simulated spectrum 2000/20000\n",
      "Simulated spectrum 2500/20000\n",
      "Simulated spectrum 3000/20000\n",
      "Simulated spectrum 3500/20000\n",
      "Simulated spectrum 4000/20000\n",
      "Simulated spectrum 4500/20000\n",
      "Simulated spectrum 5000/20000\n",
      "Simulated spectrum 5500/20000\n",
      "Simulated spectrum 6000/20000\n",
      "Simulated spectrum 6500/20000\n",
      "Simulated spectrum 7000/20000\n",
      "Simulated spectrum 7500/20000\n",
      "Simulated spectrum 8000/20000\n",
      "Simulated spectrum 8500/20000\n",
      "Simulated spectrum 9000/20000\n",
      "Simulated spectrum 9500/20000\n",
      "Simulated spectrum 10000/20000\n",
      "Simulated spectrum 10500/20000\n",
      "Simulated spectrum 11000/20000\n",
      "Simulated spectrum 11500/20000\n",
      "Simulated spectrum 12000/20000\n",
      "Simulated spectrum 12500/20000\n",
      "Simulated spectrum 13000/20000\n",
      "Simulated spectrum 13500/20000\n",
      "Simulated spectrum 14000/20000\n",
      "Simulated spectrum 14500/20000\n",
      "Simulated spectrum 15000/20000\n",
      "Simulated spectrum 15500/20000\n",
      "Simulated spectrum 16000/20000\n",
      "Simulated spectrum 16500/20000\n",
      "Simulated spectrum 17000/20000\n",
      "Simulated spectrum 17500/20000\n",
      "Simulated spectrum 18000/20000\n",
      "Simulated spectrum 18500/20000\n",
      "Simulated spectrum 19000/20000\n",
      "Simulated spectrum 19500/20000\n",
      "Simulated spectrum 20000/20000\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages/eodal/core/raster.py:1777: RasterioDeprecationWarning: Source dataset should be opened in read-only mode. Use of datasets opened in modes other than 'r' will be disallowed in a future version.\n",
      "  cog_translate(\n",
      "/mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages/openpyxl/worksheet/_reader.py:329: UserWarning: Unknown extension is not supported and will be removed\n",
      "  warn(msg)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulated spectrum 500/20000\n",
      "Simulated spectrum 1000/20000\n",
      "Simulated spectrum 1500/20000\n",
      "Simulated spectrum 2000/20000\n",
      "Simulated spectrum 2500/20000\n",
      "Simulated spectrum 3000/20000\n",
      "Simulated spectrum 3500/20000\n",
      "Simulated spectrum 4000/20000\n",
      "Simulated spectrum 4500/20000\n",
      "Simulated spectrum 5000/20000\n",
      "Simulated spectrum 5500/20000\n",
      "Simulated spectrum 6000/20000\n",
      "Simulated spectrum 6500/20000\n",
      "Simulated spectrum 7000/20000\n",
      "Simulated spectrum 7500/20000\n",
      "Simulated spectrum 8000/20000\n",
      "Simulated spectrum 8500/20000\n",
      "Simulated spectrum 9000/20000\n",
      "Simulated spectrum 9500/20000\n",
      "Simulated spectrum 10000/20000\n",
      "Simulated spectrum 10500/20000\n",
      "Simulated spectrum 11000/20000\n",
      "Simulated spectrum 11500/20000\n",
      "Simulated spectrum 12000/20000\n",
      "Simulated spectrum 12500/20000\n",
      "Simulated spectrum 13000/20000\n",
      "Simulated spectrum 13500/20000\n",
      "Simulated spectrum 14000/20000\n",
      "Simulated spectrum 14500/20000\n",
      "Simulated spectrum 15000/20000\n",
      "Simulated spectrum 15500/20000\n",
      "Simulated spectrum 16000/20000\n",
      "Simulated spectrum 16500/20000\n",
      "Simulated spectrum 17000/20000\n",
      "Simulated spectrum 17500/20000\n",
      "Simulated spectrum 18000/20000\n",
      "Simulated spectrum 18500/20000\n",
      "Simulated spectrum 19000/20000\n",
      "Simulated spectrum 19500/20000\n",
      "Simulated spectrum 20000/20000\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/mnt/ides/Lukas/venvs/GeoPython/lib/python3.11/site-packages/eodal/core/raster.py:1777: RasterioDeprecationWarning: Source dataset should be opened in read-only mode. Use of datasets opened in modes other than 'r' will be disallowed in a future version.\n",
      "  cog_translate(\n"
     ]
    }
   ],
   "source": [
    "retrieve_traits(\n",
    "    rtm_params=rtm_params,\n",
    "    output_dir=output_dir\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAFVCAYAAACJsxn/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACD90lEQVR4nO2dd3hU1db/v1OSTNqkF0IKECD0XgyoiHQQFMWCNBWvohSvInDxegVEQMTXjlhfwMLPgiKKSPOKdKQKSIcEQgpJSO/JzP79wZuRtadPJmRmsj7PM8+Tdeqek5xv9tnnu9dSCCEEGIZhGIZhGLdG2dANYBiGYRiGYeoOd+oYhmEYhmE8AO7UMQzDMAzDeADcqWMYhmEYhvEAuFPHMAzDMAzjAXCnjmEYhmEYxgPgTh3DMAzDMIwHwJ06hmEYhmEYD4A7dQzDMAzDMB4Ad+oYhmEYhmE8AO7UuTAHDhzAtGnT0L59e/j7+yM+Ph4PPPAAzp49a3L7U6dOYejQoQgICEBoaCgmTJiAnJwcss3p06cxe/ZsdOnSBYGBgWjSpAlGjBiBgwcPGh3v+++/x4MPPogWLVrAz88PSUlJmDlzJgoKCuz6Hp9++inatm0LjUaDVq1a4d133zXaZv78+VAoFEYfjUZj83kWLVqEUaNGISoqCgqFAvPnz7dpv0GDBkGhUGDatGk2n2vPnj249dZb4efnh+joaMyYMQMlJSVG21VWVmLOnDmIiYmBr68vevfuja1bt9p8HsC26wcA6enpeOCBBxAcHAytVou7774bFy9etPk89rTV1u/PMJZgjas/jfvqq6/QrVs3aDQaREREYPLkycjNzbX5XKxxbqpxgnFZ7rvvPhEdHS2mT58uPv74Y7Fw4UIRFRUl/P39xfHjx8m2aWlpIjw8XCQmJoq3335bLFq0SISEhIjOnTuLyspKw3YzZ84UwcHBYvLkyeLDDz8Ur732mkhMTBQqlUps3bqVHDMsLEx07NhR/Oc//xEff/yxmDFjhvD29hZt2rQRZWVlNn2HDz74QAAQ9913n/joo4/EhAkTBADx6quvku3mzZsnAIgVK1aIzz//3PBZs2aNzdcLgIiOjhZDhgwRAMS8efOs7vPdd98Jf39/AUBMnTrVpvMcOXJEaDQa0bVrV7FixQrx73//W/j4+IihQ4cabfvQQw8JtVotnn/+efHhhx+K5ORkoVarxc6dO206l63Xr7i4WLRq1UpERkaKpUuXijfeeEPExcWJ2NhYkZuba9O5bG2rPd+fYSzBGlc/Gvf+++8LAGLAgAFi+fLlYu7cucLPz0906tRJlJeXWz0Pa5z7ahx36lyY3bt3E7ESQoizZ88KHx8fMW7cOLL8qaeeEr6+vuLSpUuGZVu3bhUAxIcffmhYdvDgQVFcXEz2zc3NFREREaJv375k+W+//WbUptWrVwsA4uOPP7ba/rKyMhEWFiZGjBhBlo8bN074+/uLvLw8w7JawcvJybF6XHOkpKQIIYTIycmxqVNXXl4umjVrJl5++WW7OnXDhg0TTZo0EYWFhYZlH3/8sQAgNm/ebFi2f/9+AUAsW7aMnDMxMVEkJydbPY8912/p0qUCgPjjjz8My06dOiVUKpWYO3eu1XPZ01Zbvz/DWIM1zj5s0bjKykoRHBwsbr/9dqHX6w3Lf/rpJwFAvPPOO1bPwxrnvhrHnTo3pFu3bqJbt25kWWRkpLj//vuNtm3durUYMGCA1WPee++9IjQ01Op2RUVFAoB47rnnrG77888/CwDi559/Jsv37NkjAIjPP//csKxW8LKzs0VhYSERI3uxtVO3YMECER8fL8rKymzu1BUWFgq1Wi1mzZpFlldWVoqAgAAxefJkw7JZs2YJlUpFhEEIIRYvXiwAiMuXL5M2nzp1SpSWlhqW2XP9evbsKXr27GnU3sGDB4vExESy7NKlS+LUqVNkma1ttef7M4yjsMZZxpLGHTp0SAAQy5cvN1oXEBAg+vTpY/HYrHHurXHsqXMzhBC4evUqwsPDDcvS09ORnZ2NHj16GG3fq1cvHDlyxOpxs7KyyDEtbQfApm1rzyu3q3v37lAqlSbb1aJFCwQFBSEwMBDjx4/H1atXrZ7HES5fvoxXX30VS5cuha+vr837HT9+HDU1NUbfydvbG126dCHf6ciRI2jdujW0Wi3ZtlevXgCAo0ePGpa99957aNu2Lf744w+yP2D9+un1ehw7dszs7//ChQsoLi42LJs4cSLatm1LtrO1rfZ8f4ZxBNa4ulFZWQkAJnXN19cXR44cgV6vN7s/a5x7axx36tyML7/8Eunp6XjwwQcNyzIzMwEATZo0Mdq+SZMmyMvLM9zopti5cyf27t1LjmmOpUuXQqVSYcyYMVa3zczMhEqlQmRkJFnu7e2NsLAwZGRkGJaFhIRg2rRp+PDDD7F27Vo8/vjj+Prrr3HbbbehqKjI6rnsZebMmejatSseeughu/azdq1v/E6ZmZlmtwNAtjV3LluuX+3vt67nsmV/e74/wzgCa1zdaNWqFRQKBXbv3k2WnzlzBjk5OSgvL0d+fr7F7wSwxrmrxqkbugGM7Zw+fRpTp05FcnIyJk2aZFheXl4OAPDx8THap3ZmVXl5ucn12dnZePjhh9G8eXPMnj3b4vnXrFmDTz/9FLNnz0arVq2stre8vBze3t4m12k0GkO7AeCZZ54h6++77z706tUL48aNw/vvv49//etfVs9nK7/99hu+++477N+/3+59rV3rG7+TuWt+4++klvnz5xvNZLP1+tn6+69l+/btJr+XLfvb8/0Zxl5Y4+pOeHg4HnjgAaxevRpt27bF6NGjkZ6ejunTp8PLywvV1dUW71PWOPfWOB6pcxOysrIwYsQIBAUFYe3atVCpVIZ1tcPspp5UKyoqyDY3UlpairvuugvFxcVYv349AgICzJ5/586dmDx5MoYMGYJFixaRdTk5OcjKyjJ8aqd9+/r6oqqqyuTxKioqrL72fPjhhxEdHY1t27YZlt14nqysLLtvsJqaGsyYMQMTJkxAz5497doXsH6tb/xOvr6+dv9O5HPZcv0c/f3L57Jlf3u+P8PYA2vcdeqqcQDw4YcfYvjw4Xj++eeRmJiI22+/HR07dsTIkSMBwOJ1YI1zb43jTp0bUFhYiGHDhqGgoACbNm1CTEwMWV87TFw7bHwjmZmZCA0NNXrqqKqqwr333otjx45h/fr16NChg9nz//nnnxg1ahQ6dOiAtWvXQq2mA7w9e/ZEkyZNDJ/XX3/d0C6dTofs7Gyjc1+7ds3oe5giLi4OeXl55Lve+Pn666+tHuNGPvvsM5w5cwZPPvkkUlNTDR8AKC4uRmpqKsrKyszub+1a3/idmjRpYnY7AFa/v63Xr/b3W9dz2bK/Pd+fYWyFNc55GgcAQUFBWL9+PS5duoTff/8dqamp+Pzzz5GZmYmIiAgEBweb3Zc1zr01jjt1Lk5FRQVGjhyJs2fPYsOGDWjXrp3RNk2bNkVERITJ5Jp//PEHunTpQpbp9XpMnDgRv/76K9asWYN+/fqZPf+FCxcwdOhQREZGYuPGjSaf8L788kts3brV8Jk4cSIAGM4rt+vgwYPQ6/VG7ZIRQiA1NRURERGGZTeeZ+vWrRgyZIjFY8hcvnwZ1dXV6Nu3L5o3b274ANc7fM2bN8eWLVvM7t+hQweo1Wqj71RVVYWjR4+S79SlSxecPXvWyC9T+9rX2ve39foplUp07NjR5O9///79aNGiBQIDA62ey5a22vP9GcYWWOOcq3E3Eh8fj9tvvx0JCQkoKCjAoUOHMHDgQIv7sMa5ucY17ORbxhI1NTVi1KhRQq1WG035lpkyZYrw9fUlU8i3bdtmSHZ5I08//bRRbidTZGZmihYtWoiYmBhDfiR7KCsrE6GhoeKuu+4iy8ePHy/8/PzEtWvXDMuys7ON9l++fLkAIN544w27zmtpuv+pU6fEunXrjD4AxPDhw8W6detERkaGxeMPHTpUNGnSRBQVFRmWffLJJwKA+OWXXwzL9u3bZ5QXqaKiQrRs2VL07t3bqM3ydH97rt+rr74qAIgDBw4Ylp0+fVqoVCoxZ84csr+p6f72tNXW788w1mCNc77GmWPKlClCqVSSPG/mYI1zX43jTp0L88wzzwgAYuTIkSQDee3nRi5fvizCwsJEYmKieOedd8TixYtFSEiI6Nixo6ioqDBs9+abbwoAIjk52eQxS0pKDNt27txZABCzZ8822m7Lli02fYda0RozZoz4+OOPxcSJEwUAsWjRIrKdr6+veOSRR8T//M//iOXLl4uxY8cKhUIhunTpQkTAEp999plYuHChmDt3rgAg+vfvLxYuXCgWLlwoUlNTLe4LO5IPHzp0SPj4+JBs4xqNRgwePNho2/vvv9+Q8+jDDz8Uffr0EWq1Wvz+++9ku9ocVnIyVFuvX1FRkUhMTBSRkZHitddeE2+++aaIi4sTMTExRv9M+vXrJ0w9z9naVnu+P8NYgjWufjRuyZIlYty4ceKdd94R77//vhg8eLAAIF555RWbzsMa574ax506F6b2D9PcR+bEiRNi8ODBws/PTwQHB4tx48aJrKwsss2kSZMsHvPGp1VL2/Xr18/m7/HRRx+JpKQk4e3tLRITE8Wbb75plHjz8ccfF+3atROBgYHCy8tLtGzZUsyZM4c8KdXlepnKHH8j9nTqhBBi586dok+fPkKj0YiIiAgxdepUk20tLy8Xzz//vIiOjhY+Pj6iZ8+eYtOmTUbbmRM8IWy7fkJcL6M0ZswYodVqRUBAgLjrrrvEuXPnjLYzJ3i2ttWe788wlmCNqx+N27Bhg+jVq5cIDAwUfn5+4pZbbhHffPONzecRgjXOXTVOIYQQll/QMgzDMAzDMK4OT5RgGIZhGIbxALhTxzAMwzAM4wFwp45hGIZhGMYD4E4dwzAMwzCMB8CdOoZhGIZhGA+AO3UMwzAMwzAegNr6JkxDotfrkZGRgcDAQCgUioZuDsOYRAiB4uJixMTEQKn8+1mxoqLCbMHuWry9vaHRaOq7iYwLwvrGuAuOatzN1jfu1Lk4GRkZiIuLa+hmMIxNpKWlITY2FsB1sWueEICsbJ3FfaKjo5GSksIdu0YI6xvjbtircTdb37hT5+LUFilOS0uDVqtt4NY0bu4OmkhihZc3iZXxMSSujqK/r4KW9Kbu+difJH6j68d1bWKDUVRUhLi4OFJUu6qqClnZOqQcSoA20LTTo6hYj+bdL6Gqqoo7dY0Q1jfX4Z6wx0is8KH6hsR4ElYH+ZC4LIpu3/bpv0i8vPsHdWxhw+KIxjWEvnGnzsWpfSWh1WpZ9BoYtcKLxAopVqqoyAk1vYlV3jT2DqD7e8Lv19QrNP+A6x9T6LieTaOG9c11MNY3qVNnRd/U0kOudwCNPeX3a4/GOaJv6enpmDNnDn755ReUlZWhZcuWWLlyJXr06GHT/typYxgbUUVFkri8awKJL91Fn9SU5TSO3V5D4v/+0J2eQAo9BT0E9DCtbuaWMwxzc1E2pyNxuX2iSNzyidMk3ne0NYnjN9JXkPu/6UxP0LOODXRhzGmcvfqWn5+Pvn37on///vjll18QERGBc+fOISQkxOZjcKeOYZh6RQ899BbWMQzDuDPmNM5efVu6dCni4uKwcuVKw7LmzZvbdQxOacIwTL2iE8Lih2EYxp2xpm9FRUXkU1lZafI4P/74I3r06IH7778fkZGR6Nq1Kz7+2D6vNXfqGIapV2pfTZj7MAzDuDPW9C0uLg5BQUGGz5IlS0we5+LFi1ixYgVatWqFzZs346mnnsKMGTOwevVqm9vCr1+ZRsEg5f1Wt1FKs5MULajHJH1sSxI/+NivJO7ml0ri/5y+m8SlJ8NIHJjWODo0egjo2FPHMPXGYO+xdIHCeLxGGRRI4uq2VN/OP0/16ZOBH5E4QlVK4rdV1SQ+cqIjiX2zG8+9bU7javVNnt3t4+NjtC1wPW9jjx49sHjxYgBA165dceLECXzwwQeYNGmSTW3hkTqGYeoVHqljGMaTsaZvtbO7az/mOnVNmjRBu3btyLK2bdvi8uXLNreFR+oYhqlXLHnn2FPHMIy7Y07j7NW3vn374syZM2TZ2bNnkZCQYGYPY3ikjmGYekVv5cMwDOPOOEvfnn32Wezbtw+LFy/G+fPnsWbNGnz00UeYOnWqzcfgkTqmcaBUkVAdE220SUnXpiROG0qTTM7uv57EHTVpJPZS0DxNg5rSvE5fd+pDYr/0xvFMpbPgqTO3nGEY21H4+tIFical1zL7BJM4aRzVp/80+S+JW6hLSEwVFOgXTPff3rcViX3OSG3yYMxpnL361rNnT6xbtw5z587Fyy+/jObNm+Ott97CuHHjbD4Gd+oYhqlXdMJ8ZnWuKMEwjLtjTuMc0be77roLd911l8Nt4U4dwzD1Sg0UqIZxaZ3adQzDMO6MOY1rCH3jTh3DMPWKXlz/mFvHMAzjzpjTuIbQN+7UMR6BnIdOzjknencg8cUh/sbHuOsAiScHpJA43iuPxKHKChJHqOgdPFh7gsS/t6V57jK0oSRuteQNEtck0OOnPPyCUZvdAR0U0Jl5YjW3nGGYvxmkfojEqjCqHWV3tCVx3mPUDwcACzt8RuIYdT6JW6irSBykpJ44pXSvdpc8xV0TaHxZqlfa7N3/IXGXrhdI/MOt7xu12V0wp3ENoW/cqWMYpl7hTh3DMJ4Md+oYhmk06IUCemFa3MwtZxiGcRfMaVxD6Bt36hiGqVd4pI5hGE+GR+oYxsmoIiJIXNazGYkzJlC/yN1Je42OcYf2FImbqgpJHC7VOgxW0tvHS0EzOXXwLiZx/+hzJC4Kp76/4/ExJM7cS2N3RQcldGbynOtMLmUY5kbUzWjeuZzbm5C445TjJI73pf5fAIhT50kx1TM/hReJZT2TSVDTDsvAMKqfeUEBJD4aFEviAydb0APeavF0Lo05jWsIfeNOHcMw9Yqw8PpV8OtXhmHcHHMa1xD6xp06hmHqFX79yjCMJ8OvXxmGaTTohBI6Yeb1K+epYxjGzTGncQ2hb9ypY1wCOc+cjMLHh8TKls1IXNg/kcQ53egT0uCWx0gc6V1kdI42XrkkDlTSY8h5mnwV3uYbDGOPymAt9b0cLm9G4s7+NM/TEr97SNz9cZrHzruEKsber2ZabE9DoYcCejOeOj3XfmUaAdb0Tc47V902nsTZI6m/tupW6tftE3SexMGqMqNzRKmorzjAin5ZI0BJPcG3+dE2XKgOI3E3v1QSHzzfjMTN1ywm8YZbl5O4fXy6A628OZjTuIbQN+7UMQxTr/DrV4ZhPBl+/cowTKPB8utXHqljGMa9Mf/6lUfqGIbxMK6/mjCTfJhH6hiGcXPMaVxD6Bt36hiGqVf0FvLUsaeOYRh3x5zGsaeOYf4P2TisS2xK4ty2NLFlaVP6RFTjR9M+ZpRrSdw14JLROSsETbbpI+gx9NL2XqKSxBqF5WTEPbypURlIJVFODW2jLriGHj+fHk9VIbfINeHXr0xjR+FFJyWooiNJXN42msTZ3ej2pc2oFoRqqPZkVgeTOM7rmlEbyqScaSo91SOVQhpVkuRFXu8jTQRrLX1HL0UOiVOlNnZoTic+XNhEkxHPbDqGxJvo3BGXgl+/MgzTaNBDybNfGYbxWMxpHI/UMQzjceiEAjozmdXNLWcYhnEXzGlcQ+gbd+oYhqlXLNd+5ZE6hmHcG/O1X3mkjvFQ5OSbysBAEit6diRxQaI/iXVe9InHL5d6TKqCqL+jzIcaQnLLqQevQtDtASBd8rQVKyvoOZW0ALYStA3yPS177HykuIW6hMS7SpJIrCijHjqFXliMXRW9UEJvxlOnZ08d4wEM9h5LYlVEOInLhnchcWEzy/96Nbn0vqjSUi3QN6H306Vymui3q5ToFwDydDRZcJmC6le4iuqbXiGZ6uR7WBqEkj3EEVLy9sM6qsHN/Knv73RQcxJf2tqMnqAfXBZzGtcQ+sadOoZh6hUeqWMYxpPhkTqGYRoNepj3lrjH/F2GYRjzmNO4htA37tQxDFOvWJ79ano5wzCMu2B+9uvN1zfu1DE3BVXbViS+1oN6TvLb0+1rfOmwdcAl+hSk96J/unrpL1lZLuWcU1P/iCn/w7mqaKNlN9LKO4vEMT4FJPaTCmTrpee0GtC8d8FK2uhqKU+eQke/c1Ugjb3dZOKo5Tx13Klj3J/q2zuT+EovHxLruheT2NuL6lHxpSAS++RL94WS6mFxiS+JK/VUOwp01JNsaplK0qfOPjRvXBOVlFtPYd+96qf0kuJKM1tep2Vvmjs09b/N7DpfQ2I+Tx136hiG8TC4TBjDMJ4MlwljGKbRwCN1DMN4Mq40UseKyjBMvVI7M8zcxx7mz58PhUJBPm3atKmnljMMw1jHWfrmDHikjrEJOc+cjDqW1mYtb9eExFdHR5C4qkMZiX39qN8i2JfmiEsPDSaxOE49JQpqV4NQSx6USupxOVdOay8CQKVkzKuSYv9g2sY23vkkDpfyNJXpaaMqBfXRKKVais19skksP+TVaKTajVV0g+SH/ofEXiXSRQGwY8Nso2X1jV4ooDc3+9WBjOvt27fHtm3bDLFazTLG1A05z5zCm/rJFM3jSJzbg9amzp1IjxcTTf23bUOukjjKp4jEv/q0JnHOX1QvldX0PpHTn2WX07yfp8pjICPrW5C6nMTBKqrJgUqaR66Jyo/EsseuWqqVrQbVw47eVN/yAmjeOrl9p+Lp/5Q2/3mTxB2GnYHM2j4fGC27GZjTOEf0ra6wGjIMU6/oLTyxOjI7TK1WIzra8qQWhmGYm4U5jWuI2a/8+pVhmHqlNtu6uQ8AFBUVkU9lpfmZcufOnUNMTAxatGiBcePG4fLlyzfrqzAMwxhhTd9uJtypYximXtFBYfEDAHFxcQgKCjJ8lixZYvJYvXv3xqpVq7Bp0yasWLECKSkpuO2221BcXGxye4ZhmPrGmr7dTPj1K2MTcq1WJFKPSVYvmmepuIV8BGoCqSmT6qJqqZ+jVVAOiWMDCki8T0dPoM6SPDA19GbKzqV1XX+rpHnzAKBtBPW9tPTPMdrmRo5X0Vx7GgX1oESqaF4oH0G/8zU99bRopNqyUNFrJqRailJqKsilGnUa13hms1z79frytLQ0aLV//458fHxMbj9s2DDDz506dULv3r2RkJCAb775BpMnT3Ziq5nGhCqWetBK20WR+GoPmnOtIla6V6vp33d5tZRHU/JWxftQrXgw/hCJf/HuQOKzZ2n7FHp6vKvFVJ9PqIw9dYkBVM9C1KUkLtXTey5DR2N/BfU5h0geO7n2q05QQYpR+0gx9SSX6+g1Tm53nsQH8umEqONbaK1sAEAf40U3A/O1Xxv5SJ2lmW15eXmYPn06kpKS4Ovri/j4eMyYMQOFhYXkGJcvX8aIESPg5+eHyMhIzJo1CzU1fxvUt2/fbnQOhUKBrKy/ja07duzAyJEjERMTA4VCgR9++MFiu6dMmQKFQoG33nqLLM/Ly8O4ceOg1WoRHByMyZMno6SkxPRBGMZD0cHSk+x1tFot+Zjr1MkEBwejdevWOH/+vPWNGxjWN4bxTMxr3M3H7pG61NRUrF+/Hrt378bJkyeRm5sLhUKB8PBwtG3bFn379sWoUaPQvHlzhxpkbmZbRkYGMjIy8Prrr6Ndu3a4dOkSpkyZgoyMDKxduxYAoNPpMGLECERHR2PPnj3IzMzExIkT4eXlhcWLF5PznDlzhowMREb+PRuytLQUnTt3xmOPPYZ7773XYnvXrVuHffv2ISbG+Mlo3LhxyMzMxNatW1FdXY1HH30UTzzxBNasWWP/hWEYN8WWkTpHKSkpwYULFzBhwoQ6HedG6lPjWN8YxvNwpZE6mzt1GzZswOuvv45du3ZBCIHExES0aNECHTt2hBAC+fn5OHr0KL777js899xzuPXWWzFr1izcdddd9jXIzMy2Dh064LvvvjPEiYmJWLRoEcaPH4+amhqo1Wps2bIFJ0+exLZt2xAVFYUuXbpg4cKFmDNnDubPnw/vG6apR0ZGIjg42GQbhg0bRl7zmCM9PR3Tp0/H5s2bMWLECLLu1KlT2LRpEw4cOIAePXoAAN59910MHz4cr7/+ukmRBIDKykpiEi8qKjK5HcO4C85MPvz8889j5MiRSEhIQEZGBubNmweVSoWxY8da39kKN0PjWN9Y3xjPw5WSD9vUqbvlllvw559/4u6778Y333yDgQMHkqfAGykqKsLWrVuxdu1aPPDAA+jcuTP27t1rc4NqZ7ZpNBokJydjyZIliI+PN7ltYWEhtFqt4Wl379696NixI6Ki/vZDDBkyBE899RT++usvdO3a1bC8S5cuqKysRIcOHTB//nz07dvX5jYCgF6vx4QJEzBr1iy0b9/eaP3evXsRHBxsEDwAGDhwIJRKJfbv34/Ro0ebPO6SJUuwYMECu9pSVwapHiCxqpWRIQ5591CPx7XOUqIkBY2N8ipJflGFN/Vb6KXD+aqoZ+W24LMkDvehr3l+zutOYm+pdmKVoJ67Ur2xgTVNE0zPGUJf6YWq6DnlWoqHK+mX6OaTR+IgJW2Dn+RBiVRRs3+b9mkkTr3ajMReZfR8qirpIkohANx67zISK6vpRjt+cn4eO2GhTJiw00h85coVjB07FteuXUNERARuvfVW7Nu3DxEREdZ3tsDN0jjWt5uvb0OCHiNxTZeWRtukPk1zoiGW+l1rrMzDUftTvQr3ox5hpaSP/lId1E4+9F4PjKH+tSVZtBMusqk9oaCCakmN3rgzoVbSl4E9A1JIHK0uIHGxXkPik9X0HB0U9BoFKWnuUL0kQErJ7dVbQzv0GSGnSfxx6q0krgmk7Y84bPwdm71Lc3V+NepdEt+SkGq0jzMwp3H26tv8+fON7o+kpCScPn3azB7G2NSp69+/P9avX0/ExBxarRb33Xcf7rvvPmRlZeHtt9+2uTG1M9uSkpKQmZmJBQsW4LbbbsOJEycQKBn1c3NzsXDhQjzxxBOGZVlZWUZtrI1rPSVNmjTBBx98gB49eqCyshKffPIJ7rjjDuzfvx/dunWzua1Lly6FWq3GjBkzTK7PysoirzyA60/poaGhxN8iM3fuXDz33HOGuKioCHFxcWa3ZxhXx5kjdV999ZUzmmTEzdA41jfWN8YzceZIXV2Tq9u0tbn0AtaIjo62a19bZ7YVFRVhxIgRaNeuHebPn29Xm5KSkpCU9PesmT59+uDChQt488038fnnn9t0jEOHDuHtt9/G4cOHoVA4d8qyj4+PzSZxhnEHnF1Roj64GRrH+sb6xngmzqwoUdfk6i41+1XG1My24uJiDB06FIGBgVi3bh28vP6eBh0dHY2rV2laitrY0kXq1auXXbPndu7ciezsbMTHx0OtVkOtVuPSpUuYOXMmmjVrZjhfdjYti1JTU4O8vDzOhs80KpxZ+9WTYH1jGM/Amr7dzOTqDuWpKy0thZ+fn9mnuJycHJw6dQq33367I4c3IM9sKyoqwpAhQ+Dj44Mff/wRGg1955+cnIxFixYhOzvb8Gpg69at0Gq1aNeundnzHD16FE2aNDG7XmbChAkYOHAgWTZkyBBMmDABjz76qKEtBQUFOHToELp3v+73+u9//wu9Xo/evXvbfK6bQU1/+lrm7MPGv9fJvX4j8S8Z9HpmH6OvharDqcckIIx6TNpH0lc0TTUFJA7xkjwpoB48rZp6Trxiac4lfTmtK6guo9+pWklzIgFAVkUYiVMjaRyhph4Q2ReTXh1C26SgtV47etP9ZU9dhIp+h8GRJ0n8cRd6/NJKmhtQm0qvkUI2KgLQedNOVM1NGDRxh5E6mZuhcaxvN4eCu6gnsN0zJ4y2mRp2lMQ7i2gt1vWnO5M4JryAxF3CrpC4hS/NCRetpqlpAlXUj1Yk+ddk7eiWQP+xH8ynvkCFlCevJJ/mkAOAYxXUN9hem2mxTSrJE1ct5dlMqaYa39JL0mRJ32RkD14/34skPhhGZ5dnZFL9q/Y1Pn7cNuq7e6b1QyTen2CxSQ5jbaROthjMmzfP5Ci8PRYNc9jVqVu9ejVeeuklXLlyBX5+fhgzZgyWLFli9GS2ZcsWTJw4ETqdfVlaLM1sKyoqwuDBg1FWVoYvvvjC0OMFgIiICKhUKgwePBjt2rXDhAkT8NprryErKwsvvvgipk6dahjyf+utt9C8eXO0b98eFRUV+OSTT/Df//4XW7ZsMbSjpKSEPNmmpKTg6NGjCA0NRXx8PMLCwhAWRv/he3l5ITo62vDqo23bthg6dCj+8Y9/4IMPPkB1dTWmTZuGhx56yOzMMIbxRPRQmq2B2BC1ES1RnxrH+sYwnok5jatddjOTq9vcqfv111/x6KOPok2bNnj22WeRmZmJr7/+Gj/99BO+++479OvXz9ZDmcXSzLbt27dj//79AICWLelTSkpKCpo1awaVSoUNGzbgqaeeQnJyMvz9/TFp0iS8/PLLhm2rqqowc+ZMpKenw8/PD506dcK2bdvQv39/wzYHDx4kca2xd9KkSVi1apXN3+fLL7/EtGnTMGDAACiVStx333145513HLk0DOO26IQCOjMjcuaWNwT1rXGsbwzjmZjTuNpltUnV7cWR5Oo2d+peeeUVdO/eHbt37zbkQzp//jzuv/9+DB06FKtXr8YDDzxg5SiWsTSz7Y477oAQJnI0SCQkJGDjxo1m18+ePRuzZ1tO22DruW4kNTXVaFloaCgn4mQaPe7y+rW+NY71jWE8E2dOlLgRR5Kr29ypO3HiBF566SWS4LJly5bYs2cPHnzwQTz88MPIzs7GtGnT7Gs1U2cGKe8nsUKaAq1MbEbi9OE0FcFtS6nH5MEAmjMJMP7j7B99jsQ/lFF/RFfJMzc8/DiJ5anesofEW0Ffa6VVh5K4Uk+/Y6eYDBIfKKZ+DPU16qGTPXYAgHLq0ZB9NAEdqIeuqTetXShfo9OV9DWUl/Sdmqlp3rtQafT+Nj+am0+fRDf4sIbmcaoopHnzNLR5148h3fF6Vf13qoSFihKiAZJzmoM1zjUZ7DOOxMpg6iWt7ETz/F0cS/+mF8z7gcTxXjR/JABUCKoPt2npvVfZmq6/VVrfwptOGsnTUU+vXNjdX1FFYtlTJ/vX5DrUJ2OoHaA0V/LQVRrfV9XVNE/mT5do7tGQROrpbeadS+ICnZ/FWKmgHr1oFf2OAQp6DaulIlrNveg1+1fUryRu6UcnCb13bShktBvpMWtWSvkrhxjt4hTMaZy9+uaM5Oo2d+p0Op3JfCm+vr744Ycf8Pjjj+OZZ55BdnY2mVLPMEzjprYOorl1rgJrHMMwjmBO4+zVN2ckV7e5U9eqVSvs3r0bTz31lNE6pVKJ//3f/0VYWBheeeUVtG7d2sQRGIZpjOiF+dcQJiboNhiscQzDOII5jbNX35yRXN3mscHhw4fjhx9+QF6e8dB1LcuWLcPSpUtx9uxZs9swDNO4qC12be7jKrDGMQzjCK6kbzaP1E2ePBnh4eG4evUqQkNDzW43a9YstGzZEseOHXNKAxnrKLyoV0Lfi+asOjeG+jXu7UfrVHb0M/bQyRTqqV+rgy/Ny9SyHfU7BKtonrkYNTV4na6kebOu1lCfTJSc10lJcyAVKqmfI8GP/iNOjaZ/o7llNNZcNZFDSXrQqk6j1+0L0YvEA1vTenwxPrTNMnmSZ6SpVEs2SEXbFCikPFA+1KfYM47mrtrTgb4SDDluov5jubTgJrz91Fuo/WpueUPAGueaKLypF6vwDlqbumIi1ZY5ibtJHKmmhVtl/xxg/JpMo6D33l0hR0kcqKQ3UqC0fZag97JKQXNIKqXYGuFeVCtahF0j8YkSmiJDcY3+TwAAoaLDRvnpVHM/qaH1gQcmnCFxOz/qW67Q0+uYJ3nsolWWC+aWCeq585KuWaya6uUTQfRBym8E3R8A3iy6h8RRB6qNtqkPzGlcQ+ibzZ262NhYTJ061aZtR48ebbagM8MwjQt3SWnCGscwjCNYS2lyM3GoogTDMIytWHoN4UqvXxmGYRzBnMY1hL7ZdcasrCy89tpreOGFF/Drr39PN37//ffRt29fdOjQAY8++iguXLjg9IYyDOOe6KEw5HEy+rjQ61eANY5hGPsxq3Gu/Pr10qVL6NWrF3JyrufLWbp0Kd555x2UlpZi/vz56Nq1K7y9vfHFF1/g559/xsGDBxEfH2/lqIyt3NtrAdSq674JkUI9cAov6dcoJRZVVdI/rB/PdSRxUFvqD4n0ojVKAaBCWP5TkWsbyh6VLMkzJ5dUqZYSqFUI6glRKuh3ClXTnEp5NdTzV11D/RnCX8pf5G/CbyblrpOvmzKVeuy26tqSeFDbUyROkvIqyfVrKyQPiV76vXlJtRdlX06ilLtK2YOu3x2cCBnNSZpPUCpHi56T3qDnrKZt2LdmptExrSEseOqEC3XqWOMajnv6L4Fa9X/315/UOyXrm1829VJlpAeT+GtNDxI/04zmO5PvQ8BYjzSSZ0722BXr6X1UqqeeNvleLdBRffKScrRpJc+wyovuXyj51TQq2h6/AJpDs7TEWK/l+rDKCnrvlafS2qKbJX0LaEXP0VIj6Zv0naskPdNL6ysFjcsEvSYqqe6yn5L+T5ispb5uAPB+eC2JX46gFomOM98kcf/xf5D4na7/z+iYtmBO4xpC32weqVuwYAGUSiUOHjyInJwcjBgxAgsWLMA333yDEydOYPfu3Th8+DB2796NsrIyLFq0qD7bzTCMm2B2lM5CpYmGgDWOYRhHcCV9s7lTt337djzxxBPo1q0bwsLC8NJLLyE3Nxfjx49HixZ/z0bq1asXHnnkEWzdurVeGswwjHvhLilNWOMYhnEEV9I3m8+YkZGB5s3/Lr2UkJAAAGjbtq3Rtp06dUJGRobRcoZhGh/uMlLHGscwjCO4kr7Z7KkLDg5Gfv7f+YBUKhV8fHyg0WiMti0tLYWvr6/RcsZxFIUlUCive0l0ldTbYOShO36RxFFRNG9dhjf1Z2S1oH631hpaww8AotUFJJZrs8oeOtlj8ld5LImDVFKeJxX1lFRLfjM5751K8sWkVoQZtflGlF6Sp87P+E9fVSF76uh6IXk89Ffpd0yLDyFxW396HWWfTXqNlsQaBc23Faik54uQ8j519btE4kQf6nG5u/dRyPzWhnZQdmfQGrmle2h+tqCL9uXTMoW75KljjWs4VFn5UCmv3081NdQvJnT03vX6k+pbTGQbEl/ypVpwrkkUiTtraH5HwFhf5Fx2pnLb3UhaNT1nmFTXOUJFzauyh08pefb8USWtlzx2lfRvz8eL7l8eaCI/W46P8TILVBTQv/vUMvodZf2RNTtP8kmrQOt7Byrp+jLJZ1gh6PY66f+c7LEDgNEBtE0xw/6XxAvOjyTx+sNdSfwODW3GLfPUtWnTBseP/12UPTQ0FOXlcibT6/z1119o1qxZnRvHMIz7Y+mJ1ZVG6ljjGIZxBHMa59KeuvHjx9v0ZFpQUICvv/4a/fv3r1PDGIbxDNzl9StrHMMwjuBK+mbzSN3jjz9u03aBgYG4cuUK/Pz8rG/MMIzH4y4jdaxxDMM4giuN1Dm9ooRKpUJQUJD1DRm7+O7MMmi1WpPrBinvJ7G+hPo5Ao9nk1jTltZdPVVIPScd/I3z/zTzziWx7KEr1tERjkDJMyd76LwkT16c9zVpPfVTyOiUdJA5xIt6YkrLJP+IdG8pQiXDHABdOfWQKPR0J52PkNbT/U+db0riFgH0O4VKtV5l302BnnpEQpX0GjRTU5+NRkF/J2Uq6vuJMSr0CrTxpr6724No/dqXKkeRuDKP1l/s+QjNY1eLrqrC5HLAfTp1tsIa53y+O/WaeX1TP0RifQnNURn0J70PrvaKJHFuNc2/JtdgBoz9riopR6SsV/L6aK8CEuukWY/+Sqo38vHKJA+ynBfPT0nv/SsF9O9PsvsiIsy47urVCunffQ3dSVEt5eUsoR65Py4lkLiVP/2/0tmPehXl7+ivkGpdSy8KAxQ01pvIJ3gjZXrj2q8aBf2O3XwKSLyizRoSP1D0DxK3fJXq2wv3fmf4ubzE/P8kj+rU/fe//8WePXuQn5+PiIgIDBo0CN27d3dG2xiG8QAEzBuGhcmlrgVrHMMwljCncQ2hbzZ36v7zn/8gIyMDn376KQAgPz8fo0aNwp49eyBumJXy73//Gw8//DBWrVoFlUpl7nAMwzQS3GWkjjWOYRhHcKWROpsnSnz99deIi4szxE8//TT279+PxYsX49KlSygvL8eFCxcwZ84crFmzBgsXLqyXBjMM4164y0QJ1jiGYRzBlfTN5pG6K1euGOoc6vV6rFu3Di+88ALmzJlj2KZ58+ZYvHgxrl27hlWrVmH+/PlObzDDMO5FjV4J6E0/P9aYWd4QsMYxDOMI5jSuIfTN5k5dYGAgrl27bvyurKxEVVUVOnfubHLbLl26YPXq1c5pIWM/UpJGcYUmwQ09RY3EPoOoATSzOtjokNWC/qlEqGkyTTl5p1w0u7MvTQqZWhVBYnliRLCSHk8nDSrLxuImklEZCukaSJMefHyNk3NWRNBlNf701Zqykh5DXUrb5JVDr9G2lNYkbtEuh8TyRAl58km15MjwkYzEcSp6jYulIuQ6E4YOjWRebuZFTeZhWmpCz4+mJnMfOvfj73NVmX8iFUIBYeaJ1dzyhoA1zkXR079Z+W9G1rfIQ1Rbau6k901aNU2wDRgn942UEn3L5Z5CJb0LBY1zdP4klidWyPdhtaR/crLjSElvvdV0/8pqqj2mRogiYgpInF9I21hTRCdqCb2kodl0MtzP2vYkbtKCHt9PQSeHlCrts/ArYXniRIWg1wAA/ECvW5CSTn4LkRK6f9qD3sMTsp8i8eurxhh+1lVWAPjDZFvNaVxD6JvN3cgBAwbgs88+Q1VVFXx9fdGtWzesW7fO5LY//PADKbfDMEzjpTbburmPq8AaxzCMI7iSvtncdX7llVfQvXt39O/fH/PmzcOyZcvwwAMP4K677sJDDz2EqKgopKenY+XKldi1axfefffd+mw3wzBugrtMlGCNYxjGEVxpooTNnboWLVrgt99+w4QJEzB06FAoFAoIIbBx40b88ssvAAAhBHx8fLBgwQI8/fTT9dZohmHcB3d5/coaxzCMI7jS61e7XnJ36dIFx44dw+bNm/Hrr7/i3LlzKCkpga+vL2JiYtC9e3fcfffdiIqKsn4w5qahr6CJYbV/pJE4dWszEvsNM07q2F5LfSstpOLxflJyTX8FPUaUisbRUiLe1JoQEsueEjmZZ3YN9XrlSnGolnpccvPp+kA/42S5tyecJ/GBrHgSF6TQNgol9ZwopeSd1RnUs7Ivhr6u8wqhnhCNhnr6sqRC5tXSNYhQycmI6fm95IykAPL0lhN6JgTmk/iqlt7L6lIpQen/NVlYMHK4y0gdwBrnFsie4Sp6H4Tsovr2c7veJL535G6rpwjzpfoUoaJeU3/JAxelUkvrC0mcWkOTBVfpqV9XZeQXo/pXKnmIm4dQc2tqAfUJ3hZ9ATJNpUS8fxQ2I/GhS1Tv9JIPGTm0DTnpwSQ+Gkn399dK/xOk/xFXdVSjq6VXlaGSpmgUKosxAFQKy0nrVZIm3uJDj/HvAetJ/OGrow0/66rMZ51zy5G6WhQKBYYOHYqhQ4fWR3sYhvEw3GWkrhbWOIZh7MFtR+oYhmHsRVgYqXPFTh3DMIw9mNM4l579WsvFixfx888/Y//+/dCbeZVz/PhxvPzyy3VuHMMw7o/A9TdmJj8N3TgTsMYxDGMPZjWuAdpi80idXq/H5MmT8dlnnxmWtWjRAu+99x6GDBlCtj127BgWLFiAl156yXktZZyG7iotxBz3C/WKHTKRquFaM+oPA011h8GBx0kcqKT+MNnvVSnlaSrW0RxIJ6ub0jZ65ZHYW87zJKg3om0o9fwVBlCPS6xfAWTCvKhvpmUozeF2qNCPxLoa6jFRSTY9r0L6zHTwAi2IHd2e5p6SiZJy78m5+Vp4099jmJSnLlBpnMfJS3pwDJQKh3cNokW5jzanv4dyHfUm+mVeP6BFTx0UUJiZ2u9KKU1Y49wTIXlPdVdpPsiEjWEk/qpZD6Nj9Gl5ke7jTe992VPnZ+XP1tpoSZ4ugMRyXk7Zo1ygo9rTLZj6BjsGZZDYKG8ngGgvqoFlgTQvXVpoMImz87QkrvGVHnAkz92udPp/I1TSUzkXqewblONoFdVH+X+K7LkDAD8FPYbsodNJfkwfJT3IyADqRfyfuL/311Wa/6Wb07iG0DebR+o++ugjrF69Go888gjWr1+PDz74AN7e3hg+fDiWLFlSn21kGMaNqfWbmPu4CqxxDMM4givpm80jdR999BHuvfdeQ7FrAHjkkUcwdepU/Pvf/0Z6ejree++9emkkwzDui14ooHCD2a+scQzDOII5jWsIfbN5pO78+fNGryC8vLzw0UcfYenSpVixYgUefPBBVFcbl19iGKbxYtZPJ4yyU9jFq6++CoVCgX/+859OaSdrHMMwjlAf+uYoNo/UBQQEoKjItAdo1qxZiIiIwD/+8Q8MHz4c9957r9MayFhnq/5bi+sHKe8nsaih3gbleeqjitnSzugYl0aEk/irPJp3qXln6mPx8kkncbWUP0gpWUjl2rGnKmJInAHq+2unocdv73vF4vG9/Kjv5lq15BEEUKanHpMuWnrM8DbUI7LDP5HEFeV0f10F9fkpCqnfY1dGCxLrm9BnrFxv6ruRc/VlVNFropE8J7f6n4FMKzXdpkKqkVuio7USI7U0X9flAHrdhPr6dzRRhvHvbeohpcmBAwfw4YcfolOnTg7tbwrWONfEbn2TPHaqk6kkjv7ZWN/+GC3lpKyiHt8nmv5Od5A8vhVW8qMFS35Xufa1xovel3Ld1EAV3b+DpHeyX+1EeZxxGySNlY+hjKFacDqI5mIskq5JbjnVgpxCqlcni6JJLOuT7PuTv0OBVD83ULqGHSVPMQD4Sanr5N+LyorHLa2GarTm2t/XxFKeOldKaWLzSF3Hjh2xZcsWs+sfeeQRrF27Frt378asWbOc0jiGYdwfZ3vqSkpKMG7cOHz88ccICQmxvoONsMYxDOMIruSps7lTN2rUKGzbtg2nTp0yu83dd9+NTZs2Qa3m9HcMw1ynNtu6uQ8AFBUVkU9lZaXZ402dOhUjRozAwIEDndpO1jiGYRzBmr7dTGxWpsceewx33HEHIiMjLW53++2349ixY0hJSalz4xiGcX8seUtql8fF0ddF8+bNw/z58422/+qrr3D48GEcOHDAya1kjWMYxjHMaZxLe+p8fX3Rvn17m7aNj49HfHy89Q0Zl0BfQn1TwX9kGG1TEkP/6RZ1ok8g2/KoT0UTRv0TEWo55xBN6hYq1YKN96a1DYv11Osl1369XEU9fxmV1PMn56DLqaL+DwDoEkg9Jh00NBdUD5oqCrcE0lqxMleleo9fnO9F4modNYBU6untWKGn/g4vhWVfYJyG+nzkXFgAkKOk2xRIPsJOvvQ7K6OoKq0+K9WCLbu+XmHRb2LeW1IremlpadBq/86L5ePjY7RtWloannnmGWzduhUajcZofV1hjfMQ9FLeuhJ674fso35cACiV9C3lVnqMY6H0d62V9KtMyvcYJOWZ81LQHG/Rkp9MztEmlYJFWhXNtXehgj54RHpbznkJGGtsU0mTW3jR3Hz9Aqg2yPVn9dKLvg/S7yBxhc5y/W6ZYj317Mm1YuXcpGXCuPbrVR2tAyzXk41V0TZVCvp7O1BOc+1p8v/WtZpq+zWuITp1dleUMEdlZSVyc3Otb8gwTKPCFk+dVqslH1OdukOHDiE7OxvdunWDWq2GWq3G77//jnfeeQdqtRo6nYXZGk6ANY5hGFPUl6fOkRn+dnXqtmzZguHDh6NPnz546aWXUFVVhdLSUjzwwAPw9/dHVFQUmjdvjnXr1tnbdoZhPBRh5WMrAwYMwPHjx3H06FHDp0ePHhg3bhyOHj0Klcr4yd1eWOMYhrEXZ+ibjKMz/G1+/bpnzx4MHz4cERERiIiIwKJFi5Cbm4uqqiqcOnUKc+bMQVlZGb777js88MAD+P3339GnTx+7vwjDMJ6Fs1KaBAYGokOHDmSZv78/wsLCjJY7AmscwzCO4OyUJjfO8H/llVfs2tfmTt3ChQvRqVMn7N69G76+vnjhhRfwxhtvIDk5GYcPH4aX1/V31fPmzUOnTp2wbNkyfpp1EezN86RLzzTaJmY79WdVRAaT+JRkLk/0p3nmqnzoKEq0ulCKi0ncRUNz5+VINUf/LKcel3WX6NNMfj5tb2xUPolvjzL2w7X0ySJxM8n3EiyNaxd70WNeqKbpNWQPSMdIel3P5tNcVVE+1OPSWkO3l/M8Xa0OhiXSq43Tfegkj0lTFf093OZLzynnhvoyXKqbKf7P12fpkdTSI2tDVLw2A2uce2JV31QPkFiXcdVom5jfqL6caUnjS7HU0+aloDkq20p5MzVSzjW5bmkrL+oZzpFqX1+T/LA/Xu1M4rxyavCVa1331tJatoCxRy1KRY17GgXV6DJB7/2MGvodsnS0Nmy8P/XrXikLpueTas/KvkIvSPW8YXnkPa0m2GiZXEM3VCrI7aeUcokKeg0Cpe2JDdDSe01zGvd/y+T8lz4+PiYtJrXcOMPf3k6dza9fjx49ivHjx8PX9/of36RJk1BVVYVJkyYZxA4AgoOD8cgjj2Dv3r12NYRhGA/Fkt+kjp6T7du346233nJKM1njGIZxCCv6FhcXh6CgIMPHUi3p2hn+jtabtnmkrqCggCT6DA8PNzRWJiEhAfn5+UbLGYZpfNiS0sQVYI1jGMYRrKU0sWV2f+12dZ3hb3OnLiIiAllZf7+e8vb2xpAhQxAREWG0bXZ2NoKDgx1qEMMwnkV9lAmrD1jjGIZxBGueutpZ/da4cYZ/LTqdDjt27MB7772HyspKqxPCbO7Ude3alST8DAwMxC+//GJy2/3796NNmza2HppxMeTaiQCgkOrDRhymf6CZCTRn2p8FsSQuC6ReBm0g9S4A1FMXp5aLptP1IWqae6ptGK0DmOIl1ZpVWB8Skj0nKiOTBL1pvaz0Rwp01Pei9aLf2d+b5lSSc+d18qPXIEzKMyV76golX47snwOAXKkmpJcv/V2396ZtjJTOOTrpGIk37e1j9lwGLL1mdaFOHWuchyINoZjUt3OXSNz0Vzrx5nB8UxLrIyz/3cYH0FFcP2lzPyV1PlUIqgUFeur16hd+jsRyner8aqo1ck5LwNjnZw0/yWPnp6T7a/RUn5r6FJA4q0L6HyHpVTNvmh5Irk0r+6jlPHmmNKdIymcaqMwy2uZGKqXasDpZj4SZn2XMaZyd+lY7w/9GHn30UbRp0wZz5syxaYa/zZ26559/HmfPnrW6XU5ODq5cuYLJkyfbemiGYTwYd3n9yhrHMIwjOKuihDNm+Nvcqbvttttw2223Wd0uIiICBw8etPWwDMN4Om4y+5U1jmEYh7Ay+/VmwlWpGYapV9zFU8cwDOMIzs5TdyPbt2+3a3vu1DHGmBgz1pdRv0PQYepVKGhJPSd5EdTjER5KvVlyHcIyqc7pScmPJteKbepFcyINC6M+BFUY9aRkSDnbfJSyZ8/4HKWC3h6lkk0lUPKoJEp56+TaiDopDpXq0aqk+pDnK6JJfFUl+xAtU22iNqLsxZHbmK+n5ygS1CupkopS1tZG1Fmo/QrApUbkmEaOlJsMAPTl9O9ee4jWv77WXpoBPYKub+lD88TJtV5zpLrOcm1Xf0lL5Hxrvf0uGLX5Rix6Ws20KUOqPS3rWYSKtjlUSoBWraL61dU3lcRJUp5N2WMs5xoNl3KVynk5ZbxMeATLJN+d/J0L9TT3Xp7krzxb0YS2Ie/vc9TUWPEkuojGOVT79eLFi2jRogUSExPtWscwTOPDltqvrgZrHMMwtuJK+ubwSJ1er4dCYbrBltYxDNPIcBNPnQxrHMMwNuHunroWLVogNTXV7nUMwzRGFJDTwdB1rgdrHMMwtmNO49xopI5pZEg+O10arXUYs4t61s63DSVxboRUO1by0JXBfB08ANDp6c0h+98iVJIfQ8rT1M6H+jtkDx8AFEs5joql9f4KmktKo6Ln8JIey1pIvj/ZR6jXUPdDanU4iXcUJZHYV0nP39H/ColjJE9fVk0QZPykY8htOldNc939VNCVxN+f7ELi+KvXfS81NRb8L246Usd4KCZzT9B7Wa4P23QHvTeP9aa1rfsE0VrSWVLt1mCphnKF5NdVKWmbIlR0ezlnpo/UV5Dds2UmvqLsEdZLrwYliUWxnnrIgqS6qQlq6lfzU1D9qZY6NGk1NG9dRhX9n5FbQ/PSyZ68CBWtn1qsp1oFAFrp/4LsAzxfTa/U8coEEn9zqhuJE8pu+LuoMc79Z8DdR+oYhmFshjt1DMN4Mu7aqfvyyy8RHh6OIUOGAACKi4sxbtw4o+0SEhLw7rvvOqeFDMO4N25SUQJgjWMYxgGcVFHCGdjcqdu4cSMmTpyIzZs3G5ZVVVVhw4YNiIqKMhSfFULg559/xvDhwzFs2DDnt5hhGLfCXSpKsMYxDOMIzqoo4Qxs7tR99tln6NWrFwYOHGi07ssvv8Sdd95piHv16oVVq1ax4LkJW/XfWt1mkPJ+EgspZ4/6rxQSh+xtR+IDMdS7INdiTdTQ2q2RauqfkImQciQlqC3fPXrp7spTlBtvJFkm5FqJPibqKd6IXAs2Qkn9G8GSny1LyhPVwptegwJ/mtdJvmZ+Uq1aW3JVpVdSH0sbH+oj3F/aksSXyqg3Up1Ct9dcLQQA1OhoWwhu8vqVNc4zcUjfqum96nWC1obVbaB+1yOTqb4V+9H7pI3k6ZXv3Qopp6SfpDUxaqnuqaRn1ZJ4qUzk4qvQ0X00SrqP7NOT8ZJqwXopvKUtqJ8tRyfnvaOa3c6X+rJln7VGYTlPnSm9q5DycOZJl+FCdQSJj5bSXHm+f9C8nF65OYafFY5oXAPom8156vbt24fhw4fbtO2oUaOwf/9+hxvFMIwHUftqwtzHRWCNYxjGIVxI32weqcvKykJcHM2q7evri2eeeQbx8bS327RpU1y9SmcPMQzTOFGI6x9z61wF1jiGYRzBnMY1hL7Z3KnTaDQoKaHpD/z8/PDmm28abVtaWgpvb3lolmGYRombvH5ljWMYxiFc6PWrzZ26Vq1aYefOnZg2bZrVbXfu3IlWrVrVqWGMe6Erov8Mo3ZeI/HpdrSOabQ/9czFetMcR7K/wktF/RlyjqVKKc9UleQ5kR0mpnwHcm6oKkG38pbqCMoeOn8F3V4pncVLaoUc6yWPSJw3vYY5Up4nmWIdzdsk++cA4Ou/upP4J/8OJNb6Sr6YAppvK4haJ6EouX7NLPtN3GP2K2scYw5dfiGJo3fSHJQ/d+lM4pE9jpC4mXcuPaAkSColXVAN6ierEFT/yvQ6aT09XrDSWOGCJA9dtbSPnOsuUCnltZMaLe/vo6CaHSifT2pznJpewxwdzVMne+oK9FLt2DLqYwSAE8U0f+AG6f9GtZ5+y2PZtNZr8EW6vaLy7zYo9JZycbrO7FebPXX33HMP1q1bh3379lncbv/+/fj+++8xevToOjeOYRgPQFj5uAiscQzDOIQL6ZvNnboZM2YgLi4OgwcPxuuvv47MTDqbJzMzE6+//joGDx6M+Ph4zJgxw+mNZRjGDXGTTh1rHMMwDuFC+mZzpy4wMBCbN29GQkICZs+ejdjYWISFhSEhIQFhYWGIjY3F7NmzERcXh02bNiEwMND6QRmG8XzcpFPHGscwjEO4kL7ZVVGiZcuWOHr0KL755hts2LABp0+fRlFRERISEjBkyBDcddddeOCBB6BWc/WxRofklxAXL5M4fgutQ3qmWSSJm2oKSBwVSGO59mG2VFuxSFC/Q7WU90muFRunMvaAaSRPXLV0zmLJB+OnsOyXKJPa5Cd5TiJUdH2GlEcqTl0gbU+r0eolz99l0Jxyh/PpTE4A0JfQNpTq6HeoOEd/T/5X6PpDnzxrdEwAKCoqQlDQOybXuYunDmCNY8wg61sqrbscuzmYxAWdqb/1ZHlTEvfwu2jxdFk6mi+tQlRIW0i1YyWtCjVxWwUr6D4VQvbYWe6BlEmeMi9JL2VPnUbOawfqVwuV8nbKtWNlzlTT/xnnyiKNtqmRNLGiiuYLPJpCNdH/OF2/a51pfQNqNc5MFRkX8tTZrUwqlQpjx47F2LFj66M9DMN4GO6S0qQW1jiGYezBLVOaWKKqqgr79+9HZmYmkpKS0LlzZ+s7MQzTOHCTlCaWYI1jGMYsLpTSxGZP3ebNm/HYY48hN5dOzT59+jQ6dOiAO+64A2PHjkW3bt0wZswY1EhlpBiGYVwZ1jiGYdwdm0fq/vd//xfnz59HeHg4WT5u3DicP38ekyZNQo8ePbBx40asW7cO7777Lp591vz7acaz0VdSz5r/Ieqxu3o4ke4QS8PzFTSvXZKGzkSUPXOhoHnyVFJOuWI99U4UmKgrGKOm/gcv6SmrQMqdVyB57GSPiEayU1RKnpJAJU1e28JKbVm5OO0lqVMh5/bLKDLOa+cTaqLm7Q2IKlpjUl1W90dNBSy8fq3z0Z0HaxxjK/py6nHT/pFG4t37aO3rPkMukLhC0HtVrtWqk7xhBdL5/RVUa2Q1y9VVQSZKJfvw6N2XIeVwK9DTNgUrjevJUmgrZI9dlMryGJJcW/aqzrJWVeqMuy+Bavp/p7iG6plSqhHuU+CcoTRzGtcQ+mbzSN3BgweNCl0fOXIER44cwbhx47By5UpMnToVP//8M2699VZ8+eWXTm8swzBuiJvUfmWNYxjGIVxI32zu1GVlZaFly5Zk2aZNm6BQKPDII4+Q5ffccw/OnDnjlAYyDOPmuElKE9Y4hmEcwoX0zebXrwEBASgrKyPLdu3aBaVSid69e5PlwcHB0OmsvUpiGKZR4CYTJVjjGIZxCHecKNG2bVusX7/eEOfn52PHjh3o06cPAgJozrC0tDRER0fLh2AYphFSO93f3MdVYI1jGMYRXEnfbB6pmzlzJu6++24MGzYMffr0wU8//YSysjI8/fTTRttu2rQJXbt2dWpDmYZlq/5bi+sHKe+3uF53jSaWjNtaSuJtXZNIfFtTmpwzxpvuH6EuInG1VI5aq6BGZo1CKogtjP/0L9VQI3CelAC0VNCJDXLCTxVom4LVdCTHWrJOOa6WkoNWCgsFpQGU6akpuLra+Dt6e9ProNPRNukkL7RvnjVztA24yUgda1zjxW59kyc2XM0mceI3EST+oXMXEg+IPE3idpp0Egcr6YixrDV+kp5p5IlhJvWNHqNAmjxWrKcJk/2UdNJBNeQE7vScSmmMSJ74IMfW9EymSEfb66sy3l9vZWqCPp9qeGCafW0wizuO1I0cORKvvfYa9u7di3nz5uGvv/7Cf/7zHzz44INku3379mHfvn0YMWKE0xvLMIwb4iaeOtY4hmEcwoX0za7kw88//zyeffZZ5ObmIjIyEgoTZZI6d+6MnJwcBAcHO6uNDMO4Me5UUYI1jmEYe3H7ihJCCFRWVkKj0Rit8/X1ha+vr4m9GIZplLhR7ddaWOMYhrEZd6z9KoTAv//9b7z33nsoLS2FSqXCiBEj8OmnnyI0NNT6AZjGhVQcWkgzBdUnUkjs/xVN1nnlqWAS+6poQeyeAdTv0dGbJie2hilP3dEKmgG5QvLQNfXKs3jMKsnNoJR8LrJnzhplgiYQzbDiiTlZFkNivd5YUHwkn19hBW2T/1W6j+aacRJTu3ETTx1rHGMzVvRNdYJ6gvM/bU/itGepBy9URT3GEZpiEgdLic1lj52XdKurhPGNlaWjk33kBMiyj88aSmkUW2ViVNseSvRU08ukDpHs+VOaGAZTStclrTiYxL7p1NfnXWg5wbHNuJCnzuZO3apVq/Dqq68iNjYWQ4cOxYULF7B+/Xro9XoyY4xhGOZG3OX1K2scwzCO4JavX1esWIGuXbti165dhlcPzzzzDJYvX47c3Fyj0joMwzAA3GakjjWOYRiHcKGROptnv164cAETJ04kXpKnn34aer0e586dq5fGMQzjAVjK4WSn6K1YsQKdOnWCVquFVqtFcnIyfvnlF6c0kzWOYRiHcJK+OQObR+ry8/MREUFz79Q+uVZUVJjahWlE2JvnSVdCPSQh+6kn7mSXFiRuNZh6UGT/R4ya+i/yJI/LpRpa3F72uwGAViX9HUsFA8KUtM2Bks8lSKmTYuphsZcgJfWQXJLyRJ2roslvL5dS35fGxzgHk1pF21iTTc8RmE7zXykrnVA1wYkjdbGxsXj11VfRqlUrCCGwevVq3H333Thy5Ajat29v/QAWYI1jzGFV31QPkFhfSvUpbA/Vtw3J3UicdGcWiWXPXJSKjr/kSvp2sYZqjUZhfO/Ly+RzBEt56fwU9ByhKupHk/VJzqtpDdljnCF56lKrQ0j8Vyn1VdcI4zGpSh3t0mTnUt2PTKW6r86muUUdxoVG6uya/Wpqej/DMIxFnNipGzlyJIkXLVqEFStWYN++fXXu1AGscQzDOIC7dur+9a9/YcmSJYa4tvbh448/Dn9/mn1foVDgzz//dEITGYZxZ2yZKFFURJ+YfXx84OPjY2KPv9HpdPj2229RWlqK5ORkZzSVNY5hGLtx1kSJFStWYMWKFUhNTQUAtG/fHi+99BKGDRtm8zFs7tTdfvvtJp9iIyMjbT4ZwzCNEBtG6uLi4sjiefPmYf78+SZ3OX78OJKTk1FRUYGAgACsW7cO7dq1M7mtPbDGMQzjEE4aqXOGvcTmTt327dvtax3DWEKunXiFek7ittCZhqUD6KiNHpY9JhWC+j+ydYEk1iqNPVIxalpf1ktB9wmSPCdxatqGgDp66KzhBXrNSqVar+nFQSQ2pSc6PW2zqlyq11hEfYI1fvbl1jOFLSN1aWlp0Gr/9r9YGqVLSkrC0aNHUVhYiLVr12LSpEn4/fff69yxY41jHMYoL5ykbxnUM9fsR+rdrLyD3mcFUk62PF2htJ7m0JTrPnsrjf1t/pKnTil5ghMkX7Kvwo/EKql2tU44oS70DWik5yk5j55a8vh5K40fwGRPnb6Cxn7Z9BrUhFONdxRnjdQ5w17iUEUJhmEYu7AibrWzWW3B29sbLVu2BAB0794dBw4cwNtvv40PP/ywrq1kGIZxDAsadzPtJTanNHnttddw6tQpcsI//vgDJSUlRtvu27cPjz32mM2NYBjGg7FU7NoJRmK9Xo/KykrrG1qBNY5hGIewom9xcXEICgoyfG707cocP34cAQEB8PHxwZQpU+y2l9jcqfvXv/6FI0eOGOKCggIkJyfjjz/+MNr2woULWL16tc2NYBjGczGXo87Sa1lzzJ07Fzt27EBqaiqOHz+OuXPnYvv27Rg3blyd28kaxzCMI1jTt7S0NBQWFho+c+fONXusWnvJ/v378dRTT2HSpEk4efKkzW2p0+tXYaK+HMM4glw70ecgTfa66//RvE6dHk8jsQ7UX1Et1Xb1l/xw8vYAECZt00oteVAkE32A0rjYe30i5+Jr6kU9gL5etL0+appzDgCqdNRrqKUlKvHbtn/VoYVmcGJKk+zsbEycOBGZmZkICgpCp06dsHnzZgwaNKiurTTdPNY4xgmIanovag5cIPGq/x1K4ueeWEti2V8m55iTc9BVgd7nAJCgonVOvSQ981JQn55eOod8RHl9XQlV0deRsmZ7SZ46We8A4FolnaEecoj+H/ht63N1aaJ5rEyUuJn2EvbUMQxTrziz9uunn35a9wYxDMM4kfqs/WqvvYQ7dQzD1C9uUvuVYRjGIZyU0mTu3LkYNmwY4uPjUVxcjDVr1mD79u3YvHmzzcfgTh3DMPULd+oYhvFknNSpc4a9xK5O3caNG5GVdT3fTllZGRQKBb799lscPXqUbHfo0CF7Dss0RqzkddKX0dqJsT/T2q+vJd5F4tG3HCBxC98cEid60/1j1DTvEwAkelnOM1eop56UEj3NdeeloK4TnfQdKwX11YSoaB4oGbmWokZBb9cKPfXZVNTQ9UEa41x8ZVV0n8B0J9R2tYIzX7/WN6xxjFOwpm/FxSSO3ZRL4gVJd5P43h70762dXwaJO2qox1j2BwPGnjW15JIrkvSsWuqRBCmp565C0rMyPdWSSJXlPHcyci3YYilXX7mUm6+thuY2BYC9Zc1JHJRi7LurD5z1+tUZ9hK7OnVr1qzBmjVryDJz5j2uocgwDAC3GqljjWMYxm7csfZrSkpKfbaDYRhPxU06daxxDMM4hDt26hISEuqzHQzDeCju8vqVNY5hGEeoz9mv9sITJZibwlb9txbXD1LeT2I5b524dIXEzdbROqc/hXUgcb/mNA9UqzBae7Gl2nre7UpB/Rh6ySeTK3lIwlXSeiv1aEOMU0kRZE9dNWj8awHNMl5YQj0opjx1Bfk0j1P4tbpXYrCKm4zUMYyj2K1vevqHLy5eJnGzdR1JfLhFHIm1anpvt/dJJ3G4ylhcZI9vBeTcnXR9qV6+OWmtWPl4ZdLm1jx0Mvk66qM+WEqvQbWefqer1cZ53zJyg0ncvKz+PcMAXGqkzqarPmTIEOzYscPug//2228YMmSI3fsxDOM5OLOiRH3BGscwjKO4kr7Z1KlLTEzEoEGD0LZtW8yfPx87d+40WQ+xuLgY27dvx4svvoikpCQMGzbMkBmZYZhGSj3XfnUGrHEMwziMC+mbTa9f33//fcyaNQtvv/023n//fSxcuBAKhQKhoaEICQmBEAL5+fnIz8+HEAKhoaEYN24cnnnmGTRv3tz6CRiG8Vzc4PUraxzDMA7jQq9fbfbUNW/eHG+99RZef/117Ny5E3v37sXp06dx7do1AEBYWBjatGmD5ORk3HrrrfDy8rJyRIa5ATk9hOTX0EtlUjRH6EzFgN+TSJzflPrL0qrDSJyiLjJqQqya5qGrFDS3lOwxqRJ0oDulmsZyrdgIpXEt1hvRSecrkTx9x6uoh+RoTlMSK6Sx/qyiQKNz+B+n9WpV5QUW2+QMFP/3MbfOVWCNY24akh9XT+1q8D1IizKnbab6du0hmqfur0qqBUA6ZGR9KxP07gu18t4urYb+vQcraaNDlZYPIHuUZc/w8Wqa165ER/PqaaX251QZ5xX1O0R1X+91EzzDMK9xDaFvdk+UUKvV6N+/P/r3718f7WEYxtNwg5G6G2GNYxjGLtxxpI5hGMYR3CWlCcMwjCNwShOGYRoPbjZSxzAMYxc8UscwlK26byyuH6R6gMS6fFq7NXrbVRKfDGtN4iNtaJ6nkW2PG53jnuDDJNYoqAdEqaCeN38F9ch5Setl5CqEcu3YMsljcrqa5pT7tag9iROC8kh8zYdufyUnxKgNF1571mIb6w3uvDGNGKt57CR90xdQfYv7kerbVv+eJN7UhuasHNPmiNE5RgXRZVVS7dcKQfVI1rdAJVUwneQYK9BT/VMqqAeuWPIR5ulp9+NkBfUFNvGm1yBIRfPYbcjuBJkTrzeQvgEuo3HcqWMYpl7h168Mw3gy/PqVYZjGA79+ZRjGk+HXrwzDNBZ4pI5hGE/G7UbqDh8+bH0jCyQlJcHf39/6hgxjK1JON30qrQ0bt43mMLpaSOMf87obHTJxQA6JW/hQH0u1oLdLjDqfxM3UNG9TRg3d/lw1bUNnb5or75qOelT2lrYi8dVKmqcuyIt6YAqraI4mXT7N89RguMFIHWsc06BIeTnl2tf6lDQSJ2ykOd3Sy2hOyv9X0tvoFG1uo7ntotXUs5Yj1aauVlJPnJ/ksSsWNG+dN6gmByvpd6iU7vUD5c2M2ngjIepSEqdWhJP4VFq0xf1vKu42UtejRw8o5OSwdrB161bceeedDu/PMIz74g4jdaxxDMM4ituN1AHAE088gVtuucWugxcXF+OZZ56xu1EMw3gQbjBSB7DGMQzjIO42UgcAt912Gx5++GG7Dn7t2jXMmDHD7kYxDONBuEmnjjWOYRiHcLdO3bp169CzZ0/rG0potVqsW7cOnTt3tntfhiEIy3eHqJFyyv1Fa8Nqo9qSuLQp9Y8AwGcp1IcyteV2EkdI9WIrJE/J8SpaV/VcJfV8RHtRDwtAj3ephuaV25/fjMQh3tTjUi3Vnr0k5aXTZBh/x4bAHV6/ssYxLoUVj53q3GUSa5tL+pZgXId19ZU+JH46fjuJg1XUw1YqvElcoLdcT1v2GMeqC0icrpN8ztVBJO7gS32DedL2OVXUN+h9nranIXG716933323Qwf38vJyeF+GYTwENxipY41jGMZh3G2kjmEYxlEUQkBhZqTV3HKGYRh3wZzGNYS+Odypu3TpElavXo2LFy8iPz8fQmq8QqHA+vXr69xAhmHcHDcYqTMFaxzDMDbh7iN1/+///T9MmjQJNTU1CA4ORlBQkNE2dUkPwDB2I/3D1ZeUkDjgL5qDLijSOMdRWRHNg/R61SASD0w4Q+LUEuop8VbRPE5VOnp7dQiieaKqpLxQf5XHklgveeaaaKgn70RhDImri2leutBLlmvR3izcwVMnwxrHuBRSXk5dEdW34INZJC6NotoAAJdLaW3VRWXDSNw5kuqTNfxVNC9nr0Daxl2C6lGW5KELVdPv4K+sJHGZnu5fWE09yyFnXEPfADf01MnMnTsXbdq0wdq1a9G6dWvrOzAM03hxw5E61jiGYWzGhUbqjKfI2EBubi6mTJnCYscwjFVqn2LNfVwR1jiGYWzFlfTNoZG63r174/Lly9Y3ZBiGccOROtY4hmFsxt1H6t566y188cUXWLt2rbPbwzCMh+GOI3WscQzD2Ior6ZtDI3UdO3bEokWL8NBDD8Hf3x+xsbFQqajpW6FQ4M8//3RKIxlmq/5bi+sHqR6gC+SJE5fSSRxxwDhxZWGSlsQV+TT+JdZyctoaP3rOpu2ukviqhh7v7awB9PxltE1dommbZS4XBpNYWUrvwYB0amRuMNxwpI41jrmZWNU39UMW1+vS6CSHJr9pjLbJuYUmJy8rphO9dobQ9UbnCKITwZLbXiDxlSp6vMOFcST2VtIEyiPCj1k8X6k0UeLctQgSh+TShPMNiguN1DnUqXv//fcxffp0aDQaJCYmmpwZxjAMU4urjsiZgzWOYRh7cBWNc6hTt3jxYvTp0wcbNmxgsWMYxjJCmC/z5qLJh1njGIaxGXMa5y7JhwsLCzFu3DgWO4ZhrOKOeepY4xiGsRW3z1PXr18/HD9+3NltYRjHkZ+IpMSwoob6LxRnU40OofVKpPuo/EmsrKHHrNLSWFVB4/TsYBJfLaAFqXWX6fH1PjSZZkkY9ZScLo4icX42PZ5vLp335JVfAZfADT11rHGMSyElH4aC3utCR/1qivOpRocI8/Oi2+gDSOxVSo9ZGSL5kjXUU3oun3rcLhVTT15GJo2jowtIXB1Gjyd76FIq6fHLTtDjhVe7iL4BLuWpc2j264oVK/D777/jtddew7Vr15zdJoZhPAiF3vLHFWGNYxjGVlxJ3xzq1LVr1w4pKSmYO3cuIiMj4e/vD61WSz782oJhGAB/P8Wa+9jBkiVL0LNnTwQGBiIyMhL33HMPzpw5Y31HO2GNYxjGZpykb87Aodev9913H9c9ZBjGJpzpqfv9998xdepU9OzZEzU1NXjhhRcwePBgnDx5Ev7+/tYPYCOscQzD2IqzPHVLlizB999/j9OnT8PX1xd9+vTB0qVLkZSUZPMxHOrUrVq1ypHdGKbesDePnb6szGgb1UWaF84/oBmJy6K8SVwZQv/pV4RLY+2FdPtqabWK2mCgkDx7pTV0fx8VzRMFHd3+9MvP0vUvwzVw4uzXTZs2kXjVqlWIjIzEoUOHcPvttzvaQiNY4xhXYqvuG4vrjfStstJoG9W5KyQOCG1J4mIV7Q6UUQsvEEB9yblXgkms8KX6pFDSezvAm+bNlD10FXrq+btYGk7ic3Ofo+2ZC9fBSbNfnfHQ6lCnjmEYxlZsGakrKioiy318fODj42NiD0phYSEAIDQ0tE5tZBiGcRRnjdQ546HVJk/d4cOHUVBQYF/rAOh0Ohw+fBilpaV278swjIdgg6cuLi4OQUFBhs+SJUusHlav1+Of//wn+vbtiw4dOtSpiaxxDMM4jBV9KyoqIp9KEyOppnDkodWmTl3Pnj2xceNGmw9aS0FBAXr27In9+/fbvS/DMJ6BLbVf09LSUFhYaPjMnWv93crUqVNx4sQJfPXVV3VuI2scwzCOYk3fbuZDq02vX4UQOH36NHbs2GHzgYHrvUzhohnjmUaGDX+Hurx8EntLnpHKEJo3CXLqqFD69KUvk24v6RFKL3ni5OnvxZX09WMxaKwucBP3hA2eutoZpbYybdo0bNiwATt27EBsbKwTmsgax7gxtuhbIbU4+J7PJXFFKDXRqSqpYAWG0NHoYi9aqzpYS33KRaW0/mxhBY0vllM91UkCufcv6vlDMlwXK566tLQ0om+2WEtqH1p37dplV1Ns/q/wyiuvYNGiRXYdXAjBM8gYppHjzNmvQghMnz4d69atw/bt29G8efO6N/D/YI1jGMYRrHnqbuZDq02dut9++82ug8p07ty5TvszDOPGOLGixNSpU7FmzRqsX78egYGByMrKAgAEBQXB19fXyt7mYY1jGMZhnFRRwhkPrTZ16vr162f3gRmGYQDnjtStWLECAHDHHXeQ5StXrsQjjzxif+P+D9Y4hmEcxVmzX53x0OomphyGuQlIngiRlkFivzBaK7EiiOYNKqmgt5MmlNYm9NNQz13+NVq7VZ1G89LlH6GeE3+aZgrhhW7i5dKL6x9z6+yA/WsM4yB6mhhTL+lbYBStkFIW5UfikjLqA+sYS/N6RmhKSHwirwmJs3Lp8TdVtCVxxSWqhzH7QHkcros5jbNT35zx0MqdOoZh6hcnvn5lGIZxOZz4+rWucKeOYZh6RQELr19vaksYhmGcjzmNawh9404dwzD1ixPLhDEMw7gcTioT5gxs7tQlJyfj448/rnPmdoZxF+T6ieo0mtfJN4rmXdJcobUL/aILSdwjKo3EhaHU+JoVQ6e8V9TQ21NxgtZC9CqTEtu5KM6cKFGfsMYxjQlRRWuxeqVmk9g/Np7EJRnUY6dsSm/eCG/qqesffY7EqiZUr66Uh5B436FOJPbLLDfVbJfEWRMlnIFNFSUAIDU1Fd27d8cLL7yAiooK6zswDMMANpUJcwVY4xiGcQgX0jebO3VnzpzB448/jtdeew0dO3bEtm3b6rNdDMN4CAohLH5cBdY4hmEcwZX0zeZOnVarxfLly7F3715otVoMGTIEEyZMQE5OTn22j2EYd0dv5eMisMYxDOMQLqRvdk+U6NmzJw4cOIB3330X//nPf7BhwwbExcUZbadQKPDnn386pZEMU1e26r+1us0g5f10gfSUpc8vILHfZeoJ8Yumnri8HBrXRKhI7K+inpZ2IVkkrtbT7Q9ERNLz/VUNd0ChF1CYyddkbnlDwhrHuBvO0Ddd7jUSay+Gkbg4juaRO5FJ89B11NK8dxol1acEb+pJDlBRi8POgI4kVlbWwF0wp3ENoW8OzX6tqalBTk4OKisrERYWhrCwMOs7MQzTOHHD2a+scQzD2Iw7zn6tZdu2bXj66adx8eJFPP3001i0aBECAwOt78gwTKPEXWa/1sIaxzCMPbjl7NecnByMHz8eQ4YMgZ+fH/bs2YN33nmHxY5hGMvUPsWa+7gIrHEMwziEC+mbzSN1SUlJqKqqwquvvornnnsOKpXK+k4M40HIeZ1U2fkkDkintWDLLtFartcS6foo3yISx/gUkHhbVhsSq0upQKiLaXtcFYX++sfcOleBNY5p1OhobVhVGs1bF5RC89RlNKOxPonWT1BKN3egiuadO5SbQOIAmsYTqmL3SStkTuMaQt9s7tTdcssteP/999GsWbN6bA7DMB6Hm3jqWOMYhnEId/TUbdy4sT7bwTCMp2IpCafr9OlY4xiGcQxzGtcA+sa1XxmGqVcsJeF0peTDDMMwjmBO4xpC37hTxzA2IiTPiZy3zv9iAIkDImkajAt5NI5qSj11IepSEqf9FU3ilI+fs7mtLoWbvH5lmMaMkHKq6Qto7WrtyTwSF7aIIPG5rjSP5n0Rh0jsBaqff1xoRuLU5W6qb4B7vn5lGIZxCAHzmdW5T8cwjLtjTuP49SvDMJ4Gv35lGMaT4devDMM0HgQsvH69qS1hGIZxPuY0jkfqGMZ9kPPWIS2ThNpomrQ2N0vy3MVXkvhyJfXcBabYnBvctWFPHcO4HaJGqr16mdZ2DT8eTOIjvWJJPDl6B4lTq6gHT3NGU7cGuhLsqWMYptGgB6CwsI5hGMadMadxrpx8mGEYxhHYU8cwjCfDnjqGYRoP/PqVYRhPhl+/MozrsVX/rcX1g1QPkFjO6yRKSkisuUBrJwakxJG4vAetDbvpRDsSx2TQvE5uC3fqGKbBsapv6odIbKRv5bQWq9+pLBLXpFBPXUVXLxJ/cOF2EmsveZD3gjt1DMM0GrhTxzCMJ8OdOoZhGg08UYJhGE+GJ0owDNNY4IkSDMN4MjxRgmE8AWH5MUyXRT11kYdobcQ/74whceBhmrdJc416WNwWfv3KMK6PXvLwKswNr19H1reYHU1IvH1gWxJX/B5O4rD8ajsb6MLw61eGYRoNegEozIibnjt1DMO4OeY0rgH0jTt1DMPULzxSxzCMJ8MjdQzDNB4sdOq4+CvDMG6POY3jTh3DMJ4Gj9QxDOPJ8Egdw7gfW3XfWFxvlJy4qorEPn+mkrhsdSsSB5XRAtpe+R4yUUIvYPaJlT11DOMS2J18vZrqVcDBSyTe82pvEgfp6EQM3/RSe5voupjTOPbUMQzjcQi9+ZnCVmYQMwzDuDzmNK4B9I07dQzD1C/8+pVhGE+GX78yDNNo4NevDMN4Mvz6lWEaH7preSQO3XGZro8OIbGinHry3BYeqWMY90e+VwX1yNVk55I46Df66rG6TSyJFek5zmtbQ+PEkbodO3Zg2bJlOHToEDIzM7Fu3Trcc889Nu+vtPuMDMMw9iDwt+gZfRq6cQzDMHXErMbZf6jS0lJ07twZy5cvd6gp3KljGKZ+Mduhs5S/zjQ7duzAyJEjERMTA4VCgR9++KF+2swwDGMrTtI3ABg2bBheeeUVjB492qGmcKeOYZj6Ra+3/LGDuj7FMgzDOB0r+lZUVEQ+lZWV9dYU9tQxjLOQn8qsFMQWpTRP05YDbzu7Ra6BEz11w4YNw7Bhw5zQKIZh6hNRQfNs/vr7Cw3UkpuAFU9dXFwcWTxv3jzMnz+/XprCnTqGYeoXGzp1RUVFZLGPjw98fHzqu2UMwzB1x0qnLi0tDVqt1rC4PrWNX78yDFO/6IXlD64/yQYFBRk+S5YsaeBGMwzD2IgVfdNqteRTn506HqljGKZeEUIPYSazeu3ym/kkyzAM40zMaZw53atPuFPHMDcLBR0YV2g0DdSQm4wQ5pNwCvokyzCMi2LFI6xQqUis9Perz9a4FuY0zoHZryUlJTh//rwhTklJwdGjRxEaGor4+Hir+3OnjmGY+kVYqCjByYcZhnF3zGmcA/p28OBB9O/f3xA/99xzAIBJkyZh1apVVvfnTh3DMPWLXg8ozLyGsPP1RF2fYhmGYZyOOY1z4PXrHXfcAVGHh13u1DEMU784caSurk+xDMMwTseJI3V1hTt1DOMktuq/tbh+sM84usDHux5b4zoIvR7CzEidvUbiuj7FMgzjGFt131hcP8RvAolFaFB9NselMKdxPFGCYRjPgz11DMN4MjxSxzBMo0EvAAV36hiG8VDMaRx36hiG8TiEAGBuogR36hiGcXPMaVwD6JtLVZSYP38+FAoF+bRp0wYAkJeXh+nTpyMpKQm+vr6Ij4/HjBkzUFhYSI5x+fJljBgxAn5+foiMjMSsWbNQU1Nj8ny7d++GWq1Gly5djNalp6dj/PjxCAsLg6+vLzp27IiDBw+aPM6UKVOgUCjw1ltvkeV5eXkYN24ctFotgoODMXnyZJSUlNh/YRiPQOnvSz76QH/y8VSEXlj8NBZY3xhPRqHxIZ/qUD/y8WRcSd9cbqSuffv22LZtmyFWq683MSMjAxkZGXj99dfRrl07XLp0CVOmTEFGRgbWrl0LANDpdBgxYgSio6OxZ88eZGZmYuLEifDy8sLixYvJeQoKCjBx4kQMGDAAV69eJevy8/PRt29f9O/fH7/88gsiIiJw7tw5hISEGLV33bp12LdvH2JiYozWjRs3DpmZmdi6dSuqq6vx6KOP4oknnsCaNWvqfJ0Yxm0Qepgfqbv5RuKGhPWNYTwQcxrHEyWui1x0dLTR8g4dOuC7774zxImJiVi0aBHGjx+PmpoaqNVqbNmyBSdPnsS2bdsQFRWFLl26YOHChZgzZw7mz58Pb++/ZxtOmTIFDz/8MFQqFX744QdyrqVLlyIuLg4rV640LGvevLlRm9LT0zF9+nRs3rwZI0aMIOtOnTqFTZs24cCBA+jRowcA4N1338Xw4cPx+uuvmxRJhvFEhF5AmPHUNbaZrKxvDON5mNO4htA3l+vUnTt3DjExMdBoNEhOTsaSJUvMJhUtLCyEVqs1PO3u3bsXHTt2RFRUlGGbIUOG4KmnnsJff/2Frl27AgBWrlyJixcv4osvvsArr7xidNwff/wRQ4YMwf3334/ff/8dTZs2xdNPP41//OMfhm30ej0mTJiAWbNmoX379kbH2Lt3L4KDgw2CBwADBw6EUqnE/v37MXr0aJPfqbKyEpWVleQ7AkBRUZHZa8a4BzWiisRCV0lid/4d17bdlIjViEqzT6w1qK7XdrkarG+sb56KrG81NRUkdvffsSMa1xD65lKdut69e2PVqlVISkpCZmYmFixYgNtuuw0nTpxAYGAg2TY3NxcLFy7EE088YViWlZVFBA+AIc7KygJwXVT/9a9/YefOnQaxlLl48SJWrFiB5557Di+88AIOHDiAGTNmwNvbG5MmTQJw/WlXrVZjxowZJo+RlZWFyMhIskytViM0NNTQFlMsWbIECxYsMFoeFxdndh/GTSmgYVDQsgZphjMpLi5GUND1/FTe3t6Ijo7GrqyNFveJjo4mo0yeCusb61ujYi8Ng4KMHzDcEXs17mbrm0t16oYNG2b4uVOnTujduzcSEhLwzTffYPLkyYZ1RUVFGDFiBNq1a4f58+fbfHydToeHH34YCxYsQOvWrc1up9fr0aNHD4NPpWvXrjhx4gQ++OADTJo0CYcOHcLbb7+Nw4cPQ2GlyLG9zJ0715AlH7jujUlISMDly5cNf0gM42oIIVBcXExeu2k0GqSkpKCqqsrCnteFUaPR1HcTGxzWN9Y3xn1xVONutr65VKdOJjg4GK1btya1HouLizF06FAEBgZi3bp18PLyMqyLjo7GH3/8QY5RaxKOjo5GcXExDh48iCNHjmDatGkArgucEMLgWbnzzjvRpEkTtGvXjhynbdu2Bs/Lzp07kZ2dTV6b6HQ6zJw5E2+99RZSU1MRHR2N7Oxscoyamhrk5eWZ9NTU4uPjAx8fH6PlQUFB0Gq1Fq8XwzQkpv4pazSaRtFhcwTWt79hfWPcAXfQOJdKaSJTUlKCCxcuoEmTJgCuP8EOHjwY3t7e+PHHH40uZHJyMo4fP07EZuvWrdBqtWjXrh20Wi2OHz+Oo0ePGj5TpkxBUlISjh49it69ewMA+vbtizNnzpBjnz17FgkJCQCACRMm4NixY+Q4MTExmDVrFjZv3mxoS0FBAQ4dOmQ4xn//+1/o9XrDeRiGabywvjEM43SECzFz5kyxfft2kZKSInbv3i0GDhwowsPDRXZ2tigsLBS9e/cWHTt2FOfPnxeZmZmGT01NjRBCiJqaGtGhQwcxePBgcfToUbFp0yYREREh5s6da/ac8+bNE507dybL/vjjD6FWq8WiRYvEuXPnxJdffin8/PzEF198YfY4CQkJ4s033yTLhg4dKrp27Sr2798vdu3aJVq1aiXGjh1r1zUpLCwUAERhYaFd+zEM41qwvhnD+sYwzsWlOnUPPvigaNKkifD29hZNmzYVDz74oDh//rwQQojffvuttria0SclJcVwjNTUVDFs2DDh6+srwsPDxcyZM0V1dbXZc5oSPSGE+Omnn0SHDh2Ej4+PaNOmjfjoo48stt2U6F27dk2MHTtWBAQECK1WKx599FFRXFxs8/UQQoiKigoxb948UVFRYdd+DMO4FqxvxrC+MYxzUQjRyBJFMQzDMAzDeCAu7aljGIZhGIZhbIM7dQzDMAzDMB4Ad+oYhmEYhmE8AO7UMQzDMAzDeADcqasDK1asQKdOnaDVaqHVapGcnIxffvnFsP6OO+6AQqEgnylTppg81rVr1xAbGwuFQoGCggLD8l27dqFv374ICwuDr68v2rRpgzfffNNsm1599VUoFAr885//JMsrKiowdepUhIWFISAgAPfdd58hcSnDMIwpWOMYxr1w6YoSrk5sbCxeffVVtGrVCkIIrF69GnfffTeOHDliKIL9j3/8Ay+//LJhHz8/P5PHmjx5Mjp16oT09HSy3N/fH9OmTUOnTp3g7++PXbt24cknn4S/vz+pCwkABw4cwIcffohOnToZHf/ZZ5/Fzz//jG+//RZBQUGYNm0a7r33Xuzevbuul4FhGA+FNY5h3IyGzajieYSEhIhPPvlECCFEv379xDPPPGN1n/fff1/069dP/PrrrwKAyM/Pt7j96NGjxfjx48my4uJi0apVK7F161aj8xYUFAgvLy/x7bffGpadOnVKABB79+61+bsxDMOwxjGM68KvX52ETqfDV199hdLSUiQnJxuWf/nllwgPD0eHDh0wd+5clJWVkf1OnjyJl19+GZ999hmUSuu/jiNHjmDPnj3o168fWT516lSMGDECAwcONNrn0KFDqK6uJuvatGmD+Ph47N27196vyjBMI4Q1jmFcH379WkeOHz+O5ORkVFRUICAgAOvWrTMUy3744YeRkJCAmJgYHDt2DHPmzMGZM2fw/fffAwAqKysxduxYLFu2DPHx8bh48aLZ88TGxiInJwc1NTWYP38+Hn/8ccO6r776CocPH8aBAwdM7puVlQVvb28EBweT5VFRUcjKyqrjFWAYxpNhjWMY94E7dXWktlh2YWEh1q5di0mTJuH3339Hu3btiB+kY8eOaNKkCQYMGIALFy4gMTERc+fORdu2bTF+/Hir59m5cydKSkqwb98+/Otf/0LLli0xduxYpKWl4ZlnnsHWrVuNCoAzDMPUFdY4hnEjGvr9r6cxYMAA8cQTT5hcV1JSIgCITZs2CSGE6Ny5s1AqlUKlUgmVSiWUSqUAIFQqlXjppZfMnmPhwoWidevWQggh1q1bZ9in9gNAKBQKoVKpRE1NjVkfS3x8vHjjjTec88UZhmkUsMYxjOvCI3VORq/Xo7Ky0uS6o0ePAgCaNGkCAPjuu+9QXl5uWH/gwAE89thj2LlzJxITE206x4ABA3D8+HGy/tFHH0WbNm0wZ84cqFQqdO/eHV5eXvj1119x3333AQDOnDmDy5cvE28MwzCMNVjjGMZ14U5dHZg7dy6GDRuG+Ph4FBcXY82aNdi+fTs2b96MCxcuYM2aNRg+fDjCwsJw7NgxPPvss7j99tsN0/FlUcvNzQUAtG3b1uANWb58OeLj49GmTRsAwI4dO/D6669jxowZAIDAwEB06NCBHMff3x9hYWGG5UFBQZg8eTKee+45hIaGQqvVYvr06UhOTsYtt9xSb9eHYRj3hjWOYdwL7tTVgezsbEycOBGZmZkICgpCp06dsHnzZgwaNAhpaWnYtm0b3nrrLZSWliIuLg733XcfXnzxRbvOodfrMXfuXKSkpECtViMxMRFLly7Fk08+addx3nzzTSiVStx3332orKzEkCFD8P7779t1DIZhGhescQzjXiiEEKKhG8EwDMMwDMPUDc5TxzAMwzAM4wFwp45hGIZhGMYD4E4dwzAMwzCMB8CdOoZhGIZhGA+AO3UMwzAMwzAeAHfqGIZhGIZhPADu1DEMwzAMw3gA3KljGIZhGIbxALhT1whRKBSGz+uvv+604xYUFDh07G+++QahoaEoKSlxWluY+mHTpk0ICAhATk5OQzeFYUzC+sY4iifoG3fq3Izx48dDo9Hg7NmzRuteffVVKBQKbNiwwepxRo8ejc8//xwjRoxwWtv8/f3x+eef480337R5H51Oh3nz5mH69OkICAgwLG/WrBkR0Bs/Q4cONWw3f/58ss7Pzw/t2rXDiy++iKKiInKu48ePY8yYMUhISIBGo0HTpk0xaNAgvPvuu0bt0uv1+OyzzzBo0CCEh4fDy8sLkZGRGDx4MD766COzBc1vJCMjA+PHj0dSUhICAwMRHByMXr16YfXq1ZALuXz//fd48MEH0aJFC/j5+SEpKQkzZ85EQUGByWMXFxdj9uzZaN68OXx8fNC0aVOMGTMGZWVlhm127NiBUaNGIS4uDhqNBtHR0Rg6dCh2795NjlVWVobly5dj8ODBaNKkCQIDA9G1a1esWLECOp2ObDt06FC0bNkSS5Yssfr9GcZeWN9Y3wDWt7rAZcLcjOzsbLRp0wZdunTBf//7X8PylJQUtG/fHsOHD8fatWstHkOhUGDevHmYP39+vbQxNTUVzZs3x7Jly/D8889b3PaHH37Avffei7S0NDRt2tSwvFmzZggJCcHMmTON9omJicGdd94J4LroLViwACtWrEBAQABKSkqwZcsWrFu3DsnJydi9ezcUCgX27NmD/v37Iz4+HpMmTUJ0dDTS0tKwb98+XLhwAefPnzccv7y8HKNHj8bmzZvRp08fjBw5ElFRUcjLy8Pvv/+OjRs3YtKkSfj0008tfrdjx45hxowZ6Nu3L+Lj41FdXY2tW7fixx9/xNy5c7F48WLDtuHh4YiJicE999yD+Ph4HD9+HB988AFatGiBw4cPw9fX17BtYWEh+vXrhytXruCJJ55Ay5YtkZOTg507d+Lzzz9HSEgIAOCTTz7Bhg0b0LNnT0RHRyM/Px9ffPEFjh8/jp9//tnwz+PEiRPo1KkTBgwYgMGDB0Or1WLz5s1Yt24dJk6ciNWrV5PvtWLFCjz//PPIyspCYGCgxWvAMPbA+sb6xvpWRwTjdnz00UcCgFi1apVh2dChQ4VWqxVXrlyxuj8AMW/evHprX0pKigAgli1bZnXbUaNGiVtvvdVoeUJCghgxYoTV/efNmycAiJycHLL83nvvFQDEnj17hBBCDB8+XERERIj8/HyjY1y9epXETz75pAAg3nrrLZPnPHv2rFi+fLnVtpnjrrvuEv7+/qKmpsaw7LfffjPabvXq1QKA+Pjjj8nyp556SgQHB4uLFy/afe7S0lIRFRUlhgwZYliWk5MjTpw4YbTto48+KgCIc+fOkeVXr14VKpVKfPrpp3afn2Gswfr2N6xv9sH6JgS/fnVDHn/8cfTt2xfPP/88rl27hq+++gqbNm3CK6+8Qp4G7WXVqlVQKBTYtWsXZsyYgYiICAQHB+PJJ59EVVUVCgoKMHHiRISEhCAkJASzZ882Gma3h4qKCmzatAkDBw50+BjmqH3STUlJAQBcuHAB7du3R3BwsNG2kZGRhp/T0tLwySefYOjQoXjmmWdMHrtVq1Z4+umnybLMzEycPn0a1dXVVtvWrFkzlJWVoaqqyrDsjjvuMNpu9OjRAIBTp04ZlhUUFGDlypV44okn0Lx5c1RVVdn0qqQWPz8/REREkNce4eHhaN++vU3nB65fr06dOmH9+vU2n5dhbIX1zTqsb6ZhfWNPnVuiUCjw4YcforCwEE899RSeffZZ9OjRA1OnTnXK8adPn45z585hwYIFGDVqFD766CP85z//wciRI6HT6bB48WLceuutWLZsGT7//HOHz3Po0CFUVVWhW7duJtdXV1cjNzfX6FNeXm712BcuXAAAhIWFAQASEhJw6NAhnDhxwuJ+v/zyC3Q6HcaPH2/Xd5k7dy7atm2L9PR0o3Xl5eXIzc1FamoqVq9ejZUrVyI5OZm8cjBFVlYWgOuiVMuuXbtQUVGBli1bYsyYMfDz84Ovry/69u2Lo0ePmjxOUVERcnNzcfr0abzwwgs4ceIEBgwYYPU7mTp/Ld27d8eePXusHoNh7IX1jfWN9a0ONPRQIeM4c+fOFQCESqUShw4dsnk/mHk9sXLlSgFADBkyROj1esPy5ORkoVAoxJQpUwzLampqRGxsrOjXr5/RcWx9PfHJJ58IAOL48eNG6xISEgQAk58lS5YYtqt9PXHmzBmRk5MjUlJSxIcffih8fHxEVFSUKC0tFUIIsWXLFqFSqYRKpRLJycli9uzZYvPmzaKqqoqc99lnnxUAxNGjR8nyyspKkZOTY/jk5uaS9ZMmTRIAREpKitF3WbJkCWn/gAEDxOXLly1eGyGEmDx5slCpVOLs2bOGZW+88YYAIMLCwkSvXr3El19+Kd5//30RFRUlQkJCREZGhtFxhgwZYji3t7e3ePLJJ0V5ebnFc1dWVop27dqJ5s2bi+rqaqP1ixcvFgCMXu0wjLNgfWN9Y32zH3V9dRaZ+qf2CSMmJgYdOnRw2nEnT54MhUJhiHv37o29e/di8uTJhmUqlQo9evTAoUOHHD7PtWvXAMBgfJXp3bs3XnnlFaPlrVq1MlqWlJRE4vbt22P16tXw8/MDAAwaNAh79+7FkiVLsHnzZuzduxevvfYaIiIi8Mknn2DUqFEAYJhRduNMNQDYuHGjYbgeuD4T7sYUBatWrcKqVatMfo+xY8eiR48eyMnJwYYNG3D16lWrT+Nr1qzBp59+itmzZ5PvW3tOhUKBX3/91dDOrl27Ijk5GcuXLze6Zq+++ipmzpyJtLQ0rF69GlVVVaipqbF4/mnTpuHkyZP4+eefoVYby0Tt7yw3N5e83mEYZ8H69jesb6xvtsKdOjclLS0N8+bNQ4cOHXDixAm89tprePHFF51y7Pj4eBIHBQUBAOLi4oyW5+fn1/l8woxvJTw83GY/ynfffQetVgsvLy/ExsYiMTHRaJuePXvi+++/R1VVFf7880+sW7cOb775JsaMGYOjR4+iXbt2htlOck6pvn37YuvWrQCAZcuWGU2bt0RCQgISEhIAXBfAJ554AgMHDsSZM2dMvqLYuXMnJk+ejCFDhmDRokVkXe32I0eOJMJ8yy23oHnz5iZfGXTp0sXw8/jx49GtWzc88sgjZmcRLlu2DB9//DEWLlyI4cOHm9ym9nd24z9HhnEWrG8U1jfWN1thT52bMm3aNADXPRL3338/Fi1ahIsXLzrl2CqVyubl5gTLFmr9IM4Qzttvvx0DBw5Ev379TArejXh7e6Nnz55YvHgxVqxYgerqanz77bcAgDZt2gCAkTclIiICAwcOxMCBA9GkSZM6tXXMmDFIS0vDjh07jNb9+eefGDVqFDp06IC1a9caPUXGxMQAAKKiooz2jYyMtHotvb29MWrUKHz//fcmn6ZXrVqFOXPmYMqUKRb/idaex5QfhWHqCusbhfWN9c1WuFPnhqxbtw4//vgjFi5ciNjYWLz11lvw9vZ2mpH4ZlErMLUzuBqCHj16ALg+uwsAhg0bBpVKhS+//LLezlkrNoWFhWT5hQsXMHToUERGRmLjxo1Gr0iA6wZeACYNyxkZGYiIiLDp/EIIFBcXk+Xr16/H448/jnvvvRfLly+3eIyUlBSEh4fbdD6GsQfWN+fB+vY3jUXfuFPnZhQXF2PGjBno2rUrpk+fDuD6083ChQuxadMmwxOZO9C9e3d4e3vj4MGD9X6u3377zeRT98aNGwH87VmJj4/HY489hl9++QXvvfeeyWOZOo6pKf/mSs18+umnUCgUZFZcVlYWBg8eDKVSic2bN5sVk6SkJHTu3Bnr169Hbm6uYfmWLVuQlpaGQYMGGZZlZ2cb7V9QUIDvvvsOcXFxxCuyY8cOPPTQQ7j99tvx5ZdfQqm0LA2HDh1CcnKyxW0Yxl5Y3xyD9e06rG/sqXM7XnzxRWRkZOD7778nrwumTp2K1atX45///CeGDh3qFpmwNRoNBg8ejG3btuHll182Wp+eno4vvvjCaHlAQADuueceu841ffp0lJWVYfTo0WjTpg2qqqqwZ88efP3112jWrBkeffRRw7ZvvfUWUlJSMH36dHz11VcYOXIkIiMjkZubi927d+Onn34yMi7PnTsXq1evRkpKCpo1awYAWLRoEXbv3o2hQ4ciPj4eeXl5+O6773DgwAFMnz4dLVu2NOw/dOhQXLx4EbNnz8auXbuwa9cuw7qoqCgiZm+++SYGDRqEW2+9FU8++SQKCwvxxhtvoHXr1njqqacM2w0bNgyxsbHo3bs3IiMjcfnyZaxcuRIZGRn4+uuvDdtdunQJo0aNgkKhwJgxY4z+cXbq1AmdOnUyxNnZ2Th27JjbjZwwrg/rG+sb61sdaZhJt4wjHDx4UKhUKjFt2jST6//44w+hVCrFjBkzLB4HVqb8HzhwgCw3l9V80qRJwt/f3+g49mRc//7774VCoTCaAm9pyn9CQoLVtsn88ssv4rHHHhNt2rQRAQEBwtvbW7Rs2VJMnz7d5LT1mpoasXLlSnHnnXeK0NBQoVarRXh4uBgwYID44IMPjKbMm5ryv2XLFnHXXXeJmJgY4eXlJQIDA0Xfvn3FypUrSUoFIYTZ7wrAZFqFrVu3iltuuUVoNBoRGhoqJkyYIDIzM8k27733nrj11ltFeHi4UKvVIiIiQowcOVLs2LGDbPfbb79ZPL/8t7JixQrh5+cnioqKLF5zhrEH1jfWt1pY3xyHO3WNEABi1qxZIicnR5SVlTntuHq9XuTk5IjDhw/bLHo1NTWidevW4sUXX3RaO5j6pUuXLuKf//xnQzeDYUzC+sbUBXfXN/bUNVKWLVuGiIgIq4ZReygsLERERITZDOqmUKlUePnll7F8+XKjafaM67Fp0yacO3cOc+fObeimMIxZWN8YR/AEfVMIUYc524xbsm3bNsPPrVu3Nsrb5Cg1NTXYvn17vRybYRjGFljfmMYMd+oYhmEYhmE8AH79yjAMwzAM4wFwp45hGIZhGMYD4E4dwzAMwzCMB8CdOoZhGIZhGA+AO3UMwzAMwzAeAHfqGIZhGIZhPADu1DEMwzAMw3gA3KljGIZhGIbxALhTxzAMwzAM4wH8fyXb7vcBtAFBAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scoll = SceneCollection()\n",
    "trait_selected = 'lai'\n",
    "for fpath_trait in output_dir.glob('*_traits.tiff'):\n",
    "    sensing_date = datetime.strptime(\n",
    "        fpath_trait.name.split('_')[1], '%Y-%m-%d')\n",
    "    trait_ds = RasterCollection.from_multi_band_raster(\n",
    "        fpath_trait, band_names_src=[trait_selected])\n",
    "    # check if the trait is masked\n",
    "    if trait_ds[trait_selected].is_masked_array:\n",
    "        if trait_ds[trait_selected].values.mask.all():\n",
    "            continue\n",
    "    else:\n",
    "        if np.isnan(trait_ds[trait_selected].values).all():\n",
    "            continue\n",
    "    scene_properties = SceneProperties(\n",
    "        acquisition_time=sensing_date)\n",
    "    trait_ds.scene_properties = scene_properties\n",
    "    scoll.add_scene(trait_ds)\n",
    "\n",
    "# sort scenes by acquisition time\n",
    "scoll = scoll.sort()\n",
    "\n",
    "f = scoll.plot(max_scenes_in_row=len(scoll), band_selection=['lai'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The GLAI was generated for two scenes, only, as on the 9th of May all pixel values have been masked out due to clouds. The GLAI ranges between typical values for winter wheat at this time of the year in Switzerland, when the plants are close to the heading date. Around heading, we expect to see the highest GLAI values in winter wheat.\n",
    "\n",
    "The analysis could be extended using the same methodology to the entire growing period of winter wheat (in this case, between October 2021 and end of July 2022). The most time consuming factor is the generation of the lookup-tables."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "GeoPython",
   "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}