EOdal_Multitemporal_Composites Loop.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2c150542-eea8-4388-9c86-03f76fc0945e",
   "metadata": {},
   "source": [
    "# Temporal Aggregation of Single Landsat Scenes to create Multi-Temporal Landsat Composites\n",
    "\n",
    "In remote sensing, we usually work on single images (\"scenes\") acquired on a given date over a given area. While this sufficient or even necessary for many applications such as the studying of vegetation dynamics, some applications actually require information combined from many scenes, i.e., multiple \"snapshots\" of an area over time.\n",
    "\n",
    "Single scenes tend to be noisy mostly because of atmospheric disturbances and instrumentation errors. Thus single scenes might contain significant parts of no-data or poor-quality pixels due to instrument failures, clouds (in the case of optical imagery) or shadows. Taking multiple scenes acquired at different points in time and aggregating them over time by a statistical measure such as the mean or median is therefore in some cases the preferred choice.\n",
    "\n",
    "Multi-temporal composites are used for many different applications, including\n",
    "\n",
    "* inputs for classification tasks such as land use/ land cover mapping\n",
    "* generation of equidistant time series\n",
    "* generation of basemaps\n",
    "\n",
    "In this notebook, you will learn how to generate a multi-temporal mosaic\n",
    "The steps are:\n",
    "\n",
    "1. Fetch the Landsat scenes overlapping the study region\n",
    "2. Aggregate the Landsat scenes into a single multi-temporal composite using the median, mean, mnedian absolute deviation (mad) and coefficent of variance (cov) as well as some spectral indices on the basis of the median bands.\n",
    "\n",
    "This code was written by Finn Timcke and is based on a notebook by Dr. Lukas Valentin Graf (https://github.com/terensis-dev/eodal_tutorials/blob/main/01_general/04_Multitemporal_Composites.ipynb).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d810f918",
   "metadata": {},
   "source": [
    "# Imports and Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4dd12ebc-7150-4ace-bd3d-23a3ed8da625",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# we have to install the development version for the time being\n",
    "#%pip uninstall eodal -y\n",
    "#%pip install git+https://github.com/lukasValentin/eodal@landsat-dev\n",
    "#%pip install openpyxl\n",
    "\n",
    "# required imports\n",
    "import os\n",
    "from datetime import datetime\n",
    "import rasterio\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.dates import MonthLocator, DateFormatter\n",
    "from pathlib import Path\n",
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "from eodal.config import get_settings\n",
    "from eodal.core.sensors import Landsat\n",
    "from eodal.core.scene import SceneCollection\n",
    "from eodal.mapper.feature import Feature\n",
    "from eodal.mapper.filter import Filter\n",
    "from eodal.mapper.mapper import Mapper, MapperConfigs\n",
    "from eodal.core.band import Band\n",
    "from eodal.core.raster import RasterCollection\n",
    "\n",
    "\n",
    "# the function to mask clouds and shadows:\n",
    "def prepocess_landsat_scene(\n",
    "        ds: Landsat\n",
    ") -> Landsat:\n",
    "    \"\"\"\n",
    "    Mask clouds and cloud shadows in a Landsat scene based\n",
    "    on the 'qa_pixel' band.\n",
    "\n",
    "    NOTE:\n",
    "        Depending on your needs, the pre-processing function can be\n",
    "        fully customized using the full power of EOdal and its\n",
    "        interfacing libraries!\n",
    "\n",
    "    :param ds:\n",
    "        Landsat scene before cloud mask applied.\n",
    "    :return:\n",
    "        Landsat scene with clouds and cloud shadows masked.\n",
    "    \"\"\"\n",
    "    ds.mask_clouds_and_shadows(inplace=True)\n",
    "    return ds\n",
    "\n",
    "# functions to plot:\n",
    "def scatterplot_temporal_distribution(metadata_gdf, save_path, year):\n",
    "\n",
    "    # Convert lists in 'instruments' column to strings\n",
    "    metadata_gdf['instruments_str'] = metadata_gdf['instruments'].apply(lambda x: ', '.join(map(str, x)))\n",
    "    #filtered_metadata_gdf['instruments_str'] = filtered_metadata_gdf['instruments'].apply(lambda x: ', '.join(map(str, x)))\n",
    "    \n",
    "    # Convert 'sensing_date' to datetime type if it's not already\n",
    "    metadata_gdf['sensing_date'] = pd.to_datetime(metadata_gdf['sensing_date'], errors='coerce')\n",
    "    #filtered_metadata_gdf['sensing_date'] = pd.to_datetime(filtered_metadata_gdf['sensing_date'], errors='coerce')\n",
    "\n",
    "    # Count number of acquisitions \n",
    "    acquisitions = metadata_gdf['sensing_date'].nunique()\n",
    "\n",
    "    # Plotting using the scatter plot with figsize parameter\n",
    "    fig, ax = plt.subplots(figsize=(16, 4))\n",
    "\n",
    "    # Plotting the original data\n",
    "    metadata_gdf.plot.scatter(x='sensing_date', y='instruments_str', ax=ax, color='blue', label='Acquisitions')\n",
    "\n",
    "    # Plotting the filtered data with different color\n",
    "    #filtered_metadata_gdf.plot.scatter(x='sensing_date', y='instruments_str', ax=ax, color='red', label='Acquisitions of interest')\n",
    "\n",
    "    # Set x-axis limits to cover the entire year\n",
    "    start_date = pd.to_datetime(metadata_gdf['sensing_date'].min()).replace(day=1, month=1)\n",
    "    end_date = pd.to_datetime(metadata_gdf['sensing_date'].max()).replace(day=31, month=12)\n",
    "    ax.set_xlim(start_date, end_date)\n",
    "\n",
    "    # Customize x-axis ticks to show all months\n",
    "    months = pd.date_range(start=start_date, end=end_date, freq='M').strftime('%b')\n",
    "    ax.set_xticklabels(months)\n",
    "\n",
    "    # Adding labels and title\n",
    "    plt.xlabel('Sensing Date')\n",
    "    plt.ylabel('Instruments')\n",
    "    plt.title(f'Temporal distribution of all {acquisitions} acquisitions in year {year}')\n",
    "    \n",
    "    # tight layout\n",
    "    plt.tight_layout()\n",
    "\n",
    "    # create plot\n",
    "    plt.show()\n",
    "\n",
    "    # Save plot to the specified path\n",
    "    fig.savefig(save_path)\n",
    "    \n",
    "    # Close the current figure\n",
    "    plt.close()\n",
    "\n",
    "def plot_geo_barplot(metadata_gdf, save_path):\n",
    "    \n",
    "    # Convert lists in 'instruments' column to strings\n",
    "    metadata_gdf['instruments_str'] = metadata_gdf['instruments'].apply(lambda x: ', '.join(map(str, x)))\n",
    "    \n",
    "    # Convert 'sensing_date' to datetime format\n",
    "    metadata_gdf['sensing_date'] = pd.to_datetime(metadata_gdf['sensing_date'])\n",
    "    \n",
    "    # Extract the year from the 'sensing_date' column\n",
    "    metadata_gdf['year'] = metadata_gdf['sensing_date'].dt.year\n",
    "\n",
    "    # Create a complete list of years\n",
    "    all_years = range(metadata_gdf['year'].min(), metadata_gdf['year'].max() + 1)\n",
    "\n",
    "    # Group by year and instrument_str, count the number of rows, and unstack\n",
    "    grouped_data = metadata_gdf.groupby(['year', 'instruments_str']).size().unstack(fill_value=0)\n",
    "\n",
    "    # Reindex the grouped data with all years\n",
    "    grouped_data = grouped_data.reindex(all_years, fill_value=0)\n",
    "\n",
    "    # Plotting with figsize parameter\n",
    "    fig, ax = plt.subplots(figsize=(16, 9))\n",
    "    grouped_data.plot(kind='bar', stacked=True, colormap='viridis', ax=ax)\n",
    "    \n",
    "    ax.set_xlabel('Year')\n",
    "    ax.set_ylabel('Number of Scenes')\n",
    "    ax.set_title('Number of Scenes found by Year and Instrument')\n",
    "    ax.legend(title='Instrument', bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "    \n",
    "    # tight layout\n",
    "    plt.tight_layout()\n",
    "\n",
    "    # create plot\n",
    "    plt.show()\n",
    "    \n",
    "    # Save plot to the specified path\n",
    "    fig.savefig(save_path)\n",
    "    \n",
    "    # Close the current figure\n",
    "    plt.close()\n",
    "    \n",
    "# Median and Median Absolute Deviation (mad) by spectral band\n",
    "def plot_bands_median_mad(rc, colormap_dict, save_path):\n",
    "    # Create subplots\n",
    "    fig, axs = plt.subplots(6, 2, figsize=(16, 75), sharey=True)\n",
    "\n",
    "    bands = ['blue', 'green', 'red', 'nir08', 'swir16', 'swir22']\n",
    "    stats = ['median', 'mad']\n",
    "\n",
    "    for i, band in enumerate(bands):\n",
    "        for j, stat in enumerate(stats):\n",
    "            column_name = f'{band}_{stat}'\n",
    "            colormap_key = f'{column_name}'\n",
    "            \n",
    "            # Plot band_median and band_mad\n",
    "            rc.plot_band(column_name, colormap=colormap_dict[colormap_key], ax=axs[i, j])\n",
    "            axs[i, j].set_title(f'{column_name}', fontsize=15)\n",
    "\n",
    "    # Adjust layout to prevent overlapping\n",
    "    plt.tight_layout()\n",
    "\n",
    "    # create plot\n",
    "    plt.show()\n",
    "    \n",
    "    # Save plot to the specified path\n",
    "    fig.savefig(save_path)\n",
    "    \n",
    "    # Close the current figure\n",
    "    plt.close()\n",
    "\n",
    "# Mean and coefficient of variation (cov) by spectral band\n",
    "def plot_bands_mean_cov(rc, colormap_dict, save_path):\n",
    "    # Create subplots\n",
    "    fig, axs = plt.subplots(6, 2, figsize=(16, 75), sharey=True)\n",
    "\n",
    "    bands = ['blue', 'green', 'red', 'nir08', 'swir16', 'swir22']\n",
    "    stats = ['mean', 'cov']\n",
    "\n",
    "    for i, band in enumerate(bands):\n",
    "        for j, stat in enumerate(stats):\n",
    "            column_name = f'{band}_{stat}'\n",
    "            colormap_key = f'{column_name}'\n",
    "            \n",
    "            # Plot band_median and band_mad\n",
    "            rc.plot_band(column_name, colormap=colormap_dict[colormap_key], ax=axs[i, j])\n",
    "            axs[i, j].set_title(f'{column_name}', fontsize=15)\n",
    "\n",
    "    # Adjust layout to prevent overlapping\n",
    "    plt.tight_layout()\n",
    "\n",
    "    # create plot\n",
    "    plt.show()\n",
    "    \n",
    "    # Save plot to the specified path\n",
    "    fig.savefig(save_path)\n",
    "    \n",
    "    # Close the current figure\n",
    "    plt.close()\n",
    "\n",
    "# RGB of the medians and Skewness by band\n",
    "def plot_bands_skew(rc, colormap_dict, save_path):\n",
    "    # Create subplots\n",
    "    fig, axs = plt.subplots(6, 2, figsize=(16, 75), sharey=True)\n",
    "\n",
    "    skew_bands = ['blue', 'green', 'red', 'nir08', 'swir16', 'swir22']\n",
    "    stats = 'skew'\n",
    "\n",
    "    for i, skew_band in enumerate(skew_bands):\n",
    "                # Plot RGB and band_skew\n",
    "                rc.plot_multiple_bands([\"red_median\", \"green_median\", \"blue_median\"], ax=axs[i, 0])\n",
    "                axs[i, 0].set_title(\"red_median, green_median, blue_median\", fontsize=15)\n",
    "                \n",
    "                rc.plot_band(f'{skew_band}_{stats}', colormap=colormap_dict[f'{skew_band}_{stats}'], ax=axs[i, 1])\n",
    "                axs[i, 1].set_title(f'{skew_band}_{stats}', fontsize=15)\n",
    "\n",
    "    # tight layout\n",
    "    plt.tight_layout()\n",
    "\n",
    "    # create plot\n",
    "    plt.show()\n",
    "    \n",
    "    # Save plot to the specified path\n",
    "    fig.savefig(save_path)\n",
    "    \n",
    "    # Close the current figure\n",
    "    plt.close()\n",
    "\n",
    "# RGB of the medians and the Spectral indices\n",
    "def plot_rgb_and_spectral_indices(rc, rc_clean, spectral_indices, colormap_dict, save_path):\n",
    "    # Create subplots\n",
    "    fig, axs = plt.subplots(len(spectral_indices), 2, figsize=(16, 12.5*len(spectral_indices)), sharey=True)\n",
    "\n",
    "    # Iterate over spectral_indices\n",
    "    for i, index in enumerate(spectral_indices):\n",
    "        # Plot RGB\n",
    "        rc.plot_multiple_bands([\"red_median\", \"green_median\", \"blue_median\"], ax=axs[i, 0])\n",
    "        axs[i, 0].set_title(\"red_median, green_median, blue_median\", fontsize=15)\n",
    "\n",
    "        # Plot spectral index\n",
    "        rc_clean.plot_band(index, colormap=colormap_dict[index], ax=axs[i, 1])\n",
    "        axs[i, 1].set_title(index, fontsize=15)\n",
    "\n",
    "    # tight layout\n",
    "    plt.tight_layout()\n",
    "\n",
    "    # create plot\n",
    "    plt.show()\n",
    "    \n",
    "    # Save plot to the specified path\n",
    "    fig.savefig(save_path)\n",
    "    \n",
    "    # Close the current figure\n",
    "    plt.close()\n",
    "\n",
    "# RGB of the medians and numbers of observation contributing to the composite\n",
    "def plot_rgb_and_scene_count(rc, colormap_dict, save_path):\n",
    "\n",
    "    # Create subplots\n",
    "    fig, axs = plt.subplots(1, 2, figsize=(16, 20), sharey=True)\n",
    "\n",
    "    # Plot RGB\n",
    "    rc.plot_multiple_bands([\"red_median\", \"green_median\", \"blue_median\"], ax=axs[0])\n",
    "    axs[0].set_title(\"red_median, green_median, blue_median\", fontsize=15)\n",
    "\n",
    "    # Plot scene_count\n",
    "    rc.plot_band('scene_count', colormap=colormap_dict['scene_count'], ax=axs[1])\n",
    "    axs[1].set_title('scene_count', fontsize=15)\n",
    "\n",
    "    # Adjust layout to prevent overlapping\n",
    "    plt.tight_layout()\n",
    "    \n",
    "    # create plot\n",
    "    plt.show()\n",
    "    \n",
    "    # Save plot to the specified path\n",
    "    fig.savefig(save_path)\n",
    "    \n",
    "    # Close the current figure\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f38d6e7-5729-4ed7-91cc-70c859866ccf",
   "metadata": {},
   "source": [
    "# User inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "835fbaf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ------------------------- Spatial Filters ---------------------------\n",
    "shapefile_name = \"study_areas_laura_COD.shp\" #choose a shapefile from the 'shapefiles' folder with the extension '.shp' !\n",
    "\n",
    "# ------------------------- Timeinterval Filters ---------------------------\n",
    "year_start = 1984\n",
    "year_end = 2024\n",
    "\n",
    "# ------------------------- Metadata Filters ---------------------------\n",
    "excluded_sensor = ['etm+']  # all available: ['tm', 'etm+', 'oli']\n",
    "max_cc = 100 #maximal cloud coverage of the scenes in percent \n",
    "\n",
    "# ------------------------- Spectral Indices ---------------------------\n",
    "# spectral indices to be added to the rastercollection as new bands use # to unselect\n",
    "spectral_indices = [\n",
    "    'NDVI',\n",
    "    'EVI',\n",
    "    'CI_GREEN', \n",
    "    'MSAVI',\n",
    "    'NDWI',\n",
    "    'BSI'\n",
    "]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d173091f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of Landsat scenes found for the entire time periode 1984-2024: 1450\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjsAAAN6CAYAAADYfPukAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACPdUlEQVR4nOzdebRVdf0//udBZpCLIHDlowziBDjPs5YoYo6ZZjlLzopDplI5kYWaA2pOmAP20UpL+WQWiEqSaTjinIKCmopoKggoAvf8/ujH+XoD9F65cNn2eKx11nK/997v/dr7HM5yned9v9+lcrlcDgAAAAAAQEE1aewCAAAAAAAAloSwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMA+K/0l7/8JaVSKb/73e8au5Q6eeedd/Ktb30rHTt2TKlUyrBhwxq7pOXaY489lq233jpt2rRJqVTKhAkTGrukRTr33HNTKpW+8Lgdd9wx66677jKoiCSZMmVKSqVSbr755sYuBQAAqCNhBwCw1Nx8880plUpp2bJl3nzzzYX2+wG37k455ZSMHj06gwcPzq9+9avsuuuuiz125syZOeecc7LuuuumTZs26dixYzbccMOcdNJJeeutt5Zh1Y1j7ty52W+//fL+++/nsssuy69+9at07969sctarv3mN79JqVTKddddt8j9xx57bJo1a5ann356GVe2fFvaoelbb72Vc889d7kN65aGF154Ieeee26mTJnS2KUAAFAwTRu7AADgq2/OnDm54IILcuWVVzZ2KYX1wAMPZK+99sppp532ucfNnTs322+/ff7xj3/k0EMPzYknnpiZM2fm+eefz2233ZZ99tknXbt2XUZVN45XXnklr732Wq6//vp873vfa+xyCuGAAw7IiBEjcuaZZ2bvvfdOly5dKvseffTRDB8+PN///vezwQYbNGKV/33eeuutnHfeeenRo0c23HDDxi5nmXjhhRdy3nnnZccdd0yPHj0auxwAAApE2AEALHUbbrhhrr/++gwePPgr/0P7f5o1a1batGmzxP1MmzYt7du3/8LjRo4cmaeeeiq33nprvvvd79ba98knn+TTTz9d4lqWd9OmTUuSOj0v/p9rrrkmffv2zSmnnJLbbrstSTJ//vwcffTR6datW84999ylXkO5XM4nn3ySVq1aLfVrfRXNnj07rVu3buwyAACgUZjGCgBY6n74wx9m/vz5ueCCCz73uM+bJ79UKtX6sXXBWgcvv/xyDjrooFRVVaVTp04566yzUi6X88Ybb2SvvfZKu3btUl1dnUsuuWSR15w/f35++MMfprq6Om3atMmee+6ZN954Y6Hjxo8fn1133TVVVVVp3bp1dthhh/ztb3+rdcyCml544YV897vfzUorrZRtt932c+/51VdfzX777ZcOHTqkdevW2XLLLXPPPfdU9i+YCqxcLueqq65KqVT63DUeXnnllSTJNttss9C+li1bpl27drXa/vGPf2T//fdPp06d0qpVq6y99tr50Y9+VOuYN998M0cccUS6dOmSFi1apG/fvrnxxhtrHbNgOp/bb789P/3pT7PqqqumZcuW2WmnnTJp0qSFaqnL8/zoo49y8sknp0ePHmnRokU6d+6cnXfeOU8++eRi7/+www7LDjvskCTZb7/9UiqVsuOOO1b2P/DAA9luu+3Spk2btG/fPnvttVdefPHFhfpY1F+UL2p9jVKplBNOOCEjR47MuuuuW3k+o0aNWuj8hx56KJtttllatmyZXr16LXbKqM/zxBNPZOutt06rVq3Ss2fPXHvttZV9M2fOTJs2bXLSSSctdN4///nPrLDCChk6dOhi++7Ro0fOPffc/PrXv86YMWOSJFdccUUmTJiQa665Jq1bt86cOXNyzjnnZI011kiLFi2y2mqr5fTTT8+cOXNq9XXTTTfl61//ejp37pwWLVqkT58+ueaaaxZ5zd133z2jR4/OpptumlatWn3uc/nrX/+a/fbbL926datc/5RTTsnHH39c67jDDjssbdu2zZtvvpm99947bdu2TadOnXLaaadl/vz5tY798MMPc9hhh6Wqqirt27fPoYcemg8//HCxNXyRBZ+TSZMm5bDDDkv79u1TVVWVww8/PLNnz6517JgxY7Ltttumffv2adu2bdZee+388Ic/TPLvf1ObbbZZkuTwww+v/Ntf8P24YBrAJ554Ittvv31at25dOfc/vy8X6NGjRw477LDK9oLvl4ceeiiDBg1Kp06d0r59+xx99NH59NNP8+GHH+aQQw7JSiutlJVWWimnn356yuVyrT5ramoybNiw9O3bNy1btkyXLl1y9NFH54MPPljo2rvvvnseeuihbL755mnZsmVWX3313HLLLbXq2W+//ZIkX/va1yr3/Je//KXe7wMAAP99jOwAAJa6nj175pBDDsn111+fM888s0FHd3z7299O7969c8EFF+See+7J+eefnw4dOuS6667L17/+9Vx44YW59dZbc9ppp2WzzTbL9ttvX+v8n/70pymVSjnjjDMybdq0DBs2LP369cuECRMqf13+wAMPZMCAAdlkk01yzjnnpEmTJpUfc//6179m8803r9XnfvvtlzXXXDM/+9nPFvph8LPeeeedbL311pk9e3YGDRqUjh07ZsSIEdlzzz3zu9/9Lvvss0+23377/OpXv8rBBx+cnXfeOYcccsjnPo8Fa1Pccsst+fGPf/y5wcgzzzyT7bbbLs2aNctRRx2VHj165JVXXsndd9+dn/70p5Uat9xyy8qP+p06dcqf//znDBw4MDNmzMjJJ59cq88LLrggTZo0yWmnnZbp06fnoosuyoEHHpjx48dXjqnr8zzmmGPyu9/9LieccEL69OmTf/3rX3nooYfy4osvZuONN17kPR199NH5n//5n/zsZz/LoEGDstlmm1WmZLrvvvsyYMCArL766jn33HPz8ccf58orr8w222yTJ5988ktPmfPQQw/lzjvvzHHHHZcVV1wxV1xxRfbdd9+8/vrr6dixY5Lk2WefzS677JJOnTrl3HPPzbx583LOOefUmi7qi3zwwQfZbbfdsv/+++c73/lObr/99hx77LFp3rx5jjjiiLRt2zb77LNPfvvb3+bSSy/NCiusUDn317/+dcrlcg488MDPvcYpp5ySW2+9Nccee2xGjRqVs88+OwcccEB23XXX1NTUZM8998xDDz2Uo446Kr17986zzz6byy67LC+//HJGjhxZ6WfBKJE999wzTZs2zd13353jjjsuNTU1Of7442td86WXXsp3vvOdHH300TnyyCOz9tprL7a+O+64I7Nnz86xxx6bjh075tFHH82VV16Zf/7zn7njjjtqHTt//vz0798/W2yxRS6++OLcd999ueSSS9KrV68ce+yxSf49kmSvvfbKQw89lGOOOSa9e/fOXXfdlUMPPbSub8ti7b///unZs2eGDh2aJ598Mr/85S/TuXPnXHjhhUmS559/PrvvvnvWX3/9DBkyJC1atMikSZMqoV/v3r0zZMiQnH322TnqqKOy3XbbJUm23nrryjX+9a9/ZcCAATnggANy0EEH1evz9Fknnnhiqqurc9555+Xvf/97hg8fnvbt2+fhhx9Ot27d8rOf/Sx/+tOf8vOf/zzrrrture+ho48+OjfffHMOP/zwDBo0KJMnT84vfvGLPPXUU/nb3/6WZs2aVY6dNGlSvvWtb2XgwIE59NBDc+ONN+awww7LJptskr59+2b77bfPoEGDcsUVV+SHP/xhevfuXXkWAADwhcoAAEvJTTfdVE5Sfuyxx8qvvPJKuWnTpuVBgwZV9u+www7lvn37VrYnT55cTlK+6aabFuorSfmcc86pbJ9zzjnlJOWjjjqq0jZv3rzyqquuWi6VSuULLrig0v7BBx+UW7VqVT700EMrbWPHji0nKf/P//xPecaMGZX222+/vZykfPnll5fL5XK5pqamvOaaa5b79+9frqmpqRw3e/bscs+ePcs777zzQjV95zvfqdPzOfnkk8tJyn/9618rbR999FG5Z8+e5R49epTnz59f6/6PP/74L+xz9uzZ5bXXXrucpNy9e/fyYYcdVr7hhhvK77zzzkLHbr/99uUVV1yx/Nprr9Vq/+x9Dhw4sLzKKquU33vvvVrHHHDAAeWqqqry7Nmzy+Xy/3uevXv3Ls+ZM6dy3OWXX15OUn722Wcrfdf1eVZVVdXpnv/TglruuOOOWu0bbrhhuXPnzuV//etflbann3663KRJk/IhhxxSaTv00EPL3bt3X6jfBe/vZyUpN2/evDxp0qRafSYpX3nllZW2vffeu9yyZctaz/qFF14or7DCCgv1uSg77LBDOUn5kksuqbTNmTOnck+ffvppuVwul0ePHl1OUv7zn/9c6/z111+/vMMOO3zhdcrlcnn8+PHlJk2alDt06FBu3759eerUqeVyuVz+1a9+VW7SpEmtz2u5XC5fe+215STlv/3tb5W2BZ+Lz+rfv3959dVXr9XWvXv3cpLyqFGj6lTbovodOnRouVQq1Xq2hx56aDlJeciQIbWO3WijjcqbbLJJZXvkyJHlJOWLLrqo0jZv3rzydtttt9jvos9a1GdtwefkiCOOqHXsPvvsU+7YsWNl+7LLLisnKb/77ruL7f+xxx5bbB0LPhPXXnvtQvv+8/tyge7du9f6HlzwHf2f/x632mqrcqlUKh9zzDGVtgXfr5/9HP31r38tJynfeuutta4zatSohdoXvNfjxo2rtE2bNq3cokWL8ve///1K2x133FFOUh47duyiHgkAACyWaawAgGVi9dVXz8EHH5zhw4fn7bffbrB+P7sA9QorrJBNN9005XI5AwcOrLS3b98+a6+9dl599dWFzj/kkEOy4oorVra/9a1vZZVVVsmf/vSnJMmECRMyceLEfPe7382//vWvvPfee3nvvfcya9as7LTTThk3blxqampq9XnMMcfUqfY//elP2XzzzWtNddW2bdscddRRmTJlSl544YW6PYTPaNWqVcaPH58f/OAHSf49LczAgQOzyiqr5MQTT6xMN/Tuu+9m3LhxOeKII9KtW7dafSwYDVIul/P73/8+e+yxR8rlcuXe33vvvfTv3z/Tp09faEqpww8/PM2bN69sL/hr9AXPvj7Ps3379hk/fnzeeuutej+H//T2229nwoQJOeyww9KhQ4dK+/rrr5+dd9658n5/Gf369UuvXr1q9dmuXbvKPc+fPz+jR4/O3nvvXetZ9+7dO/3796/zdZo2bZqjjz66st28efMcffTRmTZtWp544olKLV27ds2tt95aOe65557LM888k4MOOqhO19l8881zzDHH5P3338/QoUMrowXuuOOO9O7dO+uss06tz8LXv/71JMnYsWMrfXx2zY3p06fnvffeyw477JBXX30106dPr3W9nj171vk5fLbfWbNm5b333svWW2+dcrmcp556aqHj//Pf4nbbbVfre+BPf/pTmjZtWhnpkfz7e+TEE0+sUz2fZ1HX/te//pUZM2Yk+X9ryvzf//3fQt8hddWiRYscfvjhS1RnkgwcOLDWKLAttthioe/RBd+vn31+d9xxR6qqqrLzzjvX+kxssskmadu2ba3PRJL06dOn8p2QJJ06dVrsdzMAANSXsAMAWGZ+/OMfZ968eV+4dkd9/OcP9VVVVWnZsmVWXnnlhdr/cw75JFlzzTVrbZdKpayxxhqZMmVKkmTixIlJkkMPPTSdOnWq9frlL3+ZOXPmLPLH27p47bXXFjllz4IpW1577bU69fOfqqqqctFFF2XKlCmZMmVKbrjhhqy99tr5xS9+kZ/85CdJ/l/4sO666y62n3fffTcffvhhhg8fvtC9L/iBdcFi4Av85/ux0korJUnl2dfneV500UV57rnnstpqq2XzzTfPueee+6V/FF3wLBf3vBcELl/Gf95z8u/7XnDP7777bj7++OOFPmuLq2dxunbtutBi92uttVaSVD6vTZo0yYEHHpiRI0dW1oe49dZb07Jly8paCHWxYK2ITTfdtNI2ceLEPP/88wu9bwtq+Oxn4W9/+1v69etXWRulU6dOlfUkvuy/lyR5/fXXK4HVgnU4FqzR8p/9tmzZMp06darV9tn3Jfn352KVVVZJ27Ztax1Xn/dlcb7o38K3v/3tbLPNNvne976XLl265IADDsjtt99er+Djf/7nf2qFiw1Va1VVVZJktdVWW6j9s89v4sSJmT59ejp37rzQ52LmzJlf+P2QLPyeAADAl2XNDgBgmVl99dVz0EEHZfjw4TnzzDMX2r+49SX+c0Hhz/rsugSf15bkc9fPWJwFPzz+/Oc/z4YbbrjIY/7zh9LP/vV5Y+vevXuOOOKI7LPPPll99dVz66235vzzz6/TuQvu/aCDDlrsGgbrr79+re0vevb1eZ77779/tttuu9x1112599578/Of/zwXXnhh7rzzzgwYMKBO9/Bl1Pdz2JCft4ZwyCGH5Oc//3lGjhyZ73znO7ntttuy++67V37A/rJqamqy3nrr5dJLL13k/gU/jL/yyivZaaedss466+TSSy/NaqutlubNm+dPf/pTLrvssoV+zK/rv5f58+dn5513zvvvv58zzjgj66yzTtq0aZM333wzhx122EL9Lu59WVa+6HPRqlWrjBs3LmPHjs0999yTUaNG5be//W2+/vWv5957761T/fX9rqnvZ3hR7Z/9XNfU1KRz5861RhJ91n+GTcvbvxUAAL5ahB0AwDL14x//OP/7v/9bWaT3sxb85fOHH35Yq/3LjnCoiwUjDRYol8uZNGlS5Uf8BdMTtWvXLv369WvQa3fv3j0vvfTSQu3/+Mc/KvsbykorrZRevXrlueeeS/Lv4ClJZXtROnXqlBVXXDHz589vsHuv7/NcZZVVctxxx+W4447LtGnTsvHGG+enP/1pvcOOBc9ycc975ZVXroyaWGmllRb6DCZf/nPYqVOntGrVaqHP2uLqWZy33nors2bNqjW64+WXX06SWourr7vuutloo41y6623ZtVVV83rr7+eK6+88kvV/lm9evXK008/nZ122ulzF76/++67M2fOnPzhD3+o9Zf8/zmlUX09++yzefnllzNixIhaC2SPGTPmS/fZvXv33H///Zk5c2at0LI+78uSaNKkSXbaaafstNNOufTSS/Ozn/0sP/rRjzJ27Nj069fvc5/z51nUZ/jTTz9t0CkEk39/Ju67775ss802DRbyftl7BgAA01gBAMtUr169ctBBB+W6667L1KlTa+1r165dVl555YwbN65W+9VXX73U6rnlllvy0UcfVbZ/97vf5e233678mL7JJpukV69eufjiizNz5syFzn/33Xe/9LV32223PProo3nkkUcqbbNmzcrw4cPTo0eP9OnTp959Pv3003nvvfcWan/ttdfywgsvVKbn6dSpU7bffvvceOONef3112sdu+CvrFdYYYXsu++++f3vf7/IUOTL3Htdn+f8+fMXmpaoc+fO6dq1a2XdkfpYZZVVsuGGG2bEiBG1fgR+7rnncu+992a33XartPXq1SvTp0/PM888U2l7++23c9ddd9X7usm/n2P//v0zcuTIWs/6xRdfzOjRo+vcz7x583LddddVtj/99NNcd9116dSpUzbZZJNaxx588MG59957M2zYsHTs2LFBRsLsv//+efPNN3P99dcvtO/jjz+uTAO24K/3P/vX+tOnT89NN920RNdfVL/lcjmXX375l+5zt912y7x583LNNddU2ubPn98g4dAXef/99xdqWzDaacFnfEGwtajw7fP06tVroe/R4cOHf+4ouS9j//33z/z58yvT433WvHnz6l138uXvGQAAjOwAAJa5H/3oR/nVr36Vl156KX379q2173vf+14uuOCCfO9738umm26acePGVf56fWno0KFDtt122xx++OF55513MmzYsKyxxho58sgjk/z7L69/+ctfZsCAAenbt28OP/zw/M///E/efPPNjB07Nu3atcvdd9/9pa595pln5te//nUGDBiQQYMGpUOHDhkxYkQmT56c3//+92nSpP5/lzJmzJicc8452XPPPbPlllumbdu2efXVV3PjjTdmzpw5OffccyvHXnHFFdl2222z8cYb56ijjkrPnj0zZcqU3HPPPZkwYUKS5IILLsjYsWOzxRZb5Mgjj0yfPn3y/vvv58knn8x99923yB9sP09dn+dHH32UVVddNd/61reywQYbpG3btrnvvvvy2GOP5ZJLLqn3c0n+PXXWgAEDstVWW2XgwIH5+OOPc+WVV6aqqqrWcznggANyxhlnZJ999smgQYMye/bsXHPNNVlrrbUWWpC9rs4777yMGjUq2223XY477rjMmzcvV155Zfr27VsrVPk8Xbt2zYUXXpgpU6ZkrbXWym9/+9tMmDAhw4cPT7NmzWod+93vfjenn3567rrrrhx77LEL7f8yDj744Nx+++055phjMnbs2GyzzTaZP39+/vGPf+T222/P6NGjs+mmm2aXXXZJ8+bNs8cee+Too4/OzJkzc/3116dz585LNLJgnXXWSa9evXLaaaflzTffTLt27fL73/9+idZ72GOPPbLNNtvkzDPPzJQpU9KnT5/ceeedCwVtS8OQIUMybty4fOMb30j37t0zbdq0XH311Vl11VWz7bbbJvl3aNG+fftce+21WXHFFdOmTZtsscUWX7jOyfe+970cc8wx2XfffbPzzjvn6aefzujRoxday2hJ7bDDDjn66KMzdOjQTJgwIbvsskuaNWuWiRMn5o477sjll1+eb33rW/Xqc8MNN8wKK6yQCy+8MNOnT0+LFi3y9a9/PZ07d27Q2gEA+OoRdgAAy9waa6yRgw46KCNGjFho39lnn5133303v/vd73L77bdnwIAB+fOf/7zUfuj64Q9/mGeeeSZDhw7NRx99lJ122ilXX311WrduXTlmxx13zCOPPJKf/OQn+cUvfpGZM2emuro6W2yxRY4++ugvfe0uXbrk4YcfzhlnnJErr7wyn3zySdZff/3cfffd+cY3vvGl+tx3333z0Ucf5d57780DDzyQ999/PyuttFI233zzfP/738/Xvva1yrEbbLBB/v73v+ess87KNddck08++STdu3fP/vvvX6vGRx99NEOGDMmdd96Zq6++Oh07dkzfvn0XORVZXdTlebZu3TrHHXdc7r333tx5552pqanJGmuskauvvjrHHnvsl7puv379MmrUqJxzzjk5++yz06xZs+ywww658MILa/143LFjx9x111059dRTc/rpp6dnz54ZOnRoJk6c+KXDjvXXXz+jR4/OqaeemrPPPjurrrpqzjvvvLz99tt1DjtWWmmljBgxIieeeGKuv/76dOnSJb/4xS8qwdxndenSJbvsskv+9Kc/5eCDD/5SNf+nJk2aZOTIkbnssstyyy235K677krr1q2z+uqr56STTqosVL722mvnd7/7XX784x/ntNNOS3V1dY499th06tQpRxxxxJe+frNmzXL33Xdn0KBBGTp0aFq2bJl99tknJ5xwQjbYYIMvfU9/+MMfcvLJJ+d///d/UyqVsueee+aSSy7JRhtt9KVrrYs999wzU6ZMyY033pj33nsvK6+8cnbYYYecd955lfVVmjVrlhEjRmTw4ME55phjMm/evNx0001fGHYceeSRmTx5cm644YZKyDZmzJjstNNODX4f1157bTbZZJNcd911+eEPf5imTZumR48eOeigg7LNNtvUu7/q6upce+21GTp0aAYOHJj58+dn7Nixwg4AAL5QqWw1OAAA+MrZZ5998uyzz2bSpEmNXQoAAMBSZ80OAAD4inn77bdzzz33NNioDgAAgOWdaawAAOArYvLkyfnb3/6WX/7yl2nWrNkSTbMGAABQJEZ2AADAV8SDDz6Ygw8+OJMnT86IESNSXV3d2CUBAAAsE9bsAAAAAAAACs3IDgAAAAAAoNCEHQAAAAAAQKFZoDxJTU1N3nrrray44ooplUqNXQ4AAAAAwFJTLpfz0UcfpWvXrmnSxN/D89Ug7Ejy1ltvZbXVVmvsMgAAAAAAlpk33ngjq666amOXAQ1C2JFkxRVXTPLvf9zt2rVr5GoAAAAAAJaeGTNmZLXVVqv8LgpfBcKOpDJ1Vbt27YQdAAAAAMB/BVP681ViQjYAAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCs2YHAAAAAAB8jvnz52fu3LmNXcZ/nebNm6dJk7qN2RB2AAAAAADAIpTL5UydOjUffvhhY5fyX6lJkybp2bNnmjdv/oXHCjsAAAAAAGARFgQdnTt3TuvWrVMqlRq7pP8aNTU1eeutt/L222+nW7duX/jshR0AAAAAAPAf5s+fXwk6Onbs2Njl/Ffq1KlT3nrrrcybNy/NmjX73GMtUA4AAAAAAP9hwRodrVu3buRK/nstmL5q/vz5X3issAMAAAAAABbD1FWNpz7PXtgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAAA0gMMOOyx77713g/TVo0ePDBs2rEH6akxTpkxJqVTKhAkTlup1hB0AAAAAAFBA8+fPT01NTWOXsVwQdgAAAAAAQAPbcccdM2jQoJx++unp0KFDqqurc+6551b2l8vlnHvuuenWrVtatGiRrl27ZtCgQZVzX3vttZxyyikplUqVhbpvvvnmtG/fPn/4wx/Sp0+ftGjRIq+//np23HHHnHzyybWuv/fee+ewww6rbPfo0SPnn39+DjnkkLRt2zbdu3fPH/7wh7z77rvZa6+90rZt26y//vp5/PHHa/Xz0EMPZbvttkurVq2y2mqrZdCgQZk1a1atfn/2s5/liCOOyIorrphu3bpl+PDhlf09e/ZMkmy00UYplUrZcccdG+DpLkzYAQAAAAAAS8GIESPSpk2bjB8/PhdddFGGDBmSMWPGJEl+//vf57LLLst1112XiRMnZuTIkVlvvfWSJHfeeWdWXXXVDBkyJG+//XbefvvtSp+zZ8/OhRdemF/+8pd5/vnn07lz5zrXc9lll2WbbbbJU089lW984xs5+OCDc8ghh+Sggw7Kk08+mV69euWQQw5JuVxOkrzyyivZdddds+++++aZZ57Jb3/72zz00EM54YQTavV7ySWXZNNNN81TTz2V4447Lscee2xeeumlJMmjjz6aJLnvvvvy9ttv58477/zyD/RzCDsAAAAAAGApWH/99XPOOedkzTXXzCGHHJJNN900999/f5Lk9ddfT3V1dfr165du3bpl8803z5FHHpkk6dChQ1ZYYYWsuOKKqa6uTnV1daXPuXPn5uqrr87WW2+dtddeO61bt65zPbvttluOPvrorLnmmjn77LMzY8aMbLbZZtlvv/2y1lpr5YwzzsiLL76Yd955J0kydOjQHHjggTn55JOz5pprZuutt84VV1yRW265JZ988kmtfo877risscYaOeOMM7Lyyitn7NixSZJOnTolSTp27Jjq6up06NBhyR7qYgg7AAAAAABgKVh//fVrba+yyiqZNm1akmS//fbLxx9/nNVXXz1HHnlk7rrrrsybN+8L+2zevPlC/X6Zerp06ZIkldEkn21bUOPTTz+dm2++OW3btq28+vfvn5qamkyePHmR/ZZKpVRXV1f6WFaEHQAAAAAAsBQ0a9as1napVKosKL7aaqvlpZdeytVXX51WrVrluOOOy/bbb5+5c+d+bp+tWrWqrOGxQJMmTSpTTy2wqH4+W8+CPhbVtqDGmTNn5uijj86ECRMqr6effjoTJ05Mr1696nSfy0rTZXo1AAAAAAAgyb+Diz322CN77LFHjj/++Kyzzjp59tlns/HGG6d58+aZP39+nfrp1KlTrXU95s+fn+eeey5f+9rXlqi+jTfeOC+88ELWWGONL91H8+bNKzUtTUZ2AAAAAADAMnbzzTfnhhtuyHPPPZdXX301//u//5tWrVqle/fuSZIePXpk3LhxefPNN/Pee+99bl9f//rXc8899+See+7JP/7xjxx77LH58MMPl7jGM844Iw8//HBOOOGETJgwIRMnTsz//d//LbRA+efp3LlzWrVqlVGjRuWdd97J9OnTl7iuRRF2AAAAAADAMta+fftcf/312WabbbL++uvnvvvuy913352OHTsmSYYMGZIpU6akV69elUW+F+eII47IoYcemkMOOSQ77LBDVl999SUe1ZH8ey2OBx98MC+//HK22267bLTRRjn77LPTtWvXOvfRtGnTXHHFFbnuuuvStWvX7LXXXktc16KUyv85kdd/oRkzZqSqqirTp09Pu3btGrscAAAAAIClxu+hdfPJJ59k8uTJ6dmzZ1q2bNnY5fxXqs97YGQHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAA+C/yl7/8JaVSKR9++GGS5Oabb0779u0bvN9lqekyvyIAAAAAABTYzk32W2bXGlNzx1K/xre//e3stttu9Tpnxx13zIYbbphhw4ZV2rbeeuu8/fbbqaqqauAKv5iwAwAAAAAA/ou1atUqrVq1WuJ+mjdvnurq6sXunz9/fkqlUpo0afhJp0xjBQAAAAAAXyFz5szJoEGD0rlz57Rs2TLbbrttHnvsscUeX99prA477LA8+OCDufzyy1MqlVIqlTJlypTFTo/1hz/8IX369EmLFi3y+uuv5y9/+Us233zztGnTJu3bt88222yT1157bYnu2cgOAAAAAID/X0NMT7Qsph2Cz3P66afn97//fUaMGJHu3bvnoosuSv/+/TNp0qR06NBhifu//PLL8/LLL2fdddfNkCFDkiSdOnXKlClTFjp29uzZufDCC/PLX/4yHTt2TIcOHbLhhhvmyCOPzK9//et8+umnefTRR1MqlZaoJmEHAAAAAAB8RcyaNSvXXHNNbr755gwYMCBJcv3112fMmDG54YYb8oMf/GCJr1FVVZXmzZundevWnzttVZLMnTs3V199dTbYYIMkyfvvv5/p06dn9913T69evZIkvXv3XuKaTGMFAAAAAABfEa+88krmzp2bbbbZptLWrFmzbL755nnxxReXeT3NmzfP+uuvX9nu0KFDDjvssPTv3z977LFHLr/88rz99ttLfB1hBwAAAAAAsFS0atVqoSmqbrrppjzyyCPZeuut89vf/jZrrbVW/v73vy/RdYQdAAAAAADwFdGrV680b948f/vb3yptc+fOzWOPPZY+ffo02HWaN2+e+fPnf+nzN9poowwePDgPP/xw1l133dx2221LVI81OwAAAAAA4CuiTZs2OfbYY/ODH/wgHTp0SLdu3XLRRRdl9uzZGThwYINdp0ePHhk/fnymTJmStm3b1nnh88mTJ2f48OHZc88907Vr17z00kuZOHFiDjnkkCWqR9gBAAAAAABfIRdccEFqampy8MEH56OPPsqmm26a0aNHZ6WVVqrT+VOmTEnPnj0zduzY7Ljjjos85rTTTsuhhx6aPn365OOPP87kyZPr1Hfr1q3zj3/8IyNGjMi//vWvrLLKKjn++ONz9NFH1/X2FqlULpfLS9TDV8CMGTNSVVWV6dOnp127do1dDgAAAADQSHZust8S9zGm5o4GqGTp8Xto3XzyySeZPHlyevbsmZYtWzZ2OcvU2LFj881vfjOvvvpqnQOSpaE+74E1OwAAAAAAgIo//elP+eEPf9ioQUd9mcYKAAAAAACo+PnPf97YJdSbkR0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0Jo2dgEAAAAAAFAkNVPXWmbXalL9cr3P2XHHHbPhhhtm2LBhDV/QcsrIDgAAAAAAoNCEHQAAAAAA8BVx2GGH5cEHH8zll1+eUqmUUqmUm2++OaVSKaNHj85GG22UVq1a5etf/3qmTZuWP//5z+ndu3fatWuX7373u5k9e3Zj38KXYhorAAAAAAD4irj88svz8ssvZ911182QIUOSJM8//3yS5Nxzz80vfvGLtG7dOvvvv3/233//tGjRIrfddltmzpyZffbZJ1deeWXOOOOMxryFL0XYAQAAAAAAXxFVVVVp3rx5Wrdunerq6iTJP/7xjyTJ+eefn2222SZJMnDgwAwePDivvPJKVl999STJt771rYwdO7aQYYdprAAAAAAA4L/A+uuvX/nvLl26pHXr1pWgY0HbtGnTGqO0JSbsAAAAAACA/wLNmjWr/HepVKq1vaCtpqZmWZfVIIQdAAAAAADwFdK8efPMnz+/sctYpho17Bg3blz22GOPdO3aNaVSKSNHjlzomBdffDF77rlnqqqq0qZNm2y22WZ5/fXXK/s/+eSTHH/88enYsWPatm2bfffdN++8884yvAsAAAAAAFh+9OjRI+PHj8+UKVPy3nvvFXa0Rn00atgxa9asbLDBBrnqqqsWuf+VV17Jtttum3XWWSd/+ctf8swzz+Sss85Ky5YtK8eccsopufvuu3PHHXfkwQcfzFtvvZVvfvOby+oWAAAAAABguXLaaadlhRVWSJ8+fdKpU6daAwi+qkrlcrnc2EUk/54L7K677sree+9daTvggAPSrFmz/OpXv1rkOdOnT0+nTp1y22235Vvf+laSf68q37t37zzyyCPZcsst63TtGTNmpKqqKtOnT0+7du2W+F4AAAAAgGLaucl+S9zHmJo7GqCSpcfvoXXzySefZPLkyenZs2etP8Bn2anPe7DcrtlRU1OTe+65J2uttVb69++fzp07Z4sttqg11dUTTzyRuXPnpl+/fpW2ddZZJ926dcsjjzzSCFUDAAAAAADL2nIbdkybNi0zZ87MBRdckF133TX33ntv9tlnn3zzm9/Mgw8+mCSZOnVqmjdvnvbt29c6t0uXLpk6depi+54zZ05mzJhR6wUAAAAAABRT08YuYHEWLJiy11575ZRTTkmSbLjhhnn44Ydz7bXXZocddvjSfQ8dOjTnnXdeg9QJAAAAAAA0ruV2ZMfKK6+cpk2bpk+fPrXae/fuXVlMpbq6Op9++mk+/PDDWse88847qa6uXmzfgwcPzvTp0yuvN954o8HrBwAAAAAAlo3lNuxo3rx5Nttss7z00ku12l9++eV07949SbLJJpukWbNmuf/++yv7X3rppbz++uvZaqutFtt3ixYt0q5du1ovAAAAAAD4T+VyubFL+K9Vn2ffqNNYzZw5M5MmTapsT548ORMmTEiHDh3SrVu3/OAHP8i3v/3tbL/99vna176WUaNG5e67785f/vKXJElVVVUGDhyYU089NR06dEi7du1y4oknZquttsqWW27ZSHcFAAAAAEDRNWvWLEkye/bstGrVqpGr+e/06aefJklWWGGFLzy2UcOOxx9/PF/72tcq26eeemqS5NBDD83NN9+cffbZJ9dee22GDh2aQYMGZe21187vf//7bLvttpVzLrvssjRp0iT77rtv5syZk/79++fqq69e5vcCAAAAAMBXxworrJD27dtn2rRpSZLWrVunVCo1clX/PWpqavLuu++mdevWadr0i6OMUtkYnMyYMSNVVVWZPn26Ka0AAAAA4L/Yzk32W+I+xtTc0QCVLD1+D627crmcqVOnLrRuNMtGkyZN0rNnzzRv3vwLj23UkR0AAAAAALC8KpVKWWWVVdK5c+fMnTu3scv5r9O8efM0aVK3pceFHQAAAAAA8DlWWGGFOq0bQeOpWyQCAAAAAACwnBJ2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAAqtaWMXAAAAAACQJDs32W+J+xhTc0cDVAIUjZEdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChNWrYMW7cuOyxxx7p2rVrSqVSRo4cudhjjznmmJRKpQwbNqxW+/vvv58DDzww7dq1S/v27TNw4MDMnDlz6RYOAAAAAAAsNxo17Jg1a1Y22GCDXHXVVZ973F133ZW///3v6dq160L7DjzwwDz//PMZM2ZM/vjHP2bcuHE56qijllbJAAAAAADAcqZpY158wIABGTBgwOce8+abb+bEE0/M6NGj841vfKPWvhdffDGjRo3KY489lk033TRJcuWVV2a33XbLxRdfvMhwBAAAAAAA+GpZrtfsqKmpycEHH5wf/OAH6du370L7H3nkkbRv374SdCRJv3790qRJk4wfP36x/c6ZMyczZsyo9QIAAAAAAIppuQ47LrzwwjRt2jSDBg1a5P6pU6emc+fOtdqaNm2aDh06ZOrUqYvtd+jQoamqqqq8VltttQatGwAAAAAAWHaW27DjiSeeyOWXX56bb745pVKpQfsePHhwpk+fXnm98cYbDdo/AAAAAACw7Cy3Ycdf//rXTJs2Ld26dUvTpk3TtGnTvPbaa/n+97+fHj16JEmqq6szbdq0WufNmzcv77//fqqrqxfbd4sWLdKuXbtaLwAAAAAAoJgadYHyz3PwwQenX79+tdr69++fgw8+OIcffniSZKuttsqHH36YJ554IptsskmS5IEHHkhNTU222GKLZV4zAAAAAACw7DVq2DFz5sxMmjSpsj158uRMmDAhHTp0SLdu3dKxY8daxzdr1izV1dVZe+21kyS9e/fOrrvumiOPPDLXXntt5s6dmxNOOCEHHHBAunbtukzvBQAAAAAAaByNOo3V448/no022igbbbRRkuTUU0/NRhttlLPPPrvOfdx6661ZZ511stNOO2W33XbLtttum+HDhy+tkgEAAAAAgOVMo47s2HHHHVMul+t8/JQpUxZq69ChQ2677bYGrAoAAAAAACiS5XaBcgAAAAAAgLoQdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBojRp2jBs3LnvssUe6du2aUqmUkSNHVvbNnTs3Z5xxRtZbb720adMmXbt2zSGHHJK33nqrVh/vv/9+DjzwwLRr1y7t27fPwIEDM3PmzGV8JwAAAAAAQGNp1LBj1qxZ2WCDDXLVVVcttG/27Nl58sknc9ZZZ+XJJ5/MnXfemZdeeil77rlnreMOPPDAPP/88xkzZkz++Mc/Zty4cTnqqKOW1S0AAAAAAACNrGljXnzAgAEZMGDAIvdVVVVlzJgxtdp+8YtfZPPNN8/rr7+ebt265cUXX8yoUaPy2GOPZdNNN02SXHnlldltt91y8cUXp2vXrkv9HgAAAAAAgMZVqDU7pk+fnlKplPbt2ydJHnnkkbRv374SdCRJv3790qRJk4wfP36x/cyZMyczZsyo9QIAAAAAAIqpMGHHJ598kjPOOCPf+c530q5duyTJ1KlT07lz51rHNW3aNB06dMjUqVMX29fQoUNTVVVVea222mpLtXYAAAAAAGDpKUTYMXfu3Oy///4pl8u55pprlri/wYMHZ/r06ZXXG2+80QBVAgAAAAAAjaFR1+yoiwVBx2uvvZYHHnigMqojSaqrqzNt2rRax8+bNy/vv/9+qqurF9tnixYt0qJFi6VWMwAAAAAAsOws1yM7FgQdEydOzH333ZeOHTvW2r/VVlvlww8/zBNPPFFpe+CBB1JTU5MttthiWZcLAAAAAAA0gkYd2TFz5sxMmjSpsj158uRMmDAhHTp0yCqrrJJvfetbefLJJ/PHP/4x8+fPr6zD0aFDhzRv3jy9e/fOrrvumiOPPDLXXntt5s6dmxNOOCEHHHBAunbt2li3BQAAAAAALEONGnY8/vjj+drXvlbZPvXUU5Mkhx56aM4999z84Q9/SJJsuOGGtc4bO3ZsdtxxxyTJrbfemhNOOCE77bRTmjRpkn333TdXXHHFMqkfAAAAAABofI0aduy4444pl8uL3f95+xbo0KFDbrvttoYsCwAAAAAAKJDles0OAAAAAACALyLsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFFq9w4433ngj//znPyvbjz76aE4++eQMHz68QQsDAAAAAACoi3qHHd/97nczduzYJMnUqVOz884759FHH82PfvSjDBkypMELBAAAAAAA+Dz1Djuee+65bL755kmS22+/Peuuu24efvjh3Hrrrbn55psbuj4AAAAAAIDPVe+wY+7cuWnRokWS5L777suee+6ZJFlnnXXy9ttv16uvcePGZY899kjXrl1TKpUycuTIWvvL5XLOPvvsrLLKKmnVqlX69euXiRMn1jrm/fffz4EHHph27dqlffv2GThwYGbOnFnf2wIAAAAAAAqq3mFH3759c+211+avf/1rxowZk1133TVJ8tZbb6Vjx4716mvWrFnZYIMNctVVVy1y/0UXXZQrrrgi1157bcaPH582bdqkf//++eSTTyrHHHjggXn++eczZsyY/PGPf8y4ceNy1FFH1fe2AAAAAACAgmpa3xMuvPDC7LPPPvn5z3+eQw89NBtssEGS5A9/+ENlequ6GjBgQAYMGLDIfeVyOcOGDcuPf/zj7LXXXkmSW265JV26dMnIkSNzwAEH5MUXX8yoUaPy2GOPZdNNN02SXHnlldltt91y8cUXp2vXrvW9PQAAAAAAoGDqHXbsuOOOee+99zJjxoystNJKlfajjjoqrVu3brDCJk+enKlTp6Zfv36VtqqqqmyxxRZ55JFHcsABB+SRRx5J+/btK0FHkvTr1y9NmjTJ+PHjs88++zRYPQAAAAAAwPKp3mFH8u9RF0888UReeeWVfPe7382KK66Y5s2bN2jYMXXq1CRJly5darV36dKlsm/q1Knp3Llzrf1NmzZNhw4dKscsypw5czJnzpzK9owZMxqqbAAAAAAAYBmrd9jx2muvZdddd83rr7+eOXPmZOedd86KK66YCy+8MHPmzMm11167NOpsUEOHDs15553X2GUAAAAAAAANoN4LlJ900knZdNNN88EHH6RVq1aV9n322Sf3339/gxVWXV2dJHnnnXdqtb/zzjuVfdXV1Zk2bVqt/fPmzcv7779fOWZRBg8enOnTp1deb7zxRoPVDQAAAAAALFv1Djv++te/5sc//nGaN29eq71Hjx558803G6ywnj17prq6ulaAMmPGjIwfPz5bbbVVkmSrrbbKhx9+mCeeeKJyzAMPPJCamppsscUWi+27RYsWadeuXa0XAAAAAABQTPWexqqmpibz589fqP2f//xnVlxxxXr1NXPmzEyaNKmyPXny5EyYMCEdOnRIt27dcvLJJ+f888/PmmuumZ49e+ass85K165ds/feeydJevfunV133TVHHnlkrr322sydOzcnnHBCDjjggHTt2rW+twYAAAAAABRQvUd27LLLLhk2bFhlu1QqZebMmTnnnHOy22671auvxx9/PBtttFE22mijJMmpp56ajTbaKGeffXaS5PTTT8+JJ56Yo446KptttllmzpyZUaNGpWXLlpU+br311qyzzjrZaaedsttuu2XbbbfN8OHD63tbAAAAAABAQZXK5XK5Pif885//TP/+/VMulzNx4sRsuummmThxYlZeeeWMGzcunTt3Xlq1LjUzZsxIVVVVpk+fbkorAAAAAGgkOzfZb4n7GFNzR+FrWNr8HspXUb2nsVp11VXz9NNP5ze/+U2eeeaZzJw5MwMHDsyBBx5Ya8FyAAAAAACAZaHeYUeSNG3aNAcddFBD1wIAAAAAAFBvXyrsmDhxYsaOHZtp06alpqam1r4F620AAAAAAAAsC/UOO66//voce+yxWXnllVNdXZ1SqVTZVyqVhB0AAAAAAMAyVe+w4/zzz89Pf/rTnHHGGUujHgAAAAAAgHppUt8TPvjgg+y3335LoxYAAAAAAIB6q3fYsd9+++Xee+9dGrUAAAAAAADUW72nsVpjjTVy1lln5e9//3vWW2+9NGvWrNb+QYMGNVhxAAAAAAAAX6TeYcfw4cPTtm3bPPjgg3nwwQdr7SuVSsIOAAAAAABgmap32DF58uSlUQcAAAAAAMCXUu81Oxb49NNP89JLL2XevHkNWQ8AAAAAAEC91DvsmD17dgYOHJjWrVunb9++ef3115MkJ554Yi644IIGLxAAAAAAAODz1DvsGDx4cJ5++un85S9/ScuWLSvt/fr1y29/+9sGLQ4AAAAAAOCL1HvNjpEjR+a3v/1tttxyy5RKpUp7375988orrzRocQAAAAAAAF+k3iM73n333XTu3Hmh9lmzZtUKPwAAAAAAAJaFeocdm266ae65557K9oKA45e//GW22mqrhqsMAAAAAACgDuo9jdXPfvazDBgwIC+88ELmzZuXyy+/PC+88EIefvjhPPjgg0ujRgAAAAAAgMWq98iObbfdNhMmTMi8efOy3nrr5d57703nzp3zyCOPZJNNNlkaNQIAAAAAACxWvUd2JEmvXr1y/fXXN3QtAAAAAAAA9VbvkR1/+tOfMnr06IXaR48enT//+c8NUhQAAAAAAEBd1TvsOPPMMzN//vyF2svlcs4888wGKQoAAAAAAKCu6h12TJw4MX369FmofZ111smkSZMapCgAAAAAAIC6qnfYUVVVlVdffXWh9kmTJqVNmzYNUhQAAAAAAEBd1Tvs2GuvvXLyySfnlVdeqbRNmjQp3//+97Pnnns2aHEAAAAAAABfpN5hx0UXXZQ2bdpknXXWSc+ePdOzZ8/07t07HTt2zMUXX7w0agQAAAAAAFispvU9oaqqKg8//HDGjBmTp59+Oq1atcr666+f7bfffmnUBwAAAAAsIzs32W+Jzh9Tc0cDVQJQP/UOO5KkVCpll112yS677NLQ9QAAAAAAANRLnaexeuSRR/LHP/6xVtstt9ySnj17pnPnzjnqqKMyZ86cBi8QAAAAAADg89Q57BgyZEief/75yvazzz6bgQMHpl+/fjnzzDNz9913Z+jQoUulSAAAAAAAgMWpc9gxYcKE7LTTTpXt3/zmN9liiy1y/fXX59RTT80VV1yR22+/fakUCQAAAAAAsDh1Djs++OCDdOnSpbL94IMPZsCAAZXtzTbbLG+88UbDVgcAAAAAAPAF6hx2dOnSJZMnT06SfPrpp3nyySez5ZZbVvZ/9NFHadasWcNXCAAAAAAA8DnqHHbstttuOfPMM/PXv/41gwcPTuvWrbPddttV9j/zzDPp1avXUikSAAAAAABgcZrW9cCf/OQn+eY3v5kddtghbdu2zYgRI9K8efPK/htvvDG77LLLUikSAAAAAABgceocdqy88soZN25cpk+fnrZt22aFFVaotf+OO+5I27ZtG7xAAAAAgK+qnZvst8R9jKm5o1FrWNLrA0BDqHPYsUBVVdUi2zt06LDExQAAAAAAANRXndfsAAAAAAAAWB4JOwAAAAAAgEKr9zRWAAAAAEtqeVirguXH8rBuiBoAiq1OIzs23njjfPDBB0mSIUOGZPbs2Uu1KAAAAAAAgLqqU9jx4osvZtasWUmS8847LzNnzlyqRQEAAAAAANRVnaax2nDDDXP44Ydn2223TblczsUXX5y2bdsu8tizzz67QQsEAAAAAAD4PHUKO26++eacc845+eMf/5hSqZQ///nPadp04VNLpZKwAwAAAAAAWKbqFHasvfba+c1vfpMkadKkSe6///507tx5qRYGAAAAAABQF3UKOz6rpqZmadQBAAAAAADwpdQ77EiSV155JcOGDcuLL76YJOnTp09OOumk9OrVq0GLAwAAAAAA+CJN6nvC6NGj06dPnzz66KNZf/31s/7662f8+PHp27dvxowZszRqBAAAAAAAWKx6j+w488wzc8opp+SCCy5YqP2MM87Izjvv3GDFAQAAAAAAfJF6j+x48cUXM3DgwIXajzjiiLzwwgsNUhQAAAAAAEBd1Tvs6NSpUyZMmLBQ+4QJE9K5c+eGqAkAAAAAAKDO6j2N1ZFHHpmjjjoqr776arbeeuskyd/+9rdceOGFOfXUUxu8QAAAAAAAgM9T77DjrLPOyoorrphLLrkkgwcPTpJ07do15557bgYNGtTgBQIAAAAAAHyeeocdpVIpp5xySk455ZR89NFHSZIVV1yxwQsDAAAAAACoi3qHHZ8l5AAAAAAAABpbvRcoBwAAAAAAWJ4IOwAAAAAAgEJbommsAAAAAIpq5yb7LXEfY2ruaIBKAIAlVa+RHXPnzs1OO+2UiRMnLq16AAAAAAAA6qVeYUezZs3yzDPPLK1aAAAAAAAA6q3ea3YcdNBBueGGG5ZGLQAAAAAAAPVW7zU75s2blxtvvDH33XdfNtlkk7Rp06bW/ksvvbTBigMAAAAAAPgi9Q47nnvuuWy88cZJkpdffrnWvlKp1DBVAQAAAAAA1FG9w46xY8cujToAAAAAAAC+lHqv2bHApEmTMnr06Hz88cdJknK53GBFAQAAAAAA1FW9w45//etf2WmnnbLWWmtlt912y9tvv50kGThwYL7//e83eIEAAAAAAACfp95hxymnnJJmzZrl9ddfT+vWrSvt3/72tzNq1KgGLQ4AAAAAAOCL1HvNjnvvvTejR4/OqquuWqt9zTXXzGuvvdZghQEAAAAAANRFvUd2zJo1q9aIjgXef//9tGjRokGKAgAAAAAAqKt6hx3bbbddbrnllsp2qVRKTU1NLrroonzta19r0OIAAAAAAAC+SL2nsbrooouy00475fHHH8+nn36a008/Pc8//3zef//9/O1vf1saNQIAAAAAACxWvUd2rLvuunn55Zez7bbbZq+99sqsWbPyzW9+M0899VR69eq1NGoEAAAAAABYrHqP7EiSqqqq/OhHP2roWgAAAAAAAOrtS4UdH3zwQW644Ya8+OKLSZI+ffrk8MMPT4cOHRq0OAAAAAAAgC9S72msxo0blx49euSKK67IBx98kA8++CBXXHFFevbsmXHjxi2NGgEAAAAAABar3iM7jj/++Hz729/ONddckxVWWCFJMn/+/Bx33HE5/vjj8+yzzzZ4kQAAAAAAAItT75EdkyZNyve///1K0JEkK6ywQk499dRMmjSpQYsDAAAAAAD4IvUOOzbeeOPKWh2f9eKLL2aDDTZokKIAAAAAAADqqk7TWD3zzDOV/x40aFBOOumkTJo0KVtuuWWS5O9//3uuuuqqXHDBBUunSgAAAAAAgMWoU9ix4YYbplQqpVwuV9pOP/30hY777ne/m29/+9sNVx0AAAAAAMAXqNM0VpMnT86rr76ayZMnf+7r1VdfbdDi5s+fn7POOis9e/ZMq1at0qtXr/zkJz+pFbqUy+WcffbZWWWVVdKqVav069cvEydObNA6AAAAAACA5VedRnZ07959adexSBdeeGGuueaajBgxIn379s3jjz+eww8/PFVVVRk0aFCS5KKLLsoVV1yRESNGpGfPnjnrrLPSv3//vPDCC2nZsmWj1A0AAAAAACw7dQo7/tNbb72Vhx56KNOmTUtNTU2tfQtCiIbw8MMPZ6+99so3vvGNJEmPHj3y61//Oo8++miSf4/qGDZsWH784x9nr732SpLccsst6dKlS0aOHJkDDjigwWoBAAAAAACWT/UOO26++eYcffTRad68eTp27JhSqVTZVyqVGjTs2HrrrTN8+PC8/PLLWWuttfL000/noYceyqWXXprk39NrTZ06Nf369aucU1VVlS222CKPPPLIYsOOOXPmZM6cOZXtGTNmNFjNAAAAAADAslXvsOOss87K2WefncGDB6dJkzot+fGlnXnmmZkxY0bWWWedrLDCCpk/f35++tOf5sADD0ySTJ06NUnSpUuXWud16dKlsm9Rhg4dmvPOO2/pFQ4AAAAAACwz9U4rZs+enQMOOGCpBx1Jcvvtt+fWW2/NbbfdlieffDIjRozIxRdfnBEjRixRv4MHD8706dMrrzfeeKOBKgYAAAAAAJa1eicWAwcOzB133LE0alnID37wg5x55pk54IADst566+Xggw/OKaeckqFDhyZJqqurkyTvvPNOrfPeeeedyr5FadGiRdq1a1frBQAAAAAAFFO9p7EaOnRodt9994waNSrrrbdemjVrVmv/gvU0GsLs2bMXGkGywgorVBZF79mzZ6qrq3P//fdnww03TPLv9TfGjx+fY489tsHqAAAAAAAAll9fKuwYPXp01l577SRZaIHyhrTHHnvkpz/9abp165a+ffvmqaeeyqWXXpojjjiicr2TTz45559/ftZcc8307NkzZ511Vrp27Zq99967QWsBAAAAAACWT/UOOy655JLceOONOeyww5ZCObVdeeWVOeuss3Lcccdl2rRp6dq1a44++uicffbZlWNOP/30zJo1K0cddVQ+/PDDbLvtthk1alRatmy51OsDAAAAAAAaX73DjhYtWmSbbbZZGrUsZMUVV8ywYcMybNiwxR5TKpUyZMiQDBkyZJnUBAAAAAAALF/qvUD5SSedlCuvvHJp1AIAAAAAAFBv9R7Z8eijj+aBBx7IH//4x/Tt23ehBcrvvPPOBisOAAAAAADgi9Q77Gjfvn2++c1vLo1aAAAAAAAA6q3eYcdNN920NOoAAAAAAAD4Uuq9ZgcAAAAAAMDypN4jO3r27JlSqbTY/a+++uoSFQQAAAAAAFAf9Q47Tj755Frbc+fOzVNPPZVRo0blBz/4QUPVBQAAAAAAUCf1DjtOOumkRbZfddVVefzxx5e4IAAAAAAAgPposDU7BgwYkN///vcN1R0AAAAAAECdNFjY8bvf/S4dOnRoqO4AAAAAAADqpN7TWG200Ua1Figvl8uZOnVq3n333Vx99dUNWhwAAAAAAMAXqXfYsffee9fabtKkSTp16pQdd9wx66yzTkPVBQAAAAAAUCf1DjvOOeecpVEHAAAAAADAl9Jga3YAAAAAAAA0hjqP7GjSpEmttToWpVQqZd68eUtcFAAAAAAAQF3VOey46667FrvvkUceyRVXXJGampoGKQoAAAAAAKCu6hx27LXXXgu1vfTSSznzzDNz991358ADD8yQIUMatDgAAAAAAIAv8qXW7Hjrrbdy5JFHZr311su8efMyYcKEjBgxIt27d2/o+gAAAAAAAD5XvcKO6dOn54wzzsgaa6yR559/Pvfff3/uvvvurLvuukurPgAAAAAAgM9V52msLrroolx44YWprq7Or3/960VOawUAAAAAALCs1TnsOPPMM9OqVausscYaGTFiREaMGLHI4+68884GKw4AAAAAAOCL1DnsOOSQQ1IqlZZmLQAAAAAAAPVW57Dj5ptvXoplAAAAAAAAfDn1WqAcAAAAAABgeSPsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEt92HHm2++mYMOOigdO3ZMq1atst566+Xxxx+v7C+Xyzn77LOzyiqrpFWrVunXr18mTpzYiBUDAAAAAADL0nIddnzwwQfZZptt0qxZs/z5z3/OCy+8kEsuuSQrrbRS5ZiLLrooV1xxRa699tqMHz8+bdq0Sf/+/fPJJ580YuUAAAAAAMCy0rSxC/g8F154YVZbbbXcdNNNlbaePXtW/rtcLmfYsGH58Y9/nL322itJcsstt6RLly4ZOXJkDjjggGVeMwAAAAAAsGwt1yM7/vCHP2TTTTfNfvvtl86dO2ejjTbK9ddfX9k/efLkTJ06Nf369au0VVVVZYsttsgjjzyy2H7nzJmTGTNm1HoBAAAAAADFtFyHHa+++mquueaarLnmmhk9enSOPfbYDBo0KCNGjEiSTJ06NUnSpUuXWud16dKlsm9Rhg4dmqqqqsprtdVWW3o3AQAAAAAALFXLddhRU1OTjTfeOD/72c+y0UYb5aijjsqRRx6Za6+9don6HTx4cKZPn155vfHGGw1UMQAAAAAAsKwt12HHKquskj59+tRq6927d15//fUkSXV1dZLknXfeqXXMO++8U9m3KC1atEi7du1qvQAAAAAAgGJarsOObbbZJi+99FKttpdffjndu3dP8u/Fyqurq3P//fdX9s+YMSPjx4/PVltttUxrBQAAAAAAGkfTxi7g85xyyinZeuut87Of/Sz7779/Hn300QwfPjzDhw9PkpRKpZx88sk5//zzs+aaa6Znz54566yz0rVr1+y9996NWzwAAAAAALBMLNdhx2abbZa77rorgwcPzpAhQ9KzZ88MGzYsBx54YOWY008/PbNmzcpRRx2VDz/8MNtuu21GjRqVli1bNmLlAAAAAADAsrJchx1Jsvvuu2f33Xdf7P5SqZQhQ4ZkyJAhy7AqAAAAAABgebFcr9kBAAAAAADwRYQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACF1rSxCwAAAACWrZ2b7LfEfYypuaMBKgEAaBhGdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFVqiw44ILLkipVMrJJ59cafvkk09y/PHHp2PHjmnbtm323XffvPPOO41XJAAAAAAAsEwVJux47LHHct1112X99dev1X7KKafk7rvvzh133JEHH3wwb731Vr75zW82UpUAAAAAAMCyVoiwY+bMmTnwwANz/fXXZ6WVVqq0T58+PTfccEMuvfTSfP3rX88mm2ySm266KQ8//HD+/ve/N2LFAAAAAADAslKIsOP444/PN77xjfTr169W+xNPPJG5c+fWal9nnXXSrVu3PPLII8u6TAAAAAAAoBE0bewCvshvfvObPPnkk3nssccW2jd16tQ0b9487du3r9XepUuXTJ06dbF9zpkzJ3PmzKlsz5gxo8HqBQAAAAAAlq3lemTHG2+8kZNOOim33nprWrZs2WD9Dh06NFVVVZXXaqut1mB9AwAAAAAAy9ZyHXY88cQTmTZtWjbeeOM0bdo0TZs2zYMPPpgrrrgiTZs2TZcuXfLpp5/mww8/rHXeO++8k+rq6sX2O3jw4EyfPr3yeuONN5bynQAAAAAAAEvLcj2N1U477ZRnn322Vtvhhx+eddZZJ2eccUZWW221NGvWLPfff3/23XffJMlLL72U119/PVtttdVi+23RokVatGixVGsHAAAAAACWjeU67FhxxRWz7rrr1mpr06ZNOnbsWGkfOHBgTj311HTo0CHt2rXLiSeemK222ipbbrllY5QMAAAAAAAsY8t12FEXl112WZo0aZJ99903c+bMSf/+/XP11Vc3dlkAAAAAAMAyUriw4y9/+Uut7ZYtW+aqq67KVVdd1TgFAQAAAAAAjWq5XqAcAAAAAADgiwg7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIXWtLELAICGUDN1rSU6v0n1yw1UCQAAAADLmpEdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACk3YAQAAAAAAFJqwAwAAAAAAKDRhBwAAAAAAUGjCDgAAAAAAoNCEHQAAAAAAQKEJOwAAAAAAgEITdgAAAAAAAIUm7AAAAAAAAApN2AEAAAAAABSasAMAAAAAACg0YQcAAAAAAFBowg4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQmva2AUAAADAslQzda0lOr9J9csNVAkAAA3FyA4AAAAAAKDQhB0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2AAAAAAAAhSbsAAAAAAAACq1pYxcAQPHVTF1ric5vUv1yA1UCAAAAwH8jIzsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKzZodFI61AQAAoLj8/zwAAEuDkR0AAAAAAEChCTsAAAAAAIBCE3YAAAAAAACFJuwAAAAAAAAKTdgBAAAAAAAUmrADAAAAAAAoNGEHAAAAAABQaMIOAAAAAACg0IQdAAAAAABAoQk7AAAAAACAQhN2/H/t3XlYVfW6wPF3bzaDgICQsFFBnAiHBscjYNqgGHmyTE+W5ZRWlnaOQ1o2aHWazO6trlOmqfVo6c3hVB6HDIcOjolhOYE5ZoqmpjgisN/7Bw87ULRusvZaq76f59nP095Leb8hLBb8WGsBAAAAAAAAAABbY7EDAAAAAAAAAADYGosdAAAAAAAAAADA1lxmB8BePHmJV/X3ne7cSioBAOthHwkAAAAAAGAOzuwAAAAAAAAAAAC2xmIHAAAAAAAAAACwNRY7AAAAAAAAAACArXHPDgAAAAAAAACwkA7Ov13121jm+aQSSgD74MwOAAAAAAAAAABgayx2AAAAAAAAAAAAW2OxAwAAAAAAAAAA2Jrl79nx2muvyfz582XHjh1SpUoVSUlJkTFjxsi1117r/TPnz5+XYcOGyezZs6WgoEA6duwoEydOlJiYGBPLAQDwLU9e4lW/Dac7txJKAABWdbVfK/g6AQAAAKuy/Jkdq1atkoEDB8q6detk2bJlUlhYKGlpaXLmzBnvnxkyZIh8/vnn8sknn8iqVavk4MGDcs8995hYDQAAAAAAAAAAfMXyZ3YsWbKk3PMZM2ZIdHS0ZGVlSdu2beXkyZPy/vvvy0cffSS33nqriIhMnz5dGjZsKOvWrZPWrVubkQ0AAAAAAAAAAHzE8md2XOzkyZMiIhIZGSkiIllZWVJYWCjt27f3/pmkpCSJj4+XtWvXmtIIAAAAAAAAAAB8x/JndpTl8Xhk8ODBkpqaKk2aNBERkby8PAkICJCIiIhyfzYmJkby8vIqfDsFBQVSUFDgfZ6fn29YMwAAAAAAAAAAMJatzuwYOHCgbNmyRWbPnn1Vb+e1116T8PBw7yMuLq6SCgEAAAAAAAAAgK/ZZrFj0KBBsnDhQlmxYoXUqlXL+7rb7ZYLFy7IiRMnyv35w4cPi9vtrvBtjRw5Uk6ePOl9/PDDD0amAwAAAAAAAAAAA1l+sUNVZdCgQbJgwQJZvny51KlTp9z25s2bi7+/v2RkZHhfy8nJkf3790tycnKFbzMwMFDCwsLKPQAAAAAAAAAAgD1Z/p4dAwcOlI8++kg+/fRTqVq1qvc+HOHh4VKlShUJDw+Xfv36ydChQyUyMlLCwsLkiSeekOTkZGndurXJ9QAAAAAAAAAAwGiWX+yYNGmSiIjcfPPN5V6fPn269OnTR0RE3nrrLXE6ndK1a1cpKCiQjh07ysSJEyu1w5OXeFV/3+nOraQSAAAAAAAAAABQluUXO1T1V/9MUFCQTJgwQSZMmOCDIgAAAAAAAAAAYCWWv2cHAAAAAAAAAADAlbDYAQAAAAAAAAAAbI3FDgAAAAAAAAAAYGssdgAAAAAAAAAAAFtjsQMAAAAAAAAAANgaix0AAAAAAAAAAMDWWOwAAAAAAAAAAAC2xmIHAAAAAAAAAACwNRY7AAAAAAAAAACArbHYAQAAAAAAAAAAbI3FDgAAAAAAAAAAYGssdgAAAAAAAAAAAFtzmR0AAAAAAH8GnrzEq34bTnduJZQAAAAAfzyc2QEAAAAAAAAAAGyNxQ4AAAAAAAAAAGBrLHYAAAAAAAAAAABb454dAACg0nA9egAAAAAAYAbO7AAAAAAAAAAAALbGYgcAAAAAAAAAALA1FjsAAAAAAAAAAICtsdgBAAAAAAAAAABsjcUOAAAAAAAAAABgayx2AAAAAAAAAAAAW2OxAwAAAAAAAAAA2BqLHQAAAAAAAAAAwNZcZgcAsCdPXuJVvw2nO7cSSgAAAAAAAAD82XFmBwAAAAAAAAAAsDUWOwAAAAAAAAAAgK2x2AEAAAAAAAAAAGyNxQ4AAAAAAAAAAGBr3KAcAACgknnyEq/q7zvduZVUAgAAAADAnwNndgAAAAAAAAAAAFtjsQMAAAAAAAAAANgaix0AAAAAAAAAAMDWuGcHYENXey14Ea4HDwAAAAAAAOCPgzM7AAAAAAAAAACArbHYAQAAAAAAAAAAbI3FDgAAAAAAAAAAYGssdgAAAAAAAAAAAFtjsQMAAAAAAAAAANgaix0AAAAAAAAAAMDWWOwAAAAAAAAAAAC2xmIHAAAAAAAAAACwNRY7AAAAAAAAAACArbHYAQAAAAAAAAAAbI3FDgAAAAAAAAAAYGssdgAAAAAAAAAAAFtjsQMAAAAAAAAAANgaix0AAAAAAAAAAMDWWOwAAAAAAAAAAAC2xmIHAAAAAAAAAACwNRY7AAAAAAAAAACArbnMDgAAAKhMnrzEq34bTnduJZQAAAAAAABf4cwOAAAAAAAAAABgayx2AAAAAAAAAAAAW2OxAwAAAAAAAAAA2BqLHQAAAAAAAAAAwNZY7AAAAAAAAAAAALbGYgcAAAAAAAAAALA1FjsAAAAAAAAAAICtsdgBAAAAAAAAAABszWV2AAAAACqfJy/xqv6+051bSSUAAAAAABiPMzsAAAAAAAAAAICtsdgBAAAAAAAAAABsjcUOAAAAAAAAAABga9yzw0a49jYAAAAAAAAAAJfizA4AAAAAAAAAAGBrLHYAAAAAAAAAAABbY7EDAAAAAAAAAADYGosdAAAAAAAAAADA1rhBOQDb8uQlXvXbcLpzK6EEAAAAAAAAgJk4swMAAAAAAAAAANgaix0AAAAAAAAAAMDWWOwAAAAAAAAAAAC2xmIHAAAAAAAAAACwNRY7AAAAAAAAAACArbHYAQAAAAAAAAAAbI3FDgAAAAAAAAAAYGssdgAAAAAAAAAAAFtzmR1QWSZMmCBjx46VvLw8ueGGG2TcuHHSqlUrs7MAwHCevMSr+vtOd24llQAAKmL2fvpq51dGAwAAAAAY7Q9xZsecOXNk6NChMnr0aNm0aZPccMMN0rFjRzly5IjZaQAAAAAAAAAAwGB/iMWO//7v/5aHH35Y+vbtK40aNZJ3331XgoODZdq0aWanAQAAAAAAAAAAg9n+MlYXLlyQrKwsGTlypPc1p9Mp7du3l7Vr11b4dwoKCqSgoMD7/OTJkyIikp+ff9k5nlPFV9XpDL782/6taLBOg9mu9n0gcvXvBxqswwqfEzRYo8EKnxM0VE7DH+HjEb8w+9/CCp8TKGGFfwuzPx5psI4iLbzqt3Gl76Fp+HM1XO18GmiozPk0/LaG0m2qetVzAKuw/WLH0aNHpbi4WGJiYsq9HhMTIzt27Kjw77z22mvy4osvXvJ6XFycIY0lwg18278VDSWs0GAFVng/0GANVngf0FCChhI0mD9fxBoNKGGFfwsrNKCE2f8WZs8XocE6wsPNfz/QYI0Gs+fTQAMNv7/h1KlTlmgFKoPtFzt+j5EjR8rQoUO9zz0ejxw/flyioqLE4XD8v99efn6+xMXFyQ8//CBhYWGVmUoDDbZsMHs+DTTQQINV59NAAw3WbDB7Pg000GDNBrPn00ADDdZsMHt+ZTWoqpw6dUpq1KhRyXWAeWy/2HHNNdeIn5+fHD58uNzrhw8fFrfbXeHfCQwMlMDAwHKvRUREXHVLWFiYaTs5GmiwYoPZ82mggQYarDqfBhposGaD2fNpoIEGazaYPZ8GGmiwZoPZ8yujgTM68Edj+xuUBwQESPPmzSUjI8P7msfjkYyMDElOTjaxDAAAAAAAAAAA+ILtz+wQERk6dKj07t1bWrRoIa1atZK3335bzpw5I3379jU7DQAAAAAAAAAAGOwPsdjRvXt3+emnn2TUqFGSl5cnN954oyxZsuSSm5YbJTAwUEaPHn3JpbF8iQYarNRg9nwaaKCBBqvOp4EGGqzZYPZ8GmigwZoNZs+ngQYarNlg9nyrNABW5FBVNTsCAAAAAAAAAADg97L9PTsAAAAAAAAAAMCfG4sdAAAAAAAAAADA1ljsAAAAAAAAAAAAtsZiBwAAAAAAAAAAsDUWOwAAAAAAAAAAgK2x2AEAAOAjBQUFUlBQYHYGLGblypVy7tw5szNMVVBQILt27frTf34cPnxY8vLyfD63uLhYDh8+LD/99JPPZ5c6efKk5OTkSE5Ojpw8edK0DpSnqlJcXGxqw4wZM/7UHxM7d+6UjIwM+f77781O8bmLP/Y2bNgg69at8/nXiv3798v69evl66+/lmPHjvl0dlkcR+JiHEMCl2KxoxL17dtXDh48aNr8wsJC2blzp2kHgidOnJApU6bI888/L1OnTjW8Iysry9C3/1sdOXJEli9f7v3/PXz4sLzxxhvy+uuvy3fffeezjt27d8uHH34oY8aMkbFjx8q8efMkPz/fZ/NFRPLy8uTTTz+VyZMny+TJk+XTTz815Rv2ipw5c0a++uorszN85uJvDNavXy9fffWVFBYWmlRk7j7yz7Z/FGEfeTEz95HLli2TO+64Q6pVqybBwcESHBws1apVkzvuuEO+/PJLnzRcyfbt26Vu3bqGz9m8ebO8/PLLMnHiRDl69Gi5bfn5+fLQQw8Z3jB16lTp3bu3TJ8+XURE5syZIw0bNpS6devK6NGjDZ9/OWlpabJ3716fzDpy5Ei559nZ2dK7d29JTU2Vbt26ycqVKw1vmDFjhqxdu1ZERM6fPy/9+vWTkJAQSUxMlNDQUBkwYIChP8y57rrr5J///Kf88MMPhs34NcePH5du3bpJfHy8PPbYY1JcXCz9+/eX2NhYqVmzpqSkpMihQ4cM7/j3v/8tbdu2lZCQEKlRo4a43W6JiIiQnj17yv79+w2fL1LyedmoUSOJjIyURo0alfvv999/3ycNV7J582bx8/MzfM6iRYukf//+MmLECNmxY0e5bT///LPceuuths4vKiqS5557Ttq1a+fdH44dO1ZCQ0MlODhYevfuLRcuXDC04XIeeeQRnxxDbtiwodwx9MKFC6Vdu3ZSs2ZNadGihXz44YeGN7z22muSkZEhIiX/7u3bt5drr71WOnToINdee62kp6fLiRMnDG2oWrWq9OvXT9asWWPonCvZt2+ftGjRQgIDAyU9PV3y8/OlQ4cO0rp1a0lJSZFGjRpJbm6u4R0TJ06U2rVrS506dSQlJUVat24t0dHR0qZNG58da3McWcLs40iOIQEbUfy/bd68ucKHv7+/LliwwPvcSGPGjNGzZ8+qqmpRUZEOGzZMAwIC1Ol0qsvl0r59++qFCxcMbejSpYt+8sknqqq6ZcsWveaaa7R69er6l7/8RWNiYtTtduu2bdsMm+9wOLRevXr6yiuv6I8//mjYnCtZsWKFhoSEqMPhULfbrdnZ2VqrVi1t0KCBXnvttRoYGKhLly41tOH06dParVs3dTgc6nA41Ol0qtvtVj8/Pw0NDdXx48cbOr+04YEHHlA/Pz91uVwaHR2t0dHR6nK51M/PTx988EE9c+aM4R1Xkp2drU6n09AZFy5c0OHDh2u9evW0ZcuW+v7775fbnpeXZ3jDwYMHNTU1Vf38/LRt27Z6/Phx7dSpk/fjIzExUQ8ePGhog9n7SPaPJdhHljB7Hzljxgx1uVx633336fTp03XRokW6aNEinT59ut5///3q7++vH374oaENv8YX+8elS5dqQECANm7cWOPj4zUqKkqXL1/u3e6L/eNbb72lISEhes8992hsbKy+/PLLGhUVpS+//LK++OKLGhYWppMnTza0oWnTphU+HA6HNmzY0PvcSE6nUw8fPqyqqqtXr1Z/f39t166dDh8+XDt06KAul0tXrVplaEOdOnV03bp1qqr65JNPakJCgs6fP1+3b9+u//rXvzQxMVGHDx9u2HyHw6FRUVHq5+enHTt21Llz52phYaFh8yry0EMPaZMmTXTcuHHarl07veuuu/T666/XzMxMXbNmjbZs2VJ79eplaMOHH36oVatW1WHDhumzzz6rbrdbn376aZ00aZK2a9dOr7nmGs3NzTW04Y033tDg4GB9+umndcWKFbpt2zbdtm2brlixQkeOHKkhISE6duxYQxt+TXZ2tjocDkNnzJo1S/38/LRTp07apk0bDQoK0pkzZ3q3+2If+dxzz2lMTIwOHTpUGzVqpAMGDNC4uDidOXOmfvDBB1qzZk0dM2aMoQ3VqlWr8OFwODQ8PNz73Chl94+fffaZOp1O7dWrl06YMEH79++vLpdL58+fb9h8VdVatWrppk2bVFW1f//+2rRpU920aZOeO3dOs7OztXXr1tqvXz9DGxwOhzZu3FgdDocmJSXpm2++qUeOHDF05sW6du2q7dq1088//1zvvfdeTU1N1ZtvvlkPHDigBw8e1I4dO+rdd99taMPYsWO1Ro0aOm7cOJ0yZYo2bNhQX3rpJV28eLH27NlTg4OD9euvvza0gePIEmYfR3IMCdgLix2/Q+kPS0p/cFL2Ufq60Tv7sgdiY8eO1WrVqum0adN069atOnPmTI2OjvbJwej27dtVVTU9PV179OihBQUFqlryQ99+/fppWlqaYfMdDoc+/PDD3h+qd+rUSRcsWKBFRUWGzbxYmzZtdODAgXrq1CkdO3as1qxZUwcOHOjd/uSTT2pKSoqhDY888oimpqbqd999pzt37tRu3brpiBEj9MyZM/r+++9rcHCwzpo1y9CGfv36aYMGDXTJkiXl3v9FRUW6dOlSTUxM1P79+xva8Gt8cRA2evRojYmJ0bFjx+qzzz6r4eHh+sgjj3i35+XlGf7Ncs+ePTUlJUU/++wz7d69u6akpOhNN92kBw4c0H379mlqamq5j1EjmL2PZP9Ygn1kCbP3kQ0aNLjigsqECRO0fv36hs1XVR0yZMgVHw8++KDh+8fk5GR95plnVFXV4/HomDFjNDQ0VBcvXqyqvvlBXlJSkvffetOmTepyuXTq1Kne7VOnTtXmzZsb2uByufT222/XF154wfsYPXq0Op1Offzxx72vGcnhcHj3kR06dNCHHnqo3PZ//OMfeuuttxraEBgYqPv27VNV1cTERO/HQalVq1ZpfHy8YfMdDof++OOPumDBAr3zzjvV5XJp9erVddiwYYYvRJeKjY3V1atXq+ovxwdffPGFd3tmZqbWrFnT0IakpCSdPXu29/nXX3+ttWrVUo/Ho6qq3bt31y5duhjaEB8fr3PmzLns9tmzZ2tcXJyhDV26dLni49ZbbzV8/3TjjTfqO++8430+Z84cDQkJ8e6jfLGPrFu3rn7++eeqqrpz5051Op3lPj7mzJmjTZo0MbQhNDRUO3XqpDNmzPA+pk+frn5+fvrKK694XzNK2f1jmzZt9Omnny63/ZVXXtHWrVsbNl+1ZP+4d+9eVVVNSEi4ZPF548aNGhsba2hD6fshOztbBw0apJGRkRoQEKD33HOPLlq0yLuPMFL16tX1m2++UVXVEydOqMPh0P/85z/e7VlZWRoTE2NoQ0JCgi5atMj7PCcnR6OioryL43//+9+1Q4cOhjZwHFnC7ONIjiEBe2Gx43e44YYbtFOnTrp9+3bdu3ev7t27V/fs2aMul0uXLVvmfc1IZQ/EmjZteskq8syZM7Vx48aGNlSpUkW///57VS35hq30N1BK5eTkaHh4uGHzS98HhYWFOnfuXL3jjjvUz89PY2JidMSIEZqTk2PY7FJhYWHe90FhYaG6XC7vQZmqam5urqHvA1XVa665Rjdu3Oh9fvz4cQ0KCvKeSTF+/Hi98cYbDW2IiIjwfsNekczMTI2IiDC04XK/CVb6CAsLM/wgrH79+t5vElVLvlGsX7++9unTRz0ej0++UY2NjdW1a9eqquqxY8fU4XDol19+6d2ekZGhdevWNbTB7H0k+8cS7CNLmL2PDAwM1B07dlx2+44dOzQoKMiw+aolC4DNmjXTm2++ucJHixYtDN83lf1YKDVr1iwNCQnRzz//3Cf7xypVqnh/wK5a8m+zZcsW7/OdO3ca/rUqMzNT69Wrp6NGjdLi4mLv6y6XS7du3Wro7FJl95Flv2aUKj0bzUi1a9f2/kZmzZo1L/nN2G3btmlISIhh88u+D1RLzop89dVXtUGDBup0OjU5OfmSszMrW3BwcLmvhf7+/vrdd995n+/evdvQ94FqyefEnj17yr3mcrm8ZwOuX7/e8M+JoKCgKy4wbd26VatUqWJog8vl0vT0dO3Tp0+Fj86dOxu+fwoJCdHdu3eXe2358uUaGhqqkyZN8sk+MigoSPfv31/ueekvbqiWfExWrVrV0IadO3d6z2o6deqU93Vf7SPL7huio6PLHT+olnzNNvpzIjExURcuXKiqJWfBXfw91jfffKNhYWGGNly8jzx//rx+9NFHetttt6nT6dRatWrp888/b2hD1apVvZ8TxcXF6nK5NDs727t9586dhn88BgcHl9tHejwedblc3rPks7OzNTQ01NAGjiNLmH0cyTEkYC/cs+N32LBhg9SvX1+6du0qx48fl9q1a0tCQoKIiNSoUUNq164ttWvXNrzD4XCISMnNslJSUsptS0lJkT179hg6//rrr5fly5eLiIjb7ZZ9+/aV275v3z6pUqWKoQ0iIi6XS7p27Sr//ve/Zd++fTJw4ECZO3euNGzYUNq2bWvo7ICAADl//ryIiFy4cEE8Ho/3uYjIuXPnxN/f39CGoqIiCQsL8z4PDQ2VoqIiOXPmjIiUXMPx4uv+VjaPxyMBAQGX3R4QECAej8fQhoKCAnnooYfkrbfeqvAxbNgwQ+eLiPz444/SpEkT7/P69evLypUrZc2aNdKzZ0+f3Nzx559/lpo1a4qISGRkpAQHB5fbH9WvX9/wa4BbYR/J/vEX7CPN3Uc2btz4itecnzZtmjRq1Miw+SIln/dDhgyRFStWVPiYMmWKofNFRAIDAy+5xniPHj1k6tSp0r17d1mwYIHhDcHBwd5/dxGR6tWrS2hoaLk/U1RUZGhDamqqZGVlSW5urqSkpMiuXbsMnXc5p06dkvz8fAkKCpLAwMBy24KCguTs2bOGzn/ggQfk2WeflRMnTkjPnj3lpZdektOnT4uIyNmzZ+WFF16Q1NRUw+aXfo0oFRsbKyNHjpTc3FzJyMiQevXqyd///nfD5ouINGjQQBYuXCgiIosXL5agoCD54osvvNuXLl0qderUMbQhISFBNm7c6H2+adMmcTqdEhMTIyIlxxFG3+urZcuW8vrrr1f4uVdcXCxjxoyRli1bGtrQsGFD6dq1q0yfPr3Cx4svvmjofBGRsLAwOXz4cLnXbrnlFlm4cKEMHz5cxo0bZ3hDeHh4uf10s2bNpGrVqt7nBQUFl3zuVLb69evLmjVrxO12y4033iirV682dF5Ftm3bJt9++61UqVKlwu9fjP468fDDD8vw4cPl+++/l0GDBsmTTz7p/VqxZ88eGTJkiKSlpRnacPG/c2BgoNx///3y5Zdfyq5du6RPnz4yY8YMQxsaN24s06ZNExGRDz74QKKiomT27Nne7R9//LEkJiYa2pCYmCjLli3zPl+xYoUEBASI2+0WkZKvl0Z/TnAcWcLs40iOIQGbMXu1xc4WLVqktWrV0ldffdX72wa+/K28V155Rd955x2NjY295PTWzZs3G3o9U1XVhQsXamRkpE6fPl2nT5+uCQkJOnXqVF29erVOmzZN4+LiDL3ectlL1VTkyy+/1B49ehg2X1X1rrvu0r/+9a+amZmpjzzyiLZo0UI7deqkp0+f1jNnzmi3bt309ttvN7ShQ4cO5S4LM3bs2HKnNm/atMnw39Ds0aOH93qyF9u0aZM2b95cH3jgAUMbUlJS9O23377sdl9cxqpOnTrlzqIo9eOPP2piYqJ26NDB8Ib4+Hhdv3699/lTTz2lx44d8z7Pzs42/OOhlFn7SPaPJdhHljB7H1l635LrrrtOhwwZoq+//rq+/vrrOmTIEL3++us1NDTU8Psj9OjRQwcPHnzZ7b64Hn2HDh0ue939jz76SP39/Q3fP6amppa7JMvFPv/8c8Mv0VLWtGnT1O126+TJk9Xf39+nx5BOp9N7WcH33nuv3PZPP/3U8EtiFBQUaOfOnbVatWraoUMHDQoK0uDgYG3QoIGGhIRofHy8oWefXfxbyxU5efKkYfNVS84y9PPz0/r162tgYKB+8sknWqNGDb333nv1vvvu04CAAMPvKTR+/HgNDw/XESNG6KhRo7RGjRrl7gUwc+ZMw6//vXnzZnW73RoVFaVdunTRAQMG6IABA7RLly4aFRWlsbGx5c54MUKfPn308ccfv+z2bdu2aUJCgqENd911l44aNarCbaVfR4zeR95yyy1XvETU//7v/xp+mZayMjIyND4+XkeOHOmzfeTFl2J96623ym3/+OOPtVGjRoZ3PPHEE+rv769JSUkaFBSkTqfTe/+5Fi1a6KFDhwyd/1v2kUZfymrJkiUaFBSkAQEBGhQUpKtWrdLExERt1aqVtm7dWv38/K54CbzKMGfOHPX399d7771Xe/XqpaGhoeUubfbuu+9qcnKyoQ0cR5Yw+ziSY0jAXljsuEp5eXmanp6uN910k08XO2rXrq0JCQnex8UHYm+//bbh1xNVVZ07d67WqlXrkuvzBwUF6eDBgw29NvxvOQgzWm5urjZo0MB7U6gDBw5o586d1eVyea//nJWVZWhDVlaWRkZGqtvt1vj4eA0ICNCPP/7Yu338+PGG3+Dy+PHjevvtt6vD4dDIyEhNSkrSpKQkjYyMVKfTqenp6frzzz8b2vDKK69c8RqV+/fv1z59+hja0K9fv0uufV7qwIEDWr9+fcO/Ue3cufMVF33Gjx9v+LXYyzJjH8n+sQT7yBJW2Efu2bNHR4wYoW3bttXExERNTEzUtm3b6lNPPXXJJWSMcOjQIcMvr/lr5s+ff8VvlGfNmqU333yzoQ2ZmZnlLqN2sQkTJui4ceMMbbhYbm6utmzZUh0Oh8+OIVeuXFnucfGiwttvv61vvPGGT1oWL16sjz/+uN5+++2alpamvXv31vfee09Pnz5t6Nw+ffpofn6+oTN+i8zMTH3zzTe9l6nZunWr9uzZU7t27WrofQnKmjhxoqakpGjz5s31mWee0XPnznm35ebmlruMkVHy8/N14sSJ2qtXL01LS9O0tDTt1auXTpo0yfBFJ9WSS/SUXtrQLCtXrtRXX331stuXL19u+HFsTk7OJZfSKmvWrFmG/3D5YkePHtUuXbpoRETEFS/lU1lKL7Va+jh69Gi57R988IF+8MEHhneoliyyvfHGGzpgwAB95JFHdPTo0frFF1/45H4ZL7zwgumfE6olx09z5871Hivl5eXp888/r8OGDSt3c2ojLVq0SHv06KFdu3a95JcDjh49esnHiBE4jjT/OJJjSMBeHKqqZp9d8kfwP//zP7JixQoZN26c1KpVy+wcWbdunQQGBkrTpk0Nn1VcXCybNm2S3bt3i8fjkdjYWGnevHm5056NsGrVKklNTRWXy2XonN/i2LFjEhUV5X2ekZEh586dk+Tk5HKvG+XQoUOycOFCKSgokFtvvdXwU1kvZ/v27bJu3TrJy8sTkZLL9yQnJ0tSUpIpPb62b98+2bFjh3Ts2LHC7QcPHpRly5ZJ7969fVz2iw0bNkhwcHC5y235gpX2kb7eP2ZlZcmePXt8un8UYR9ZllX2kUBFPB6PnDp1SsLCwgy/HAYAAAD+GDiGBCrGYgcAAIDBioqKZOvWrd7F4NjYWGnYsKHh9yy5UoPb7ZZGjRrRQMOftsHs+TT8usLCQjl06JDEx8eb1lBUVCQHDx6kgQbT54vwOUEDAODXmP/rnja1efNmycrKkptvvlnq1q0rW7dulQkTJojH45EuXbpc9je7afhjzbdawy233CJ16tQxpUFEZPny5ZKZmSmHDh0Sp9MpdevWlc6dO0uDBg18Mp8GazfUq1dP7rzzTp81WPF9QMOfr8Hj8cioUaNkwoQJcvLkyXLbwsPDZdCgQfLiiy+K0+mkgQYafNRg9nwafrtt27ZJs2bNpLi42LSGrVu30kCDJeaL8DnxZ2yYOHGizJ8/XyIjI+XRRx+V2267zbvt6NGj0qpVK9m9ezcNBjeYPd8qDYAtmHsVLXuaN2+e+vn5aVRUlIaGhuqyZcs0IiJC27dvrx07dlQ/Pz+dNWsWDQY3mD2fhl8cPnxYW7VqpU6nU10ulzqdTm3evLm63W718/Mz/EbMNNBgpfk00FDW8OHDtXr16vruu+/qnj179OzZs3r27Fnds2ePTp48WaOjo3XEiBE00ECDDxvMnk/Db5ednW34/c5ooMEu82n48zW88847GhwcrAMHDtQHH3xQAwICyt3fJy8vjwYfNJg93yoNgF2w2PE7NGvWTF9++WVVVf344481IiJCX3rpJe/2N998U2+88UYaDG4wez4Nv+jevbvefffdevLkST1//rwOGjTIe8PfjIwMjYqKuuJNs2mg4Y80nwYayoqJidElS5ZcdvuSJUs0OjrasPk00ECD9ebT8IumTZte8ZGUlGT4D29ooMEq82mg4WKNGjUq94uLq1ev1urVq+vzzz+vqr75ATcN5s+3SgNgF9yz43cIDQ2VLVu2SEJCgqiqBAYGSlZWllx33XUiIrJ792654YYb5NSpUzQY2GD2fBp+ER4eLmvWrJHGjRuLiMiZM2ekWrVqcvToUQkLC5OZM2fKyy+/LDt27KCBBsMbzJ5PAw1lhYSEyLp167z75It9++23kpKSIqdPnzZkPg000GC9+TT8IigoSO677z6pU6dOhdsPHTokU6ZMMfQyMTTQYJX5NNBwseDgYNm2bZskJCR4X9uyZYu0b99e+vbtK4MHD5YaNWrQYHCD2fOt0gDYhnnrLPbldrt148aNqqp6/PhxdTgcumLFCu/2DRs2qNvtpsHgBrPn0/CL6tWr69atW73Pz549q06nU48dO6aqqrt27dLAwEAaaPBJg9nzaaChrDvuuEPT0tL0p59+umTbTz/9pLfffrt26tTJsPk00ECD9ebT8IvmzZvrxIkTL7v9m2++Mfw3VWmgwSrzaaDhYnFxcfrVV19d8vrWrVs1JiZGe/XqRYMPGsyeb5UGwC7Mu9ucjbVv314GDhwos2bNkt69e0taWpqMHDlSduzYITk5OTJ8+HBp06YNDQY3mD2fhl+0adNGRo0aJWfOnJHCwkJ55plnpG7duhIZGSkiIj/99JNUq1aNBhp80mD2fBpoKOvdd9+VgwcPSmxsrDRr1kzS09MlPT1dmjVrJrGxsXLw4EGZNGmSYfNpoIEG682n4RepqamSk5Nz2e1Vq1aVtm3b0kCDTxrMnk8DDRdr06aNzJ8//5LXGzVqJBkZGbJ48WJD59NgjflWaQDsgstY/Q6HDx+Wnj17ytq1ayU1NVXmzJkjzz33nEyYMEEcDofUq1dPFi9eLPXq1aPBwAaz59Pwi927d0taWprs27dPHA6HhISEyCeffCLt27cXEZEZM2ZITk6OvPbaazTQYHiD2fNpoOFiHo9Hli5dKuvWrZO8vDwREXG73ZKcnCxpaWnidBr/uyc00ECDtebTAAD4Nd9++61kZWVJ3759K9y+ZcsWmTdvnowePZoGAxvMnm+VBsAuWOyoRLt375azZ89KUlKSuFwuGkxqMHv+n7Xh7NmzkpmZKRcuXJDWrVvLNddcY/hMGmiw6nwaaAAAAAAAAL7FYgcAAIDBNmzYIGvXri33m9MpKSnSsmVLGmigwaQGs+fTcOWG5ORkadWqFQ00+LzB7Pk00ECDNRvMnm+VBsDqWOz4nc6dOycff/yxZGZmyqFDh8TpdErdunXl7rvvlttuu40GHzWYPZ8GGmiwZoPZ82mgodSRI0eka9eusnr1aomPj5eYmBgRKbkE4f79+yU1NVXmzZsn0dHRNNBAg48azJ5PAw00WLPB7Pk00FBRwz333CNr1qyhwcQGs+dbpQGwDfPujW5fO3fu1Nq1a2t0dLTGxcWpw+HQTp066V/+8hf18/PTv/3tb1pYWEiDwQ1mz6eBBhqs2WD2fBpoKKtr166anJysO3bsuGTbjh07NCUlRbt162bYfBpooMF682mggQZrNpg9nwYaaLBmg9nzrdIA2AWLHb9Denq6Pvroo+rxeFRV9fXXX9f09HRVVc3NzdWEhAQdPXo0DQY3mD2fBhposGaD2fNpoKGs0NBQ3bRp02W3b9y4UUNDQw2bTwMNNFhvPg000GDNBrPn00ADDdZsMHu+VRoAu2Cx43cIDg7W3Nxc7/OCggL19/fXo0ePqqrqv/71L01ISKDB4Aaz59NAAw3WbDB7Pg00lBUVFaUrV6687PYVK1ZoVFSUYfNpoIEG682ngQYarNlg9nwaaKDBmg1mz7dKA2AXTrMvo2VHERERcurUKe/zs2fPSlFRkQQEBIiIyPXXXy+HDh2iweAGs+fTQAMN1mwwez4NNJTVvXt36d27tyxYsEDy8/O9r+fn58uCBQukb9++cv/99xs2nwYaaLDefBpooMGaDWbPp4EGGqzZYPZ8qzQAtmH2aosd9e7dW9u1a6fbt2/X3bt3a/fu3bVp06be7StXrtS4uDgaDG4wez4NNNBgzQaz59NAQ1nnz5/XAQMGaEBAgDqdTg0KCtKgoCB1Op0aEBCgjz32mJ4/f96w+TTQQIP15tNAAw3WbDB7Pg000GDNBrPnW6UBsAuHqqrZCy52c+TIEbnrrrtk/fr14nA4JC4uThYsWCBNmzYVEZG5c+fKoUOH5IknnqDBwAaz59NAAw3WbDB7Pg00VCQ/P1+ysrIkLy9PRETcbrc0b95cwsLCDJ1LAw00WHc+DTTQYM0Gs+fTQAMN1mwwe75VGgCrY7HjKuzcuVMKCgokKSlJXC4XDSY1mD2fBhposGaD2fNpoAEAAAAAAPgO9+y4Cg0aNJAmTZpc8kOTH374QR566CEafNRg9nwaaKDBmg1mz6eBhlLnzp2TzMxM2bZt2yXbzp8/Lx9++KGh82mggQbrzaeBBhqs2WD2fBpooMGaDWbPt0oDYAvmXkXrjyk7O1udTicNJjeYPZ8GGmiwZoPZ82n4czXk5ORo7dq11eFwqNPp1LZt2+qPP/7o3Z6Xl2f4+4AGGmiw1nwaaKDBmg1mz6eBBhqs2WD2fKs0AHbBdRx+h88+++yK23fv3k2DDxrMnk8DDTRYs8Hs+TTQUNZTTz0lTZo0kY0bN8qJEydk8ODB0qZNG1m5cqXEx8cbOpsGGmiw5nwaaKDBmg1mz6eBBhqs2WD2fKs0ALZh9mqLHZWupDocjss+jF5RpcH8+TTQQIM1G8yeTwMNZUVHR+u3337rfe7xeHTAgAEaHx+vu3bt8slvYdFAAw3Wmk8DDTRYs8Hs+TTQQIM1G8yeb5UGwC64Z8fvEBsbK/PnzxePx1PhY9OmTTT4oMHs+TTQQIM1G8yeTwMNZZ07d67cfUIcDodMmjRJ7rzzTmnXrp3k5uYaOp8GGmiw3nwaaKDBmg1mz6eBBhqs2WD2fKs0AHbBYsfv0Lx5c8nKyrrsdofDIapKg8ENZs+ngQYarNlg9nwaaCgrKSlJNm7ceMnr48ePl7vuuks6d+5s2GwaaKDBmvNpoIEGazaYPZ8GGmiwZoPZ863SANiG708msb+vvvpKFy9efNntp0+f1pUrV9JgcIPZ82mggQZrNpg9nwYaynr11Vc1PT39stsfe+wxdTgchs2ngQYarDefBhposGaD2fNpoIEGazaYPd8qDYBdOFQN/pVKAAAAAAAAAAAAA3EZKwAAAAAAAAAAYGssdgAAAAAAAAAAAFtjsQMAAAAAAAAAANgaix0AAAAAAAAAAMDWWOwAAAAAAAAAAAC2xmIHAAAAcJVUVdq3by8dO3a8ZNvEiRMlIiJCDhw4YEIZAAAAAPw5sNgBAAAAXCWHwyHTp0+X9evXy+TJk72v79mzR0aMGCHjxo2TWrVqVerMwsLCSn17AAAAAGBnLHYAAAAAlSAuLk7eeecdefLJJ2XPnj2iqtKvXz9JS0uTpk2bSnp6uoSGhkpMTIz07NlTjh496v27S5YskTZt2khERIRERUXJX//6V9m1a5d3+969e8XhcMicOXOkXbt2EhQUJLNmzTLjfxMAAAAALMmhqmp2BAAAAPBHcffdd8vJkyflnnvukX/+85+ydetWady4sfTv31969eol586dk6eeekqKiopk+fLlIiIyb948cTgccv3118vp06dl1KhRsnfvXsnOzhan0yl79+6VOnXqSEJCgvzXf/2XNG3aVIKCgiQ2Ntbk/1sAAAAAsAYWOwAAAIBKdOTIEWncuLEcP35c5s2bJ1u2bJH//Oc/snTpUu+fOXDggMTFxUlOTo4kJiZe8jaOHj0q1atXl++++06aNGniXex4++235R//+Icv/3cAAAAAwBa4jBUAAABQiaKjo+XRRx+Vhg0byt133y2bN2+WFStWSGhoqPeRlJQkIuK9VNXOnTvl/vvvl7p160pYWJgkJCSIiMj+/fvLve0WLVr49P8FAAAAAOzCZXYAAAAA8EfjcrnE5So51D59+rTceeedMmbMmEv+XOllqO68806pXbu2TJkyRWrUqCEej0eaNGkiFy5cKPfnQ0JCjI8HAAAAABtisQMAAAAwULNmzWTevHmSkJDgXQAp69ixY5KTkyNTpkyRm266SUREMjMzfZ0JAAAAALbGZawAAAAAAw0cOFCOHz8u999/v3z99deya9cuWbp0qfTt21eKi4ulWrVqEhUVJe+99558//33snz5chk6dKjZ2QAAAABgKyx2AAAAAAaqUaOGrF69WoqLiyUtLU2uu+46GTx4sERERIjT6RSn0ymzZ8+WrKwsadKkiQwZMkTGjh1rdjYAAAAA2IpDVdXsCAAAAAAAAAAAgN+LMzsAAAAAAAAAAICtsdgBAAAAAAAAAABsjcUOAAAAAAAAAABgayx2AAAAAAAAAAAAW2OxAwAAAAAAAAAA2BqLHQAAAAAAAAAAwNZY7AAAAAAAAAAAALbGYgcAAAAAAAAAALA1FjsAAAAAAAAAAICtsdgBAAAAAAAAAABsjcUOAAAAAAAAAABgayx2AAAAAAAAAAAAW/s/v/JOcZVOwnsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1600x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Area of interest:     study_areas_laura_COD.shp\n",
      "Start year:           1984\n",
      "End year:             2024\n",
      "\n",
      "Directories\n",
      "base directory:       \\\\ites-sr.ethz.ch\\Satellites\\Hiwi_23_TimckeFinn\\Landsat-composites\n",
      "shapefile directory:  \\\\ites-sr.ethz.ch\\Satellites\\Hiwi_23_TimckeFinn\\Landsat-composites\\shapefiles\n",
      "\n",
      "Outputs:\n",
      "geotiff:               \\\\ites-sr.ethz.ch\\Satellites\\Hiwi_23_TimckeFinn\\Landsat-composites\\outputs\\study_areas_laura_COD\\geotiff\n",
      "metadata:              \\\\ites-sr.ethz.ch\\Satellites\\Hiwi_23_TimckeFinn\\Landsat-composites\\outputs\\study_areas_laura_COD\\metadata\n",
      "HTML reports:          \\\\ites-sr.ethz.ch\\Satellites\\Hiwi_23_TimckeFinn\\Landsat-composites\\outputs\\study_areas_laura_COD\\html\n",
      "plots:                 \\\\ites-sr.ethz.ch\\Satellites\\Hiwi_23_TimckeFinn\\Landsat-composites\\outputs\\study_areas_laura_COD\\plots\n",
      "\n",
      "EOdal will extract data for 34 time intervals and aggragate the scenes to one composite\n",
      "\n"
     ]
    }
   ],
   "source": [
    "## prepare the EOdal Mapper \n",
    "metadata_filters = [Filter('eo:cloud_cover', '<=', max_cc)]\n",
    "\n",
    "for s in excluded_sensor:\n",
    "    metadata_filters = metadata_filters + [Filter('instruments', '!=', s)]\n",
    "\n",
    "# set directories and create new folders\n",
    "base_dir: Path = Path('./')\n",
    "shapefiles_dir: Path = Path(base_dir.joinpath(\"shapefiles\")) # replace \"shapefiles\" with the actual folder where the shapefiles are located\n",
    "\n",
    "aoi_folder = shapefile_name[:-4]\n",
    "\n",
    "# Check if the folder already exists\n",
    "if not os.path.exists(f'outputs/{aoi_folder}'):\n",
    "    \n",
    "    # If it doesn't exist, create the folder\n",
    "    os.makedirs(f'outputs/{aoi_folder}')\n",
    "    os.makedirs(f'outputs/{aoi_folder}/geotiff')\n",
    "    os.makedirs(f'outputs/{aoi_folder}/metadata')\n",
    "    os.makedirs(f'outputs/{aoi_folder}/html')\n",
    "    os.makedirs(f'outputs/{aoi_folder}/plots')\n",
    "    \n",
    "output_dir: Path = Path(base_dir.joinpath(f'outputs/{aoi_folder}')) # replace  with the actual folder where you want to save all the different outputs\n",
    "output_geotiff_dir: Path = Path(output_dir.joinpath(\"geotiff\")) # replace with the actual folder where you want to save the geotiff\n",
    "output_metadata_dir: Path = Path(output_dir.joinpath(\"metadata\")) # replace with the actual folder where you want to save the metadata\n",
    "output_html_dir: Path = Path(output_dir.joinpath(\"html\")) # replace with the actual folder where you want to save the HTML\n",
    "output_plots_dir: Path = Path(output_dir.joinpath(\"plots\")) # replace with the actual folder where you want to save the geotiff\n",
    "\n",
    "# setting up the EOdal Mapper\n",
    "Settings = get_settings()\n",
    "\n",
    "# Increase the maximum item number if the limit is reached: mapper.metadata.shape[0]\n",
    "Settings.MAX_ITEMS  = 3000\n",
    "\n",
    "# we use STAC, i.e., Microsoft Planetary Computer\n",
    "Settings.USE_STAC = True\n",
    "\n",
    "# -------------------------- Collection -------------------------------\n",
    "collection = 'landsat-c2-l2'\n",
    "\n",
    "# ---------------------- Spatial Feature  ------------------------------\n",
    "geom: Path = Path(shapefiles_dir.joinpath(shapefile_name))\n",
    "feature = Feature.from_geoseries(gpd.read_file(geom).dissolve().geometry)\n",
    "\n",
    "# ------------------------- Time Range ---------------------------------\n",
    "time_start = datetime(year_start, 1, 1)\n",
    "time_end = datetime(year_end, 12, 31)\n",
    "\n",
    "# set up the Mapper configuration\n",
    "mapper_configs = MapperConfigs(\n",
    "    metadata_filters=metadata_filters,\n",
    "    collection=collection,\n",
    "    feature=feature,\n",
    "    time_start=time_start,\n",
    "    time_end=time_end)\n",
    "\n",
    "# get a new mapper instance\n",
    "mapper = Mapper(mapper_configs)\n",
    "\n",
    "# query the collection for all scenes in the time range\n",
    "mapper.query_scenes()\n",
    "\n",
    "print(f'Number of Landsat scenes found for the entire time periode {year_start}-{year_end}: {mapper.metadata.shape[0]}')\n",
    "\n",
    "# the metadata is stored as a GeoDataFrame \"mapper.metadata\"\n",
    "# create copy to plot it\n",
    "metadata_gdf = mapper.metadata\n",
    "\n",
    "# Convert the 'timestamp' column to timezone-unaware datetimes\n",
    "metadata_gdf['sensing_time'] = metadata_gdf['sensing_time'].dt.tz_localize(None)\n",
    "\n",
    "# Save metadata as xlsx-file\n",
    "output_file_path = output_metadata_dir.joinpath(f'scenes_metadata_{shapefile_name[:-4]}_{year_start}-{year_end}.xlsx')\n",
    "metadata_gdf.to_excel(output_file_path, index=True)\n",
    "\n",
    "# Assuming metadata_gdf is defined somewhere in your code\n",
    "# Call the function with the appropriate parameters\n",
    "plot_geo_barplot(metadata_gdf, save_path=f'{output_plots_dir}\\scenes_{year_start}_{year_end}.png')\n",
    "\n",
    "# prints the directories and output locations\n",
    "print('Area of interest:    ', shapefile_name)\n",
    "print('Start year:          ', year_start)\n",
    "print('End year:            ', year_end)\n",
    "print()\n",
    "print('Directories')\n",
    "print('base directory:      ', base_dir.resolve())\n",
    "print('shapefile directory: ', shapefiles_dir.resolve())\n",
    "print()\n",
    "print('Outputs:' )\n",
    "print('geotiff:              ', output_geotiff_dir.resolve())\n",
    "print('metadata:             ', output_metadata_dir.resolve())\n",
    "print('HTML reports:         ', output_html_dir.resolve())\n",
    "print('plots:                ', output_plots_dir.resolve())\n",
    "print()\n",
    "\n",
    "# Plot the temporal distribution of all acquisitions and highlight those within the months of interest\n",
    "# Extracting year from the 'sensing_date' column\n",
    "metadata_gdf['year'] = pd.to_datetime(metadata_gdf['sensing_date']).dt.year\n",
    "\n",
    "# Get a list of all unique years in the dataset\n",
    "years_with_scenes = list(metadata_gdf['year'].unique())\n",
    "#print(f'years_with_scenes: {years_with_scenes}')\n",
    "\n",
    "# Get a list of years without scenes\n",
    "years_without_scenes = [year for year in range(min(years_with_scenes), max(years_with_scenes) + 1) if year not in years_with_scenes]\n",
    "#print(f'years_without_scenes: {years_without_scenes}')\n",
    "\n",
    "print(f'EOdal will extract data for {len(years_with_scenes)} time intervals and aggragate the scenes to one composite')\n",
    "print()\n",
    "\n",
    "#reset metadata_gdf\n",
    "metadata_gdf = gpd.GeoDataFrame()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5661562d",
   "metadata": {},
   "source": [
    "# Start Extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "39aed757",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time interval 1: Number of Landsat scenes found in year 1984: 19 \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\ftimcke\\AppData\\Local\\Temp\\ipykernel_34128\\3540126039.py:77: UserWarning: FixedFormatter should only be used together with FixedLocator\n",
      "  ax.set_xticklabels(months)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjUAAAGGCAYAAAAzegNcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABbiElEQVR4nO3deXxM9/7H8fdk3ySCkNBIbLEmFK1dLEGspWqrElsX2qKtLrRKlS5ulZZbtNXQXruiqpaiQqmtaq+qJVG9QjS1pdYk5/dHfplrJIiYmBxez8cjj4c56+dsM+O853y/FsMwDAEAAAAAAAAAAORzTo4uAAAAAAAAAAAAICcINQAAAAAAAAAAgCkQagAAAAAAAAAAAFMg1AAAAAAAAAAAAKZAqAEAAAAAAAAAAEyBUAMAAAAAAAAAAJgCoQYAAAAAAAAAADAFQg0AAAAAAAAAAGAKhBoAAAAAAAAAAMAUCDUAAABwx0JDQ9WrV69czZuQkCCLxaLp06dbh40cOVIWi8U+xd1Co0aN1KhRI+vruLg4WSwWLViw4K6sv1evXgoNDb0r68qtlJQU9evXT4GBgbJYLBo8ePAdL/P6/Z7deYCby811crv72WKxaOTIkbdf3F1w/TkEAACA+wOhBgAAwC1YLJYc/cXFxTm61Pva8ePHNXLkSO3cudPRpWSRn2vLiXfeeUfTp09X//799dVXX6lHjx4OrWfMmDFq166dihUrdtOb7gsXLlSXLl1UunRpeXl5qXz58nrppZd05syZu1pvfrds2bJ8G1wge3PnztUTTzyhcuXKyWKx3DTc2b59u6Kjo+Xr66sCBQqoefPm2b4Xpaena8qUKapWrZp8fHxUrFgxtWzZUj/99NNNaxkzZowsFouqVKlyh1sFAACQMy6OLgAAACC/++qrr2xef/nll1q1alWW4RUrVrybZd3T3njjDb322mu3Nc/x48f11ltvKTQ0VNWqVcvxfN9///1tVnf7blbbZ599pvT09Dyv4U788MMPql27tkaMGOHoUiRlnB+BgYF68MEHtXLlyhtO99RTT6l48eJ64oknVLJkSe3Zs0eTJk3SsmXL9Msvv8jT0/MuVm1/ublOQkJCdPHiRbm6ulqHLVu2TP/+97+zDTYuXrwoF5f8+d/Gu3Ht5leTJ0/W9u3b9dBDDyk5OfmG0/3yyy+qX7++goODNWLECKWnp+uTTz5RZGSktm7dqvLly1unffnll/Xhhx/qiSee0IABA3TmzBlNnTpVkZGR2rhxox5++OEsy//zzz/1zjvvyNvbO0+2EwAAIDv589spAABAPvLEE0/YvN68ebNWrVqVZfi9Ij09XVeuXJGHh4fDanBxccnzG6kXLlyQl5eX3Nzc8nQ9t3LtzeX8KikpSZUqVXJ0GVbx8fEKDQ3VX3/9pYCAgBtOt2DBgiy/YK9Ro4ZiYmI0c+ZM9evXL48rzVu5uU4sFsttXduOfB+4FUdfu3np0qVLcnNzk5NT9o0rfPXVVypRooScnJxu+oTE8OHD5enpqU2bNqlw4cKSMj7TwsLCNGzYMH399deSpNTUVE2ePFmPPfaYTWDfqVMnlS5dWjNnzsw21BgyZIhq166ttLQ0/fXXX3eyyQAAADlG81MAAAB2kJ6ergkTJqhy5cry8PBQsWLF9PTTT+v06dM204WGhqpNmzaKi4tTzZo15enpqfDwcGvTVQsXLlR4eLg8PDxUo0YN7dixw2b+Xr16ycfHR0eOHFGLFi3k7e2t4sWLa9SoUTIMw2baf/75Ry+99JKCg4Pl7u6u8uXL64MPPsgyncVi0XPPPaeZM2eqcuXKcnd314oVKyRJH3zwgerWravChQvL09NTNWrUuKO+Js6cOaNevXrJz89PBQsWVExMTLZNAWXXV8CqVatUv359FSxYUD4+PipfvryGDRsmKaMfjIceekiS1Lt3b2uTYJn9BjRq1EhVqlTR9u3b1bBhQ3l5eVnnvVG7/GlpaRo2bJgCAwPl7e2tdu3a6dixYzbT3KgvkWuXeavasutT43aP3eLFi1WlShW5u7urcuXK1uN3K0lJSerbt6+KFSsmDw8PVa1aVTNmzLCOz+xfJD4+Xt9995219oSEhBsuMzY2Vk2aNFHRokXl7u6uSpUqafLkyTmqJ6dy2gdJdse1Q4cOkqT9+/ffcv5vvvlGrVu3VvHixeXu7q4yZcro7bffVlpaWpZpt2zZolatWsnf31/e3t6KiIjQRx99ZDNN5nHy8PBQlSpVtGjRoizHP3OfX9+cXU77nrnZdZLdcnr16qV///vfkmyb2suUXfNeO3bsUMuWLeXr6ysfHx81bdpUmzdvtplm+vTpslgs2rhxo1588UUFBATI29tbHTp00KlTp2ym/fnnn9WiRQsVKVJEnp6eKlWqlPr06ZNlH1/vRv3hzJs3T2PGjNEDDzwgDw8PNW3aVIcOHbrpstauXSuLxaJFixZlGTdr1ixZLBZt2rTJOuy3337TY489pkKFCsnDw0M1a9bUkiVLbOb7+++/NWTIEIWHh8vHx0e+vr5q2bKldu3aZTNdZt1z5szRG2+8oRIlSsjLy0vnzp27Yb3BwcE3DDyu9eOPPyoqKsoaaEhSUFCQIiMjtXTpUqWkpEiSrl69qosXL6pYsWI28xctWlROTk7ZPtW0fv16LViwQBMmTLhlHQAAAPbEkxoAAAB28PTTT2v69Onq3bu3Bg4cqPj4eE2aNEk7duzQxo0bbX6Nf+jQIT3++ON6+umn9cQTT+iDDz5Q27ZtNWXKFA0bNkwDBgyQJL377rvq3LmzDhw4YHPzKi0tTdHR0apdu7bGjh2rFStWaMSIEUpNTdWoUaMkSYZhqF27dlq7dq369u2ratWqaeXKlXr55Zf13//+V+PHj7ep/4cfftC8efP03HPPqUiRItabrB999JHatWun7t2768qVK5ozZ446deqkpUuXqnXr1re1jwzD0COPPKINGzbomWeeUcWKFbVo0SLFxMTcct59+/apTZs2ioiI0KhRo+Tu7q5Dhw5p48aNkjKa/ho1apTefPNNPfXUU2rQoIEkqW7dutZlJCcnq2XLlurataueeOKJLDfvrpfZTvyrr76qpKQkTZgwQVFRUdq5c+dtNVuUk9qudbvHbsOGDVq4cKEGDBigAgUK6OOPP1bHjh31xx9/2NzIvN7FixfVqFEjHTp0SM8995xKlSql+fPnq1evXjpz5owGDRqkihUr6quvvtILL7ygBx54QC+99JIk3fTpiMmTJ6ty5cpq166dXFxc9O2332rAgAFKT0/Xs88+m+P9lldOnDghSSpSpMgtp50+fbp8fHz04osvysfHRz/88IPefPNNnTt3Tv/617+s061atUpt2rRRUFCQBg0apMDAQO3fv19Lly7VoEGDJGU0ldSxY0dVqlRJ7777rpKTk9W7d2898MADdtu2W10n2Xn66ad1/PjxbJvUu9E6GjRoIF9fX73yyitydXXV1KlT1ahRI61bt061atWymf7555+Xv7+/RowYoYSEBE2YMEHPPfec5s6dKykjWGvevLkCAgL02muvqWDBgkpISNDChQtzvR/ee+89OTk5aciQITp79qzGjh2r7t27a8uWLTecp1GjRgoODtbMmTOtwVemmTNnqkyZMqpTp451H9SrV08lSpTQa6+9Jm9vb82bN0/t27fX119/bZ3/yJEjWrx4sTp16qRSpUrp5MmT1uacfv31VxUvXtxmPW+//bbc3Nw0ZMgQXb582S5Poly+fDnb9ysvLy9duXJFe/fuVe3ateXp6alatWpp+vTpqlOnjho0aKAzZ87o7bfflr+/v5566imb+dPS0vT888+rX79+Cg8Pv+M6AQAAbosBAACA2/Lss88a136N+vHHHw1JxsyZM22mW7FiRZbhISEhhiTjp59+sg5buXKlIcnw9PQ0jh49ah0+depUQ5Kxdu1a67CYmBhDkvH8889bh6WnpxutW7c23NzcjFOnThmGYRiLFy82JBmjR4+2qemxxx4zLBaLcejQIeswSYaTk5Oxb9++LNt64cIFm9dXrlwxqlSpYjRp0sRmeEhIiBETE5Nl/mtl1jR27FjrsNTUVKNBgwaGJCM2NtY6fMSIETb7ePz48YYk6/ZlZ9u2bVmWkykyMtKQZEyZMiXbcZGRkdbXa9euNSQZJUqUMM6dO2cdPm/ePEOS8dFHH91yu69f5s1qi4mJMUJCQqyvb/fYubm52QzbtWuXIcmYOHFilnVda8KECYYk4z//+Y912JUrV4w6deoYPj4+NtseEhJitG7d+qbLy3T9OWMYhtGiRQujdOnSNsOu30fx8fE33Ec3curUKUOSMWLEiBzP07dvX8PZ2dn4/fffbzltdtvy9NNPG15eXsalS5cMw8g4h0uVKmWEhIQYp0+ftpk2PT3d+u9q1aoZQUFBxpkzZ6zDvv/+e0OSzfHPPP+uve4NI/v9k5vrJLvlXP+edq3r92/79u0NNzc34/Dhw9Zhx48fNwoUKGA0bNjQOiw2NtaQZERFRdnshxdeeMFwdna27odFixYZkoxt27bdsOYbudG1W7FiRePy5cvW4R999JEhydizZ89Nlzd06FDD3d3d5hglJSUZLi4uNvugadOmRnh4uPUcMIyMY123bl2jXLly1mGXLl0y0tLSbNYRHx9vuLu7G6NGjcpSd+nSpbM9526lcuXKNvvhWuHh4UZYWJiRmppqHXb58mWjZMmShiRjwYIF1uEHDx40qlevbkiy/pUuXdr47bffsix30qRJhp+fn5GUlGQYRsaxqFy58m3XDgAAkBs0PwUAAHCH5s+fLz8/PzVr1kx//fWX9a9GjRry8fHR2rVrbaavVKmS9Re/kqy/bG7SpIlKliyZZfiRI0eyrPO5556z/juzCaIrV65o9erVkjI6/nV2dtbAgQNt5nvppZdkGIaWL19uMzwyMjLbPhOu/YXv6dOndfbsWTVo0EC//PLLzXdKNpYtWyYXFxf179/fOszZ2VnPP//8LectWLCgpIzmgHLbqba7u7t69+6d4+l79uypAgUKWF8/9thjCgoK0rJly3K1/py63WMXFRWlMmXKWF9HRETI19c32/Pm+vUEBgaqW7du1mGurq4aOHCgUlJStG7dulzVf+05c/bsWf3111+KjIzUkSNHdPbs2Vwt015mzZqladOm6aWXXlK5cuVuOf2123L+/Hn99ddfatCggS5cuKDffvtNUkZTTPHx8Ro8eLD1PM2U2YxTYmKidu7cqZiYGPn5+VnHN2vWzK59ldjjOrmZtLQ0ff/992rfvr1Kly5tHR4UFKTHH39cGzZsyNJk0lNPPWXTnFWDBg2Ulpamo0eP2tS8dOlSXb161S519u7d2+Yph8yno251TfTs2VOXL1+2aWJv7ty5Sk1Ntfah9Pfff+uHH35Q586drefEX3/9peTkZLVo0UIHDx7Uf//7X0kZ7zmZT9mlpaUpOTnZ2iRYdu+hMTExdu+8fsCAAfr999/Vt29f/frrr9q7d6969uypxMRESRlPbGUqUKCAKleurGeffVYLFy7UJ598otTUVLVv396mv4zk5GS9+eabGj58+E2f2gIAAMgrhBoAAAB36ODBgzp79qyKFi2qgIAAm7+UlBQlJSXZTH9tcCHJepMzODg42+HX98vh5ORkc0NRksLCwiTJ2tfB0aNHVbx4cZub8lJGU0iZ469VqlSpbLdt6dKlql27tjw8PFSoUCEFBARo8uTJubo5ffToUQUFBcnHx8dmePny5W85b5cuXVSvXj3169dPxYoVU9euXTVv3rzbunFbokSJ22rO5fqb3haLRWXLlr1pfxL2cLvH7vrzSZL8/f2znDfZradcuXJZ2uW/0XpyauPGjYqKipK3t7cKFiyogIAAa58Ojgw1fvzxR/Xt21ctWrTQmDFjcjTPvn371KFDB/n5+cnX11cBAQHWm9uZ23L48GFJumlnzZn7MrsgJSfnf07Z4zq5mVOnTunChQvZ1lyxYkWlp6dn6Xfm+vPT399f0v/e1yIjI9WxY0e99dZbKlKkiB555BHFxsbq8uXLua7zVuu8kQoVKuihhx7SzJkzrcNmzpyp2rVrq2zZspIymg80DMN6Q//avxEjRkiS9T0/PT1d48ePV7ly5eTu7q4iRYooICBAu3fvzvZauNH78J145plnNGzYMM2aNUuVK1dWeHi4Dh8+rFdeeUWSrO/HqampioqKkp+fnyZNmqQOHTqof//+Wr16tQ4fPmzT3Nobb7yhQoUK5SiQBgAAyAv0qQEAAHCH0tPTVbRoUZsbYde6/peszs7O2U53o+HGdZ1D54Xsfh38448/ql27dmrYsKE++eQTBQUFydXVVbGxsZo1a1ae13R9fevXr9fatWv13XffacWKFZo7d66aNGmi77///ob77vpl2Nv1nTRnSktLy1FN9uDI8+Z6hw8fVtOmTVWhQgV9+OGHCg4Olpubm5YtW6bx48fnydMDObFr1y61a9dOVapU0YIFC+Ticuv/Bp05c0aRkZHy9fXVqFGjVKZMGXl4eOiXX37Rq6++mmfbcrNz6lbscZ3Y263OT4vFogULFmjz5s369ttvtXLlSvXp00fjxo3T5s2bs4Sg9ljnzfTs2VODBg3Sn3/+qcuXL2vz5s2aNGmSdXzmcR8yZIhatGiR7TIyA5B33nlHw4cPV58+ffT222+rUKFCcnJy0uDBg7M9f/LiPUrK6B9oyJAh2rdvn/z8/BQeHm4NGjMD8fXr12vv3r368MMPbeYtV66cKlasaO2X5eDBg/r00081YcIEHT9+3DrdpUuXdPXqVSUkJMjX11eFChXKk20BAACQCDUAAADuWJkyZbR69WrVq1cvz25KXSs9PV1Hjhyx3oySpN9//12SrB18h4SEaPXq1Tp//rzNL/4zm8wJCQm55Xq+/vpreXh4aOXKlXJ3d7cOj42NzVXdISEhWrNmjVJSUmxuVB44cCBH8zs5Oalp06Zq2rSpPvzwQ73zzjt6/fXXtXbtWkVFRd3wZnBuHTx40Oa1YRg6dOiQIiIirMP8/f115syZLPMePXrU5mma26nNHscup+vZvXu30tPTbZ7WuJP1fPvtt7p8+bKWLFli82v565tgu5sOHz6s6OhoFS1aVMuWLcvxTfK4uDglJydr4cKFatiwoXV4fHy8zXSZTX/t3btXUVFR2S4rc19ef05JWc//zKcKrj+vcvrkzK2uk+zk9PwMCAiQl5dXttfsb7/9JicnpyxPnOVU7dq1Vbt2bY0ZM0azZs1S9+7dNWfOHPXr1y9Xy8utrl276sUXX9Ts2bN18eJFubq6qkuXLtbxmde1q6vrDfdnpgULFqhx48aaNm2azfAzZ87kqKN6e/L391f9+vWtr1evXq0HHnhAFSpUkCSdPHlSUvbh2dWrV5WamipJ+u9//6v09HQNHDgwSxN5UsbTJoMGDdKECRPyYCsAAAAy0PwUAADAHercubPS0tL09ttvZxmXmpqa7U3vO3XtL4cNw9CkSZPk6uqqpk2bSpJatWqltLQ0m+kkafz48bJYLGrZsuUt1+Hs7CyLxWJzkyshIUGLFy/OVc2tWrVSamqqJk+ebB2WlpamiRMn3nLev//+O8uwatWqSZK1mRpvb29JWW8G59aXX36p8+fPW18vWLBAiYmJNvuuTJky2rx5s65cuWIdtnTp0ixN8NxObfY4djnRqlUrnThxQnPnzrUOS01N1cSJE+Xj46PIyMjbXmbmL+Sv/UX82bNncx2E3akTJ06oefPmcnJy0sqVK2+r/f/stuXKlSv65JNPbKarXr26SpUqpQkTJmQ5vpnzBgUFqVq1apoxY4ZNs0OrVq3Sr7/+ajNPSEiInJ2dtX79epvh1683Ozm5TrKT0/PT2dlZzZs31zfffGPTDNvJkyc1a9Ys1a9fX76+vres81qnT5/O8gRFTmrOK0WKFFHLli31n//8RzNnzlR0dLRNAFG0aFE1atRIU6dOtfZLca1Tp05Z/+3s7Jxl2+bPn2/tc8NR5s6dq23btmnw4MHWQDMzJJ8zZ47NtL/88osOHDigBx98UFJGM2uLFi3K8le5cmWVLFlSixYtUt++fe/uBgEAgPsOT2oAAADcocjISD399NN69913tXPnTjVv3lyurq46ePCg5s+fr48++kiPPfaY3dbn4eGhFStWKCYmRrVq1dLy5cv13XffadiwYdabtm3btlXjxo31+uuvKyEhQVWrVtX333+vb775RoMHD7bpWPpGWrdurQ8//FDR0dF6/PHHlZSUpH//+98qW7asdu/efdt1t23bVvXq1dNrr72mhIQEVapUSQsXLsxRPwujRo3S+vXr1bp1a4WEhCgpKUmffPKJHnjgAeuvj8uUKaOCBQtqypQpKlCggLy9vVWrVq1ct1NfqFAh1a9fX71799bJkyc1YcIElS1bVk8++aR1mn79+mnBggWKjo5W586ddfjwYf3nP//Jsn9vpzZ7HLuceOqppzR16lT16tVL27dvV2hoqBYsWKCNGzdqwoQJWfr0yInmzZvLzc1Nbdu21dNPP62UlBR99tlnKlq0aLY3gHPrq6++0tGjR3XhwgVJGU3njB49WpLUo0cP65MR0dHROnLkiF555RVt2LBBGzZssC6jWLFiatas2Q3XUbduXfn7+ysmJkYDBw6UxWLRV199leUmtZOTkyZPnqy2bduqWrVq6t27t4KCgvTbb79p3759WrlypSTp3XffVevWrVW/fn316dNHf//9tyZOnKjKlSsrJSXFujw/Pz916tRJEydOlMViUZkyZbR06dIsffNkJyfXSXZq1KghSRo4cKBatGghZ2dnde3aNdtpR48erVWrVql+/foaMGCAXFxcNHXqVF2+fFljx469ZY3XmzFjhj755BN16NBBZcqU0fnz5/XZZ5/J19dXrVq1uu3l2UPPnj2t79nZhdX//ve/Vb9+fYWHh+vJJ59U6dKldfLkSW3atEl//vmndu3aJUlq06aNRo0apd69e6tu3bras2ePZs6cmaVPpNxYv369Nfg6deqU/vnnH+s10LBhQ+vTRevXr9eoUaPUvHlzFS5cWJs3b1ZsbKyio6M1aNAg6/Jq1KihZs2aacaMGTp37pyaN2+uxMRETZw4UZ6enho8eLCkjNCnffv2WerJfDIju3EAAAB2ZwAAAOC2PPvss0Z2X6M+/fRTo0aNGoanp6dRoEABIzw83HjllVeM48ePW6cJCQkxWrdunWVeScazzz5rMyw+Pt6QZPzrX/+yDouJiTG8vb2Nw4cPG82bNze8vLyMYsWKGSNGjDDS0tJs5j9//rzxwgsvGMWLFzdcXV2NcuXKGf/617+M9PT0W64707Rp04xy5coZ7u7uRoUKFYzY2FhjxIgRWbY/JCTEiImJyX6HXSM5Odno0aOH4evra/j5+Rk9evQwduzYYUgyYmNjrdNdv441a9YYjzzyiFG8eHHDzc3NKF68uNGtWzfj999/t1n+N998Y1SqVMlwcXGxWWZkZKRRuXLlbGuKjIw0IiMjra/Xrl1rSDJmz55tDB061ChatKjh6elptG7d2jh69GiW+ceNG2eUKFHCcHd3N+rVq2f8/PPPWZZ5s9piYmKMkJAQm2nv9Njl9HicPHnS6N27t1GkSBHDzc3NCA8PtzkO1y4vu/M2O0uWLDEiIiIMDw8PIzQ01Hj//feNL774wpBkxMfHW6e7fh9lnu/Zrf96kZGRhqRs/9auXWud7kbTSMpyfLKzceNGo3bt2oanp6dRvHhx45VXXjFWrlyZZT2GYRgbNmwwmjVrZhQoUMDw9vY2IiIijIkTJ9pM8/XXXxsVK1Y03N3djUqVKhkLFy7M9vifOnXK6Nixo+Hl5WX4+/sbTz/9tLF37167XCfZ7efU1FTj+eefNwICAgyLxWKzTEnGiBEjbOr75ZdfjBYtWhg+Pj6Gl5eX0bhxY+Onn36ymSY2NtaQZGzbts1meOb1lbn/fvnlF6Nbt25GyZIlDXd3d6No0aJGmzZtjJ9//jm7Q2LjRtfu/Pnzbaa7nXPLMAzj8uXLhr+/v+Hn52dcvHgx22kOHz5s9OzZ0wgMDDRcXV2NEiVKGG3atDEWLFhgnebSpUvGSy+9ZAQFBRmenp5GvXr1jE2bNuW47pvJPPbZ/V17vA4dOmQ0b97cKFKkiPV9/N133zUuX76cZZkXLlwwRo0aZVSqVMnw9PQ0/Pz8jDZt2hg7duy4ZT03e48FAACwN4thOKAHQQAAAORKr169tGDBAptfdgMwr169eikuLs6mOSc4VmpqqooXL662bdtm6Q8DAAAAjkefGgAAAAAA/L/Fixfr1KlT6tmzp6NLAQAAQDboUwMAAAAAcN/bsmWLdu/erbffflsPPvigIiMjHV0SAAAAssGTGgAAAACA+97kyZPVv39/FS1aVF9++aWjywEAAMAN0KcGAAAAAAAAAAAwBZ7UAAAAAAAAAAAApkCoAQAAAAAAAAAATOG+6Cg8PT1dx48fV4ECBWSxWBxdDgAAAAAAAAAAuIZhGDp//ryKFy8uJ6cbP49xX4Qax48fV3BwsKPLAAAAAAAAAAAAN3Hs2DE98MADNxx/X4QaBQoUkJSxM3x9fR1cDQAAAAAAAAAAuNa5c+cUHBxsvZ9/I/dFqJHZ5JSvry+hBgAAAAAAAAAA+dStupCgo3AAAAAAAAAAAGAKhBoAAAAAAAAAAMAUCDUAAAAAAAAAAIAp3Bd9agAAAAAAAAAAHCctLU1Xr151dBlwIFdXVzk7O9/xcgg1AAAAAAAAAAB5wjAMnThxQmfOnHF0KcgHChYsqMDAwFt2Bn4zhBoAAAAAAAAAgDyRGWgULVpUXl5ed3QzG+ZlGIYuXLigpKQkSVJQUFCul0WoAQAAAAAAAACwu7S0NGugUbhwYUeXAwfz9PSUJCUlJalo0aK5boqKjsIBAAAAAAAAAHaX2YeGl5eXgytBfpF5LtxJ/yqEGgAAAAAAAACAPEOTU8hkj3OBUAMAAAAAAAAAAJgCoQYAAAAAAAAAAPlIXFycLBaLzpw5Y7dpQ0NDNWHCBLvU50iEGgAAAAAAAAAAXGfTpk1ydnZW69at7/q669atq8TERPn5+d32tNOnT1fBggWzTLdt2zY99dRT9i71riPUAAAAAAAAAADgOtOmTdPzzz+v9evX6/jx43d13W5ubgoMDMxRHxQ5nTYgIOCe6LSdUAMAAAAAAAAAkK/9/ru0fLl08ODdWV9KSormzp2r/v37q3Xr1po+fbrN+G+//VYPPfSQPDw8VKRIEXXo0ME6LikpSW3btpWnp6dKlSqlmTNn2jT9lJCQIIvFop07d1rnOXPmjCwWi+Li4iRlbVLq6NGjatu2rfz9/eXt7a3KlStr2bJlWaaNi4tT7969dfbsWVksFlksFo0cOVJS1uan/vjjDz3yyCPy8fGRr6+vOnfurJMnT1rHjxw5UtWqVdNXX32l0NBQ+fn5qWvXrjp//rx1mgULFig8PFyenp4qXLiwoqKi9M8//9zZzr8FQg0AAAAAAAAAQL70999SdLRUvrzUqpUUFpbx+vTpvF3vvHnzVKFCBZUvX15PPPGEvvjiCxmGIUn67rvv1KFDB7Vq1Uo7duzQmjVr9PDDD1vn7dWrl44dO6a1a9dqwYIF+uSTT5SUlHRH9Tz77LO6fPmy1q9frz179uj999+Xj49Plunq1q2rCRMmyNfXV4mJiUpMTNSQIUOyTJeenq5HHnlEf//9t9atW6dVq1bpyJEj6tKli810hw8f1uLFi7V06VItXbpU69at03vvvSdJSkxMVLdu3dSnTx/t379fcXFxevTRR637Ka+45OnSAQAAAAAAAADIpccfl1avth22erXUrZu0YkXerXfatGl64oknJEnR0dE6e/as1q1bp0aNGmnMmDHq2rWr3nrrLev0VatWlST9/vvvWr58ubZu3aqHHnrIuqyKFSveUT1//PGHOnbsqPDwcElS6dKls53Ozc1Nfn5+slgsCgwMvOHy1qxZoz179ig+Pl7BwcGSpC+//FKVK1fWtm3brLWnp6dr+vTpKlCggCSpR48eWrNmjcaMGaPExESlpqbq0UcfVUhIiCRZ68tLPKkBAAAAAAAAAMh3fv9dWrlSSkuzHZ6WljE8r5qiOnDggLZu3apu3bpJklxcXNSlSxdNmzZNkrRz5041bdo023n3798vFxcX1ahRwzqsQoUK2XbcfTsGDhyo0aNHq169ehoxYoR27959R8vbv3+/goODrYGGJFWqVEkFCxbU/v37rcNCQ0OtgYYkBQUFWZ86qVq1qpo2barw8HB16tRJn332mU7n9SM0ItQAAAAAAAAAAORDhw/ffPyhQ3mz3mnTpik1NVXFixeXi4uLXFxcNHnyZH399dc6e/asPD0972j5Tk4Zt+Wvbabp6tWrN52nX79+OnLkiHr06KE9e/aoZs2amjhx4h3VkROurq42ry0Wi9LT0yVJzs7OWrVqlZYvX65KlSpp4sSJKl++vOLj4/O0JkINAAAAAAAAAEC+U6bMzceXLWv/daampurLL7/UuHHjtHPnTuvfrl27VLx4cc2ePVsRERFas2ZNtvNXqFBBqamp2r59u3XYgQMHrB1+S1JAQICkjD4pMl3bafiNBAcH65lnntHChQv10ksv6bPPPst2Ojc3N6Vd/3jLdSpWrKhjx47p2LFj1mG//vqrzpw5o0qVKt2ylkwWi0X16tXTW2+9pR07dsjNzU2LFi3K8fy5QZ8aAAAAAAAAAIB8JyxMatEiow+Na+/ROztLUVFSuXL2X+fSpUt1+vRp9e3bV35+fjbjOnbsqGnTpulf//qXmjZtqjJlyqhr165KTU3VsmXL9Oqrr6p8+fKKjo7W008/rcmTJ8vFxUWDBw+2ebrD09NTtWvX1nvvvadSpUopKSlJb7zxxk3rGjx4sFq2bKmwsDCdPn1aa9euvWE/HaGhoUpJSdGaNWtUtWpVeXl5ycvLy2aaqKgohYeHq3v37powYYJSU1M1YMAARUZGqmbNmjnaV1u2bNGaNWvUvHlzFS1aVFu2bNGpU6fuuP+QW+FJDQAAAAAAAABAvjR7dkaAca2oqIzheWHatGmKiorKEmhIGaHGzz//rEKFCmn+/PlasmSJqlWrpiZNmmjr1q3W6WJjY1W8eHFFRkbq0Ucf1VNPPaWiRYvaLOuLL75QamqqatSoocGDB2v06NE3rSstLU3PPvusKlasqOjoaIWFhemTTz7Jdtq6devqmWeeUZcuXRQQEKCxY8dmmcZiseibb76Rv7+/GjZsqKioKJUuXVpz587NyW6SJPn6+mr9+vVq1aqVwsLC9MYbb2jcuHFq2bJljpeRGxbj2oa77lHnzp2Tn5+fzp49K19fX0eXAwAAAAAAAAD3vEuXLik+Pl6lSpWSh4fHHS3r4MGMPjTKls2bJzTyWmhoqAYPHqzBgwc7uhSHutk5kdP7+DQ/BQAAAAAAAADI18qVM2eYAfuj+SkAAAAAAAAAAGAKPKkBAAAAAAAAAEAeSkhIcHQJ9wye1AAAAAAAAAAAAKZAqAEAAAAAAAAAAEyBUAMAAAAAAAAAkGfS09MdXQLyCXucC/SpAQAAAAAAAACwOzc3Nzk5Oen48eMKCAiQm5ubLBaLo8uCAxiGoStXrujUqVNycnKSm5tbrpdFqAEAAAAAAAAAsDsnJyeVKlVKiYmJOn78uKPLQT7g5eWlkiVLyskp941IEWoAAAAAAAAAAPKEm5ubSpYsqdTUVKWlpTm6HDiQs7OzXFxc7vhpHUINAAAAAAAAAECesVgscnV1laurq6NLwT2AjsIBAAAAAAAAAIApEGoAAAAAAAAAAABTINQAAAAAAAAAAACmQKgBAAAAAAAAAABMgVADAAAAAAAAAACYAqEGAAAAAAAAAAAwBUINAAAAAAAAAABgCoQaAAAAAAAAAADAFAg1AAAAAAAAAACAKRBqAAAAAAAAAAAAUyDUAAAAAAAAAAAApkCoAQAAAAAAAAAATIFQAwAAAAAAAAAAmAKhBgAAAAAAAAAAMAVCDQAAAAAAAAAAYAqEGgAAAAAAAAAAwBQINQAAAAAAAAAAgCkQagAAAAAAAAAAAFMg1AAAAAAAAAAAAKZAqAEAAAAAAAAAAEyBUAMAAAAAAAAAAJgCoQYAAAAAAAAAADAFQg0AAAAAAAAAAGAKhBoAAAAAAAAAAMAUCDUAAAAAAAAAAIApEGoAAAAAAAAAAABTINQAAAAAAAAAAACmQKgBAAAAAAAAAABMgVADAAAAAAAAAACYAqEGAAAAAAAAAAAwBUINAAAAAAAAAABgCoQaAAAAAAAAAADAFAg1AAAAAAAAAACAKRBqAAAAAAAAAAAAUyDUAAAAAAAAAAAApkCoAQAAAAAAAAAATIFQAwAAAAAAAAAAmAKhBgAAAAAAAAAAMAVCDQAAAAAAAAAAYAqEGgAAAAAAAAAAwBQINQAAAAAAAAAAgCkQagAAAAAAAAAAAFMg1AAAAAAAAAAAAKZAqAEAAAAAAAAAAEyBUAMAAAAAAAAAAJgCoQYAAAAAAAAAADAFQg0AAAAAAAAAAGAKhBoAAAAAAAAAAMAUCDUAAAAAAAAAAIApEGoAAAAAAAAAAABTINQAAAAAAAAAAACmQKgBAAAAAAAAAABMgVADAAAAAAAAAACYAqEGAAAAAAAAAAAwBUINAAAAAAAAAABgCoQaAAAAAAAAAADAFAg1AAAAAAAAAACAKRBqAAAAAAAAAAAAUyDUAAAAAAAAAAAApkCoAQAAAAAAAAAATIFQAwAAAAAAAAAAmAKhBgAAAAAAAAAAMAVCDQAAAAAAAAAAYAqEGgAAAAAAAAAAwBQINQAAAAAAAAAAgCkQagAAAAAAAAAAAFMg1AAAAAAAAAAAAKZAqAEAAAAAAAAAAEyBUAMAAAAAAAAAAJgCoQYAAAAAAAAAADAFQg0AAAAAAAAAAGAKhBoAAAAAAAAAAMAUCDUAAAAAAAAAAIApEGoAAAAAAAAAAABTINQAAAAAAAAAAACmQKgBAAAAAAAAAABMgVADAAAAAAAAAACYAqEGAAAAAAAAAAAwBUINAAAAAAAAAABgCrkKNY4dO6Y///zT+nrr1q0aPHiwPv30U7sVBgAAAAAAAAAAcK1chRqPP/641q5dK0k6ceKEmjVrpq1bt+r111/XqFGj7FogAAAAAAAAAACAlMtQY+/evXr44YclSfPmzVOVKlX0008/aebMmZo+fbo96wMAAAAAAAAAAJCUy1Dj6tWrcnd3lyStXr1a7dq1kyRVqFBBiYmJ9qsOAAAAAAAAAADg/+Uq1KhcubKmTJmiH3/8UatWrVJ0dLQk6fjx4ypcuLBdCwQAAAAAAAAAAJByGWq8//77mjp1qho1aqRu3bqpatWqkqQlS5ZYm6UCAAAAAAAAAACwJ4thGEZuZkxLS9O5c+fk7+9vHZaQkCBvb28FBATYrUB7OHfunPz8/HT27Fn5+vo6uhwAAAAAAAAAAHCNnN7Hz9WTGk2aNNH58+dtAg1JKlSokLp06ZKbRQIAAAAAAAAAANxUrkKNuLg4XblyJcvwS5cu6ccff7zjogAAAAAAAAAAAK7ncjsT79692/rvX3/9VSdOnLC+TktL04oVK1SiRAn7VQcAAAAAAAAAAPD/bivUqFatmiwWiywWi5o0aZJlvKenpyZOnGi34gAAAAAAAAAAADLdVqgRHx8vwzBUunRpbd261aZDcDc3NxUtWlTOzs52LxIAAAAAAAAAAOC2Qo2QkBBJUnp6ep4UAwAAAAAAAAAAcCO3FWpc6+DBg1q7dq2SkpKyhBxvvvnmHRcGAAAAAAAAAABwrVyFGp999pn69++vIkWKKDAwUBaLxTrOYrEQagAAAAAAAAAAALvLVagxevRojRkzRq+++qq96wEAAAAAAAAAAMiWU25mOn36tDp16mTvWgAAAAAAAAAAAG4oV6FGp06d9P3339u7FgAAAAAAAAAAgBvKVfNTZcuW1fDhw7V582aFh4fL1dXVZvzAgQPtUhwAAAAAAAAAAEAmi2EYxu3OVKpUqRsv0GLRkSNH7qgoezt37pz8/Px09uxZ+fr6OrocAAAAAAAAAABwjZzex8/Vkxrx8fG5LgwAAAAAAAAAACA3ctWnRqYrV67owIEDSk1NtVc9AAAAAAAAAAAA2cpVqHHhwgX17dtXXl5eqly5sv744w9J0vPPP6/33nvPrgUCAAAAAAAAAABIuQw1hg4dql27dikuLk4eHh7W4VFRUZo7d67digMAAAAAAAAAAMiUqz41Fi9erLlz56p27dqyWCzW4ZUrV9bhw4ftVhwAAAAAAAAAAECmXD2pcerUKRUtWjTL8H/++ccm5AAAAAAAAAAAALCXXIUaNWvW1HfffWd9nRlkfP7556pTp459KgMAAAAAAAAAALhGrpqfeuedd9SyZUv9+uuvSk1N1UcffaRff/1VP/30k9atW2fvGgEAAAAAAAAAAHL3pEb9+vW1c+dOpaamKjw8XN9//72KFi2qTZs2qUaNGvauEQAAAAAAAAAAQBbDMAxHF5HXzp07Jz8/P509e1a+vr6OLgcAAAAAAAAAAFwjp/fxc9X8VKakpCQlJSUpPT3dZnhERMSdLBYAAAAAAAAAACCLXIUa27dvV0xMjPbv36/rH/SwWCxKS0uzS3EAAAAAAAAAAACZchVq9OnTR2FhYZo2bZqKFSsmi8Vi77oAAAAAAAAAAABs5CrUOHLkiL7++muVLVvW3vUAAAAAAAAAAABkyyk3MzVt2lS7du2645U3atRIgwcPvuPlAAAAAAAAAACAe1+untT4/PPPFRMTo71796pKlSpydXW1Gd+uXTu7FAcAAADzWrlS2rJFqlNHatbs7qzz99+lw4elsmWlcuXuzjqR/3FeAObCNfs/jvgsza/YF47FdZkhv56HjqgrP+4LR52n+XFf3OssxvU9fefAt99+qx49eujcuXNZF5jDjsJ79eqlGTNm2AyLjY1V7969tWLFCr322mv67bffVKdOHc2ZM0fbt2/Xiy++qP/+979q06aNPv/8c3l5eeWo3nPnzsnPz09nz56Vr69vzjYSAAAAuXL4sFSrlpSc/L9hhQtL27ZJpUrlzTr//lt6/PGM/1BkatFCmj1b8vfPm3Ui/+O8AMyFa/Z/HPFZml+xLxyL6zJDfj0PHVFXftwXjjpP8+O+MLuc3sfPVagRGhqqNm3aaPjw4SpWrFiuCjx79qxatmypKlWqaNSoUZKkffv2KSoqSrVr19YHH3wgLy8vde7cWSVKlJC7u7vee+89paSkqEOHDnr55Zf16quv5mhdhBoAAAB3T5Eitl/sMxUuLP31V96sMzpaWr1auva3Nc7OUlSUtGJF3qwT+R/nBWAuXLP/44jP0vyKfeFYXJcZ8ut56Ii68uO+cNR5mh/3hdnl9D5+rvrUSE5O1gsvvJDrQEOS/Pz85ObmJi8vLwUGBiowMFDOzs6SpNGjR6tevXp68MEH1bdvX61bt06TJ0/Wgw8+qAYNGuixxx7T2rVrb7jsy5cv69y5czZ/AAAAyHsrV2b/xV7KGL5qlf3X+fvvGeu9/mHhtLSM4QcP2n+dyP84LwBz4Zr9H0d8luZX7AvH4rrMkF/PQ0fUlR/3haPO0/y4L+4nuQo1Hn300ZuGCncqIiLC+u9ixYrJy8tLpUuXthmWlJR0w/nfffdd+fn5Wf+Cg4PzrFYAAAD8z5YtNx+/aZP913n48M3HHzpk/3Ui/+O8AMyFa/Z/HPFZml+xLxyL6zJDfj0PHVFXftwXjjpP8+O+uJ/kqqPwsLAwDR06VBs2bFB4eHiWjsIHDhx4R0VduzyLxZJl+RaLRenp6Tecf+jQoXrxxRetr8+dO0ewAQAAcBfUqnXz8XXq2H+dZcrcfHzZsvZfJ/I/zgvAXLhm/8cRn6X5FfvCsbguM+TX89ARdeXHfeGo8zQ/7ov7Sa5Cjc8//1w+Pj5at26d1q1bZzPOYrHkONRwc3PLUafit8vd3V3u7u52Xy4AAABurkWLjDZkb9S2bLNm9l9nWFjGem/Ujm65cvZfJ/I/zgvAXLhm/8cRn6X5FfvCsbguM+TX89ARdeXHfeGo8zQ/7ov7Sa6an4qPj7/h35EjR3K8nNDQUG3ZskUJCQn666+/bvr0BQAAAMxh27aML/LXKlw4Y3hemT074z8t14qKyhiO+xfnBWAuXLP/44jP0vyKfeFYXJcZ8ut56Ii68uO+cNR5mh/3xf0iV09q2MuQIUMUExOjSpUq6eLFi4qNjXVkOQAAALCDUqWkv/7K6Bxv06aMR6/z+pdK/v7SihUZHQEeOpTxmPn98utB3BjnBWAuXLP/44jP0vyKfeFYXJcZ8ut56Ii68uO+cNR5mh/3xf3CYhiGcbsz9enT56bjv/jii1wXlBfOnTsnPz8/nT17Vr6+vo4uBwAAAAAAAAAAXCOn9/Fz9aTG6dOnbV5fvXpVe/fu1ZkzZ9SkSZPcLBIAAAAAAAAAAOCmchVqLFq0KMuw9PR09e/fX2Vu1eU8AAAAAAAAAABALuSqo/BsF+TkpBdffFHjx4+31yIBAAAAAAAAAACs7BZqSNLhw4eVmppqz0UCAAAAAAAAAABIymXzUy+++KLNa8MwlJiYqO+++04xMTF2KQwAAAAAAAAAAOBauQo1duzYYfPayclJAQEBGjdunPr06WOXwgAAAAAAAAAAAK5126GGYRiaMWOGAgIC5OnpmRc1AQAAAAAAAAAAZHHbfWoYhqGyZcvqzz//zIt6AAAAAAAAAAAAsnXboYaTk5PKlSun5OTkvKgHAAAAAAAAAAAgW7cdakjSe++9p5dffll79+61dz0AAAAAAAAAAADZshiGYdzuTP7+/rpw4YJSU1Pl5uaWpW+Nv//+224F2sO5c+fk5+ens2fPytfX19HlAAAAAAAAAACAa+T0Pv5tdxQuSePHj5fFYsl1cQAAAAAAAAAAALcrV6FGr1697FwGAAAAAAAAAADAzeWqTw1nZ2clJSVlGZ6cnCxnZ+c7LgoAAAAAAAAAAOB6uQo1btQNx+XLl+Xm5nZHBQEAAAAAAAAAAGTntpqf+vjjjyVJFotFn3/+uXx8fKzj0tLStH79elWoUMG+FQIAAAAAAAAAAOg2Q43x48dLynhSY8qUKTZNTbm5uSk0NFRTpkyxb4UAAAAAAAAAAAC6zVAjPj5ektS4cWMtXLhQ/v7+eVIUAAAAAAAAAADA9XLVp8batWttAo20tDTt3LlTp0+ftlthAAAAAAAAAAAA18pVqDF48GBNmzZNUkag0bBhQ1WvXl3BwcGKi4uzZ30AAAAAAAAAAACSchlqzJ8/X1WrVpUkffvtt0pISNBvv/2mF154Qa+//rpdCwQAAAAAAAAAAJByGWokJycrMDBQkrRs2TJ16tRJYWFh6tOnj/bs2WPXAgEAAAAAAAAAAKRchhrFihXTr7/+qrS0NK1YsULNmjWTJF24cEHOzs52LRAAAAAAAAAAAECSXHIzU+/evdW5c2cFBQXJYrEoKipKkrRlyxZVqFDBrgUCAAAAAAAAAABIuQw1Ro4cqSpVqujYsWPq1KmT3N3dJUnOzs567bXX7FogAAAAAAAAAACAJFkMwzAcXUReO3funPz8/HT27Fn5+vo6uhwAAAAAAAAAAHCNnN7Hz9WTGpK0Zs0arVmzRklJSUpPT7cZ98UXX+R2sQAAAAAAAAAAANnKVajx1ltvadSoUapZs6a1Xw0AAAAAAAAAAIC8lKtQY8qUKZo+fbp69Ohh73oAAAAAAAAAAACy5ZSbma5cuaK6devauxYAAAAAAAAAAIAbylWo0a9fP82aNcvetQAAAAAAAAAAANxQrpqfunTpkj799FOtXr1aERERcnV1tRn/4Ycf2qU4AAAAAAAAAACATLkKNXbv3q1q1apJkvbu3WvPegAAAAAAAAAAALKVq1Bj7dq19q4DAAAAAAAAAADgpm4r1Hj00UdvOY3FYtHXX3+d64IAAAAAAAAAAACyc1uhhp+fX17VAQAAAAAAAAAAcFO3FWrExsbmVR0AAAAAAAAAAAA35eToAgAAAAAAAAAAAHKCUAMAAAAAAAAAAJgCoQYAAAAAAAAAADAFQg0AAAAAAAAAAGAKhBoAAAAAAAAAAMAUCDUAAAAAAAAAAIApEGoAAAAAAAAAAABTINQAAAAAAAAAAACmQKgBAAAAAAAAAABMgVADAAAAAAAAAACYAqEGAAAAAAAAAAAwBUINAAAAAAAAAABgCoQaAAAAAAAAAADAFAg1AAAAAAAAAACAKRBqAAAAAAAAAAAAUyDUAAAAAAAAAAAApkCoAQAAAAAAAAAATIFQAwAAAAAAAAAAmAKhBgAAAAAAAAAAMAVCDQAAAAAAAAAAYAqEGgAAAAAAAAAAwBQINQAAAAAAAAAAgCkQagAAAAAAAAAAAFMg1AAAAAAAAAAAAKZAqAEAAAAAAAAAAEyBUAMAAAAAAAAAAJgCoQYAAAAAAAAAADAFQg0AAAAAAAAAAGAKhBoAAAAAAAAAAMAUCDUAAAAAAAAAAIApEGoAAAAAAAAAAABTINQAAAAAAAAAAACmQKgBAAAAAAAAAABMgVADAAAAAAAAAACYAqEGAAAAAAAAAAAwBUINAAAAAAAAAABgCoQaAAAAAAAAAADAFAg1AAAAAAAAAACAKRBqAAAAAAAAAAAAUyDUAAAAAAAAAAAApkCoAQAAAAAAAAAATIFQAwAAAAAAAAAAmAKhBgAAAAAAAAAAMAVCDQAAAAAAAAAAYAqEGgAAAAAAAAAAwBQINQAAAAAAAAAAgCkQagAAAAAAAAAAAFMg1AAAAAAAAAAAAKZAqAEAAAAAAAAAAEyBUAMAAAAAAAAAAJgCoQYAAAAAAAAAADAFQg0AAAAAAAAAAGAKhBoAAAAAAAAAAMAUCDUAAAAAAAAAAIApEGoAAAAAAAAAAABTINQAAAAAAAAAAACmQKgBAAAAAAAAAABMgVADAAAAAAAAAACYAqEGAAAAAAAAAAAwBUINAAAAAAAAAABgCoQaAAAAAAAAAADAFAg1AAAAAAAAAACAKRBqAAAAAAAAAAAAUyDUAAAAAAAAAAAApkCoAQAAAAAAAAAATIFQAwAAAAAAAAAAmAKhBgAAAAAAAAAAMAVCDQAAAAAAAAAAYAqEGgAAAAAAAAAAwBQINQAAAAAAAAAAgCkQagAAAAAAAAAAAFMg1AAAAAAAAAAAAKZAqAEAAAAAAAAAAEyBUAMAAAAAAAAAAJgCoQYAAAAAAAAAADAFQg0AAAAAAAAAAGAKhBoAAAAAAAAAAMAUCDUAAAAAAAAAAIApEGoAAAAAAAAAAABTINQAAAAAAAAAAACmQKgBAAAAAAAAAABMgVADAAAAAAAAAACYAqEGAAAAAAAAAAAwBUINAAAAAAAAAABgCoQaAAAAAAAAAADAFAg1AAAAAAAAAACAKRBqAAAAAAAAAAAAUyDUAAAAAAAAAAAApkCoAQAAAAAAAAAATIFQAwAAAAAAAAAAmAKhBgAAAAAAAAAAMAVCDQAAAAAAAAAAYAqEGgAAAAAAAAAAwBQINQAAAAAAAAAAgCkQagAAAAAAAAAAAFMg1AAAAAAAAAAAAKZAqAEAAAAAAAAAAEyBUAMAAAAAAAAAAJgCoQYAAAAAAAAAADAFQg0AAAAAAAAAAGAKhBoAAAAAAAAAAMAUCDUAAAAAAAAAAIApEGoAAAAAAAAAAABTINQAAAAAAAAAAACmQKgBAAAAAAAAAABMgVADAAAAAAAAAACYAqEGAAAAAAAAAAAwBUINAAAAAAAAAABgCi6OLuBuMAxDknTu3DkHVwIAAAAAAAAAAK6Xef8+837+jdwXoUZycrIkKTg42MGVAAAAAAAAAACAGzl//rz8/PxuOP6+CDUKFSokSfrjjz9uujNwbzt37pyCg4N17Ngx+fr6Oroc2BnHFxLnwf2IYw7OgfsPxxycA5A4D+5HHHNwDtx/OOb3H8MwdP78eRUvXvym090XoYaTU0bXIX5+flwAkK+vL+fBPYzjC4nz4H7EMQfnwP2HYw7OAUicB/cjjjk4B+4/HPP7S04eSqCjcAAAAAAAAAAAYAqEGgAAAAAAAAAAwBTui1DD3d1dI0aMkLu7u6NLgQNxHtzbOL6QOA/uRxxzcA7cfzjm4ByAxHlwP+KYg3Pg/sMxx41YDMMwHF0EAAAAAAAAAADArdwXT2oAAAAAAAAAAADzI9QAAAAAAAAAAACmQKiB+5rFYtHixYsdXQYAAADuEN/rAAAAgPvDPRFq9OrVS+3bt3d0GXCQXr16yWKxZPk7dOiQo0uDHWQe32eeeSbLuGeffVYWi0W9evW6+4XBITZt2iRnZ2e1bt3a0aUgj3DN43p8z7t/cezvT3zW399OnTql/v37q2TJknJ3d1dgYKBatGihjRs3Oro05LFjx46pT58+Kl68uNzc3BQSEqJBgwYpOTk5R/PHxcXJYrHozJkzeVso7C7z+/97771nM3zx4sWyWCwOqgp56dr7eK6uripWrJiaNWumL774Qunp6Y4uDyZxT4QaQHR0tBITE23+SpUq5eiyYCfBwcGaM2eOLl68aB126dIlzZo1SyVLlryjZV+9evVOy8NdNG3aND3//PNav369jh8/fkfLSktL4wtTPpWX1zwAIH+z52c9zKdjx47asWOHZsyYod9//11LlixRo0aNcnxjG+Z05MgR1axZUwcPHtTs2bN16NAhTZkyRWvWrFGdOnX0999/O7pE5DEPDw+9//77On36tKNLwV2SeR8vISFBy5cvV+PGjTVo0CC1adNGqampji4PJnDPhRorVqxQ/fr1VbBgQRUuXFht2rTR4cOHreMTEhJksVi0cOFCNW7cWF5eXqpatao2bdrkwKpxpzJ/xXPtn7Ozs7755htVr15dHh4eKl26tN56660sb46JiYlq2bKlPD09Vbp0aS1YsMBBW4EbqV69uoKDg7Vw4ULrsIULF6pkyZJ68MEHrcNyev3PnTtXkZGR8vDw0MyZM+/qtiD3UlJSNHfuXPXv31+tW7fW9OnTreMyf5n13XffKSIiQh4eHqpdu7b27t1rnWb69OkqWLCglixZokqVKsnd3V1//PGHA7YEt2Kva75JkyZ67rnnbJZ96tQpubm5ac2aNXm/IbC70NBQTZgwwWZYtWrVNHLkSOtri8Wizz//XB06dJCXl5fKlSunJUuW3N1CYXc5OfYwv5t91md+jl8ru1/xjh49WkWLFlWBAgXUr18/vfbaa6pWrVreF487dubMGf344496//331bhxY4WEhOjhhx/W0KFD1a5dO+s0/fr1U0BAgHx9fdWkSRPt2rXLuoyRI0eqWrVqmjp1qoKDg+Xl5aXOnTvr7Nmzjtos5MCzzz4rNzc3ff/994qMjFTJkiXVsmVLrV69Wv/973/1+uuvS5IuX76sV199VcHBwXJ3d1fZsmU1bdo0JSQkqHHjxpIkf39/nuw1oaioKAUGBurdd9+94TRff/21KleuLHd3d4WGhmrcuHHWccOGDVOtWrWyzFO1alWNGjUqT2rGncm8j1eiRAlVr15dw4YN0zfffKPly5dbP/9v9Z4vSd9++60eeugheXh4qEiRIurQoYMDtgaOcM+FGv/8849efPFF/fzzz1qzZo2cnJzUoUOHLL/Gff311zVkyBDt3LlTYWFh6tatG0ngPebHH39Uz549NWjQIP3666+aOnWqpk+frjFjxthMN3z4cHXs2FG7du1S9+7d1bVrV+3fv99BVeNG+vTpo9jYWOvrL774Qr1797aZJqfX/2uvvaZBgwZp//79atGixV2pH3du3rx5qlChgsqXL68nnnhCX3zxhQzDsJnm5Zdf1rhx47Rt2zYFBASobdu2Nk/jXLhwQe+//74+//xz7du3T0WLFr3bm4Ecssc1369fP82aNUuXL1+2zvOf//xHJUqUUJMmTe7OhsAh3nrrLXXu3Fm7d+9Wq1at1L17d37lCZhATj7rb2bmzJkaM2aM3n//fW3fvl0lS5bU5MmT87Bi2JOPj498fHy0ePFim8/ua3Xq1ElJSUlavny5tm/frurVq6tp06Y27/GHDh3SvHnz9O2332rFihXasWOHBgwYcLc2A7fp77//1sqVKzVgwAB5enrajAsMDFT37t01d+5cGYahnj17avbs2fr444+1f/9+TZ06VT4+PgoODtbXX38tSTpw4IASExP10UcfOWJzkEvOzs565513NHHiRP35559Zxm/fvl2dO3dW165dtWfPHo0cOVLDhw+33vzu3r27tm7davMDp3379mn37t16/PHH79Zm4A41adJEVatWtf647Vbv+d999506dOigVq1aaceOHVqzZo0efvhhR24C7ibjHhATE2M88sgj2Y47deqUIcnYs2ePYRiGER8fb0gyPv/8c+s0+/btMyQZ+/fvvxvlws5iYmIMZ2dnw9vb2/r32GOPGU2bNjXeeecdm2m/+uorIygoyPpakvHMM8/YTFOrVi2jf//+d6V23Frm9Z2UlGS4u7sbCQkJRkJCguHh4WGcOnXKeOSRR4yYmJhs573R9T9hwoS7uAWwl7p161qP3dWrV40iRYoYa9euNQzDMNauXWtIMubMmWOdPjk52fD09DTmzp1rGIZhxMbGGpKMnTt33vXakXP2vOYvXrxo+Pv7W88BwzCMiIgIY+TIkXdjU2An137PCwkJMcaPH28zvmrVqsaIESOsryUZb7zxhvV1SkqKIclYvnz5XagW9pSbY79o0aK7Vh/s72af9bGxsYafn5/N9IsWLTKu/S9trVq1jGeffdZmmnr16hlVq1bNy7JhRwsWLDD8/f0NDw8Po27dusbQoUONXbt2GYZhGD/++KPh6+trXLp0yWaeMmXKGFOnTjUMwzBGjBhhODs7G3/++ad1/PLlyw0nJycjMTHx7m0Icmzz5s03ff/+8MMPDUnGli1bDEnGqlWrsp0u8/8Dp0+fzrtikSeu/byvXbu20adPH8MwbN/jH3/8caNZs2Y287388stGpUqVrK+rVq1qjBo1yvp66NChRq1atfK4euTGze7jdunSxahYsWKO3vPr1KljdO/ePa/LRT51zz2pcfDgQXXr1k2lS5eWr6+vQkNDJSlLEyMRERHWfwcFBUmSkpKS7lqdsK/GjRtr586d1r+PP/5Yu3bt0qhRo6y/+PHx8dGTTz6pxMREXbhwwTpvnTp1bJZVp04dntTIhwICAqzNEMTGxqp169YqUqSIzTQ5vf5r1qx5t8qGnRw4cEBbt25Vt27dJEkuLi7q0qWLpk2bZjPdtddzoUKFVL58eZvr2c3Nzeb9H/mXPa55Dw8P9ejRQ1988YUk6ZdfftHevXtpjuA+cO117u3tLV9fX77nAflcTj/rb7WM63+hyS82zaVjx446fvy4lixZoujoaMXFxal69eqaPn26du3apZSUFBUuXNjm/3jx8fE2v84uWbKkSpQoYX1dp04dpaen68CBA47YJOSQcYunshISEuTs7KzIyMi7VBEc4f3339eMGTOy3JPZv3+/6tWrZzOsXr16OnjwoNLS0iRlPK0xa9YsSRnn0+zZs9W9e/e7UzjsxjAMWSyWHL3n79y5U02bNnVwxXAUF0cXYG9t27ZVSEiIPvvsMxUvXlzp6emqUqWKrly5YjOdq6ur9d+Z7bDSYax5eXt7q2zZsjbDUlJS9NZbb+nRRx/NMr2Hh8fdKg121KdPH2v7+P/+97+zjM/p9e/t7X1X6oX9TJs2TampqSpevLh1mGEYcnd316RJk3K8HE9PzyxtbyP/ssc1369fP1WrVk1//vmnYmNj1aRJE4WEhNy1bYB9OTk5ZbnpcW0Tc5mu/Z4nZXzX43ueueX02MO8bvVZzzlw//Dw8FCzZs3UrFkzDR8+XP369dOIESM0YMAABQUFKS4uLss81/e3AvMoW7asLBaL9u/fn21b+Pv375e/v3+Wpqlwb2rYsKFatGihoUOH3vYPkbp166ZXX31Vv/zyiy5evKhjx46pS5cueVMo8sz+/ftVqlQppaSk3PI9n/eF+9s9FWokJyfrwIED+uyzz9SgQQNJ0oYNGxxcFRylevXqOnDgQJaw43qbN29Wz549bV5f2xEt8o/o6GhduXJFFoslS18YXP/3rtTUVH355ZcaN26cmjdvbjOuffv2mj17tipUqCAp4/otWbKkJOn06dP6/fffVbFixbteM+zDHtd8eHi4atasqc8++0yzZs26rRAM+U9AQIASExOtr8+dO6f4+HgHVoS7hWN/b8vJZ31ISIjOnz+vf/75x/oDlZ07d9pMW758eW3bts3mu/22bdvyvH7krUqVKmnx4sWqXr26Tpw4IRcXF+vTmdn5448/dPz4cWtAtnnzZjk5Oal8+fJ3qWLcjsKFC6tZs2b65JNP9MILL9jcpDxx4oRmzpypnj17Kjw8XOnp6Vq3bp2ioqKyLMfNzU2SrL/ah3m99957qlatms01W7FiRW3cuNFmuo0bNyosLEzOzs6SpAceeECRkZGaOXOmLl68qGbNmtGHosn88MMP2rNnj1544QU98MADt3zPj4iI0Jo1a7L0vYj7wz0Vavj7+6tw4cL69NNPFRQUpD/++EOvvfaao8uCg7z55ptq06aNSpYsqccee0xOTk7atWuX9u7dq9GjR1unmz9/vmrWrKn69etr5syZ2rp162095o67x9nZ2foYauYXl0xc//eupUuX6vTp0+rbt6/8/PxsxnXs2FHTpk3Tv/71L0nSqFGjVLhwYRUrVkyvv/66ihQpovbt2zugatiDva75fv366bnnnpO3t3e2vwCEeTRp0kTTp09X27ZtVbBgQb355ptZzg3cmzj297acfNavXLlSXl5eGjZsmAYOHKgtW7ZYO4jN9Pzzz+vJJ59UzZo1VbduXc2dO1e7d+9W6dKl7+LWILeSk5PVqVMn9enTRxERESpQoIB+/vlnjR07Vo888oiioqJUp04dtW/fXmPHjlVYWJiOHz9u7Sg2s4lZDw8PxcTE6IMPPtC5c+c0cOBAde7cWYGBgQ7eQtzIpEmTVLduXbVo0UKjR49WqVKltG/fPr388ssqUaKExowZo0KFCikmJkZ9+vTRxx9/rKpVq+ro0aNKSkpS586dFRISIovFoqVLl6pVq1by9PSUj4+PozcNuRAeHq7u3bvr448/tg576aWX9NBDD+ntt99Wly5dtGnTJk2aNEmffPKJzbzdu3fXiBEjdOXKFY0fP/5ul47bcPnyZZ04cUJpaWk6efKkVqxYoXfffVdt2rRRz5495eTkdMv3/BEjRqhp06YqU6aMunbtqtTUVC1btkyvvvqqozcPd8E90adGenq6XFxc5OTkpDlz5mj79u2qUqWKXnjhBeuNLtx/WrRooaVLl+r777/XQw89pNq1a2v8+PFZmh156623NGfOHEVEROjLL7/U7NmzValSJQdVjVvx9fWVr69vluFc//euadOmKSoqKstNDinjRsfPP/+s3bt3S8r4Vc+gQYNUo0YNnThxQt9++631V1swJ3tc8926dZOLi4u6detG84MmlPk9T5KGDh2qyMhItWnTRq1bt1b79u1VpkwZB1eIvMKxv3/k5LP+zz//1H/+8x8tW7ZM4eHhmj17tkaOHGkzbffu3TV06FANGTJE1atXV3x8vHr16sV7v0n4+PioVq1aGj9+vBo2bKgqVapo+PDhevLJJzVp0iRZLBYtW7ZMDRs2VO/evRUWFqauXbvq6NGjKlasmHU5ZcuW1aOPPqpWrVqpefPmioiIyHLjE/lLuXLl9PPPP6t06dLq3LmzypQpo6eeekqNGzfWpk2bVKhQIUnS5MmT9dhjj2nAgAGqUKGCnnzySf3zzz+SpBIlSuitt97Sa6+9pmLFilmbMIU5jRo1yqbp0OrVq2vevHmaM2eOqlSpojfffFOjRo3K0kTVY489puTkZF24cIEft+VzK1asUFBQkEJDQxUdHa21a9fq448/1jfffCNnZ+ccvec3atRI8+fP15IlS1StWjU1adJEW7dudfCW4W6xGLfqjckEoqOjVbZsWZqUAID7VFxcnBo3bqzTp0/TpjKySEhIUJkyZbRt2zZVr17d0eXgNvE97/7FsYc9NGvWTIGBgfrqq68cXQrugpEjR2rx4sVZmiYDAAD3FlM3P3X69Glt3LhRcXFxeuaZZxxdDgAAyEeuXr2q5ORkvfHGG6pduzaBhsnwPe/+xbFHbl24cEFTpkxRixYt5OzsrNmzZ2v16tVatWqVo0sDAACAHZk61OjTp4+2bduml156SY888oijywEAAPnIxo0b1bhxY4WFhWnBggWOLge3ie959y+OPXIrs6mKMWPG6NKlSypfvry+/vrrbDsVBgAAgHndE81PAQAAAAAAAACAe9890VE4AAAAAAAAAAC49xFqAAAAAAAAAAAAUyDUAAAAAAAAAAAApkCoAQAAAAAAAAAATIFQAwAAAAAAAAAAmAKhBgAAAIB8oVevXmrfvr2jywAAAACQjxFqAAAAANCpU6fUv39/lSxZUu7u7goMDFSLFi20cePGu1bDRx99pOnTp+f5enr16iWLxSKLxSJXV1cVK1ZMzZo10xdffKH09PTbWtb06dNVsGDBvCkUAAAAQBYuji4AAAAAgON17NhRV65c0YwZM1S6dGmdPHlSa9asUXJy8l2rwc/P766tKzo6WrGxsUpLS9PJkye1YsUKDRo0SAsWLNCSJUvk4sJ/lQAAAID8iCc1AAAAgPvcmTNn9OOPP+r9999X48aNFRISoocfflhDhw5Vu3btbKbr16+fAgIC5OvrqyZNmmjXrl3W8SNHjlS1atX01VdfKTQ0VH5+furatavOnz9vnWbBggUKDw+Xp6enChcurKioKP3zzz+SsjY/1ahRIw0cOFCvvPKKChUqpMDAQI0cOdKm9t9++03169eXh4eHKlWqpNWrV8tisWjx4sU33ebMp1FKlCih6tWra9iwYfrmm2+0fPlym6dFPvzwQ4WHh8vb21vBwcEaMGCAUlJSJElxcXHq3bu3zp49a33yI7O+y5cva8iQISpRooS8vb1Vq1YtxcXF5fygAAAAAMgWoQYAAABwn/Px8ZGPj48WL16sy5cv33C6Tp06KSkpScuXL9f27dtVvXp1NW3aVH///bd1msOHD2vx4sVaunSpli5dqnXr1um9996TJCUmJqpbt27q06eP9u/fr7i4OD366KMyDOOG65wxY4a8vb21ZcsWjR07VqNGjdKqVaskSWlpaWrfvr28vLy0ZcsWffrpp3r99ddzvR+aNGmiqlWrauHChdZhTk5O+vjjj7Vv3z7NmDFDP/zwg1555RVJUt26dTVhwgT5+voqMTFRiYmJGjJkiCTpueee06ZNmzRnzhzt3r1bnTp1UnR0tA4ePJjr+gAAAAAQagAAAAD3PRcXF02fPl0zZsxQwYIFVa9ePQ0bNky7d++2TrNhwwZt3bpV8+fPV82aNVWuXDl98MEHKliwoBYsWGCdLj09XdOnT1eVKlXUoEED9ejRQ2vWrJGUEWqkpqbq0UcfVWhoqMLDwzVgwAD5+PjcsLaIiAiNGDFC5cqVU8+ePVWzZk3r8latWqXDhw/ryy+/VNWqVVW/fn2NGTPmjvZFhQoVlJCQYH09ePBgNW7cWKGhoWrSpIlGjx6tefPmSZLc3Nzk5+cni8WiwMBABQYGysfHR3/88YdiY2M1f/58NWjQQGXKlNGQIUNUv359xcbG3lF9AAAAwP2OUAMAAACAOnbsqOPHj2vJkiWKjo5WXFycqlevbm2KadeuXUpJSVHhwoWtT3b4+PgoPj5ehw8fti4nNDRUBQoUsL4OCgpSUlKSJKlq1apq2rSpwsPD1alTJ3322Wc6ffr0TeuKiIiweX3t8g4cOKDg4GAFBgZaxz/88MN3tB8Mw5DFYrG+Xr16tZo2baoSJUqoQIEC6tGjh5KTk3XhwoUbLmPPnj1KS0tTWFiYzb5at26dzb4CAAAAcPvo/Q4AAACAJMnDw0PNmjVTs2bNNHz4cPXr108jRoxQr169lJKSoqCgoGz7hShYsKD1366urjbjLBaL0tPTJUnOzs5atWqVfvrpJ33//feaOHGiXn/9dW3ZskWlSpXKtqabLS8v7N+/31pLQkKC2rRpo/79+2vMmDEqVKiQNmzYoL59++rKlSvy8vLKdhkpKSlydnbW9u3b5ezsbDPuZk+lAAAAALg1Qg0AAAAA2apUqZK1w+3q1avrxIkTcnFxUWhoaK6XabFYVK9ePdWrV09vvvmmQkJCtGjRIr344ou3vazy5cvr2LFjOnnypIoVKyZJ2rZtW65r++GHH7Rnzx698MILkqTt27crPT1d48aNk5NTxkPumU1PZXJzc1NaWprNsAcffFBpaWlKSkpSgwYNcl0PAAAAgKwINQAAAID7XHJysjp16qQ+ffooIiJCBQoU0M8//6yxY8fqkUcekSRFRUWpTp06at++vcaOHauwsDAdP35c3333nTp06KCaNWvecj1btmzRmjVr1Lx5cxUtWlRbtmzRqVOnVLFixVzV3axZM5UpU0YxMTEaO3aszp8/rzfeeEOSbJqQys7ly5d14sQJpaWl6eTJk1qxYoXeffddtWnTRj179pQklS1bVlevXtXEiRPVtm1bbdy4UVOmTLFZTmhoqFJSUrRmzRpVrVpVXl5eCgsLU/fu3dWzZ0+NGzdODz74oE6dOqU1a9YoIiJCrVu3ztX2AgAAAKBPDQAAAOC+5+Pjo1q1amn8+PFq2LChqlSpouHDh+vJJ5/UpEmTJGWEBMuWLVPDhg3Vu3dvhYWFqWvXrjp69Kj1KYlb8fX11fr169WqVSuFhYXpjTfe0Lhx49SyZctc1e3s7KzFixcrJSVFDz30kPr166fXX39dUkZTWjezYsUKBQUFKTQ0VNHR0Vq7dq0+/vhjffPNN9Ymo6pWraoPP/xQ77//vqpUqaKZM2fq3XfftVlO3bp19cwzz6hLly4KCAjQ2LFjJUmxsbHq2bOnXnrpJZUvX17t27fXtm3bVLJkyVxtKwAAAIAMFsMwDEcXAQAAAAD2sHHjRtWvX1+HDh1SmTJlHF0OAAAAADsj1AAAAABgWosWLZKPj4/KlSunQ4cOadCgQfL399eGDRscXRoAAACAPECfGgAAAABM6/z583r11Vf1xx9/qEiRIoqKitK4ceMcXRYAAACAPMKTGgAAAAAAAAAAwBToKBwAAAAAAAAAAJgCoQYAAAAAAAAAADAFQg0AAAAAAAAAAGAKhBoAAAAAAAAAAMAUCDUAAAAAAAAAAIApEGoAAAAAAAAAAABTINQAAAAAAAAAAACmQKgBAAAAAAAAAABMgVADAAAAAAAAAACYwv8BOGvrqqpljEoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-02-26 10:13:19,997 eodal        INFO     Starting extraction of landsat scenes\n",
      "2024-02-26 10:17:44,264 eodal        INFO     Finished extraction of landsat scenes\n"
     ]
    }
   ],
   "source": [
    "# Extraction of scenes per time interval and Composite generation\n",
    "# make sure that there is at least one scene within time_start and time_end, otherwise the code will crash.\n",
    "\n",
    "for i, year in enumerate(years_with_scenes):\n",
    "    \n",
    "    time_start = datetime(year, 1, 1)\n",
    "    time_end = datetime(year, 12, 31)\n",
    "\n",
    "    # set up the Mapper configuration\n",
    "    mapper_configs = MapperConfigs(\n",
    "        metadata_filters=metadata_filters,\n",
    "        collection=collection,\n",
    "        feature=feature,\n",
    "        time_start=time_start,\n",
    "        time_end=time_end)\n",
    "\n",
    "    # get a new mapper instance\n",
    "    mapper = Mapper(mapper_configs)\n",
    "    mapper.query_scenes()\n",
    "    \n",
    "    print(f'Time interval {i+1}: Number of Landsat scenes found in year {year}: {mapper.metadata.shape[0]} ')\n",
    "    \n",
    "    # the metadata is stored as a GeoDataFrame \"mapper.metadata\"\n",
    "    # create copy to plot it\n",
    "    metadata_gdf = mapper.metadata\n",
    "\n",
    "    # Call the function with the appropriate parameters\n",
    "    scatterplot_temporal_distribution(metadata_gdf, save_path=f'{output_plots_dir}\\Aquisitions_{year}.png', year=year)\n",
    "    \n",
    "    # Convert the 'timestamp' column to timezone-unaware datetimes\n",
    "    metadata_gdf['sensing_time'] = metadata_gdf['sensing_time'].dt.tz_localize(None)\n",
    "\n",
    "    # Save metadata as an Excel file\n",
    "    metadata_gdf.to_excel(output_metadata_dir.joinpath(f'scenes_{year}_metadata.xlsx'), index=True)\n",
    "\n",
    "    # We tell EOdal how to load the Landsat scenes using `Landsat.from_usgs`and pass on some kwargs, e.g., the selection of bands we want to read. \n",
    "    # In addition, we tell EOdal to mask out clouds and shadows and the fly while reading the data using the qa_pixel band (therefore, we set the `read_qa` flag to True.\n",
    "    scene_kwargs = {\n",
    "        'scene_constructor': Landsat.from_usgs,\n",
    "        'scene_constructor_kwargs': {\n",
    "            'band_selection': ['blue', 'green', 'red', 'nir08', 'swir16', 'swir22'],\n",
    "            'read_qa': True},\n",
    "        'scene_modifier': prepocess_landsat_scene,\n",
    "        'scene_modifier_kwargs': {}\n",
    "        }\n",
    "\n",
    "    # now we load the scenes\n",
    "    mapper.load_scenes(scene_kwargs=scene_kwargs)\n",
    "\n",
    "    # the scenes are now loaded into a EOdal SceneCollection object named \"mapper.data\"\n",
    "\n",
    "    # all scenes have the same shape, i.e., the same number of bands, rows and columns\n",
    "    shapes = [{timestamp: scene.get_values().shape} for timestamp, scene in mapper.data]\n",
    "\n",
    "    #The mapper returns the single scenes. As we told the EOdal to mask out clouds, a significant share if the pixels is masked out. We will aggregate them in the next step.\n",
    "    #f = mapper.data.plot(band_selection=['red', 'green', 'blue'], figsize=(10, 20), max_scenes_in_row=3)\n",
    "\n",
    "    # open arrays for storing the data per band\n",
    "    shape = (len(mapper.data), list(shapes[0].values())[0][1], list(shapes[0].values())[0][2])\n",
    "    blue = np.ma.masked_array(data=np.ndarray(shape, dtype=float), mask=False)\n",
    "    red = np.zeros_like(blue)\n",
    "    green = np.zeros_like(blue)\n",
    "    nir08 = np.zeros_like(blue)\n",
    "    swir16 = np.zeros_like(blue)\n",
    "    swir22 = np.zeros_like(blue)\n",
    "\n",
    "    scene_count = np.zeros_like(blue)\n",
    "\n",
    "    # Next, we loop over the scenes.\n",
    "\n",
    "    # loop over scenes and store values\n",
    "    idx = 0\n",
    "    for _, scene in mapper.data:\n",
    "        blue[idx, :, :] = scene['blue'].values\n",
    "        red[idx, :, :] = scene['red'].values\n",
    "        green[idx, :, :] = scene['green'].values\n",
    "        nir08[idx, :, :] = scene['nir08'].values\n",
    "        swir16[idx, :, :] = scene['swir16'].values\n",
    "        swir22[idx, :, :] = scene['swir22'].values\n",
    "        idx += 1\n",
    "        \n",
    "    # Finally, we aggregate the data using the median reflectance. \n",
    "    # Since we deal with masked arrays, we have to call [numpy.ma.median](https://numpy.org/doc/stable/reference/generated/numpy.ma.median.html) instead of the \"normal\" numpy.median function.\n",
    "\n",
    "    # Calculate Spectral-temporal metrics:\n",
    "\n",
    "    # Calculate the median reflectance per spectral band\n",
    "    blue_median = np.ma.median(blue, axis=0)\n",
    "    green_median = np.ma.median(green, axis=0)\n",
    "    red_median = np.ma.median(red, axis=0)\n",
    "    nir08_median = np.ma.median(nir08, axis=0)\n",
    "    swir16_median = np.ma.median(swir16, axis=0)\n",
    "    swir22_median = np.ma.median(swir22, axis=0)\n",
    "\n",
    "    # Calculate the median absolute deviation (mad) per spectral band\n",
    "    # mad = np.median(abs(data - np.median(data))) (https://datagy.io/mean-absolute-deviation-python/)\n",
    "    blue_mad = ((np.ma.median(np.abs(blue - np.ma.median(blue, axis=0)), axis=0)) / np.ma.median(blue, axis=0)) * 100\n",
    "    green_mad = ((np.ma.median(np.abs(green - np.ma.median(green, axis=0)), axis=0)) / np.ma.median(green, axis=0)) * 100\n",
    "    red_mad = ((np.ma.median(np.abs(red - np.ma.median(red, axis=0)), axis=0)) / np.ma.median(red, axis=0)) * 100\n",
    "    nir08_mad = ((np.ma.median(np.abs(nir08 - np.ma.median(nir08, axis=0)), axis=0)) / np.ma.median(nir08, axis=0)) * 100\n",
    "    swir16_mad = ((np.ma.median(np.abs(swir16 - np.ma.median(swir16, axis=0)), axis=0)) / np.ma.median(swir16, axis=0)) * 100\n",
    "    swir22_mad = ((np.ma.median(np.abs(swir22 - np.ma.median(swir22, axis=0)), axis=0)) / np.ma.median(swir22, axis=0)) * 100\n",
    "\n",
    "    # Calculate the mean reflectance per spectral band\n",
    "    blue_mean = np.ma.mean(blue, axis=0)\n",
    "    green_mean = np.ma.mean(green, axis=0)\n",
    "    red_mean = np.ma.mean(red, axis=0)\n",
    "    nir08_mean = np.ma.mean(nir08, axis=0)\n",
    "    swir16_mean = np.ma.mean(swir16, axis=0)\n",
    "    swir22_mean = np.ma.mean(swir22, axis=0)\n",
    "\n",
    "    # Calculate the coefficient of variation (cov) per spectral band\n",
    "    # cov = np.std(data) / np.mean(data) (https://datagy.io/python-coefficient-of-variation/)  \n",
    "    blue_cov = (np.ma.std(blue, axis=0) / np.ma.mean(blue, axis=0)) * 100\n",
    "    green_cov = (np.ma.std(green, axis=0) / np.ma.mean(green, axis=0)) * 100\n",
    "    red_cov = (np.ma.std(red, axis=0) / np.ma.mean(red, axis=0)) * 100\n",
    "    nir08_cov = (np.ma.std(nir08, axis=0) / np.ma.mean(nir08, axis=0)) * 100\n",
    "    swir16_cov = (np.ma.std(swir16, axis=0) / np.ma.mean(swir16, axis=0)) * 100\n",
    "    swir22_cov = (np.ma.std(swir22, axis=0) / np.ma.mean(swir22, axis=0)) * 100\n",
    "\n",
    "    # Calculate the difference between mean and median, also known as skewness or skew per spectral band\n",
    "    # skew = np.mean(data) - np.median(data) https://www.turing.com/kb/calculating-skewness-and-kurtosis-in-python#skewness  \n",
    "    blue_skew = np.ma.mean(blue, axis=0) - np.ma.median(blue, axis=0)\n",
    "    green_skew = np.ma.mean(green, axis=0) - np.ma.median(green, axis=0)\n",
    "    red_skew = np.ma.mean(red, axis=0) - np.ma.median(red, axis=0)\n",
    "    nir08_skew = np.ma.mean(nir08, axis=0) - np.ma.median(nir08, axis=0)\n",
    "    swir16_skew = np.ma.mean(swir16, axis=0) - np.ma.median(swir16, axis=0)\n",
    "    swir22_skew = np.ma.mean(swir22, axis=0) - np.ma.median(swir22, axis=0)\n",
    "\n",
    "    # calculate the pixel count (e.g. for the blue band)\n",
    "    scene_count = np.ma.count(blue, axis=0)\n",
    "\n",
    "    # We store the results in a new RasterCollection.\n",
    "\n",
    "    ## Create new dictionaries for the metrics bands:\n",
    "    # median bands (changed the order of the first 3 bands from Landsats default order to match RGB when displayed in ArcGIS Pro)\n",
    "    median_bands = {'red': red_median, \n",
    "                    'green': green_median,\n",
    "                    'blue': blue_median, \n",
    "                    'nir08': nir08_median,\n",
    "                    'swir16' : swir16_median,\n",
    "                    'swir22' : swir22_median}\n",
    "\n",
    "    # mad bands\n",
    "    mad_bands = {'blue': blue_mad, \n",
    "                'green': green_mad, \n",
    "                'red': red_mad, \n",
    "                'nir08': nir08_mad,\n",
    "                'swir16' : swir16_mad,\n",
    "                'swir22' : swir22_mad}\n",
    "\n",
    "    # mean bands\n",
    "    mean_bands = {'blue': blue_mean, \n",
    "                    'green': green_mean, \n",
    "                    'red': red_mean, \n",
    "                    'nir08': nir08_mean,\n",
    "                    'swir16' : swir16_mean,\n",
    "                    'swir22' : swir22_mean}\n",
    "\n",
    "    # cov bands\n",
    "    cov_bands = {'blue': blue_cov, \n",
    "                    'green': green_cov, \n",
    "                    'red': red_cov, \n",
    "                    'nir08': nir08_cov,\n",
    "                    'swir16' : swir16_cov,\n",
    "                    'swir22' : swir22_cov}\n",
    "\n",
    "    # skew bands\n",
    "    skew_bands = {'blue': blue_skew, \n",
    "                    'green': green_skew, \n",
    "                    'red': red_skew, \n",
    "                    'nir08': nir08_skew,\n",
    "                    'swir16' : swir16_skew,\n",
    "                    'swir22' : swir22_skew}\n",
    "\n",
    "    ## Create new RasterCollection rc:\n",
    "    rc = RasterCollection()\n",
    "\n",
    "    ### add bands to RasterCollection rc as new bands:\n",
    "    # add median bands to rc\n",
    "    for band_name, band_value in median_bands.items():\n",
    "        rc.add_band(\n",
    "            band_constructor=Band,\n",
    "            values=band_value,\n",
    "            band_name=f'{band_name}_median',\n",
    "            geo_info=scene[band_name].geo_info)\n",
    "        \n",
    "    # add mad bands to rc\n",
    "    for band_name, band_value in mad_bands.items():\n",
    "        rc.add_band(\n",
    "            band_constructor=Band,\n",
    "            values=band_value,\n",
    "            band_name=f'{band_name}_mad',\n",
    "            geo_info=scene[band_name].geo_info)\n",
    "\n",
    "    # add mean bands to rc\n",
    "    for band_name, band_value in mean_bands.items():\n",
    "        rc.add_band(\n",
    "            band_constructor=Band,\n",
    "            values=band_value,\n",
    "            band_name=f'{band_name}_mean',\n",
    "            geo_info=scene[band_name].geo_info)\n",
    "        \n",
    "    # add cov bands to rc\n",
    "    for band_name, band_value in cov_bands.items():\n",
    "        rc.add_band(\n",
    "            band_constructor=Band,\n",
    "            values=band_value,\n",
    "            band_name=f'{band_name}_cov',\n",
    "            geo_info=scene[band_name].geo_info)\n",
    "        \n",
    "    # add skew bands to rc\n",
    "    for band_name, band_value in skew_bands.items():\n",
    "        rc.add_band(\n",
    "            band_constructor=Band,\n",
    "            values=band_value,\n",
    "            band_name=f'{band_name}_skew',\n",
    "            geo_info=scene[band_name].geo_info)\n",
    "        \n",
    "    # add pixel count to rc as new band\n",
    "    rc.add_band(\n",
    "            band_constructor=Band,\n",
    "            values=scene_count,\n",
    "            band_name='scene_count',\n",
    "            geo_info=scene['blue'].geo_info)\n",
    "\n",
    "    # add band aliases to the band names, which are used in the calculation of the spectral indices\n",
    "    rc._band_aliases = ['red', 'green', 'blue', 'nir_1', 'swir_1', 'swir_2',\n",
    "                        'blue_mad', 'green_mad', 'red_mad', 'nir08_mad', 'swir16_mad', 'swir22_mad', \n",
    "                        'blue_mean', 'green_mean', 'red_mean', 'nir08_mean', 'swir16_mean', 'swir22_mean', \n",
    "                        'blue_cov', 'green_cov', 'red_cov', 'nir08_cov', 'swir16_cov', 'swir22_cov', \n",
    "                        'blue_skew', 'green_skew', 'red_skew', 'nir08_skew', 'swir16_skew', 'swir22_skew', \n",
    "                        'scene_count'] + spectral_indices\n",
    "\n",
    "    # calculate some SI\n",
    "    for index in spectral_indices:\n",
    "        rc.calc_si(index, inplace = True)\n",
    "\n",
    "    # reset to original band aliases\n",
    "    rc._band_aliases = ['red_median', 'green_median', 'blue_median', 'nir08_median', 'swir16_median', 'swir22_median',\n",
    "                        'blue_mad', 'green_mad', 'red_mad', 'nir08_mad', 'swir16_mad', 'swir22_mad', \n",
    "                        'blue_mean', 'green_mean', 'red_mean', 'nir08_mean', 'swir16_mean', 'swir22_mean', \n",
    "                        'blue_cov', 'green_cov', 'red_cov', 'nir08_cov', 'swir16_cov', 'swir22_cov', \n",
    "                        'blue_skew', 'green_skew', 'red_skew', 'nir08_skew', 'swir16_skew', 'swir22_skew', \n",
    "                        'scene_count'] + spectral_indices\n",
    "\n",
    "\n",
    "    # save as GeoTiff \n",
    "    rc.to_rasterio(f'S:\\Hiwi_23_TimckeFinn\\Landsat-composites\\{output_geotiff_dir}\\Landsat_Composite_{year}_temp.tif')\n",
    "\n",
    "    ### clean the GeoTiff \n",
    "    # replace high fillvalues with NaN in the raster to make it readable by ArcGIS Pro\n",
    "\n",
    "\n",
    "    # Define input and output paths for TIFF files\n",
    "    input_tif_path = f'S:\\Hiwi_23_TimckeFinn\\Landsat-composites\\{output_geotiff_dir}\\Landsat_Composite_{year}_temp.tif'\n",
    "    output_tif_path = f'S:\\Hiwi_23_TimckeFinn\\Landsat-composites\\{output_geotiff_dir}\\Landsat_Composite_{year}_temp2.tif'\n",
    "\n",
    "    # Open the TIFF file using rasterio\n",
    "    with rasterio.open(input_tif_path) as src:\n",
    "        # Read all bands of the raster data\n",
    "        raster_data = src.read()\n",
    "\n",
    "        # Replace specific values (e.g., 1e+20) with NaN for each band\n",
    "        replaced_data = np.where(raster_data == 1e+20, np.nan, raster_data)\n",
    "\n",
    "        # Create a new raster file with NaN values\n",
    "        with rasterio.open(\n",
    "            output_tif_path,\n",
    "            'w',\n",
    "            driver='GTiff',\n",
    "            height=src.height,\n",
    "            width=src.width,\n",
    "            count=src.count,  # Number of bands\n",
    "            dtype=replaced_data.dtype,\n",
    "            crs=src.crs,\n",
    "            transform=src.transform,\n",
    "        ) as dst:\n",
    "            dst.write(replaced_data)\n",
    "\n",
    "    # Remove the original input TIFF file\n",
    "    os.remove(input_tif_path)\n",
    "\n",
    "    # Initialize a RasterCollection from the cleaned TIFF file\n",
    "    rc_clean = RasterCollection.from_multi_band_raster(fpath_raster=output_tif_path)\n",
    "\n",
    "    # Set band aliases for the RasterCollection\n",
    "    rc_clean._band_aliases = ['red_median', 'green_median', 'blue_median', 'nir08_median', 'swir16_median', 'swir22_median',\n",
    "                        'blue_mad', 'green_mad', 'red_mad', 'nir08_mad', 'swir16_mad', 'swir22_mad', \n",
    "                        'blue_mean', 'green_mean', 'red_mean', 'nir08_mean', 'swir16_mean', 'swir22_mean', \n",
    "                        'blue_cov', 'green_cov', 'red_cov', 'nir08_cov', 'swir16_cov', 'swir22_cov', \n",
    "                        'blue_skew', 'green_skew', 'red_skew', 'nir08_skew', 'swir16_skew', 'swir22_skew', \n",
    "                        'scene_count'] + spectral_indices\n",
    "\n",
    "    # Save the RasterCollection to a new GeoTIFF file\n",
    "    rc_clean.to_rasterio(f'S:\\Hiwi_23_TimckeFinn\\Landsat-composites\\{output_geotiff_dir}\\Landsat_Composite_{year}.tif',\n",
    "                        use_band_aliases=True)\n",
    "\n",
    "    # Remove the temporary cleaned TIFF file\n",
    "    os.remove(output_tif_path)\n",
    "\n",
    "    ### Generate plots of the metric bands and spectral indices ##\n",
    "    # We plot the bands. As you can see, after the compositing there are still some no-data pixels left (white areas). This is because for some locations there were not enough cloud-free observations to compute the statistical measure. This is most likely due to Landsat orbit design and overpass frequency.\n",
    "\n",
    "\n",
    "\n",
    "    # Define colormap dictionary\n",
    "    colormap_dict = {\n",
    "        'scene_count': \"Greys_r\",\n",
    "\n",
    "        'blue_median': \"Blues\",\n",
    "        'green_median': \"Greens\",\n",
    "        'red_median': \"Reds\",\n",
    "        'nir08_median': \"Greys\",\n",
    "        'swir16_median': \"Greys\",\n",
    "        'swir22_median': \"Greys\",\n",
    "\n",
    "        'blue_mad': \"Blues\",\n",
    "        'green_mad': \"Greens\",\n",
    "        'red_mad': \"Reds\",\n",
    "        'nir08_mad': \"Greys\",\n",
    "        'swir16_mad': \"Greys\",\n",
    "        'swir22_mad': \"Greys\",\n",
    "\n",
    "        'blue_mean': \"Blues\",\n",
    "        'green_mean': \"Greens\",\n",
    "        'red_mean': \"Reds\",\n",
    "        'nir08_mean': \"Greys\",\n",
    "        'swir16_mean': \"Greys\",\n",
    "        'swir22_mean': \"Greys\",\n",
    "\n",
    "        'blue_cov': \"Blues\",\n",
    "        'green_cov': \"Greens\",\n",
    "        'red_cov': \"Reds\",\n",
    "        'nir08_cov': \"Greys\",\n",
    "        'swir16_cov': \"Greys\",\n",
    "        'swir22_cov': \"Greys\",\n",
    "\n",
    "        'blue_skew': \"RdBu\",\n",
    "        'green_skew': \"RdBu\",\n",
    "        'red_skew': \"RdBu\",\n",
    "        'nir08_skew': \"RdBu\",\n",
    "        'swir16_skew': \"RdBu\",\n",
    "        'swir22_skew': \"RdBu\",\n",
    "\n",
    "        'NDVI' : 'YlGn',\n",
    "        'EVI' : 'YlGn',\n",
    "        'CI_GREEN' : 'YlGn', \n",
    "        'MSAVI' : 'YlGn',\n",
    "        'NDWI' : 'GnBu',\n",
    "        'BSI' : 'pink_r'\n",
    "    }\n",
    "\n",
    "    # Assuming rc, rc_clean, spectral_inidices and colormap_dict are defined somewhere in your code above:\n",
    "    # Call the function with the appropriate parameters, the plots will be saved in the plots folder\n",
    "    plot_bands_median_mad(rc, colormap_dict, save_path=f'{output_plots_dir}\\median_mad_{year}.png')\n",
    "    plot_bands_mean_cov(rc, colormap_dict, save_path=f'{output_plots_dir}\\mean_cov_{year}.png')\n",
    "    plot_bands_skew(rc, colormap_dict, save_path=f'{output_plots_dir}\\skew_{year}.png')\n",
    "    plot_rgb_and_spectral_indices(rc, rc_clean, spectral_indices, colormap_dict, save_path=f'{output_plots_dir}\\spectral_indices_{year}.png')\n",
    "    plot_rgb_and_scene_count(rc, colormap_dict, save_path=f'{output_plots_dir}\\scene_count_{year}.png')\n",
    "\n",
    "    # HTML Report generation for a single time interval\n",
    "    from jinja2 import Environment, FileSystemLoader\n",
    "\n",
    "    # 2. Create a template Environment\n",
    "    env = Environment(loader=FileSystemLoader('templates'))\n",
    "\n",
    "    # 3. Load the template from the Environment\n",
    "    template = env.get_template('report_template.html')\n",
    "\n",
    "    # 4. Render the template with variables\n",
    "    try:\n",
    "        html = template.render(page_title_text=f'Landsat Composite Report',\n",
    "                            title_text=f'Landsat Composite Report {year}',\n",
    "                            timeinterval=f'{str(time_start)[:-9]} - {str(time_end)[:-9]}',\n",
    "                            geom=geom,\n",
    "                            max_cc=max_cc,\n",
    "                            excluded_sensor=excluded_sensor,\n",
    "                            year=year,\n",
    "                            aoi_folder=aoi_folder\n",
    "                            )\n",
    "        \n",
    "    except Exception as e:\n",
    "        print(f\"Error rendering template: {e}\")\n",
    "        # Handle the exception or exit gracefully\n",
    "\n",
    "    # 5. Write the HTML file at the specified location\n",
    "    output_path = output_html_dir / f'Landsat_Composite_Report_{year}.html'\n",
    "    with open(output_path, 'w') as f:\n",
    "        f.write(html)\n"
   ]
  }
 ],
 "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
}