UC3 - Crop type feature engineering using random forest.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b48bd388",
   "metadata": {},
   "source": [
    "<img src=\"https://avatars.githubusercontent.com/u/74911464?s=200&v=4\"\n",
    "     alt=\"OpenEO Platform logo\"\n",
    "     style=\"float: left; margin-right: 10px;\" />\n",
    "# OpenEO Platform - UC3\n",
    "Crop type feature engineering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "947f7d3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import openeo\n",
    "from openeo.processes import ProcessBuilder, array_modify, quantiles, sd, array_concat\n",
    "from eo_utils import *\n",
    "from shapely.geometry import box\n",
    "import pandas as pd\n",
    "from helper import *\n",
    "import hvplot.xarray\n",
    "from sklearn.model_selection import train_test_split, GridSearchCV\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.metrics import accuracy_score\n",
    "import pickle\n",
    "import urllib.request"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97c842e9",
   "metadata": {},
   "source": [
    "## Objectives and approach\n",
    "\n",
    "The goal of this section is to show how openEO functionality can be integrated into a basic feature engineering pipeline. We will do this using crop classification as an example.\n",
    "\n",
    "The constantly increasing demand of food has resulted in a highly intensified agricultural production. This intensification on the one hand requires more planning and management, and on the other, threatens ecosystem services that need to be monitored by scientists and decision makers who rely on detailed spatial information of crop cover in agricultural areas.\n",
    "\n",
    "Crop classification on a large scale is a challenging task, but with the recent advances in satellite sensor technology and the push of a.o. the ESA for higher resolution open satellite data with a frequent revisit time this task has become possible.\n",
    "\n",
    "There are various approaches to crop classification. One can use basic rule-based classification as shown in the rule-based notebook, or use more sophisticated methods such as one of various machine learning models. In this example, we will have a look at this machine learning approach and more specifically, look at implementing a random forest model.\n",
    "\n",
    "Generally, any classification task will contain the following steps:\n",
    "\n",
    "* Preprocessing & feature engineering\n",
    "* Training\n",
    "* Classification & model evaluation\n",
    "\n",
    "We will have a more detailed look at all three of these steps, and provide code examples along the way."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43902a21",
   "metadata": {},
   "source": [
    "First, make a connection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "24531079",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Authenticated using refresh token.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Connection to 'https://openeo.vito.be/openeo/1.1/' with OidcBearerAuth>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "connection = openeo.connect(\"openeo.vito.be\")\n",
    "connection.authenticate_oidc()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67bfd1a1",
   "metadata": {},
   "source": [
    "Next, we'll set some variables we will need later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dc306130",
   "metadata": {},
   "outputs": [],
   "source": [
    "year = 2020  # The selected worldcover connection will only work for 2020\n",
    "temp_ext_s2 = [str(year-1)+\"-09-01\", str(year+1)+\"-04-30\"]\n",
    "temp_ext_s1 = [str(year)+\"-01-01\",str(year)+\"-12-31\"]\n",
    "temp_ext_wc = str(year)\n",
    "\n",
    "job_options = {\n",
    "        \"driver-memory\": \"2G\",\n",
    "        \"driver-memoryOverhead\": \"1G\",\n",
    "        \"driver-cores\": \"2\",\n",
    "        \"executor-memory\": \"2G\",\n",
    "        \"executor-memoryOverhead\": \"2G\",\n",
    "        \"executor-cores\": \"1\",\n",
    "        \"max-executors\": \"100\"\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b06c447",
   "metadata": {},
   "source": [
    "## Load the dataset\n",
    "\n",
    "Connect to the vito backend and load a large number of S2 bands (they will be used to calculate a variety of vegetation indices) and the scene classification layer (SCL) for the year of 2020.\n",
    "We will use the scene classification band to perform cloud masking and subsequently drop the SCL layer as we will not need it anymore in the rest of this notebook.\n",
    "\n",
    "We use a temporal range that is slightly bigger than the actual temporal range we need to improve the interpolation that we will do later on. For the S1 imagery, no interpolation is needed as no cloud masking needs to be performed, hence the temporal extent is kept to the range we need to avoid loading in unnecessary data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4f71234f",
   "metadata": {},
   "outputs": [],
   "source": [
    "s2 = connection.load_collection(\"TERRASCOPE_S2_TOC_V2\",\n",
    "                                temporal_extent=temp_ext_s2,\n",
    "                                bands=[\"B03\",\"B04\",\"B05\",\"B06\",\"B07\",\"B08\",\"B11\",\"B12\",\"SCL\"])\n",
    "s2 = s2.process(\"mask_scl_dilation\", data=s2, scl_band_name=\"SCL\").filter_bands(s2.metadata.band_names[:-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "da8c62e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = connection.load_collection(\"S1_GRD_SIGMA0_ASCENDING\",\n",
    "                                temporal_extent=temp_ext_s1,\n",
    "                                bands=[\"VH\",\"VV\"]\n",
    "                               )\n",
    "s1 = s1.apply_dimension(dimension=\"bands\",process=lambda x: array_modify(data=x, values=x.array_element(0)/x.array_element(1), index=0))\n",
    "s1 = s1.rename_labels(\"bands\", [\"ratio\"] + s1.metadata.band_names)\n",
    "#TODO: this scaling does not seem appropriate for sentinel-1\n",
    "s1 = s1.apply(process=lambda x: x.linear_scale_range(0,1,0,250))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f520871c",
   "metadata": {},
   "source": [
    "In addition, we want to load in a crop mask to ensure that we only classify areas that are actually cropland."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "bfa2c98e",
   "metadata": {},
   "outputs": [],
   "source": [
    "cropland_mask = connection.load_collection(\"ESA_WORLDCOVER_10M_2020_V1\",\n",
    "                                temporal_extent=temp_ext_wc,\n",
    "                                bands=[\"MAP\"]\n",
    "                    )\n",
    "cropland_mask = cropland_mask.band(\"MAP\") != 40\n",
    "\n",
    "s2 = s2.mask(cropland_mask.max_time().resample_cube_spatial(s2))\n",
    "s1 = s1.mask(cropland_mask.max_time().resample_cube_spatial(s1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fc10ebc",
   "metadata": {},
   "source": [
    "## Feature selection for crop type mapping"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f05f3e98",
   "metadata": {},
   "source": [
    "For this use case, the following crop types will be mapped:\n",
    "\n",
    "* Wheat\n",
    "* Maize \n",
    "* Potato \n",
    "* Sugar beet\n",
    "* Soy\n",
    "\n",
    "A separability analysis was conducted based on the reference data currently available in the WorldCereal Common Input Baseline. We only focused on European reference data (i.e. LPIS and SIGPAC data). All classification features currently supported by the WorldCereal system were tested. Separability was first computed for each crop type separately in a binary one-vs-all case. Afterwards, an average separability of each feature was computed and features were ranked accordingly. Based on the ranking of features for individual crop types, the overall ranking for all types and ease of implementation, a proposal is made of features to be included for the openEO CCN."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0705e52e",
   "metadata": {},
   "source": [
    "**Sentinel-2 features**\n",
    "For the final feature selection, the following timeseries should be computed/derived (see satio.rsindices for definitions of the band combinations mentioned below) from Sentinel-2:\n",
    "* bands B06 and B12\n",
    "* vegetation indices NDVI, NDMI, NDGI, ANIR, NDRE1, NDRE2 and NDRE5\n",
    "\n",
    "**Sentinel-1 features**\n",
    "* bands VV and VH\n",
    "* VH/VV ratio\n",
    "\n",
    "Based on all these timeseries, the following features should be computed: \n",
    "* Percentiles (p10, p50, p90) \n",
    "* Standard deviation \n",
    "* Tsteps. This results in a set of 6 features. Ts0, ts1, ts2, ts3, ts4 and ts5, of which ts0, ts1 and ts5 are the most important ones. The others can be dropped to limit the total number of features to be used. \n",
    "\n",
    "**Other features**\n",
    "* Altitude derived from Copernicus 30m DEM\n",
    "* Derived from AgERA5 meteo data: Sum of precipitation over growing season\n",
    "* P10, p50 and p90 of mean temperature \n",
    "\n",
    "In the end, we decided to stick with the Sentinel-1 and Sentinel-2 features (percentiles, std. dev. and t-steps), and leave the other features out for now as the accuracy was already very decent. With that we got to an accuracy of about 97% on the test set, but bear in mind that the test set was made up of samples from fields known to be one of the 5 crop types, and not any other crop classes / land cover classes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "228851e8",
   "metadata": {},
   "source": [
    "## Feature computation\n",
    "\n",
    "Features are generally computed over the temporal dimension, and result in multiple bands being added. In the openEO API, the 'apply_dimension' process allows us to work on the temporal dimension, and store the result in the bands dimension."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "905558ab",
   "metadata": {},
   "source": [
    "We calculate a number of vegetation indices, aggregate to 10-day composites, interpolate linearly to fill up any potential missing data values and crop back to the the final temporal range (01-01-2020 to 31-12-2020). Finally, we merge all the cubes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3ddf596a",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_list = [\"NDVI\",\"NDMI\",\"NDGI\",\"ANIR\",\"NDRE1\",\"NDRE2\",\"NDRE5\"]\n",
    "s2_list = [\"B06\",\"B12\"]\n",
    "\n",
    "s1_dekad = s1.aggregate_temporal_period(period=\"dekad\", reducer=\"mean\")\n",
    "s1_dekad = s1_dekad.apply_dimension(dimension=\"t\", process=\"array_interpolate_linear\")\n",
    "#s1_dekad = s1_dekad.resample_cube_spatial(s2)\n",
    "\n",
    "indices = compute_indices(s2, idx_list).filter_bands(s2_list+idx_list)\n",
    "idx_dekad = indices.aggregate_temporal_period(period=\"dekad\", reducer=\"mean\")\n",
    "idx_dekad = idx_dekad.apply_dimension(dimension=\"t\", process=\"array_interpolate_linear\").filter_temporal([str(year)+\"-01-01\", str(year)+\"-12-31\"])\n",
    "\n",
    "base_features = idx_dekad.merge_cubes(s1_dekad)\n",
    "base_features = base_features.rename_labels(\"bands\",s2_list+idx_list+[\"ratio\",\"VV\",\"VH\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cf72f88",
   "metadata": {},
   "source": [
    "From these temporal cubes we will calculate the features that we need: quantiles, standard deviation and t-steps of all indices and band 6 and 12. We use linear scaling to ensure that all values are int8s for the sake of computational efficiency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b88db03f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def computeStats(input_timeseries:ProcessBuilder):\n",
    "    tsteps = list([input_timeseries.array_element(6*index) for index in range(0,6)])\n",
    "    return array_concat(array_concat(input_timeseries.quantiles(probabilities=[0.1,0.5,0.9]),input_timeseries.sd()),tsteps)\n",
    "\n",
    "features = base_features.apply_dimension(dimension='t',target_dimension='bands', process=computeStats).apply(lambda x: x.linear_scale_range(0,250,0,250))\n",
    "\n",
    "tstep_labels = [ \"t\"+ str(6*index) for index in range(0,6) ]\n",
    "all_bands = [band + \"_\" + stat for band in base_features.metadata.band_names for stat in [\"p10\",\"p50\",\"p90\",\"sd\"] + tstep_labels ]\n",
    "features = features.rename_labels('bands',all_bands)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8bb061d",
   "metadata": {},
   "source": [
    "## Dataset sampling\n",
    "\n",
    "Once we have our features defined, we will want to generate a set of input features for model calibration. Instead of calculating features\n",
    "for a very large area, openEO can sample the datacube at selected locations.\n",
    "\n",
    "Sampling a datacube in openEO currently requires polygons as sampling features. Other types of geometry, like lines\n",
    "and points, should be converted into polygons first by applying a buffering operation. Using the spatial resolution\n",
    "of the datacube as buffer size can be a way to approximate sampling at a point.\n",
    "\n",
    "To indicate to openEO that we only want to compute the datacube for certain polygon features, we use the\n",
    "`~openeo.rest.datacube.DataCube.filter_spatial` method.\n",
    "\n",
    "Next to that, we will also indicate that we want to write multiple output files. This is more convenient, as we will\n",
    "want to have one or more raster outputs per sampling feature, for convenient further processing. To do this, we set\n",
    "the 'sample_by_feature' output format property, which is available for the netCDF and GTiff output formats.\n",
    "\n",
    "Sampling only works for batch jobs, because it results in multiple output files, which can not be conveniently transferred\n",
    "in a synchronous call.\n",
    "\n",
    "Combining all of this, results in the following sample code:\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1610f593",
   "metadata": {},
   "outputs": [],
   "source": [
    "#crops = {\"maize\":1200,\"potatos\":5100,\"sugarbeet\":8100,\"barley\":1500,\"soy\":4100}\n",
    "crops = [\"maize\", \"potatos\", \"sugarbeet\", \"barley\", \"soy\"]\n",
    "crop_samples = {name:gpd.read_file(\"UC3_resources/\"+ name + \"_2019.geojson\", driver='GeoJSON') for name in crops}\n",
    "\n",
    "points_per_type = point_sample_fields(crop_samples, 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79be9600",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for i in points_per_type.keys():\n",
    "    sampled_features = features.filter_spatial(eval(points_per_type[i]))\n",
    "    job = sampled_features.execute_batch(\n",
    "        title=\"Point feature extraction\",\n",
    "        description=\"Feature extraction for p10,p50,p90,sd and tsteps\",\n",
    "        out_format=\"netCDF\",\n",
    "        sample_by_feature=True,\n",
    "        job_options=job_options)\n",
    "    results = job.get_results()\n",
    "    results.download_files(\"./data/rf_300_\"+i)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f73029e7",
   "metadata": {},
   "source": [
    "We can then load this data into a dataframe which can be used to train the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "264c0833",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: total: 45 s\n",
      "Wall time: 45.2 s\n"
     ]
    },
    {
     "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>Crop type</th>\n",
       "      <th>Iteration nr</th>\n",
       "      <th>B06_p10</th>\n",
       "      <th>B06_p50</th>\n",
       "      <th>B06_p90</th>\n",
       "      <th>B06_sd</th>\n",
       "      <th>B06_t0</th>\n",
       "      <th>B06_t6</th>\n",
       "      <th>B06_t12</th>\n",
       "      <th>B06_t18</th>\n",
       "      <th>...</th>\n",
       "      <th>VH_p10</th>\n",
       "      <th>VH_p50</th>\n",
       "      <th>VH_p90</th>\n",
       "      <th>VH_sd</th>\n",
       "      <th>VH_t0</th>\n",
       "      <th>VH_t6</th>\n",
       "      <th>VH_t12</th>\n",
       "      <th>VH_t18</th>\n",
       "      <th>VH_t24</th>\n",
       "      <th>VH_t30</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>barley</td>\n",
       "      <td>openEO_0</td>\n",
       "      <td>46.889999</td>\n",
       "      <td>53.099998</td>\n",
       "      <td>84.099998</td>\n",
       "      <td>12.156557</td>\n",
       "      <td>62.50</td>\n",
       "      <td>53.200001</td>\n",
       "      <td>79.5</td>\n",
       "      <td>45.0</td>\n",
       "      <td>...</td>\n",
       "      <td>13.7</td>\n",
       "      <td>26.875</td>\n",
       "      <td>52.049999</td>\n",
       "      <td>17.337147</td>\n",
       "      <td>34.0</td>\n",
       "      <td>75.5</td>\n",
       "      <td>29.0</td>\n",
       "      <td>23.0</td>\n",
       "      <td>20.25</td>\n",
       "      <td>18.75</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>barley</td>\n",
       "      <td>openEO_1</td>\n",
       "      <td>42.750000</td>\n",
       "      <td>66.777779</td>\n",
       "      <td>99.811111</td>\n",
       "      <td>20.900290</td>\n",
       "      <td>39.00</td>\n",
       "      <td>61.333332</td>\n",
       "      <td>83.0</td>\n",
       "      <td>55.0</td>\n",
       "      <td>...</td>\n",
       "      <td>9.9</td>\n",
       "      <td>22.250</td>\n",
       "      <td>69.199997</td>\n",
       "      <td>20.769791</td>\n",
       "      <td>28.0</td>\n",
       "      <td>47.5</td>\n",
       "      <td>27.5</td>\n",
       "      <td>10.0</td>\n",
       "      <td>13.00</td>\n",
       "      <td>52.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>barley</td>\n",
       "      <td>openEO_10</td>\n",
       "      <td>54.466667</td>\n",
       "      <td>63.472221</td>\n",
       "      <td>88.733337</td>\n",
       "      <td>12.992085</td>\n",
       "      <td>67.75</td>\n",
       "      <td>55.333332</td>\n",
       "      <td>86.0</td>\n",
       "      <td>43.0</td>\n",
       "      <td>...</td>\n",
       "      <td>12.7</td>\n",
       "      <td>29.000</td>\n",
       "      <td>50.000000</td>\n",
       "      <td>14.174424</td>\n",
       "      <td>29.0</td>\n",
       "      <td>56.0</td>\n",
       "      <td>50.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>17.00</td>\n",
       "      <td>24.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>barley</td>\n",
       "      <td>openEO_100</td>\n",
       "      <td>55.000000</td>\n",
       "      <td>62.000000</td>\n",
       "      <td>90.099998</td>\n",
       "      <td>13.240349</td>\n",
       "      <td>61.75</td>\n",
       "      <td>57.000000</td>\n",
       "      <td>86.0</td>\n",
       "      <td>44.0</td>\n",
       "      <td>...</td>\n",
       "      <td>8.7</td>\n",
       "      <td>32.000</td>\n",
       "      <td>84.599998</td>\n",
       "      <td>31.028559</td>\n",
       "      <td>64.0</td>\n",
       "      <td>76.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>10.00</td>\n",
       "      <td>35.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>barley</td>\n",
       "      <td>openEO_101</td>\n",
       "      <td>33.000000</td>\n",
       "      <td>52.000000</td>\n",
       "      <td>88.833336</td>\n",
       "      <td>19.758512</td>\n",
       "      <td>31.75</td>\n",
       "      <td>45.000000</td>\n",
       "      <td>76.5</td>\n",
       "      <td>50.0</td>\n",
       "      <td>...</td>\n",
       "      <td>10.7</td>\n",
       "      <td>25.500</td>\n",
       "      <td>60.000000</td>\n",
       "      <td>18.545933</td>\n",
       "      <td>20.0</td>\n",
       "      <td>73.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>31.00</td>\n",
       "      <td>36.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 122 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "  Crop type Iteration nr    B06_p10    B06_p50    B06_p90     B06_sd  B06_t0  \\\n",
       "0    barley     openEO_0  46.889999  53.099998  84.099998  12.156557   62.50   \n",
       "1    barley     openEO_1  42.750000  66.777779  99.811111  20.900290   39.00   \n",
       "2    barley    openEO_10  54.466667  63.472221  88.733337  12.992085   67.75   \n",
       "3    barley   openEO_100  55.000000  62.000000  90.099998  13.240349   61.75   \n",
       "4    barley   openEO_101  33.000000  52.000000  88.833336  19.758512   31.75   \n",
       "\n",
       "      B06_t6  B06_t12  B06_t18  ...  VH_p10  VH_p50     VH_p90      VH_sd  \\\n",
       "0  53.200001     79.5     45.0  ...    13.7  26.875  52.049999  17.337147   \n",
       "1  61.333332     83.0     55.0  ...     9.9  22.250  69.199997  20.769791   \n",
       "2  55.333332     86.0     43.0  ...    12.7  29.000  50.000000  14.174424   \n",
       "3  57.000000     86.0     44.0  ...     8.7  32.000  84.599998  31.028559   \n",
       "4  45.000000     76.5     50.0  ...    10.7  25.500  60.000000  18.545933   \n",
       "\n",
       "   VH_t0  VH_t6  VH_t12  VH_t18  VH_t24  VH_t30  \n",
       "0   34.0   75.5    29.0    23.0   20.25   18.75  \n",
       "1   28.0   47.5    27.5    10.0   13.00   52.00  \n",
       "2   29.0   56.0    50.0    16.0   17.00   24.00  \n",
       "3   64.0   76.0    60.0    16.0   10.00   35.00  \n",
       "4   20.0   73.0    18.0    17.0   31.00   36.00  \n",
       "\n",
       "[5 rows x 122 columns]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time df = prep_df(2020,all_bands)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cdf7ad8",
   "metadata": {},
   "source": [
    "## Training a random forest model\n",
    "Next, we'll start training a model. We first split our dataset into a training and test set, and then perform grid search to find the optimal parameters for our model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "548dbe6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = df[all_bands]\n",
    "y = df[\"Crop type\"]\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3bd1aef5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting 3 folds for each of 48 candidates, totalling 144 fits\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'max_depth': 100,\n",
       " 'max_features': 'log2',\n",
       " 'min_samples_leaf': 1,\n",
       " 'min_samples_split': 2,\n",
       " 'n_estimators': 300}"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "param_grid = {\n",
    "    'max_depth': [100, None],\n",
    "    'max_features': [4, 'log2'],\n",
    "    'min_samples_leaf': [1, 2],\n",
    "    'min_samples_split': [2, 3],\n",
    "    'n_estimators': [100, 200, 300]\n",
    "}\n",
    "rf = RandomForestClassifier()\n",
    "grid_search = GridSearchCV(estimator = rf, param_grid = param_grid, cv = 3, n_jobs = -1, verbose = 2)\n",
    "\n",
    "grid_search.fit(X_train, y_train)\n",
    "grid_search.best_params_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f5e46907",
   "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>mean_fit_time</th>\n",
       "      <th>std_fit_time</th>\n",
       "      <th>mean_score_time</th>\n",
       "      <th>std_score_time</th>\n",
       "      <th>param_max_depth</th>\n",
       "      <th>param_max_features</th>\n",
       "      <th>param_min_samples_leaf</th>\n",
       "      <th>param_min_samples_split</th>\n",
       "      <th>param_n_estimators</th>\n",
       "      <th>params</th>\n",
       "      <th>split0_test_score</th>\n",
       "      <th>split1_test_score</th>\n",
       "      <th>split2_test_score</th>\n",
       "      <th>mean_test_score</th>\n",
       "      <th>std_test_score</th>\n",
       "      <th>rank_test_score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>1.704370</td>\n",
       "      <td>0.023810</td>\n",
       "      <td>0.048929</td>\n",
       "      <td>0.003428</td>\n",
       "      <td>100</td>\n",
       "      <td>log2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>300</td>\n",
       "      <td>{'max_depth': 100, 'max_features': 'log2', 'mi...</td>\n",
       "      <td>0.982639</td>\n",
       "      <td>0.979167</td>\n",
       "      <td>0.993031</td>\n",
       "      <td>0.984946</td>\n",
       "      <td>0.005891</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.481911</td>\n",
       "      <td>0.005506</td>\n",
       "      <td>0.025906</td>\n",
       "      <td>0.003910</td>\n",
       "      <td>100</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>100</td>\n",
       "      <td>{'max_depth': 100, 'max_features': 4, 'min_sam...</td>\n",
       "      <td>0.986111</td>\n",
       "      <td>0.979167</td>\n",
       "      <td>0.989547</td>\n",
       "      <td>0.984942</td>\n",
       "      <td>0.004318</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>0.906782</td>\n",
       "      <td>0.040823</td>\n",
       "      <td>0.029147</td>\n",
       "      <td>0.004681</td>\n",
       "      <td>None</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>200</td>\n",
       "      <td>{'max_depth': None, 'max_features': 4, 'min_sa...</td>\n",
       "      <td>0.979167</td>\n",
       "      <td>0.986111</td>\n",
       "      <td>0.986063</td>\n",
       "      <td>0.983780</td>\n",
       "      <td>0.003262</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.884070</td>\n",
       "      <td>0.035430</td>\n",
       "      <td>0.033424</td>\n",
       "      <td>0.005056</td>\n",
       "      <td>100</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>200</td>\n",
       "      <td>{'max_depth': 100, 'max_features': 4, 'min_sam...</td>\n",
       "      <td>0.975694</td>\n",
       "      <td>0.979167</td>\n",
       "      <td>0.993031</td>\n",
       "      <td>0.982631</td>\n",
       "      <td>0.007490</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>1.361126</td>\n",
       "      <td>0.010292</td>\n",
       "      <td>0.047401</td>\n",
       "      <td>0.010305</td>\n",
       "      <td>None</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>300</td>\n",
       "      <td>{'max_depth': None, 'max_features': 4, 'min_sa...</td>\n",
       "      <td>0.982639</td>\n",
       "      <td>0.979167</td>\n",
       "      <td>0.986063</td>\n",
       "      <td>0.982623</td>\n",
       "      <td>0.002815</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    mean_fit_time  std_fit_time  mean_score_time  std_score_time  \\\n",
       "14       1.704370      0.023810         0.048929        0.003428   \n",
       "0        0.481911      0.005506         0.025906        0.003910   \n",
       "25       0.906782      0.040823         0.029147        0.004681   \n",
       "1        0.884070      0.035430         0.033424        0.005056   \n",
       "26       1.361126      0.010292         0.047401        0.010305   \n",
       "\n",
       "   param_max_depth param_max_features param_min_samples_leaf  \\\n",
       "14             100               log2                      1   \n",
       "0              100                  4                      1   \n",
       "25            None                  4                      1   \n",
       "1              100                  4                      1   \n",
       "26            None                  4                      1   \n",
       "\n",
       "   param_min_samples_split param_n_estimators  \\\n",
       "14                       2                300   \n",
       "0                        2                100   \n",
       "25                       2                200   \n",
       "1                        2                200   \n",
       "26                       2                300   \n",
       "\n",
       "                                               params  split0_test_score  \\\n",
       "14  {'max_depth': 100, 'max_features': 'log2', 'mi...           0.982639   \n",
       "0   {'max_depth': 100, 'max_features': 4, 'min_sam...           0.986111   \n",
       "25  {'max_depth': None, 'max_features': 4, 'min_sa...           0.979167   \n",
       "1   {'max_depth': 100, 'max_features': 4, 'min_sam...           0.975694   \n",
       "26  {'max_depth': None, 'max_features': 4, 'min_sa...           0.982639   \n",
       "\n",
       "    split1_test_score  split2_test_score  mean_test_score  std_test_score  \\\n",
       "14           0.979167           0.993031         0.984946        0.005891   \n",
       "0            0.979167           0.989547         0.984942        0.004318   \n",
       "25           0.986111           0.986063         0.983780        0.003262   \n",
       "1            0.979167           0.993031         0.982631        0.007490   \n",
       "26           0.979167           0.986063         0.982623        0.002815   \n",
       "\n",
       "    rank_test_score  \n",
       "14                1  \n",
       "0                 2  \n",
       "25                3  \n",
       "1                 4  \n",
       "26                5  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gs_results = pd.DataFrame(grid_search.cv_results_).sort_values(by=[\"mean_test_score\"],axis=0,ascending=False)\n",
    "gs_results.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b47edea4",
   "metadata": {},
   "source": [
    "From playing around with gridsearch we can concluse that the most important model parameters are:\n",
    "* max_features: low (around 3-5)\n",
    "* min_samples_leaf: low (1)\n",
    "* min_samples_split: low (2)\n",
    "\n",
    "And we can run the best model on our test set and calculate the accuracy (note that our test set consists of either one of the 5 crop types, so the results on a set consisting of other crop types or even land cover types will be lower."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1715dad7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy on test set: 0.972\n"
     ]
    }
   ],
   "source": [
    "y_pred = grid_search.predict(X_test)\n",
    "print(\"Accuracy on test set: \"+str(accuracy_score(y_test,y_pred))[0:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e9d3aa6",
   "metadata": {},
   "source": [
    "We can then save this model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8391b55e",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = './data/rf_model.pkl'\n",
    "pickle.dump(grid_search, open(filename, 'wb'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc578b18",
   "metadata": {},
   "source": [
    "After saving the model, we can upload it to any public repository. The model on this repository can then be accessed from the UDF that will be used for performing inference.\n",
    "\n",
    "We select an area for which we want to perform model inference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "d02640ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "44f253012dd7414c90dbb17b8af2c187",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[50.8, 4.75], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_title', 'zoom_out…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "center = [50.8, 4.75]\n",
    "zoom = 13\n",
    "\n",
    "eoMap = openeoMap(center,zoom)\n",
    "addS2Tiles(eoMap)\n",
    "eoMap.map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b82e9f2",
   "metadata": {},
   "source": [
    "We can then use the model to do inference and store the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "9c7d434d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Authenticated using refresh token.\n",
      "0:00:00 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': send 'start'\n",
      "0:00:16 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': queued (progress N/A)\n",
      "0:00:22 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': queued (progress N/A)\n",
      "0:00:28 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': queued (progress N/A)\n",
      "0:00:36 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': queued (progress N/A)\n",
      "0:00:46 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': queued (progress N/A)\n",
      "0:00:59 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': queued (progress N/A)\n",
      "0:01:14 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': queued (progress N/A)\n",
      "0:01:34 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:01:58 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:02:28 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:03:05 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:03:52 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:04:51 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:05:51 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:06:51 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:07:52 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:08:52 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:09:54 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:10:55 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:11:55 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:12:56 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:13:56 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:14:57 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:15:58 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:16:59 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': running (progress N/A)\n",
      "0:17:59 Job 'j-10aeb2573c0a4327a07d7edc5cc19cbf': finished (progress N/A)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[WindowsPath('data/test_klein_gebied/openEO.tif'),\n",
       " WindowsPath('data/test_klein_gebied/job-results.json')]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "terrascope_settings = {\n",
    "    \"SENTINEL2\":\"TERRASCOPE_S2_TOC_V2\",\n",
    "    \"SENTINEL1\":\"S1_GRD_SIGMA0_ASCENDING\",\n",
    "    \"mask\":True,\n",
    "    \"url\":\"openeo.vito.be\",\n",
    "    \"props\": None\n",
    "}\n",
    "creo_settings = {\n",
    "    \"SENTINEL2\":\"SENTINEL2_L2A\",\n",
    "    \"SENTINEL1\":\"SENTINEL1_GRD\",\n",
    "    \"mask\":False,\n",
    "    \"url\":\"openeo.creo.vito.be\",\n",
    "    \"props\": None\n",
    "}\n",
    "\n",
    "s = terrascope_settings\n",
    "connection = openeo.connect(s[\"url\"]).authenticate_oidc()\n",
    "bbox = eoMap.getBbox()\n",
    "\n",
    "udf_rf = \"\"\"\n",
    "from typing import Dict\n",
    "from openeo_udf.api.datacube import DataCube\n",
    "import pickle\n",
    "import urllib.request\n",
    "import xarray\n",
    "import sklearn\n",
    "from openeo.udf.xarraydatacube import XarrayDataCube\n",
    "import functools\n",
    "\n",
    "@functools.lru_cache(maxsize=25)\n",
    "def load_model():\n",
    "    return pickle.load(urllib.request.urlopen(\"https://artifactory.vgt.vito.be:443/auxdata-public/openeo/rf_model.pkl\"))\n",
    "\n",
    "def impute(a):\n",
    "    a[np.isnan(a)] = np.int16(np.nanmean(a))\n",
    "    return a\n",
    "\n",
    "def apply_datacube(cube: XarrayDataCube, context: dict) -> XarrayDataCube:\n",
    "    array = cube.get_array()\n",
    "    stacked_array = array.stack(pixel=(\"x\",\"y\"))\n",
    "    stacked_array = stacked_array.transpose()\n",
    "    stacked_array_f = stacked_array[~np.isnan(stacked_array).all(axis=1)]\n",
    "    \n",
    "    if len(stacked_array_f) == 0:\n",
    "        result = np.full(np.multiply(*array.shape[1:]),np.nan)\n",
    "    else:\n",
    "        nr_feats = 10\n",
    "        int_ars = []\n",
    "        for rng in [range(nr_feats*i,nr_feats*i+nr_feats) for i in range(len(stacked_array_f[0])//nr_feats)]:\n",
    "            ar = stacked_array_f[:,rng]\n",
    "            np.apply_along_axis(impute, 1, ar)\n",
    "            int_ars.append(ar)\n",
    "        stacked_array_filtered = np.concatenate(int_ars,axis=1)\n",
    "\n",
    "        clf = load_model()\n",
    "        probs = clf.predict_proba(stacked_array_filtered)\n",
    "        pred_array = np.argmax(probs,axis=1)\n",
    "        none_indices = np.where(np.amax(probs,axis=1)<0.5)\n",
    "        pred_array[none_indices] = 99\n",
    "        result = np.full(len(stacked_array),np.nan)\n",
    "        result[~np.isnan(stacked_array).all(axis=1)] = pred_array\n",
    "    da = xarray.DataArray(np.transpose(result.reshape(array.shape[1:])),coords=[array.coords[\"x\"],array.coords[\"y\"]], dims=[\"x\",\"y\"])\n",
    "    return DataCube(da)\n",
    "\"\"\"\n",
    "\n",
    "def computeStats(input_timeseries:ProcessBuilder):\n",
    "    tsteps = list([input_timeseries.array_element(6*index) for index in range(0,6)])\n",
    "    return array_concat(array_concat(input_timeseries.quantiles(probabilities=[0.1,0.5,0.9]),input_timeseries.sd()),tsteps)\n",
    "\n",
    "def rf_classification(bbox, udf, year=2020, rel_orbit_filtering=False):\n",
    "    temp_ext_s2 = [str(year-1)+\"-09-01\", str(year+1)+\"-04-30\"]\n",
    "    temp_ext_s1 = [str(year)+\"-01-01\",str(year)+\"-12-31\"]\n",
    "    spat_ext = {'west': bbox[0],'east': bbox[2],'south': bbox[1], 'north': bbox[3],'crs': 'EPSG:4326'}\n",
    "    \n",
    "    if rel_orbit_filtering:\n",
    "        from eodag import EODataAccessGateway\n",
    "        dag = EODataAccessGateway()\n",
    "        search_criteria = {\n",
    "            \"productType\": \"S1_SAR_GRD\",\n",
    "            \"start\": temp_ext_s1[0],\n",
    "            \"end\": temp_ext_s1[1],\n",
    "            \"geom\": {\"lonmin\": bbox[0], \"latmin\": bbox[1], \"lonmax\": bbox[2], \"latmax\": bbox[3]}\n",
    "        }\n",
    "        all_prods = dag.search_all(**search_criteria)\n",
    "        rel_orbits = [prod.properties[\"relativeOrbitNumber\"] for prod in all_prods]\n",
    "        (uni_vals, counts) = np.unique(rel_orbits, return_counts=True)\n",
    "        creo_settings[\"props\"] = {\"relativeOrbitNumber\": uni_vals[np.argmax(counts)]}\n",
    "\n",
    "    ### Sentinel 2 data\n",
    "    s2 = connection.load_collection(s[\"SENTINEL2\"],                                 \n",
    "                                temporal_extent=temp_ext_s2,\n",
    "                                spatial_extent=spat_ext,\n",
    "                                bands=[\"B03\",\"B04\",\"B05\",\"B06\",\"B07\",\"B08\",\"B11\",\"B12\",\"SCL\"]\n",
    "                    )\n",
    "    #s2._pg.arguments['featureflags'] = temporal_partition_options\n",
    "    s2 = s2.process(\"mask_scl_dilation\", data=s2, scl_band_name=\"SCL\").filter_bands(s2.metadata.band_names[:-1])\n",
    "\n",
    "    ### Sentinel 1 data\n",
    "    s1 = connection.load_collection(s[\"SENTINEL1\"],                                 \n",
    "                                    temporal_extent=temp_ext_s1,\n",
    "                                    spatial_extent=spat_ext,\n",
    "                                    bands=[\"VH\",\"VV\"],\n",
    "                                    properties=s[\"props\"]\n",
    "                                   )\n",
    "    #s1._pg.arguments['featureflags'] = temporal_partition_options\n",
    "    if s[\"SENTINEL1\"] == \"SENTINEL1_GRD\":\n",
    "        s1 = s1.sar_backscatter(coefficient=\"sigma0-ellipsoid\",options={\"implementation_version\":\"2\",\"debug\":True,\"tile_size\":256})\n",
    "\n",
    "    s1 = s1.apply_dimension(dimension=\"bands\",process=lambda x: array_modify(data=x, values=x.array_element(0)/x.array_element(1), index=0)).rename_labels(\"bands\",[\"ratio\",\"VV\",\"VH\"])\n",
    "    s1 = s1.apply(lambda x: x.linear_scale_range(0,1,0,250))\n",
    "    \n",
    "    if s[\"mask\"]:\n",
    "        ### Cropland mask\n",
    "        cropland_mask = connection.load_collection(\"ESA_WORLDCOVER_10M_2020_V1\",\n",
    "                                        temporal_extent=temp_ext_s2,\n",
    "                                        spatial_extent=spat_ext,\n",
    "                                        bands=[\"MAP\"]\n",
    "                                        )\n",
    "        #cropland_mask._pg.arguments['featureflags'] = temporal_partition_options\n",
    "        cropland_mask = cropland_mask.band(\"MAP\") != 40\n",
    "\n",
    "        s2 = s2.mask(cropland_mask.resample_cube_spatial(s2).max_time())\n",
    "        s1 = s1.mask(cropland_mask.resample_cube_spatial(s1).max_time())\n",
    "\n",
    "    ### Base feature calculation\n",
    "    idx_list = [\"NDVI\",\"NDMI\",\"NDGI\",\"ANIR\",\"NDRE1\",\"NDRE2\",\"NDRE5\"]\n",
    "    s2_list = [\"B06\",\"B12\"]\n",
    "\n",
    "    s1_dekad = s1.aggregate_temporal_period(period=\"dekad\", reducer=\"mean\")\n",
    "    s1_dekad = s1_dekad.apply_dimension(dimension=\"t\", process=\"array_interpolate_linear\")\n",
    "    s1_dekad = s1_dekad.resample_cube_spatial(s2)\n",
    "\n",
    "    indices = compute_indices(s2, idx_list, 250).filter_bands(s2_list+idx_list)\n",
    "    idx_dekad = indices.aggregate_temporal_period(period=\"dekad\", reducer=\"mean\")\n",
    "    idx_dekad = idx_dekad.apply_dimension(dimension=\"t\", process=\"array_interpolate_linear\").filter_temporal([str(year)+\"-01-01\", str(year)+\"-12-31\"])\n",
    "\n",
    "    base_features = idx_dekad.merge_cubes(s1_dekad)\n",
    "    base_features = base_features.rename_labels(\"bands\",s2_list+idx_list+[\"ratio\",\"VV\",\"VH\"])\n",
    "\n",
    "    ### Advanced feature calculation\n",
    "    features = base_features.apply_dimension(dimension='t',target_dimension='bands', process=computeStats)\n",
    "\n",
    "    tstep_labels = [ \"t\"+ str(6*index) for index in range(0,6) ]\n",
    "    features = features.rename_labels('bands',[band + \"_\" + stat for band in base_features.metadata.band_names for stat in [\"p10\",\"p50\",\"p90\",\"sd\"] + tstep_labels ])\n",
    "\n",
    "    clf_results = features.apply_dimension(code=udf_rf, runtime=\"Python\", dimension=\"bands\").rename_labels(\"bands\",[\"pixel\"]).apply(lambda x: x.linear_scale_range(0,250,0,250))\n",
    "    return clf_results\n",
    "\n",
    "\n",
    "job = rf_classification(bbox, udf_rf).execute_batch(\n",
    "        title=\"RF model inference on small user-selected area, GTiff\",\n",
    "        out_format=\"GTiff\",\n",
    "        job_options=job_options)\n",
    "results = job.get_results()\n",
    "results.download_files(\"./data/test_klein_gebied\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0b70181",
   "metadata": {},
   "source": [
    "Finally, we can have a look at the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "25ddacbf",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAELCAYAAADz3X2AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABrlUlEQVR4nO3dd5hcZdn48e85Z/rM7mzL9k2y6b0XCC2URIKhd1SKvgpSNAKiiL4Gf5gIauSVpigSBOlNkCKBkFBCSSW9kZ7sZnufnXLO8/tjNptMtmQ32b7357rmurKnzT1nJzv3POV+NKWUQgghhBCiC9E7OwAhhBBCiKNJgiKEEEKILkcSFCGEEEJ0OZKgCCGEEKLLkQRFCCGEEF2OJChCCCGE6HIkQRFCCCFElyMJihBCCCG6HElQhBBCCNHlSIIihBBCiC6nUxOURx99lNzcXFwuFxMnTuTjjz/uzHCEEEII0UV0WoLywgsvMGfOHO655x5Wr17NaaedxqxZs9izZ09nhSSEEEKILkLrrMUCp06dyoQJE3jsscfqtw0fPpyLLrqI+fPnN3uuZVkcOHCAuLg4NE1r71CFEEII0QaUUlRWVpKZmYmuN99GYuugmGKEQiFWrlzJz3/+85jtM2fOZNmyZQ2ODwaDBIPB+p/379/PiBEj2j1OIYQQQrS9vXv3kp2d3ewxnZKgFBUVYZomaWlpMdvT0tLIz89vcPz8+fO59957G2wfdOP/YjhdxxWDvUaR9vJmrKoaAAIzxlI6OHo7vPkW/ldXo6zjb1wyBvRl/6xUVDdr4EncFsH17iqqL5hIRT8DAN9+i/h/r271tapmj6Oir3FccSRtCeP6YC0Aenwc+ZcOIeI+rksdkzffwv/aV6CsVp+ruV0UXjYCR7Ui/qhrGOmp7L8gB8sGmoKs94qwvt7d4Bq6x03BZcMJ+aJvlvQva9C+WA+ALTOd/Rf1xTLqrvFOAdauvcf5ShsxdggHTolrk0tlLS5Fbd4Rs6327DGUDbCT9douzKLiNnme46U5HRReMQp7tTqu33fNN8ZSmW0j47UdWCWlrXvykYM5cEZ8zCZfnkX8a43/vzKnjODgpNg3fMKOCJ53v2rxU5onjeDghMb/02R9WIba9HWD7ZrDTtFlowkmtO4PV9xek7g31zR9gKZTduk4alJjr5v1filq644mTjrqEk4HA98O80DGqlbFJrqWiiqLfhN2ERd37L87nZKgHHJ094xSqtEum7vvvpvbb7+9/ueKigpycnIwnK7jTlCMiMKmObC0MAA2uwvDGb0dZooicu4UPNtLMbftPL7rG04Mp6vbJSjBbIX9vJMJp9swnNHgbXYLm2Zv9bWi9/T4EpTavk5s504BwLJraD4bRju9WyN9IPKNSXi2l2J9vQsA3R9PzZSBuA9UoTZsa/JcTXNE34chVXePjkhQ9Oh7QKtLUGyGE6ux+xhSpK8JUZvhoWS4HZvNQre5ME8eTXWqA91tcOi/RdWkbNz90nB8vgkVCp3wazfyKslZ7qZkuJfa5BN7s9oMJ6ru9RnpaVSPzaS27n1k0x1ox/EeajOjhxLI9kK8jYin4e+7OUZWBtWj0vFUm8SvDqPVhNFb+1rq/h4cKdKn7u/MthKsHbGJq2Zr+Let6BSFdc2EFj9lSmIlkxL3N7rvq7OzqK5u5FqaYlTWAeIdgRY/D0BejZ+dlzcf25DMg6S4qup/tpTO2rMzqGksjkZommJu338Q75LJpz1BS4ZndEqCkpKSgmEYDVpLCgoKGrSqADidTpxOZ7vGpJkKPQLKgLBXo3i4HVtNPPbdsX+IlGmhTLNdY+lMwQSNYIId3QQ9mruhR46vJUkzqb+nrU3UDsXRGN0E1PFdt8G1ImA5oHi4HSPgx7HXgQpH0Fwuygfa0Swvzg3NX0OLRF9rAyr6nuLQ7TOb+MZumrBuC57aXEqG90HpGprbTdkgV4OkoTJHJ5jgIG3t4W/GKhw5rtYfALOwCK2wCHf6BELxsX8OlAGqBZ8FmhV9/Zqp6l+qivdSOij6+9MjxxVaPd3jgRP8/182NI7yAYdfTGGCnT4qEU9J+THPDfVNoXCsnYxlJvqufPB40DyeVj1/ON6B5YjdFnRoFCbaSTWTcJdWxOyr9dkaHO8cXs66qc+26nmb1K9tLhNj+HGc07+1J0hy0pt0SoLicDiYOHEiixYt4uKLL67fvmjRIi688MLOCAnPZ9vxboon7xuZhOua2otHuDAGjI85LnFrLfqnazsjxA6VurwGx84CAFQwyPF8/MV9tI24RD9538ho0+6ZlK+CuHYUkv+NLEJxx5+h6GHIercAy+fkwHQ/JcOd2PqNJXXRHszCIjJfsVC1zb92VRsk7a0dqEgE66gkwSoqJvMVE+oGglmlZS2Kq2iMG2PoKMJNvLawV6Pg4iH1iU/qp0VYx9nSd4jv4+3ELY/9xl4xNYeyAcduAYvfY+L/bC9W2eEPe7V7P9kv131btizM1naJHGHLvSOYM/Od4z4fINP+BQlGdcy2wkg8hZH4Js44LM74ihx7MfuvTaTc9B7X8/uN1WTZG78H+dclUBLxxWxLsq0i3VYWsy3TqARalxgJ0Z11WhfP7bffzne+8x0mTZrEySefzOOPP86ePXu46aabOiUeq7ISLRRCszLrt0XcEHHHfkgEUh34+2VjFRRhBWobXEfTNYyMdMIpvgb7ujLdBGeJwnRCKF7DVhnEPFhwQte0yivQLQUqo01iNILgLFc4D1ZhFRajWVkndD1NASXlGNUOPAfjCMVpBBM1IplJ2AAz/+CxL6KspsdW6Dq4jvjmb7Ssu+vI953lgECmCc3kYcpIwTcwsUXXbo3SIQY1OcduLbTsBkZtNtDcgLe+xx1H/9EHuC2x4didE1da92ghdxFQ1A5xFNc9jkWSE9G7dFqCcuWVV1JcXMxvfvMb8vLyGDVqFG+//Tb9+rVH22PbqehvUNk3i6y37dDIt1bN6eTguTmE4rRuNf7EVq1I+fdGzCF9OXB62wyabGvefJO4N9ZgmSZaG3b5RfYfIOH5g4Smj6VwnJO8U+NwlfhIfr4YFTn+/gk9OYn952cdHiT7th3VgjEPRwqkWXx54QLidEeTx4QvbJ8uRx0dowX9xKZSWMfVxtYyTs0OHN9YJiFE99Wpg2Rvvvlmbr755s4MoVH2KkXC9hCaAmVolA5yEKn78qK0aN981fAknOkNP8gjNp2IS2tR331HcFQo/Duir+Vo1Rl2atJ0/DtM3AVBVG0Q28EyUlc70IvKae3Hnu7xEB4zAGWLvnjnlgNY5RWkrAsQSHNQ3r/1HzKB8QEyU8oAKKzwcXBKdECdMhQJg4tJdIRbfc1DwpbO5kEDMGoHAWBl1pLVJw+Aylon20aPA+v4s0zLY5E5MB+97uZvGtoHo7JPk8ebcSY5uXkx28bHF+PXXdi1pu+dszMHn0KzrTtCCHG8OjVB6aoclQrHx+ujgxcNA0f6BMwjunqUBqWDbLTr7TuRBEdRPz7BUalwLtsEViPfcKeNpDrDSdyWUtSOaAVf62AhjoOFWIBma93r0/zxFEz2YNV9XmaVJ6OVlmFbvom4YQMoH5DQ6pdy/ejP+GXK5sMbTmn1JZo3tpl9k9v4uUYf74nSeiCE6H0kQWlEbZJO2WXjSdxQgVq/lZSl+1CO6KduzeBkike07zdW0wWzLv6cwe4WjIFoxP0fn4dvWzTGwomQe6uz/lv8kbLcK8l1FrL2ihwqwsc3XftIbqOU8+K+wqhr7t94RSYloeigwiTHNi70HGj1Nc/zbQG613geIYQQJ67XJihK19ASEzASEzCToh+AptNGIE0j7FVUZ2nYauOIz0vCOlhYPxbBFe/GkX3skf+H1CZrhPytm6Zreiy+m/wpIx3HN/XliaxyykqTAXDmVvJs7ofNn5DQeK2EE9Ym15XkRAgheqNem6CYLth/fhblY0L855yHANBRxOlWfYP6ecP/h+r0QWS+tQ9zX/Tbv7bha9K2NT1g8Wib7h/Gsm8uaFVshqaRrB//vNz3xz1FzZjoCBKXpiOj/4UQQnQ3vTZBURqYDrAX2bly1f9w14j/colvH7M3XkOKu4oXBrxHdcBB+r4I1AZBiw4KUaaFamR6cZN0RYatY1oBvgyG+WfRqdyQ8jETne1UF14IIYToAF1krknnSV1ukX3FFuatncUB08R1Txx5fx5IBBO1x4vn7TWYxcdfZKoj/bPoVJa+MJHnS6d2dihCCCHECem1LShGENKWV2ErqMAE0p90cdXiO0nf8zX0z+3s8IQQQoherdcmKHpEYWzejVkVLX/tWrQal6ajpaYQ9Lddw5JeY7A2VMsgm46nmWJbJ8JUFnlmDUHLRsQLHuPEF5HrLEEVZkvYJEmPkN1BXWNCCCG6nl7fxXMkPSmBxJdqWPCLx9qs+NWwBfu48+obeaTseFbSapkN4RBnvHonnx/oxzPXP8icpOXt9lzt7YXKDO749k1Mf/HOzg5FCCFEJ+q1LSjK0FD9MtFr6yqRFpWg6TrnJG3k9LqSIJF4E31AX9SBg6ja4HE9j5lfgK2yktLw8S0y1hK1ysBRohPONJjodADt01LTESotN/Yt+3GWDO7sUIQQQnSiXpugRFyw/5zDC6xlfuTG2FcYc0xW/yL2z0on+99hzH15R19CCCGEEO2kVyYoidsj2Gosioc7cFQqErbVYuSV1O//Mhjm9i1XcmB7H5pY8b5L6W8LMfYbm5nob48VXztHxme1DIn/IQCROIsPZv+RXLuMSRFCiN6iVyYovm3lcLAIffBQnBUWxpcbsWw29LQ+2DWTVYH+lC9Ox9sWC7TqGthsjZaabyuphpfncxe32/U7mmYzsH2yngEfR38B+oB+7J3lI7eT18QTQgjRcXr9INmqDIOSqybw9d8HMem1bZzn2dum17cmDqPg4iEMcUkXUUv0dxSSd2Eu5kkjOzsUIYQQnahHt6BoCuyVCj26jA4Rl0bEA5EENzYrCaWD6YRAHw2vJ4jfCKBrbdunE46zE0jR8OodO/U3rEw+rbUTp9fWDZztHrxaqG79IjtuQE9MJNInrm4Bwl6fTwshRK/RoxMUFKQtKUDVDXCNTBhC/kkeDk7xgPJgHfHqU37v4v2C8aS9Wd5JwbatPDPA91+5FSujlq/PfrKzwzluFacPoHiUQR8jALTfTCghhBBdS49IULx5Fs5yi/JcG6Yzus1dqPAURAhnxKP1iQOgJiO60zIaXsMImmiV1YSVwTjXHlynF2EpCEVsmB/HQTebxROn6fSdsJ/+vpJjH9yF5Ngq8J96kD0DE3FNmkJoYIC+6SUkSeOJEEL0Kj0iQfF/XYO+YSfVGaMwndEumvg9QeyfbaTkivHUpLau2+Ykl8HyCS8CUG4FuDD1R7jaPOr2lWh4+GDEG50dRqsNtPv4bOwrMPboPdJ6IoQQvUmPSFCKR3qx544k7DmciJQNcuJIG0fQ3w3mCQshhBAiRo9IUIKJGsHE2EQk6NdalZyYLht6nBe7ZrZJTJqho/m8hPw2wn6FVz++SrRCCCFEb9QjEpS2UDjejel0k2MvbpPraTmZ5M1Ix3X+Qd4Z8RTZNicghTyEEEKIlpAEpY7pANMFDkygkVG0raXrmG7o46lmoFRAFUIIIVpF5kYIIYQQostp8xaU+fPn8+qrr7J582bcbjfTpk3j/vvvZ+jQofXHKKW49957efzxxyktLWXq1Kk88sgjjBzZuuqhiVsjVIxu61cQy0CjaIydZGNcq84LeyX3E0IIIY5XmycoS5cu5ZZbbmHy5MlEIhHuueceZs6cycaNG/F6o1NFH3jgARYsWMDChQsZMmQI9913HzNmzGDLli3ExcW1+Lm8uyuoGna4+0TZQLXxpB0dnZpsEy0ivWFCCCFER2nzT91333035ucnn3yS1NRUVq5cyemnn45SigcffJB77rmHSy65BICnnnqKtLQ0nn32WW688cYWP5favY/sFwKgaeCwk39uNqHusPywEEIIIZrV7v0Q5eXR0vFJSUkA7Ny5k/z8fGbOnFl/jNPp5IwzzmDZsmWNXiMYDFJRURHzALBqQ0QKiogcLMQqLEY7wdWHNRNeKZvEQ6X9uHX/VN6rsWNoGp7MKqoGRKgaECGY0n6rEgshhBAiql37LZRS3H777Zx66qmMGjUKgPz8fADS0tJijk1LS2P37t2NXmf+/Pnce++97RkqAHoI3nnhZDz5iuRX13Pj729g5wWPs3rqP7GIZj/f2jGLLf8e0u6xCCGEEL1Zu7ag3Hrrraxdu5bnnnuuwT7tqFWDlVINth1y9913U15eXv/Yu3dvu8QL0VaU6iyNbXNHcs74DQDYNQOnZsep2bHpJ9hMI4QQQohjarcWlNtuu4033niDjz76iOzs7Prt6enpQLQlJSMjo357QUFBg1aVQ5xOJ06n89hPqkBTJzZQVlNQ0zfC1vMfQ0cjqEycmhRYE0IIITpSm7egKKW49dZbefXVV1m8eDG5ubkx+3Nzc0lPT2fRokX120KhEEuXLmXatGnH/7yhEGmL9pH+WTXacQ4T0cOQubicpNUGFhYjPrme6XfexiNlOccdlxBCCCFar81bUG655RaeffZZ/v3vfxMXF1c/5sTv9+N2u9E0jTlz5jBv3jwGDx7M4MGDmTdvHh6Ph2uuuea4n1dZisjeA9gNA5QXjqMVRVNg5BXhKolOXbZ2e0n4z3p2/TQFaF23kh4BW42iMtSClh8hhBBCxGjzBOWxxx4DYPr06THbn3zySa6//noA7rrrLgKBADfffHN9obb33nuvVTVQujp3oUXi6+vYPGkEtK7+nBBCCNHrtXmCotSx+1c0TWPu3LnMnTu3zZ5XMwy0oQMIZPqOq/WkMRlj89lz62hu8P2LfZEqHig4k1V7cmhJm4hmgRWoBUvqsgghhBCt1WPKo2p2G4UnJxL0t11C8NHo16CulP67Ncl88MpknOE2u7wQQgghmtAjFozRJgyn8rzRhD0nlpxYNig9I5eiMTr6UbdmgrOEC6/4hPCUyhZdq2SURuHrg7jz9HdOKCYhhBCiN+rWLSi6142uOajO9lDe3zjh6ykdKnN0QqkNm0lSDS/z0tayoyaFDQw75rVCqRE2TnwOQ2s+BzSVRZ5Zg0vTSDG8xx27EEII0ZN06wSl4IoRGA4XprP7jvPYHA5ywet3kjSohOUTXuzscIQQQoguoVsnKGGvhtWNkxMAl2aipdaSEVfR2aEIIYQQXUa3TlB6goF2H1vO+Mcxu4KEEEKI3kQ+FRvh2WVn5FO38u1d0zvk+SQ5EUIIIWLJJ2Mj4vZYDPj1Sj5dI6sWCyGEEJ1BEhQhhBBCdDmSoLRSWdCNq1hhhDo7EiGEEKLnkgSllbbuTSP52VV4CqzODkUIIYTosWQWzzEEVZjT1lxDOGLw6cR/cvawLXz8i4lEBgXQdYu0l1yYDo3K7GiuZ9nBmlDJadm7e9zg18fLM3m3cCRr92WR5K/m03HPY9dOvECeEEIIcTRJUJqia6BBrYrgeTgBW8Ck/J8h/pbzKXz/UwDeqPbw17vPwUrwUZmTAEQTlKcm/YMpTnsnBt8+/rLtNMzFyQx9ZQ9VYzMJ/iUsCYoQQoh2IQlKI0pGa/i/24f7sl9p9riTXIX88+kaXEYF0+OWA+DUwwy1R4Cel6D8fuTLfNF/ECsu6Uuu+yvcmqOzQxJCCNFDSYLSiHBihP8M+zd2zaDcgup0G7ZaA0OLrVqbanh5eeD7jVzB3TGBdpAqq5bdEcUEh+Js9xZI2VK3p2d1YQkhhOg6JEE5Br/u5olf/wmIJiS90d/Lh/GXV2Zxxrlr+Gv2Z50djhBCiF5AEpRG2IttnL72CnRNYWiKh4c+xxiHq7PD6nA1VoiHSkfyzPbJJK+z2HJSGmR3dlRCCCF6A2mjb0TSBoX/gr3Ezd6D7/Jini45ubND6hQHzRD/eOMc4p6PJ+61lezJT+rskIQQQvQS0oJyDEqpzg6hU3xvz6ks2T6YfkvD1KTa2PHMCH4y+oPODksIIUQvIQmKiBFUYSqtEEu2DSbuSzeuHXmUDspg02kLe1xdFyGEEF2XJCgixryicfzrvdNxlmmEEmD0Szv5UdzbkpwIIYToUL02QdFNcBUrTAcEEzSc5QojCLVJ2rFP7oGqrFr+U53Bpsp0TJ9JwAeaJ8JNSZ+Qa/d1dnht7vNak13hFGZ78/DpvW8AtBBCdHW9NkExahXJr29A9c9k38wkUtZUoW3cQcG3xnR2aJ1ifcjOvc9ejT6mnM0XPAKAjo5d63nJCcD3v/oOkTUJ9L/uIU6S/EQIIbqcdm+3nz9/PpqmMWfOnPptSinmzp1LZmYmbreb6dOns2HDhvYOJYZl16idMpjKwX4AKvt7CJ00DNOpEfZphE8bjdG3d8yp/dGByfxq50Wkn7qfqwevxKnZcWr2HlXGvsis5pLtM5ix6XxmbDqf1LgqcqfvItsW6OzQhBBCNKJdW1CWL1/O448/zpgxsa0SDzzwAAsWLGDhwoUMGTKE++67jxkzZrBlyxbi4uLaM6R6pgMKxznrf67oZ1DRL/qBbDo1CiY4STcTsRUWd0g8ncFUFkEV4c3V4zDKbCy76g89thjdQVNn44eDsVcBGky7bHVd0bme2UIkhBDdXbu1oFRVVfGtb32Lv/3tbyQmJtZvV0rx4IMPcs8993DJJZcwatQonnrqKWpqanj22WfbKxzRiLvyJzHyrVuYNW4d8y58lsReMBajalCYOTe8yh1pjS1RIIQQoqtotwTllltu4Zvf/CbnnHNOzPadO3eSn5/PzJkz67c5nU7OOOMMli1b1ui1gsEgFRUVMY+OEHHb0FOScOqRDnm+jlYW9mBUGpwev4UrfOU9qkvnSJtCNXwaGIhmAXbFYGc++aaXz2tNgirc2eEJIYRoRLt08Tz//POsWrWK5cuXN9iXn58PQFpaWsz2tLQ0du/e3ej15s+fz7333tv2gR5DyXA7xSOzuMb9SYc/d0d4KHsxVVf+t67lpGcmJwAXLLsZ+wYPthB4tzi4acfNAJhuxdPX/JkpzmNcQAghRIdr8xaUvXv38uMf/5hnnnkGl6vpLgPtqJWBlVINth1y9913U15eXv/Yu3dvm8bcFMsAyw4GVoc8X0fz6A5SDW+PbTk5ZHK/3YRH1mDZQDPBqI0+bFUa/7vzIn5fMrCzQxRCCHGUNk9QVq5cSUFBARMnTsRms2Gz2Vi6dCl//vOfsdls9S0nh1pSDikoKGjQqnKI0+kkPj4+5iFESz2b+yGvT3sM86h8WQ/D/rf78einZxFWZucEJ4QQolFt3sVz9tlns27duphtN9xwA8OGDeNnP/sZAwYMID09nUWLFjF+/HgAQqEQS5cu5f7772/rcBrVZ00QZ0F1k/tDKV4KxztRvbNmW4+UbYNvX7yYf6w6haEPBigdGU91hlTHFUKIrqrNE5S4uDhGjRoVs83r9ZKcnFy/fc6cOcybN4/BgwczePBg5s2bh8fj4ZprrmnrcBrlLKyBr2O7iVQwiIpEB8O6BvSH8Zn1+zbVZvKlowCAdCNIX5tMTe1u/LqbX6Zs5sPsIejlGvZAHLba6D4tLImKEEJ0NZ1SSfauu+4iEAhw8803U1paytSpU3nvvfc6rAZK/qmJaCcnxGxLW1aKWru5wbF6CJbdOoUvS0cDsPV/Evj6ir90RJiiHTwz5FmWv5fKrx69nqyntwDgOG8I1sUWPXmgsBBCdDcdkqAsWbIk5mdN05g7dy5z585t1+fVFLgLFZqlqEnTsVUrXKUq5piwR6M2WaM2Mw53bXSwZCg1LuYajr3FqEAtpWcNQE+tbdeYe4J3a5w8uv9MfpbzNqe4ulbrRIbNxwW2Gn48Joj37MEkvL8VW6BnDoIWQojurGevxaMg5aP9UBskcOVA4ndHcL4TO/VZHzOMfTOTKB5ph5GHB+kePf4kPCyLF+7/A1mGhw5YIaBb+/HKK8m9dgt3v3EpH41+rbPDadT2mX9j85lB7rjs+50dihBCiEb07ARFg4oJGdgrIqQtr0bZdcIzJ+H6ag+qvILwySOoTHUADROSQywbbLozg7jsClJ0B4Ymycmx3DjyEx569CyMtXaGrP4hrmFlnJq1k0ezPu/s0OoZmk62Dfb8HBJ9hdike0cIIbqUHv1pqzQoG2BQ0c+BsWYbetCkeIQDkhPQfF5KhjqpztDRI0SrjDZ2DQPuOfvffDbpKTy6o2NfQDd1e9IOds76O3G7NQY/vAdtSSLvbRvW2WE14NfdbJz2DJ+OeVUSTyGE6GI0pZQ69mFdS0VFBX6/n6E/mofhPPb6MXoEHOWKuP1h3BsOUDUhm5BPJ+mLfAiGACibllO/WGAMDWr7KMgOsPr0v+DrBevVtJUXq/xsrc1guOsAOfZipjjtnR2SEEKITlRRaZE4ZAfl5eXHrGnWs7t4jmLURIjsP4A+KpNDjUeqJoBZWootkN34SQpcBRo1dhcm3S6X61RX+MrBV173kyQnQgghWq5XJCj2SkXKC19h1QYBcCxajTspgbzLhuDLT8X9+pedHKEQQgghjtQrEhTTqRGZMATNiraA2LbshUgEZWgEEg1c08ZSm9T8IElblca3tl/CBWlf8QP/gY4IWwghhOi1esXIwIgH8k/ykDfNS940L6QmR3coCCZp5E3zUpOqoTXTg+Mo19jzZi5/2XZaxwQthBBC9GK9IkE5WtHUFCrPGkbG0mL6rIl2+6SsD5P13xJsNY2f46hUZL9TRPn2xA6MVAghhOidenyCYguArYaY1pFAikZVuoFWVom9pBZ7lcJ1oBK1bSfOcgtboOF19DBY23bhKOvxt0wIIYTodD3601azIOP9AjLf2Y9mxu4z3XDg4v7U9POS9vRa1MbtqGCQxNfWkvFhUbPdPUIIIYRoXz1+kGw4PQ4touCoSrGaCZ6DJq6DQazqamzZWah4L9a2Xei1oQbXMR2gjRyEHoGfHRzHqpIcdBTPD3mRRMPTQa9GCCGEOHFvVHtYWjGM5UX9yPSV83T/Rdi1rlVRu0cnKEqH/MnuRvcZQUX8uxuxKisBqBqXRfkAG5kFxY0eH/Zp7DsnAT0E7zwzjex3ilC6zrY37UzpWr9TIYQQollzN52PuTiZ7Jd2kT9uIMG/vi0JSldhOjSqzhmBHo725TjKw2R8VIVVWYURH9fsuUqD0rFJhOI04rQwUoRMCCGEaFu9NkGx7FA28HC2mPVhNdaajS0+vypLJ5SgOGDGkWxWkGp42yPMTlNgVlNZVzfGo0GGzdfJEQkhhGgrTnuEShegd92hqL02QWkL7oMav7v+WvZ8w8WW7z7W2eG0qdOfvJN+b1cDsPdsHxtvebSTIxJCCNFWnh3xFF8MyOHJJed3dihNkgSllezVCkdlXUXagMK2ZjuesaM7Oaq2F45X1GS46/8thBCi58i1+4jT9zFvmo+IG/QuOKlXEpRWit8VwfnuivqfLa3r/VLbwsbLH8K8LJqYGJqGjLMRQoieJcXwsvyO/wPAqTk6OZqGJEFpJU0BKvrBrbtcRCYNoyqn57UwODV7g6nZQgghehan1nW/fPbMr/9tQLNoUKxNs6hPTgA0t5uicR5UToCgCmMqq2ODFEIIIXooSVAaYR3IJ/u1fSR8Hanf5qhQZL+Rh3vFjgbH27d6GPbGLSwoHdyRYQohhBA9lnTxHMVITgLdQFUH6mukQLQ1RQsEUbqB0adPdKPfh9JAi4BRZVBjOjspaiGEEKJnkQTlCJrTSeEFQwnHRQdfqCOK6oXiNfZd1j+m20dpoGwQHBbgizMeJk53IINJhRBCiBPXLl08+/fv59vf/jbJycl4PB7GjRvHypUr6/crpZg7dy6ZmZm43W6mT5/Ohg0bWv08jipF0uYIyZvC9Y/4PdZxLfRXk+3FmjCMiEdDi4B/R4SEbRGSNkcfiVsiJHwdwb8j+jjUupLwdQQKnaQY3i492EgIIYToTtq8BaW0tJRTTjmFM888k3feeYfU1FS+/vprEhIS6o954IEHWLBgAQsXLmTIkCHcd999zJgxgy1bthAX13yZ+SO5Sixc//kyZps3tx+V2VmtnoFSOtgGg6O3w3fAwvX2Smhm0Gv4iqnUJoL3g03EZ/a8OihCCCFEZ2rzBOX+++8nJyeHJ598sn5b//796/+tlOLBBx/knnvu4ZJLLgHgqaeeIi0tjWeffZYbb7zxhJ7fKigi+78ulKaBoVE4JYGgv+lsxVOgSFpVGo3NbnDwZH+jxxkJCVSdMRjP7mrUmo0kfpmHcjmwAgEylhQxwfohY69fz5N9Pz6h+IUQQgjRDl08b7zxBpMmTeLyyy8nNTWV8ePH87e//a1+/86dO8nPz2fmzJn125xOJ2eccQbLli1r9JrBYJCKioqYB4DSNXSvN+YBoHbtg5170XbuRw81H68tYMHOvbBzL/qePOzVCiPUsI9Iczkp728jkOVF9/mwikpQew6gTBNz0zb6/PVzlmwawspgiKAKt/a2ieOwL1JFXqSqs8MQQgjRDtq8BWXHjh089thj3H777fziF7/gyy+/5Ec/+hFOp5Nrr72W/Px8ANLS0mLOS0tLY/fu3Y1ec/78+dx7770Ntlen64S/M6bZeEx38309VVkGgbpr6BFIe38/VnEpVhPdO2UDbVRmjwLAXgNJr6zFqqkBYMQv8/h57g+45cmXucgrH5ztqdSs4bR3foLhi7Bp+t+73DLhQgghTkybJyiWZTFp0iTmzZsHwPjx49mwYQOPPfYY1157bf1xmhabOCilGmw75O677+b222+v/7miooKcnBwsG4SdJ1bu1LKB5dOoybSw3BZpixRWVWxyYeuXQyQ9AbTo+jvOsmgLi73GAtOsPy5yIA97bS13rbyEjwev548Zq04oNtE0XdOIS6sizhVEl5K3QgjR47R5gpKRkcGIESNitg0fPpxXXnkFgPT0dADy8/PJyMioP6agoKBBq8ohTqcTp7N9a4ycdco6bk9bxE/+/n3Yc3i7ZhgUnZ5FoI+G0iBuz1Fr8Rx1HbOklNxryvns8imEFyyXb/btxK+7WTvlubqfpN6gEEL0NG3+l/2UU05hy5YtMdu2bt1Kv379AMjNzSU9PZ1FixbV7w+FQixdupRp06a1dTgt9v7aEVz85Y0YpdUx25WlSFpTSvpn1WQsq8a7pfDYF1MW3v21TF15DX8py2qniIUQQoieq81bUH7yk58wbdo05s2bxxVXXMGXX37J448/zuOPPw5Eu3bmzJnDvHnzGDx4MIMHD2bevHl4PB6uueaaE35+TYF2uNcFZUQLqh1JN4GjxsHGb7Cjh+yoigOxO5SFtW5LfSdChJaxH6wgsjSdt/2juSlhfytegRBCCCHaPEGZPHkyr732GnfffTe/+c1vyM3N5cEHH+Rb3/pW/TF33XUXgUCAm2++mdLSUqZOncp7773XqhooTXGWKFLf3wuWBbpOwYwcapNiM5S0L6qx726kJcSyMMvKTzgGIYQQQpyYdil1P3v2bGbPnt3kfk3TmDt3LnPnzj2h5zGC4K5UhH0aEdcRO3Qt2myiazjLLZSmE0rQYltSHHaC/ZJRxtEDLA93ybh2FhPZ2fjMIgDd7SYyeSiWLdpT5tqSj3mwEGvKCEr7uajqa9HfW3xCr1EIIYTojbr1WjyeApPEt1YSOHccJUOjLyWYpLHvkhwg2tWT9eZ+/KEw+y7vj6p7tfkneQn5vSz8zkMMstc2ef1pT99J7j1NJyhavyzmPPEcU53RJOTsP/yU7BcifPPxD/l2/Cbsml5X/l4GygohhBCt0a0TlLBPx5oygtpEA2UDa1IF8Z7DCYdp6Wzun4FmaviHFWE3Ds+5yXDWMtQeJNHwNnn9IdN2se3/ndT08/stxjqKSDF8ADhmFLElJ5ffet4g0fAc12v6fclAvqrIYUHOW6Q2E5sQQgjRk3XrBKU2USMv3Q0aWHb4y4RnON111EETmrtC80nEf4a8A0OOFYWv/l/LJ7xY93zHnhJtKgtDaziJ6u8bTkHf6GPX9YtI1hsvFtfYeUIIIURP0q0TlEPspxXzg4GfM9YRANydHc4xnb7uYoLPpHPunR9xb5/YVZz/b+ILfDB4BD/4049xlTRMUMoH6iz//gJ8+tGZmBBCCNFz9IgEZWzqAW5L3M2JJCdrgkG+DvdpsN3QLM50F+LXTzzxCaowW8Ime3f0Yfh7X/PGZaM5xbuVM9219QXdzvUE6Wf7lC/2TMa7o6LBNSxbIubRc6SFEEKIHqZHJCht4eqnfsKAvzccEKvcTla8tI37Uted8HMsDXi47YX/wVuhkXf5IFLnV/E737Wk//0RxjgOD6QdYnfx2IP/R61qOLjWo0Xw68c3vkUIIYToLnplgvK9PaeyuiC2wmvKWpPIvoYF1XSXi4Bpb/Z6b1R7uHdz09OqdQ3G9dnP/ho/9koNV4nCXWxh21eMzW7jss9u5LwhG3gwI1pC39B0RjqaarFxNP/ihBBCiB6gVyYo6x8dTZ9/ftZm17t382z6XLyjyf2azcZXV40nkKKhAXF7wjg/WMOhgrcDvrOfD2+ZAj9b0eQ1hBBCiN6kVyYopSPAO3MSjqXrUMFgs8daoTCf/N9URvQ7CWt0JbrecPyHtjIeiE1QjLRUqiZko3RAg7DvcEG4ymw71qwJ+FbtwzxY0BYvSQghhOhRemWCEk4JU9HXTmqcD8tSqHAIzWYDw4gmLJqGXrd6slKKhKc/JyU9jX1XDsBq5I75d5loR622bCXGUzrE1mAdIIBgokYowYZndzxaWfkxkyQhhBCit+mVCUriSjtpi/ZR+M0hePPDOP67gvDpY6ns6yDl35vRfF4OfiMHNMCCtLd3Nnu9yiyD2mvGxGyzbA0XKTyS0uDgKYk4RieQ+MaGpg8UQggheqFemaDYq0GVVxD2aQRSbHhy+1HZx07Qr6H6phPxOgj5o2v3aBaY2X2iiw82kXBYdgj6m8lGmhD2amhWdG2izvRujZP94UQAhjkPcIpLCsEJ0d2tDdWyPNCf87zbybD5jn2CEF1Mr0xQjlSTrlNzcRbUrS+4/+zE+n8DKB0OnB5dZbm5FpHuylQWP/zwWnzbojOV1EnlrD/pX50clRDiRP1i18XsejcX87rX+YH/QGeHI0Sr9foExVGmSNweXb/HMnRKRjhjV0YmNjHxFCjid9dSOsR1XK0mRzOdGoFpQ0GH6esvavF5F2WtYU7irmMe92mtxS+/vhjVRHZlKQ17kQ3qitYGdscx3XcRvxjwFjM94RbHI4ToGt6rsfPD/3wPe7mO3Yr+HxeiO+r1CYq71EJfuho0DZvHgz5oDLia/g/tKYigL12NI/2ktklQHFA0Ktp6UfxeZovPe+Fso0UJyhc1gyj8IAvNhKYK0B45vNdzQKcoP5OPUocx03PixemEEB1rSeVwhv72a4Kj+1Iw6djrggnRVclggzrWKWMpuXQMEW/P+rZxdfxa7rr+RQLja1p0fNWIIHfd8CI/SGq7OjFCiI5zU/IyBr1Txq7vNr7YqBDdRa9MUKqyNcJjo1OGTYeGkZZKIM1JTaqG1bC6fKMclSaOCoXWScviVNS4WBLQWRLQ+bzWxFSN/zHKsPm4Nr6IYZkHCaRbBNItgsmqyQG/7rgg18YX0VcG1TWpyKxmSUCnwKzu7FCEaKCvzcefM5czfdA2AmkWSbaqzg5JdHEFZjV/Kcviy2DX6tbXlFLdbuW5iooK/H4/Q380D8PpYtJl63iy78ctPr/ArGZ5MJk7n/ku9grQTFBGdEDssaSsD+N4dzmazYaenETe5YMajFnpCMpGfTIVTLH45NI/NDtSv8YKUaOib75Xqwbz56cuQos0PM48qYKN055pj5B7jF8WjObVl05j5iVf1i9PIERXE1RhKq0Qft1VvxipEI352cFxrL8oh00/zWTHJX9t1+eqqLRIHLKD8vJy4uPjmz22V45BSTW89LeVojSwVyni9kc/qZWuUZ5rwzxiuRtNQdweE3tNtIXCCJpok0ahb98Loc7LNrUIGHUJhqNc57bdFzEzZWOTo/U9ugNP3To+CUbDb/6mC9TYSmb039JuMfcUQcuGEYRQY1X7hOginJodp9H8OmJCAEz07uTN70xj4PC9nR1KjF7/F9ZTZOF8azkQXRiw+jvjMR2H+z80CxI+309k7z4A1MljOXBGPNmVfaCguFNiPpqtCrb8ewjrp2Twg1OfPq5rhOMUn570F1INbxtHJ4QQoiu7wlfOFTc/2tlhNNAjEpSlX4xk2J7+vDblrwx3eFp1bnWajrp0KhDt4gl7Gh+cYSQmUnXGYAKJXbep1NruY1D4BgAczjDvTvlLo2NJprn2c+nVSzGP6NNKslXj12Wl5Jb4TtJnOL8d4TL/CmLnQAkhhGgrPSJB8e7VMQvjKJzkYXgrzw17NcpzDYxQ3ViUI+6IHgYjqMCy0FxOyvvbsGzRbh/lsqN7W5cMtTdniQYlbgAiHjerxqZjJ58Mm48is5pKS5Ftc5Nt83Fvn8bK60tzcEuMcbgYk7oOSU6EEKL99MpZPI1JXVlDxsvbsVcdHjOcsj5I+rMbiOQdjDlWaZB/WiIHLuiH2UU/o2wB+Plz3+HMz38IwJVbrmbm8z/l01pJQoQQQnR9bZ6gRCIRfvnLX5Kbm4vb7WbAgAH85je/wbIOT4NVSjF37lwyMzNxu91Mnz6dDRs6bsG8F6v8PFJ4JvoRs1hC8XZIScR70CR+j0X8HgtnfhVmRQW2vllE+qXGTM2NuCDi6cLl7xXYKzRC+R7uKxrGzv0pOEs1/pI/nYUVqU1OSxZCCCG6gjbv4rn//vv5y1/+wlNPPcXIkSNZsWIFN9xwA36/nx//+McAPPDAAyxYsICFCxcyZMgQ7rvvPmbMmMGWLVuIi4tr65Aa+Plnl+Jd78I4YoJ18Ug7ZYOSyXxuC2ZRdPCrCaAbFJ2eRU1qV81EmufdY/DC02fhrXutG14bxpf9B3PFRQ/h0WTMiRBCiK6pzVtQPvvsMy688EK++c1v0r9/fy677DJmzpzJihXRehFKKR588EHuueceLrnkEkaNGsVTTz1FTU0Nzz77bFuH0zilNV72vXvmIMd21Gt1FBucsvJa7i8e3DnxCCGEEMfQ5gnKqaeeygcffMDWrVsB+Oqrr/jkk08477zzANi5cyf5+fnMnDmz/hyn08kZZ5zBsmXLGr1mMBikoqIi5tFeNJsN9LqZOrqB7rCjtOjAWD0SnXYMoJvEdBF1J/ZKiHycxGt7x1JuBaS7RwghRJfT5gnKz372M66++mqGDRuG3W5n/PjxzJkzh6uvvhqA/Px8ANLS0mLOS0tLq993tPnz5+P3++sfOTk5bR02AJYdDs4eQHDWBNA0mDKSomvGE0zUcFQosl/ZQ8raIJqC9M+qyHr7IEawXULpEGXLUxn/yk94vTqhs0MRQgghYrR5gvLCCy/wzDPP8Oyzz7Jq1Sqeeuop/vCHP/DUU0/FHKdpsf0pSqkG2w65++67KS8vr3/s3ds+1e6UBsEEjZoUG7bsLGr7uKhN0qJTi00w8wtwFlTjKlLY8kqhsKTT1uJpC/ZK8O7T+XfxON6rsRNWZmeHJIQQQgDtMEj2pz/9KT//+c+56qqrABg9ejS7d+9m/vz5XHfddaSnpwPRlpSMjIz68woKChq0qhzidDpxOjtuPm9tisa+y/o2OkPHWreFpI0GkUgEIzGxw2JqNwpW/XsUy/qMYNkVf5RKskIIIbqENm9BqampQddjL2sYRv0049zcXNLT01m0aFH9/lAoxNKlS5k2bVpbh3NclBZdiM9TYJGyLoSt5vA+W2YG1tRRqFPGERw/oMWrH3dlWgTsFTrf2XolC0oGdHY4QgghRNu3oJx//vn89re/pW/fvowcOZLVq1ezYMECvvvd7wLRrp05c+Ywb948Bg8ezODBg5k3bx4ej4drrrnmhJ7bUjqmio5cNbTW5V6HumoOtZpoCuJ21aIvW4cjbUp9IhLu24e8k7tWBdm2YKuGA+/05e8nJfDjk7e3+v4JIYQQbanNE5SHHnqIX/3qV9x8880UFBSQmZnJjTfeyP/+7//WH3PXXXcRCAS4+eabKS0tZerUqbz33nsnVAPFCMEN//0fcERbai4Zv4o/Zqxq8fnJ60O491aQf0YyRi30+bQAiksxLZOUD/diJcVRfskEQnE9dS5ylLkpjkFlP+CXp/yH7/kbH7QshBBCtDdNKdXthnlWVFTg9/sZ9b3fYjhcmE4NZQNbjULpGhE3uM8s5PERzzDEruE5ahG83Dd+QOIag4hXQzPBFlCkrK7CtreQ/AtysdVA8jtbwTxi+m2fJPafl4bVI1YvOrbRF2/i51nvMNxux66deD9WuRVgezh6HQPFcIeOU5Oy+0II0ZtUVFokDtlBeXk58fHxzR7brROUsxO+g01zUDV9KOX9bGS/eQDldrLvG8mYTjDdit9e9S8u9cXWTZn0qx+S9vZODlw6AM9Bk/i31hGYPpKygXYsB6CiLTJHUjqYvajwqumEUILizcv/2OoVohvz+5KB/OOlb9RdW/HkVY9wiku6kYQQojdpTYLSvdsDkhKx9hWhmYAGoaxELIcOWnQlYt3UeHTPdHakb2BO4tb6loDiqRFM9wBMJ0TcOgzuhxG08B0wqcw2sOzRtXZ6MyMIjnKNe/fPJstVBsAliStalVQsrEhlZ7APP01eyWDnQczhVaA07I4ISXot0HZjed6tcfJu+WhWFPYl3lnLy4Nfb9ByJoQQovvo1i0oky7+fyS+tYnAueMoGdp0rhXIsFh+xQISjcMfiBtCAS596g7sldGfsz4sR9vwNYXfGU/Q37PHmRyv0Rdv4tncD1t8/MDFN2Df4eK1a9umFaY509dfRMm7meS8uo9wVhL/eO5hsm2+dn1OIYQQrdN7WlDq+DYU4CzyUzjBixFSJH9VhaaiM3KKxvkw3TpTP70Jo251QGurD0epxtEjIKxQmNRlJdRm+Cga6yR+l4lvR1XMMdX9fZQN7AFzi4/D56uHMPxANgD9kkt4c+gbMeNTtoaruXz1/3Ba1g4ezvqiQ2J6ojydB5+8hOpRtcy4ajUrzswh3lVKkrSeCCFEt9atE5RDa+GYew9gFBZjHzoKI6Qwtu5BhcJoloVt2HiMABgrDn+bzvqgHKOwjAMX9avvylH2unV3vt6Ny8qBsX1wF4RQqzfEPKfbNZaygb2zmJl3lwG7ojOttvb1sG9gAIBapdPPZqPEdFG7IYHPjH6Q9QV2RwTTBUY7ltv9qjqHvi/uY+OgNP6a/Rlkf1a3RxIUIYTozrp1F8/ZSdejVx0ezaon+FEZyRROTSRpfQ3G6i2UXDGemtTYLhtbDWiWIuLV6uue2AJghOpaWGwQ9mqkfxFA/3h1zLnq5LHkndI7E5QjWXYI+RX+bZC0sYbxj37Fb1KX81UI+uhBcu0+toarqbTsjHEYbTITqDEFZjXvVvdjnGsfYxy9fOCQEEJ0cb2mi0dzu+CIBMUqK0dPiCMUp1Gb5sSXm9PozJuIByA2aYm4IeKO3Rby2/EN6B/9IRwhsm9/276AbkwPg6tIw11iYiurYfH+IYz37OZSXxF2LdpatSucwJZgJo8ezCFYNz/72tRPmekJt1kcqYaXa+OLWBnUebw8k4t826RcvxBC9ADdOkEpn5xFwtsVje4rG2ijbGCfRtfTaanSoTZKh2YC4CxRJD+dT7drbmpn5QMMygf0QS2FX26+knOu/CMpdQnCzV98G9/nbjJf3I4qj/6ebv7Ht9g+fWGbx3Hn9sspWJyF9/oX+FZccZtfXwghRMfq1oUoggk6oTNGY2SlN9inNE4oOTnyGkqDiEcjNH0slf3dJ3bRHqb+PiuwV+rMXncd9xaOAOBbo77EOrsUzeFARSKoSOTEfylNxaG0br2ytBBCiFjdOkEJ+TSKxjiJpPqPeaxmcdwfYJoVLVxWNMZJdWa3vmXtylYN1YtTeenr8QDc22cDj435F8rjAk1HczrR2imLMHQLpUNYGQRV23UhCSGE6Bzduosn+608bLoTVVqO1cxxjkpF2pJCArmJFI1qXXl1zYLMjysxiqJdFDVDU1t9DQHVs8cz/p5V/DLpH7RHXvynQS+yKrsv81bPYp6mWHXq4/h0GTQrhBDdVbdOUFRFFWZlC8YbKNCCIbRIc2lMQ7ZAdH0f40AxkbyDANgzEqBBBRVxpECNgw8CBmMdFXg1GxWjkykeafDnzOW0V6PdGIeLofb9zNMUkVC3flsLIYSgm08znnLB/8P/7qaYfXr/bPadn4464nNQU6BFAB2sVsx2TVtZi/2zjdHxE1b0NqmTRpE3TWaJNEcZ0ana11y5mLuTN5Jn1uDStPrBs+2p3IrWZvHrMlZICCG6ml4zzbi5ZMNdpHAXm1T0tRFxgzpGo4c338JVasZscxyowAzFjmewFVSQvKnxImBVmTYpkw9oJhgmPL99Il/X9AFgbNxe5iTuAqDKquWsr76D1xHivRGvtmmNFElMhOj6PggYPFM4DYA4Wy3/L/0j+b/bhEfKcthUncn/S/8wZrmW3qBbJyjNid8dxP7ZRmqvGN+gvkljErZUo1ZujNlmNnKcuWMPzh17Gr2Gdf5Egv4ee0tbTVvmZwWjAVg6ZjBzzvkHACVWhITfeqjtk0j4EbPdirgJIbqmZwqnseLl6N+GkF9x67VL8Mv8g0Y9uOZsbFs9/M91H5HYy/5U9thP07KBTpx9xnVoi0bcqgN49vg5OM1PRL4MxLDvdjLwgxu4ZfxSbkrYTNLv9+I1Qji1HvsWFEKI43bt7tPZ+MRIPF4N0w0XL76Fobl5vDvsrc4OrcP0uJxVi5jYAgrTqVHTR8fqwPGskf15aFt24qhU2GqiJfV1mfEKgKNcw/OVm8VFQ8kzQzzdfxFP9P0EQ+txb0EhRCP2RarYGq5ma7iawtrDa6Nplsa2cDJ5kapmzu49girMp7UWn24fSOorm/HlWWgR8G1ysGV3w5pfPVmP+/pq5R0k/cVK0DQwDIpmDSSQ0nGtKFZtkOTXN4ARbYurmTaI4uEy6+eQHR/kcm7cT3n6soc5SWYBC9ErhJXJGUtvw7Yr+p/eCGkc6q2wV8Mdz96AbVQF60/6V+cF2UW8U5PIo9dfRn+nTt7Vw7F68bqnPSJBMdJTsfyHM/L6aUl6dFE7PQzegxYhr0Ywsf2TFbOquv7frvwa4r3RFYBNB9Sk6e1VTLVbMAKgRTT+nHcO78flAXC6bzOnS7IiRI+W0aecA1b0j19jDcuDkko6NqAuqtayY99ThJ7qJ+JxdnY4napHJChVYzMpGdb0S3GWKvxvrsMcO6jDpwir1Zvx1S2IbMvKIHBpX1SPuOvHTw/DuteGs47hAHx47hA+GPFGJ0clhGgvds3gkzGvdnYYopvpGQMAjtEiEfFo1Jw5kooBnTtyVZVXkP55FfG7G5sf1Hvt3JTBxJVX8GVQBuwIIXq3JKOKklOzKR8aV7/NckD6rL18f+InnRhZx+ve3+U1Dc1hjynK1hjLAWWDbMc8rr2ZVdXw5Qa8+iiqsrwoo93WzutWvLsNgvtS+GLAIEbYt8SUqK+xQoTrJnzbMfDovbhDVgjR4yUb1ZSM1LDVHP5wCPsUfxz4EmMcvasvvFtXkh1x0zxcYSdhr4bZTFedo0KRvugAgcF9KBzT+Adc1tLKBnVQ2ovh86LFx1F0Vl9qUiVDOaQ2RWFl1bLqzEfrizZNXnUFJduSAEgeUsyX41/qzBCFEKJdBVWYT2pdhI8YC2DXIpzqqsWpdf8JF62pJNvqNoWPPvqI888/n8zMTDRN4/XXX4/Zr5Ri7ty5ZGZm4na7mT59Ohs2bIg5JhgMctttt5GSkoLX6+WCCy5g3759rQ0FI6ioTWo+OQFAA2UzukxrhVlVTeRAPu7CCK4Shda6JYJ6LFeRhrHPxQuVA1kTDAJgM0yUTaFsCrshXWNCiJ7Nqdk5221yridY/zjbbfaI5KS1Wp2gVFdXM3bsWB5++OFG9z/wwAMsWLCAhx9+mOXLl5Oens6MGTOorKysP2bOnDm89tprPP/883zyySdUVVUxe/ZsTLN1H0B9Vla06LiwT2P/7HSKR3at7gHX+1+R+vYOqZVyBGeJxp+fuojvrrsWgCVjnmf9xQ+x/uKHWDz6hU6OTgghREdp9RiUWbNmMWvWrEb3KaV48MEHueeee7jkkksAeOqpp0hLS+PZZ5/lxhtvpLy8nCeeeIKnn36ac845B4BnnnmGnJwc3n//fb7xjW+0OJZwgpM+a0NUZtupTWq6eURp0QXsuhplmqhIBK3bdbK1o7qFHSNWNHd2anb2mFVMX/Ij+qaXsGTU650bnxBC9GCmsrgrfxJbKtMAOCtlC7cn7eiUWNp02OjOnTvJz89n5syZ9ducTidnnHEGy5YtA2DlypWEw+GYYzIzMxk1alT9MUcLBoNUVFTEPAACyXYcS9fhKZSm/1bTj3p0ke6vQ5TSMFW072tHOJ5h8yupfCWDoAoTViamsggrk7CS370QQrQVC8WrKyay++1cdr+dy5PbTuq0WNp0Fk9+fj4AaWlpMdvT0tLYvXt3/TEOh4PExMQGxxw6/2jz58/n3nvvbbA9kGqgXzqe2oSeMVu6o9SmKn5w3nu4tEj9ti016Xz42kT0UCcGdoTQBj+Dim7kvtNfZbZ3H1kLD7Blq48z77iNkuE6wbQIycsNKnNh7Q1/7pX9s0II0dbsmsGfznqOvacmAzDatbfTYmmXacaaFvt1XCnVYNvRmjvm7rvv5vbbb6//uaKigpycHFBQldnNkxPTwl6lAI1IO84gi3jBdEf7kqyMWn6UuDnmQ32NZyvv9BmPo1zH1gWWxHCUaTjK7Lw/egSjnft5NPsjLq/1UbslBT3sI+g36PPmdtynDsC8XnW5FqDOUm4F2B42yLGFSTU6tiihEKJnuMhbBd7oB0GRWc3KoJ1BdrN+dmVHadNP9/T06EJGR7eEFBQU1LeqpKenEwqFKC0tbfKYozmdTuLj42MeAKkrWjZItiszyytIfmktqV9WHvvgE9DvjN0svub3LL7m93w2/eEGLQ6jHXY+ufQPjJ2xuV3jaK0v3hzNJS/PYU8kwFMDX+eB156gfKBBykvrscrKOzu8Luef5cP49sI5/Gz/uZ0dihCiB/jNwel8e+EcHikZ1+HP3aYJSm5uLunp6SxatKh+WygUYunSpUybNg2AiRMnYrfbY47Jy8tj/fr19ce0VG1K51aGbStWoBZbYQUJX5s4KttmxKyyQc3YAMFJVQQnVXFu2gaybT6ybb5Gv1kbmk6GzceMpI3159SMC3ToatCNMYLgKNf51b7zebkylzEOF+aUCvK/MxojORFPXoCpy69nQcmAzg20ixjgPIg1ooqxcZ3XLCuE6DlGefdjDq9iRVlfflkwmiqrtsOeu9WF2qqqqti+fTsA48ePZ8GCBZx55pkkJSXRt29f7r//fubPn8+TTz7J4MGDmTdvHkuWLGHLli3ExUVL9/7whz/kP//5DwsXLiQpKYk777yT4uJiVq5ciWEce7rNoUJtQ380D8PZNv0iHVmorTk1F06ibOCJTzmK+OC56/7EOOfxLza1M1zFuc/8FEdp1+g/qRkX4OuzngTg63AVN195M9qK6O9s658msuPSv3ZmeEII0WMNXHwD9h0uXrv2jwx3eI77Oq0p1NbqMSgrVqzgzDPPrP/50NiQ6667joULF3LXXXcRCAS4+eabKS0tZerUqbz33nv1yQnAn/70J2w2G1dccQWBQICzzz6bhQsXtig5EU2rGhDhm5O/AsBrC9LP1rNmuBh7XAz75DvcPvoDronbQej/lbNv9WQG/mpVZ4cmhBA92s8m/pddI1PItHXcF9ZuXepeWlBihadUsvnUp9sslj2RKs566U6cxXqXmd0DMPTCrTzU73VSDQ8/yZvK15dnsvHuVHbO/htFZjUAKTJAVHSgneEqauqKLfUxLBmgLEQT2rUFRfQeWYaHly/5Px46eDZfvjIGukgqu+GDIZzqv5MXL/oz/5u2hPffzWa08wBBZeOkj27FDBqsmflQh484F71TuRXgint/Ssqq6KDtzT/wsfOixzs5KiG6v24+R1ccKVjh5JGyHL4Ot808YUPTGed0MtBT2CbXayu2KnAV6DxRdDovVA6j2PRRqwx0dOLjarC5IjxRNpJPa2WRI9G+Pq81eahkPIlbAlhfbcL6ahNJawyu2Xkma0MdN5hQiJ5IWlB6EN8WB3/Zej57rk7m/rQ1nR1Ou9LD8PGLE/hYmwDAm+fu5t1hb7F84nMsCri5feH3eGZKMasmyfo9ov3csfVyKhelk7NnD4fKHqb8/UtKnrIxd9H5vDpoUbPnCyGaJi0oXZBlB+9ZBUSmVraqAFlNlkXOrF3Mil/bpvGc5dtIznm7yPnmLlLP3Uc47tjndAgFWNHHlq1ZTF9/ERvCIUY4ihky42uSvTWcvu5iPpIvsqKN/bMihcHP/JDgS2mkfVGDKj+iJpOyUOEIlpI/r0KcCGlB6WI0S6EMeGDYy3xWPZhnvpzR4nNtGTW8NuTf2DBoy9zzJJfB20PfBqDUrGHyV7djr+xaf3x9220U7cpkXb8svhVXzOuD/8u1u09n1b9H8UH6SCY4VuLTXZjKIsLh2U1SIr93CKowtSqCHQO7Fn0c7zUAXsyfzMBfLEfVrcB+5Hw5zTDQ3G5ses+aRSdER5MEpYvxLduJb1si/+P9IZYN3K0YmKpv8jG86FZ+Mv1dbkvc3X5BdhO/znybz2/YwG+fvpK3D5zOwv9dwJKaIfzxw/MAUA6L/8z4MyMdMpi2pxv271vI+S/kTTPwDStlxaRnMbTWJdnDX72Vvm9HxzU5S4NgNr522MEbp3DGd7/k5pR3AZnNI8TxkgSlTsRnx5magllcWv+tqHMCiaAFgvj2QNjbuvnmtiqwVRnsOTkJ6IUJioLlVbkMtBcw2akx0O5joL2Y/9unSPmymDLLyf5gIt6dBiG/IpwEJhqlZg2rQnGYjTTJ97eXMcQuHzLdVYFZzZpgAr6dNuK+2kdVRhZlHj+LRrrROTyIOl6vZbJTazRp2Rqu5j+Vo0lepeN898uGT6Lp6GOGYnodAJSNDfNgxgokORHixEiCUqdwnAtj2ADSX96OWVzSaXFUnzqYsoE2lNSsazXNhEUvT+GNzImsvvhB/FrTLSOTz97EQzlvE6+7WFiRyR//dQlaI5N+4k87yGdjX2nHqEV7eqRkCi+9eAaGgn0XZqNs4N1rcMc/vhdzXCDLZPWFjb9nrtt4LYnfLiW5alWjM+11t4vMv+zh15nvAuDXDUBa5YQ4UZKg1LEM0Owa6O1XJc+Wk8nuq3Owmkk+AgNCxCeXonF4fGzVTj+e/V1jzIdTszFo/F62pabiWevuMrVRDtHDoIc0rCPqD5bMDFAyMplMo6Z+25d7+vJjLbqg3oaidIwgjb6W/D1JXJtwOndmvMcYRzsuNy3azOyts9ial8q4nH1sOFj3u4XD60pZNCg8aC/V+f6u2biMcP22moiDr/Zm417pIb5kR8zxxpCB7L4kFbTodX+YvJC+Nl87vioheh9JUDpQ7cBUFt/y+1ZXmTzJcRnV+1PbKarW8egO3h32Fs9nJXLfhqvRIsc+p7Ntn76w7l+HP0CcK3ysXjHqmOf6ttlZvX0U/71uP2OSvm6fAEWbMZVFwZP9GfTGZrZdM5yW1upzlGlsen1ozDZ7lWLQv9ZjVTWsK1Q2LoU1tz50XINthRAtIwnKESwblJ49AE9+DrYlaxo9xhg6gPJRSY3u868pxPw6OvbDPH0sobtLMfTD/QY5vm0k6q3/Fv6bIf/ms6zBLT7+Mv9KenMTs6NEZ+IHt4Gu0DTFA1Nf4VJfxbFPFN3a70sG8tiqM4iP1yj7xlAsR+uvEbfXImF9KUWTk4g08l/VSEhg68P9OXfIaklOhGhnvSpB0SNHNO1qYLpAHdGjo3SoytSxDAfJiX6oGyxrVQdAWeg+L6G0OCr6GUQ8YNlj+wRstSl4AbX3ADWJdu4Z+DYuLdpknGTU4NJMdh5Hi0N/O/RPWA6AS1MtaEpu/+TEpYcJxSu0uvHE9mqty7SmGEHwrj+8ivOWsRlQl6D4jCAhv8JWo6GHm7qC6E5qrBArQg7eOjAa71oXphMqs4+vS9ReY8H+fIyxiURcsd29trRUIv3T+MOUl7jI2zbVmoUQTetVCYp/VwTfx9sB0OJ8HJid3ei3pNoUjbwrh0WPsxQZb+2FSIT883MxXRpoMOOC5dzaZ0nMeSXfcrGytj//ueoUfO9v5M/nzY7usBnsOzeFYOKJD9gIZYbZ8I1H8ejH8fWwDc3ylDLiW38AoFYZXPTWj/Hu7PrfKH+U9BWXXbuK2ctuxrlSZln0BM9U9ufVa86EdA+MPLFrlQ20UZEzEssJRiD2/+vG3+bw1Jl/Z6ozDEj9HCHaW69KUPSQwqqr+KgDKKLTTfsFmz7J0jDdfcGCinFBNCP6R6so6OO5sskAHAzFs6ksDYCSag9ZNbWY1TWwfRcQLdyUnBNPKP7EP8Br8h3MyLqKq3JWdmqtE6dmZ4g9+kc6rEyU0frkK+KB4MBa9AIn7oPtNzj53bwReIwgP/Bvxae7GKJDuNpOyl6L6nT98OBJ0a2YymJhRSbP7ZuM09BQRsveQ7Za8Bw0CSbohOI1PActjGD0/Vs+0KByQrT0sKo10NSo+pa2k4dt5nQXSHIiRMfoVQkKcLjGiRX9g2QODLCjfhBlE2bG/mgqi0HvfZ91G4YDkLDDxPP6CgDSia0qeeg5HYtW0xZtHj6Ax+Cheedx27WPtcEVO09teoTN5/yV8zZdSsG72e32PKWLMngsMZ3zv72egXUtT76tDvxvryN49ShC9vZLjkT7iWBy34cX4Nthg+ktP89VYuF7ZQWuM8dRNMZJ8sf7iezZB0DJ/57MjnP+cfjg2W0bsxCi5XpMgmLZIfWs/WR5y5s8Zuc5Sey5ajwAhs1ibM5WTk3a3uTxMzadz/4lOYSG1xDnPbygi6U07Hmd28WS816IEVU3A2C6FPYRFZzTb0tdgaiOZdcMrj/5E5YMHEzB4qxmx3YEJtRwUv9dAAz0FtaV5e944y9ezycjhuDeqhHyKyaeuoUvNg3At7mR36uCx1ZM59nEyTGbfzbsv1wVV9pBEYsj3V88mH99PQln0XGUrPdrBM+dgK3GJP3TSqySsrYPUAhxwrp1gqJHQDeig1uVDj/pv4iz3CX4jpgpU2XVElZHVOBqpI+6tInCsbuWZzNowVcUXzaGQB9PzD7nEf9WGugOOyoSQVkdUxjEtmQNOUui/zaSkzhw1WDePmUEv079iDCKUF0dEJemkdLKac3H49d9NnJO3Hq+99ktaOqIAbNadHbUoaIuZw7cxt9yPj3iTB2nEcFyRH+fNFIsrU0oKLMc1FghPLqDf/b7iILsdzi59A6shDBP91/EmZWXU745vdHTveudROp+65oVfXyYMYyr4j5rp4BFY0xlUaWCvLF/NOqTxCY7WzQVLdyHFq1xpJvUv7dMh0bpYDtpK8KwfD1aXBxGQgIAlqPzCvsc+bfKo9tlnSjR62lKqS5WauvYKioq8Pv9nJN1EzbdSflJ2ZTnGtSmKGwDqlg7bSF2zSCowoz55HtE9nuOfdFGuAt0XIWKcJwW/ZBtghECW7Ui9fMSzE1Nt8i0l0MJSigBgokWfVZB3O7ouJr9Z7hZ98OHW73uyPGosUJ8FnQzb+c367tsIl649co3GeyIrlsywlFK9lGzkHaGq1gXSuUnb12LZ1/7xKkMqE1RpIws5PNxLwPRD7vPg+DSIkx0Ojh17SWUv994gnKkuH0WiZ/sJe9RH6smvdAu8YrGvV7t4463vo2jTMfWzEQaWwAy3zlAOCOBvJM9pK2sxbn5QMwxVlU1hMPkv9CfawZGZ8lNdu9kuru9suSm1Vghxnz0A1R+9MvV+MnbeXng+x0ehxDtraLSInHIDsrLy4mPj2/22G7dgoLdBroNpWugwFWoUePw8lxlGg7NpNpyom/wkbIrmoNF3BCKa914g2DisY83HdFvZZark25nJIK7yMIIahgBncR1pVjrtgKQGj+Bn+RNRdcUBhYjPAfw6rGDgs/y7Gt18bjGeHQHZ7tNPkvdyj9yM3AeNFCG4izPFoY7DiWJDadI59p9ZNrKsZwWbbkKcwwNLKeFy3Z4LrSh6ZziAupGBw1OKGRp/z64DxgNKo0eKRivUT02k/S4/e0Tq2jAVBbvBjy8WjQBzwH9mC1tmqWwCosxfC7Ag60yRCT/IABGSgqBif1RGlhOjesGvcecxF3tGv+eSBXLAjkA2DWTWZ6i+pl4r1f7eKd0NPpuN87i6N+btfszeT41EQCPHmSWp1Lqrohep1snKPu/mYnhdHHkGm+eAzoPLLwCiDbD576ah7U7+kFiTR5O3rSeN7XULK/A//IqGDuE/dNjM1LXf1ezbUld14TDwZLLpxJMOpx0KQ0qrn2dH/hjv12eiLuTN3L7hWsZ/eGN2Hd2jfLwYZ/izQseZIjdAU2Me3k8ZwnlWe8y5fXb8e5u+sOg/IxaFp/xKB6t6WuJtlVqBfjRGzfjPqg3umZSa1SeNoBX/m8BzrpWRZ/mpN0S4zr/V3g6/33pJAAiHsXA7zzIuLp+4rte/g4D563HfznU1tWAdK7ycd+aqwGoTVZMvvL3ZEgpfdHLdOsExTKgwZcKRf34B01BzdAU7FkJANSkOWlPVbk+3L4xANjKautbMTqCMk1sB8vos8aJXlRe/wVTWQoViA7w1UJhktfVEPEe8WvXNRa4LmJ+VohBuQeZmbaRn9aVdM+LVHHq0ttI71POp2NebXEshqbj0RxcNOIrNmWn06cFU5BtGJw1fiNbcqMl/fPXpdH3vyEqch0E/W0wy0YDl2Y2+y3Urhl4NPvhRZCaOm6bm3HlP+bm0z+ov1ei/SwoGcCbeaNxVOj1hQEb4y5SxO0NUTbQUb/Ypl5YRupqL8bBMiyPhz0/GktkXBXJurtDuj0POcu/kS+m9wfAaw+RZoRZGdT49e4LsVVphCcOji0MVzfOCcBeqXH15m8zK2MDP0ve1mExC9HZunWCcixKg+LhdjqqbkF5f4Py/tEqrr4DTvwbtA4bNAsQ2Z+HfX8eTRV0VaaJ9sX6Bnej3/tg9EnmwBWDeOK0BH50ymZ0dHZEPAz9XQ1Fk9OoGRXCqdla9Uf9jxmrIAOOtey8qSwimDzR9xPoG902qOR6HEvW40ya1HSCcuSKisegWpPjaBz+Qq1osIhgyloT7xur+Mdz05gzbas0vbcjU1n8Y+vJaMv8x/xj5SmIYPtwDa7kyQSSo7/wSP5BbPkHiWg6tow0fnbdi1wbX0R7t5gc7ZueWr55RJIfVm7+WT6M3W/lYig4OKnplkZbDRS/n8nCyXH87BRJUETv0aMTFNFyVlk5me/kEVkRz1kv30bRGI1wokXyVB1bQDHj9h/h+H4eH478d5s/97gvv43vpXgu/8V73J6049gnEJ1WfsYlqxjpbdk4kDg9QLbt2C1oTs3G72c8R37ED8DjW0/F+iSx0WNz/qRz8ru38odf/LVTBlb2dAsrUrn3owtxFNhaVEOoop8N51WTqU3S0czYrHLXvVOYNmM953p3c6yEub3tDFdx8R/uwlmmIO3Yx4f8imsuWMoZvs3tH5wQXUirE5SPPvqI3//+96xcuZK8vDxee+01LrroIgDC4TC//OUvefvtt9mxY0d0ps055/C73/2OzMzM+msEg0HuvPNOnnvuOQKBAGeffTaPPvoo2dntV6yro1mGhp6chKqqxgrUHvuETqbCEcwde9B2gO8LsFdOoDLHhmUDR6WF+91VbJ09pslS4ptCNVSqaNtMHz1Irj3aXx5WJhtCEeL0MAPtjfehV5e7ydxRQ0nk8AdHvC8AY4dSlaETTGnYCmXZFX0clQxzHuBsdzPt/jGO3ZJmaHrdwoLRisPvJJSRV5lATZqG6YnGUVVr4B09FD0Qwb8zxCfVQ8i0rWSIveeNb+poX4erKLaiieT7JSPwbbVjq41WgW6JkE9DDyn0MKhh/dHD0fdG4sRCnuz7MV0hOVkaGEDa8iq0oEkg1d9s617IrwilRvhh0vI2GcguRHfS6gSlurqasWPHcsMNN3DppZfG7KupqWHVqlX86le/YuzYsZSWljJnzhwuuOACVqw4XEBszpw5vPnmmzz//PMkJydzxx13MHv2bFauXIlh9Izm8kAfjQNXDCJ1RTXaF+s7O5xWc3+4Dveh34VpYlmqQVfHIWFl8s0lt+HcEf1g8U0pYvmEFwHYFwlw2St3YuTUsOW0fzZ6/udn/5nC6Tq5NoNDM2o+GP8kB17R8GgmjRV63R3x8Msf/oB30k/n43l/breaEZv2pjP0hY0EXkrg6aH/AiCsoObO6L3JN3384I3v82TKyWw964kOHdfQE533+c3oG6KJrB6JdsT0WVWNsaaVXRuD+nLxCx9xkjvaIpdjWMDxlRtoSxeu/j6RlYnYp0Z/brbrUYOrLviIW5K+lORE9EqtTlBmzZrFrFmzGt3n9/tZtGhRzLaHHnqIKVOmsGfPHvr27Ut5eTlPPPEETz/9NOeccw4AzzzzDDk5Obz//vt84xvfOI6X0fUoHUwnKJve0mESXYoVCgOxJWHjVrkYHX9Ng2OV0kj4woEvLzr6pSSUwuhw9LjagIOU1RDY62O0o+G5AN8etLx+8N/WcDVXr72BBHctIxPyWFOcTVXQweT0PTj1w6NrioI+TLeOvcZi4hfX0zexlCHxBfw4ZUl9601z/lWZzNqaHH6a8mmzheyUqWNVB/DaPcTpBvcXTqXKPNxVVBZ2Y6/QiVhOfnzgZC5OWtGKFh1x9O87nO/BWxN7jB62UKEQ+pABWM5oImrkF9dPG7alpXLg0oExtYpCiXC2Z2uTrXYd7cUqP/9vwzep3erHFaDJ9Z8iHjBHVKFpoGmK032bJTkRvVa7j0EpLy9H0zQS6io1rly5knA4zMyZhxe4yczMZNSoUSxbtqzRBCUYDBIMHq7dUVFR0d5hi0akP/QFPHTs47LeaLjND02e+/iCs/nZFdEEZVH1MNJuKMYckMmHp2eT804xvv35fHnNeCJHfwEeBN58i+wrtxE8fTTvT+7P5Bt2kmsvPmaM89afi1rj5+IbVpLSwka7nWGd1/59aoMCYXbAXqmz9IWJ7Dk/kbMHv9uyC4ro7/v6IsyBWXx4enaTHTCazUbBKcn1A6bTlzvQ6xKU8MAM3vzZAw0KADZWc6ez3LfxPDIv3kjgoimUDmr6z24wxWT9qX/v9NXKhegK2jVBqa2t5ec//znXXHNNfcW4/Px8HA4HiYmxAw/T0tLIz89v9Drz58/n3nvvbc9Quy3rlDFEPAaupevRM9OpGJtG/Oo8Inu6TxGxAa8EGbcjuq6Qo1yRXLUaY/dBsj6wIL8QK1BL+kclKPsRmYRN5+DkOIJ+nerzxxNIbH3Xih6Caz7+PjbnEcXbDMVTk/7BlLpv6j+f+g6PvHQGoUIbl+69CVczC1+LlguqMMNfuZU+K3QSq1Y1eoy7UJH8VTna/kLMUJjUT4oIZsRTON7J15c5iPvJMADS44pI6mIf6KayGPTmTfi+tlEzNoC+30XgoilUpzWeDSsDsmbu4Zy0TTg1mbsgBLRjghIOh7nqqquwLItHH330mMcrpdC0xjtD7r77bm6//fb6nysqKsjJyWmzWLs6zW5Ds9lQwWCDactVOU5CcRruVXGYyXFUZhv4tnV+X3tr6J+uJe2I5XksgMJiKCyuXxlabdyGZhhorrqic04H+vghhOI1ygbU/dFXcDDsp8DcE9MsXmRWU3vUig7hsIHTBO86F5pF/QKHygafDB9Kf9s6Ug0vP/Af4AdTniP3P9/Hv96OZW9+3EB12MG+SBVphlumHzcjrExyFilc//m8qaFN2AMKtu3Gqms91UorsPndRDxOxo3dwauDjuxO7vwEJaxMDpoBgnXrPmUt0ohfvIm8bw0n4qJBy4kyot3AaNFB3z/pt4hzPUE6egq0EF1VuyQo4XCYK664gp07d7J48eKYevvp6emEQiFKS0tjWlEKCgqYNm1ao9dzOp04ne1bZK0rM08aQelgF2lv7SRysLDB/rBXI//SQcTvjZD5/DZUZWUnRNn+9AF9OTi9DxBNEiKe2ExBM+Gvr32DhzLPYs03HsKvR2vSTF/+A0JbYyvs2ioPn+suUiT/92uwFJqhs1A7l4dGn8m6GY/g012YyiJphY20xXkcmJWB2cxb8cBH2Uxf/VMevvgfdR824nhVZerUfnsMqZ+Xoe87iPflCN9Oe5l0Wzk5RpCu1IUD8FGtgxtfu4XEjRopX5YSt3sj2Jr+E1udG+G5mY9hr6s+N9wOXSHREqKraPME5VBysm3bNj788EOSk5Nj9k+cOBG73c6iRYu44opoSfq8vDzWr1/PAw88cMLPb4TAXWQRitNave5Oewgl2PEO7Ifan49Ve3wfWHoggqNSoayGtTbcBWEOTZ+1l4cxi0tOJNwuTdltTf5OHZUKR6WCA1BR4yA802JDKMDi6mGEN8STvOno7+mHf3YXRTALo+NWNF0jaVNfTJeT8DmH73f5EIWzIg13sUXYrTW5RpOtBoygxjMFJ1OW/BWX+4plZk8LGJW1+PY3PgKlfLif8IQEfpR+KOnruOKLrVFtOXEW6zjLTbTaIKSloHzumKU4INpyUt0/Qlb/Ik5yGchyCUI0rtUJSlVVFdu3H16xd+fOnaxZs4akpCQyMzO57LLLWLVqFf/5z38wTbN+XElSUhIOhwO/38/3vvc97rjjDpKTk0lKSuLOO+9k9OjR9bN6ToSzVOF/ZTXhU0dRML7zW11Kh9goG5hO9ishrP15x3eRFRuJX6VhNlKV1r54DQ49+mHZkVVru5rErUFsn0Snc/tOGkntDxXzD8ziq9dGMPDtYqwtTZekP/K+KUvh/O8qMgJjqb0pmqAYms6Wqx/l9QsS+Me5Z2H6vew/29/k9TQTvnptBMuyhzH70gfxaV1jPaKuzNy8nfgtjRfp+/qZ0Ww644lu02VW0d+gon9G/c9HdwmaTnjh3EeY6DCQ7hwhmtbqBGXFihWceeaZ9T8fGhty3XXXMXfuXN54IzqFY9y4cTHnffjhh0yfPh2AP/3pT9hsNq644or6Qm0LFy5skxoomlLQhT6olQaaBugn1prTXPLRmxOTo6mJw6js6+TCr75LSb4fnwKUavU9shfVcN7q/+Hy3NX8ImULhqYzznmATb9OQi9w4j54rEDAUaZzztpvc2nf1bJmz1Gcmp2y71VScd5kAFKWGyQ9+Xmjx2rQLZKTkY4Chs7aRshq/s+qxxaivy2Eocn0YSGa0+oEZfr06SjVzIdlM/sOcblcPPTQQzz0UAvmrLaS0jQ0u611a6+Ibu/Q771kqJdAHw2WphCn6hZca8F78mh6aQXq/f48d9ZEbklag4FGts3Jjhn/4GcHx/HWC4fHS2lNLH5kr4Tqxam8OH0ityRsavVaRj2ZXTP4aspz9T8PSbiOlOfdjR6rG91jGYGBdt9RA3ebI8mJEMeiqZZkFF1MRUUFfr+foT+ah+GMbT7Xw+AsU4S9WsO6GZ1ENyH71b1E9h7o7FC6NX3EYPadm9ToPnu1wlYDwQStvgiWf5eJf0UeVmFxq5cb0Ow2jPRUzKR4atM9HJxiJziwlk+n/xkT+CQQnUVWGInn4dfOw1HWdEYcjleEki3mnfMiV8WVtiqO3uLzWpOXSic3uu/yxOV1YzWEEN1dRaVF4pAdlJeXx0ygaUyPm3Bv2aNl5rsaK8GHUZtcPxizOZquYaQkoyIRzJKy9g+uC9K9HvT4OKzSMlQohJGSTNjf9FiOsFcjfNSXUluN1ep6MJquweihmD5H/RRne1UEI2hHhaKtH9k2X32iUWTuY0HWOUQ8zQ/a1EIaISUfsk05yWVwUkbj9VBkEKkQvVOPS1C6IsuAA2cl4ipOIPGFclS4iT6BOprTycHZA3BWKHwvL++gKLsWNaw/+6bFkbWoCO3AQfIvHNBgWnF70JxO0h/bw7zMd2K2uzQdu6bj02OntqYYXtbNeASzyWoeh/k0JzIoUgghWqZHJyiagvjdJpoF5f2NBtP9OpJlQNgL1qTh2A+UEtm9r+mDTZOEr4MYgeYTmfagpo5i74zO7x8Pppok9C1mS78kbDXJ2IZX4LRHWvVdevewOFyTprbqeZVNcWPy02Q0KJveNJ8us3SEEKKt9egEBQX+lfkQClPRt3+nJigQ7YbIO9lDyno7jmYSFCsUxlj6VQdGdljeqV423XTsyr8dZlJnByCEEKIz9OwEpRVsp5Vwae6aDnmu3YFkdtyZ1SHP1VpXpi7u7BCEEEKI7pmgHJp4ZAabn5mhWRCxgmCFMYO1WGbTx07wfc2PnBvbMsymOYGEjnmq41HRMyvlCyGE6GQVVdGyAS2ZQNwtpxnv27evVy0WKIQQQvQke/fuJTs7u9ljumWCYlkWW7ZsYcSIEezdu/eYc6lFyxxaJVruaduRe9r25J62PbmnbU/uaeOUUlRWVpKZmYmuNz8wtFt28ei6TlZWdAxHfHy8/PLbmNzTtif3tO3JPW17ck/bntzThvz+ptcyO5IUZRBCCCFElyMJihBCCCG6nG6boDidTn7961/jdDo7O5QeQ+5p25N72vbknrY9uadtT+7pieuWg2SFEEII0bN12xYUIYQQQvRckqAIIYQQosuRBEUIIYQQXY4kKEIIIYTocrptgvLoo4+Sm5uLy+Vi4sSJfPzxx50dUrcwd+5cNE2LeaSnp9fvV0oxd+5cMjMzcbvdTJ8+nQ0bNnRixF3PRx99xPnnn09mZiaapvH666/H7G/JPQwGg9x2222kpKTg9Xq54IIL2Lev6RWue7pj3dPrr7++wfv2pJNOijlG7mms+fPnM3nyZOLi4khNTeWiiy5iy5YtMcfIe7V1WnJP5b3adrplgvLCCy8wZ84c7rnnHlavXs1pp53GrFmz2LNnT2eH1i2MHDmSvLy8+se6devq9z3wwAMsWLCAhx9+mOXLl5Oens6MGTOorJQVBA+prq5m7NixPPzww43ub8k9nDNnDq+99hrPP/88n3zyCVVVVcyePRvTbGZFyx7sWPcU4Nxzz41537799tsx++Wexlq6dCm33HILn3/+OYsWLSISiTBz5kyqq6vrj5H3auu05J6CvFfbjOqGpkyZom666aaYbcOGDVM///nPOymi7uPXv/61Gjt2bKP7LMtS6enp6ne/+139ttraWuX3+9Vf/vKXDoqwewHUa6+9Vv9zS+5hWVmZstvt6vnnn68/Zv/+/UrXdfXuu+92WOxd1dH3VCmlrrvuOnXhhRc2eY7c02MrKChQgFq6dKlSSt6rbeHoe6qUvFfbUrdrQQmFQqxcuZKZM2fGbJ85cybLli3rpKi6l23btpGZmUlubi5XXXUVO3bsAGDnzp3k5+fH3Fun08kZZ5wh97aFWnIPV65cSTgcjjkmMzOTUaNGyX1uxpIlS0hNTWXIkCF8//vfp6CgoH6f3NNjKy8vByApKQmQ92pbOPqeHiLv1bbR7RKUoqIiTNMkLS0tZntaWhr5+fmdFFX3MXXqVP75z3/y3//+l7/97W/k5+czbdo0iouL6++f3Nvj15J7mJ+fj8PhIDExscljRKxZs2bxr3/9i8WLF/PHP/6R5cuXc9ZZZxEMBgG5p8eilOL222/n1FNPZdSoUYC8V09UY/cU5L3alrrlasYAmqbF/KyUarBNNDRr1qz6f48ePZqTTz6ZgQMH8tRTT9UP5JJ7e+KO5x7KfW7alVdeWf/vUaNGMWnSJPr168dbb73FJZdc0uR5ck+jbr31VtauXcsnn3zSYJ+8V49PU/dU3qttp9u1oKSkpGAYRoNMs6CgoME3AXFsXq+X0aNHs23btvrZPHJvj19L7mF6ejqhUIjS0tImjxHNy8jIoF+/fmzbtg2Qe9qc2267jTfeeIMPP/yQ7Ozs+u3yXj1+Td3Txsh79fh1uwTF4XAwceJEFi1aFLN90aJFTJs2rZOi6r6CwSCbNm0iIyOD3Nxc0tPTY+5tKBRi6dKlcm9bqCX3cOLEidjt9phj8vLyWL9+vdznFiouLmbv3r1kZGQAck8bo5Ti1ltv5dVXX2Xx4sXk5ubG7Jf3ausd6542Rt6rJ6BzxuaemOeff17Z7Xb1xBNPqI0bN6o5c+Yor9erdu3a1dmhdXl33HGHWrJkidqxY4f6/PPP1ezZs1VcXFz9vfvd736n/H6/evXVV9W6devU1VdfrTIyMlRFRUUnR951VFZWqtWrV6vVq1crQC1YsECtXr1a7d69WynVsnt40003qezsbPX++++rVatWqbPOOkuNHTtWRSKRznpZnaq5e1pZWanuuOMOtWzZMrVz50714YcfqpNPPlllZWXJPW3GD3/4Q+X3+9WSJUtUXl5e/aOmpqb+GHmvts6x7qm8V9tWt0xQlFLqkUceUf369VMOh0NNmDAhZpqXaNqVV16pMjIylN1uV5mZmeqSSy5RGzZsqN9vWZb69a9/rdLT05XT6VSnn366WrduXSdG3PV8+OGHCmjwuO6665RSLbuHgUBA3XrrrSopKUm53W41e/ZstWfPnk54NV1Dc/e0pqZGzZw5U/Xp00fZ7XbVt29fdd111zW4X3JPYzV2PwH15JNP1h8j79XWOdY9lfdq29KUUqrj2muEEEIIIY6t241BEUIIIUTPJwmKEEIIIbocSVCEEEII0eVIgiKEEEKILkcSFCGEEEJ0OZKgCCGEEKLLkQRFCCGEEF2OJChCCCGE6HIkQRFCCCFElyMJihBCCCG6HElQhBBCCNHlSIIihBBCiC7n/wNv/u9zB35tswAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds_xr = rasterio.open(\"./data/test_klein_gebied/openEO.tif\").read(1)\n",
    "plt.imshow(ds_xr)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}