EX-Landsat-download-reanalysis-2020.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Folium version: 0.11.0\n"
     ]
    }
   ],
   "source": [
    "from src.tools_dev import Mapdisplay, random_point_in_shp\n",
    "import ee\n",
    "import tensorflow as tf\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Point, Polygon\n",
    "from folium.plugins import MarkerCluster\n",
    "import numpy as np\n",
    "from src.tools_dev import get_bound_points, generate_box_buffer, \\\n",
    "                            generate_ee_points, get_landsat_data\n",
    "import geopy.distance\n",
    "import rasterio\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "import folium\n",
    "print('Folium version: ' + folium.__version__)\n",
    "\n",
    "# Define the URL format used for Earth Engine generated map tiles.\n",
    "EE_TILES = 'https://earthengine.googleapis.com/map/{mapid}/{{z}}/{{x}}/{{y}}?token={token}'\n",
    "ee.Initialize()\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_features(shpfile):\n",
    "    df = gpd.read_file(shp).to_crs(\"EPSG:4326\")\n",
    "    df = df[df.area==np.max(df.area)].reset_index()\n",
    "    xmin, ymin, xmax, ymax = df.bounds.values[0]\n",
    "\n",
    "    x,y = df.geometry[0].exterior.coords.xy\n",
    "    cords = np.dstack((x,y)).tolist()\n",
    "    g = ee.Geometry.Polygon(cords)\n",
    "    feature = ee.Feature(g)\n",
    "    fc = ee.FeatureCollection([feature])\n",
    "    return(fc,g)\n",
    "\n",
    "def normalize(img):\n",
    "    return (img-np.nanmin(img))/(np.nanmax(img)-(np.nanmin(img)))\n",
    "\n",
    "def clipper(img):\n",
    "    return img.clip(fc)\n",
    "  \n",
    "def addTime(image): \n",
    "    constant = 1000 * 60 * 60 * 24 * 365\n",
    "    return image.addBands(image.metadata('system:time_start').divide(constant))\n",
    "\n",
    "def fill(image):\n",
    "    filled1a = image.focal_mean(2, 'square', 'pixels', 1).blend(image)\n",
    "    return filled1a.set('system:time_start',image.get('system:time_start'))\n",
    "\n",
    "def l8remap(image):\n",
    "    return image.rename(\\\n",
    "                        ['B0','B1','B2','B3','B4','B5','B6','B7','B8','B9','B10','B11'])\n",
    "\n",
    "def addNDVI(image): \n",
    "    ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI')\n",
    "    return image.addBands(ndvi)\n",
    "\n",
    "def addEVI(image):\n",
    "    evi = image.expression('2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))',\\\n",
    "                           {'NIR': image.select('B4'),\\\n",
    "                            'RED': image.select('B3'),\n",
    "                            'BLUE': image.select('B1')}).rename('EVI')\n",
    "    return image.addBands(evi)\n",
    "\n",
    "def srtm_image_generator(bounds):\n",
    "    elev = ee.Image(\"USGS/SRTMGL1_003\").clip(bounds)\n",
    "    slope = ee.Terrain.slope(elev)\n",
    "    aspect = ee.Terrain.aspect(elev)\n",
    "    out = elev.addBands(ee.Image(slope)).addBands(ee.Image(aspect)).float()\n",
    "    return out\n",
    "    \n",
    "def landsat_collection_generator(bounds, cloud_thr={'l5':20,'l7':20,'l8':10}):\n",
    "    L5coll = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')\\\n",
    "        .filter(ee.Filter.lt('CLOUD_COVER',cloud_thr['l5']))\\\n",
    "        .filterBounds(bounds)\n",
    "\n",
    "    L7coll = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')\\\n",
    "        .filter(ee.Filter.lt('CLOUD_COVER',cloud_thr['l7']))\\\n",
    "        .filterBounds(bounds)\\\n",
    "        .map(fill)\\\n",
    "        .filter(ee.Filter.calendarRange(2012,2013,'year'))\n",
    "    # add filter to only use 1 year of L7\n",
    "\n",
    "    L8coll = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')\\\n",
    "        .filter(ee.Filter.lt('CLOUD_COVER',cloud_thr['l8']))\\\n",
    "        .filterBounds(bounds)\\\n",
    "        .map(l8remap)\n",
    "\n",
    "    #Build landsat collection\n",
    "    landsat = ee.ImageCollection(L5coll.merge(L7coll.merge(L8coll)))\\\n",
    "        .map(addNDVI)\\\n",
    "        .map(addEVI)\n",
    "        \n",
    "    return landsat\n",
    " \n",
    "def image_collection_export(ic, outbands, BUCKET, PREFIX, DESC, REGION):\n",
    "    task = ee.batch.Export.image.toCloudStorage(\\\n",
    "                                         image = ic,\\\n",
    "                                         description = DESC,\\\n",
    "                                         bucket = BUCKET,\\\n",
    "                                         fileNamePrefix = PREFIX,\\\n",
    "                                         scale = 30,\\\n",
    "                                         region = REGION)\n",
    "    task.start()    \n",
    "    return \n",
    "\n",
    "def landsat_aggregate(landsat_ic, outbands, years, month_range):\n",
    "    landsat_images = []\n",
    "    for year in years:\n",
    "        renamed = [f'{b}_{year}' for b in outbands]\n",
    "        landsat_annual = landsat_ic.filter(ee.Filter.calendarRange(int(year),int(year),'year'))\\\n",
    "            .filter(ee.Filter.calendarRange(int(month_range[0]),int(month_range[1]),'month'))\\\n",
    "            .select(outbands)\\\n",
    "            .median()\\\n",
    "            .rename(renamed)\n",
    "        landsat_images.append(landsat_annual)\n",
    "\n",
    "    out = ee.Image.cat(landsat_images)    \n",
    "    return out\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:3: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.\n",
      "\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n",
      "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:3: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.\n",
      "\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbIAAAD4CAYAAACJx2OiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOxdZ5gUxdZ+a5clS5IoiIiAWVkXUVCyYECuSFAM1wCKmFjSrphAL/gJKMl0RVHBgImLiiSVJMGAu7IiqCiCJFFJooAsMHO+H6eL7unpnumJ3TNT7/PU093VVV2nZqr67ao6dY4gIigoKCgoKKQqstwWQEFBQUFBIRYoIlNQUFBQSGkoIlNQUFBQSGkoIlNQUFBQSGkoIlNQUFBQSGmUcVsAM2rWrEmNGjVyW4yMQnFx8S4iquW2HAoKCgrRwHNE1qhRIxQVFbktRkZBCLHZbRkUFBQUooWaWlRQUFBQSGkoIlNQUFBQSGkoIlNQUFBQSGkoIlNQUFBQSGkoIlNQUFBQSGkoIlNQUFBQSGkoIlNQUFBQSGlkNpFt2AAsXAgcOuS2JAoKCgoKUSIziWzOHKBdO+CMM4DOnYEKFYC9ewEi4L77gAUL3JZQQUFBQcEhPGfZI+H48kugb19g587A+FdeAbZuBSZN4uvLLku+bAoKCgoKESNzRmR//w0MHAi0asUkVqUKx3frxsf//IdJbNAgYMwY9+RUUFBQUIgImUFkH37I04hPPw1kZQGFhcCvvzJpffghp9m3D+jeHZgwARDCXXkVFBQUFBwjvYlsxw6gd2/gX/8Ctm0DWrQAioqAsWOBSpWYtIzYs8cdORUUFBQUokZ6EpnfD0yZApx+OjBzJpPWpEnAF18AzZtzGiJgyJDAfMuWAR99lHx5FRQUFBSiRvoR2XffAW3bAgMG8HRh164cl58PZGdzGklick3M7wcuvpjvXXMNsGWLe/IrKCgoKESE9CGyQ4eAkSN5xLVyJVCnDvDOO7wG1rBhYNrhw3USk2tiH30E1K3LSiHnnsvTkgoKCgoKnkf6qN/36gXMncvn/fuz5mH16tZpO3Tg45gxumJHxYrAunU8HfnHH8AllwBLlwK1lONkBQUFBS8jPUZk06cziVWsyOtcU6YwidltcL7sMlb4MGsn1qjBZHbGGTwd2aULb5RWUFBQUPAs0oPIpk3jY/36+lqXXAcbNw5YsoTjFixgYiMKzC8Jb/Ro4IkngE8+AZo2BUpKgMsvB/76S1n8UFBQUPAo0mNq8eBBPv70E5PXhAmByhxyg/OSJUxshw/ra2NGxY+WLYFVq/j+woWsNPLllzxC276dn6EsfigoKCh4CqlPZL/8wuRTsSJw661MSNLMlFGZA2BCO3xYv28mvM6dgREj9PuLFrHyyPbtwHHHAa1bJ716CgoKCgphQESeCnl5eRQRnniCCCC65hoiv5/PZfD7g9P7/USDBgWmGzSI4wsL+To3N/B+mTJ8POUUotLSyORLAQAoIg/89yqooIIK0YTUXyN7910+9uoVvMF5yJDg9TAhgi16yFHbmDE8Mlu9OvD+0aNA+fLAzz8DN9zA1woKCgoKnkBqE9nmzfq04qefBm5wHjSIr81kRhYWPWQaIYDx44Hc3MD7ubmsDVmlClsKef75xNdNQUFBQcERUpvIZs7k44knAs8+y+TVpQtveB4/Xiez4cOZqAoLgR49Agmve3e+HjyYr4cODR6RrV4NzJgBTJ3K11JhREFBQUHBdaS2soecVuzZk6f7xoxh0ho3DvjxR6BJE6CgAGjfXlfqAIC8PCY8IVjNHgAmTwY+/5xHeI0bAxs36uXk5nLeMmWA449nv2UzZgC33JLM2irEEUKI8gCWASgH7gcziWikEGIagHYA9mlJbyGiEov84wB0BX8MfgIgH0BlAMsNyRoAeJ2IBiWqHgoKCqlMZJs3s2p8xYrAgw/yEQjWTMzPZ/NTkyfz6GvLFqC4mK12XH45b4w+fJjvZ2UxyRUXc97mzVkF/6mnOP6nn4Ddu/Vy/v1v3X6jQqqhFEBHItovhMgBsEIIMV+7V0BEM+0yCiFaA7gIwDla1AoA7YhoKYDmhnTFAGaFE6RmzZrUqFGjqCqhoJBJKC4u3kVEQeaWUpfI5LRi1646iQG6MgcRk9PkyRyfn8/H998P3FsmBDBxIh8l+QFMYiUlQLt2nHfyZCa4e+8FZs8G1q/nZ/Xsmfi6KsQdREQA9muXOVog+xyB2QGUB1AWgNDy/m5MIIRoBqA2AkdolmjUqBGKioocFq2gkLkQQmy2vOG22qQ5OFa/v+ACIoDonXes75tV8c2q9uHS+3zBavrdu3O6Z5/l6/POs35WigEZqn4PIBtACZjQxmpx0wCsB7AGwEQA5WzyPgngT/AU5GMW90cAeDJE2f0BFAEoatiwYVT/m4JCpsHuXZWayh5btujTildcEXyfLDQTJaw8QFulHzqUFUaMmDWL8956K1C7NvD118D//V/09VBwFUTkI6Lm4LWslkKIswDcD+A0AOcDqAHgPnM+IUQTAKdr+eoD6CiEaGNK1gfAmyHKfoGIWhBRi1rKMLWCQkxwTGRCiGwhxGohxBztepoQYpMQokQLzW3yjRNCrBNCfC+EeEoIM4tEAeO0YqVKgffmzwfOP5+nCfPz9SlFiauvDlTHN6Y3q+7XqxeYd/Bgzlu+vK4kMmVKzNVRcBdE9CeAJQAuI6Id2sdfKYBXALS0yHI1gC+IaD8R7QcwH0AreVMIcS6AMkRUnATxFRQyHpGMyPIBfG+KKyCi5lqw0uwyLoqfBf7KbRetsMfwzjt87N07+N4jj/BaVvPm+jrZwIG6Z+gPPgjcWybT5+byCEwI4Mkn2X3LH3+wXzNJhpMnM5kNHsw+z8qWZQ3GFStirpJCciGEqCWEqKadVwDQGcAPQoh6WpwA0B3AWovsWwC0E0KU0RRF2iGwb1yHEKMxBQWF+MIRkQkhGoBVjadG+Hzjong5WCyKRww5rVihAk8rmi3aP/IIaxiWlLC24VVXMTmVlHD8VVfpe8uM6Vev5ulEIuCCC4CdO3n68OWXWRnESGaTJ/Oo7T5t1mnYMB7JKaQS6gFYIoRYA+ArAJ8Q0RwAbwghvgXwLYCaAEYDgBCihRBCtv+ZAH7W0nwD4Bsi+tDw7GugiCw0iot1Y98KCrHCauHMHMAdNw9AewBzKI6L4hTpwvf48axo0bs3X0v7iEYlDp+P7SJaKXlIm4rz5+vPtLK/mJfHzzGmMdtx/Osvorp1+fr118MuVHoVyFBlD6+EiO2LegWlpUTt2hHVqUPUrx/RnDlE//wTPt/bb3Of6dAh4SIqpBfs3lVOSOxKAM9p50YiqwdWPS4HYDqAERZ5mwCYC94oWhnA5wDahCovbKe+8EIK0Fb0+1mbECDKz7cmJUk88+cziZk1Df1+ooKC4PTG+3aGhl95ha/r1yfavz+Kv8Z9KCJTRBYVRo8O7meVK7MB7zffJNq3LzhPaan+8Se1gxUUHCIWInscwDYAvwD4DcBBsLUCY5pjBGeKLwDwsOF6BIDCUOWF7NRbtrDIFSoEkoaZhACi5s0piHhkOuPoze9nArRT0zeSmIwzXh89ymr4ANEjj8T8R7kBRWSKyCLG+vVE5crp/eDii4O9RpQtyzMjDz2ke4145hn9fnY20aZNrlZDIbUQNZEFJDaNyLSjADAJwBiL9NcCWAjeeJ0DYBGAbqHKCNmpJ0xgkXv1Coy3IyM5SpP3jOeSlIz5zCO6UOQn0xQWEi1bRscIdsuWmP4oN6CITBFZEDZuJNq71/qez8dTigDRLbcETu9v3Eg0cSJR27ZEQuh9q2ZNvl+7Nl/Le23bpsVeTIXkIBFEthi82L0WwOsAKmvxLQBM1c6zAUwBa3R9B2BCuDJCdmo5rfj224HxdkQmiclIPAUF1lOPMq18nkzTq5f9dKRxra13b05/ww0R/jXuQxGZIrIAvPoqUVYWr32tXh18/8UXua3Xrk20a1foWYu2bYnOOiuwr11wAdEff/A0JKBGZQqOERciS0aw7dR204rGTmMmMzM5SeIxK24UFIQnqnDYtEmfavn8c2d5PAJFZB5t827gxRcDR1JVqhAtXarf37GDqFo1vvfmm3p8qHVkv5/oq6+I7ryTlaiKizmPJLI9e5JbR4WUReoTmZxWPPXUQNKR0xp5eURXXRXYkSRJGRGqw8WKBx6gY1+cKTRdoojMo20+2TCuX40erc8ylCvH61x79uhxV1xh/fFnpzBlxpEjOmEeOMBx//zDH4GTJxONHKlGagpBSH0ia9WKLIln3jwmMfM9OTqT9hGJwituxEo+KaqOr4jMo20+mZAfigCvcRGxItM99+jxlSrpx82bA/NH84F49tmc7oQTuA+XKROYv2xZoqFDiXbvTly9FVIKqU1kxmnFu+8mSyIKtc5VWMhxVnvOrNLFgpdf5melkDq+IjIPtvlk4v/+T+9Dzz0XfH/lSqJLLtHTTJ4ceD/aD8Q5c4iaNtWfKwSvp/XtS9Snjx5frRrRuHHO9qgppDVSm8jk12LPntZfft27O1vnCrWPrLCQaNQoZ4odoeDzpZw6viIyD7b5ZMDv5yk8SSJTp4ZOv3w5K4KY937F8oHo9/N04tKlPKNhRHExUadOej9v2JBo+nS19yyDkdpEJqcVpbZiJHPxkSBeI7ZPP+X0NWrER64EQxGZB9t8ouH3Ew0fzu00K4votdeif1a4D0SnClN2ci5YQHTOOXp/79cv+ucppDRSl8jM2oqJVNaI1xraX3/RsbWEFIAiMo+1+UTD7ycaPJjbaHZ28HYWL+LoUaJp0/g9APAUvkLGIXWJbOJEspxWTISyBlGgySsrNWInX5g7d3K+44+PXZ4kQBGZx9p8IuHzEd11F7fPnByi995LXtnxgDQJJwTRv/5lbQZLIW1h964qA6/j3Xf52Ls3W6yXfsOkg8wJE/j+pEnsVmXs2NjKE0L3NSYhyxgyhMsBgMsus3/GoUN8LFcuNlkUFOKNJ58EnnuOz99/39oxrZdxyy3Ajz+yy6XZs9nJ7cyZwc5yFTIK3iayrVuBzz5jly1duwJVq3L8mDF6w5VkVrYs0KFD7GUSAYcPB8YNHszHyZOB7t31sZqx8xAx0XbooBNh+fKxy6OgEE9s3szHOnVSj8Qk/u//mMDOP5+9tj/5JFBQ4LZUCi4iEseaycfGjXysUweoWJFHQWPHBn99CcHxoUZJTkDEo67JkwO9S0sfZPn5TFJPPBHonFPmGzcOWLJEH5EpIlPwGmrU4OP11wffM/v2kyDi+AULrJ8ZTb5Y8zRtCrz6KscPHw707Mn3N25kUps1S/k7yyRYzTe6GQLWC0pLiU48kcc/iZ7Lnz+fqGVLCloTM1sJcbJOV1TE17m5iZU5ToBaI/NOm080OnbktjlrVvC9aLV2o8kXrzwjRuj98/rriRo31q+rViWaNEmp66cR7N5Vrndicwjq1JMns5gtWybW7JPR1JXPZ22IWFoJCac5uXIlx7VqlTh54whFZB5r84lErVrcNlu0IJo5k7UBJaJVpoomX7zy2Hm9MKrrX3op0a+/xuf3U3AVqUtkBw6wCwiAaPHiuP4oATAbHza7dzG7fwm1l23RIo5r3z5x8sYRisg81uYTiSlTWJtWtttmzYheeIHo0CG+H+32lmjyxSvPvfcSde6sX191Fad9/329rjVrEn3wQWS/lYLnkLpERkT0n/+wqJ07R1d7pxs2rTqJU1cwxg44dy7HXXZZdPImGYrIPNjmE4kDB4iefpropJP09luvHtHYsUR//hm9wYFo8sUrz759+vXo0XraX38NJLkBA3QjxQoph9Qmsj17dJcPRUWR1z6S+fhQHcvvZxKT+8zspkRmzaJjU5EpAEVkHmzzycDhw2zcWhrvBdhtS4sWgX3A6yMymef334meeILot98C8/h8bOaubFlOf+qpuisZhZRCahMZEdGwYSyu2Ts0UfgR17x5zubjnXQsJ6Q4YwafX3ut7R/iJSgi82ibTxb8fu4j0uszwBY/bruN6OabrfuBOb9ba2SRGEQoKSE64wxOn5PDhoiVIkhKIfWJbPt2/qISguiHHwLvGcll3jy+9vn0Rt6yJS9qG929SMWOefP4GcaFY6nUYdVJnExTSgv4N9/s5L9xHZlIZADKA1gF4BsA6wA8qsVPA7AJQIkWmtvkH6fl+x7AUwCEFl8WwAsAfgTwA4Ce4WRxncgkZD865RTdV5gQuoV6L2otRmoL9eDBQNc0HTsSbd0aPp+CJ5D6REZEdPvtLHLfvoHxxsYsySo31/raHKTjzYKC4LhoDAYTEf33v5znjjuc53ERGUpkAkBl7TwHwJcALtSIrFeYvK0BrASQrYXPAbTX7j0KYLR2ngWgZjhZPENkxo+0778nuuUWvU9ccYW9abZojAYnK48d5s4lql2b61a9OmtwKnge6UFkP/3ElrpzcoLVaa2mBXNzeSQmSUyqHhvDwIGB9hWNPs3kcyPtJNI+5MCBzvO4iEwkMmMAUBHA1wAucEhkrQAUA6ig5S0CcLp2byuASpGU7xkiM2PSJG7H9eoR7d3rtjTxx2+/MUHLd0HfvkR//+22VAohkB5ERqRv6JwzJ/ieWVHDGMwk1rx5cJp4GR4eM4YiHsW5iEwlMm00VQJgP4CxWtw0AOsBrAEwEUA5m7xPAvgTwD4Aj2lx1TQim6AR47sA6tjk768RYFHDhg2j/u8Shp9/JqpYkdtxOqut+/1EzzxDVL4817VJE6Ivv3RbKgUb2L2rvG2iygoNG/Jxx47AeCI2E2WHnTsDr4uLdRNUEl26sLkbftEEPjsSEz1Go8Gh8im4CiLyEVFzAA0AtBRCnAXgfgCnATgfQA0A95nzCSGaADhdy1cfQEchRBuw7dIGAD4jovPAU45P2pT9AhG1IKIWtWrVin/lYgERcPvtbOKpTx/gX/9yW6LEQQjg7ruBoiLgnHOADRuA1q2Bxx4DfD63pVNwiNQjsnr1+Ggksvnz2YCotIzv8wG5uYH5atcOvDbaSpR4+GG2lxjKjqIVliwJzCeJbOHC0PkUPAEi+hPAEgCXEdEO7eOvFMArAFpaZLkawBdEtJ+I9gOYD55u3A3gIIBZWrp3AZyX8ArEG1OnAosXAzVrAk895bY0ycGZZwKrVnEf9vmAhx4C2rfXjSwreBtWwzQ3Q9ipxaee4imAO+/U46SNRLkmJtfK5HSinCLJzWVtRrPVDuP6mlxPi8VEzwsv6M+P13RlAoEMnFoEUAtANe28AoDlAK4EUE+LEwAmARhjkfdaAAvBI7AcAIsAdNPuvQWgo3Z+C4B3w8niqTWyrVt5LxlA9OabbkvjDj7+mNcFAf4tZsxwWyIFDXbvKtdfKOYQtlO/+y4dU5GXmDdP106sU4cCCEleSw2lwsLwGopmNf1oN4QKoSulRKM0kiRkKJGdA2A1eC1sLYARWvxiAN9qca9D12xsAWCqdp4NYApY9f47ABMMzz0JwDLtuYsANAwni2eIbN483ejuv/4VqOIequ16XWsxmnxvvcXrZbIv33ijbvXEo/04E5A+RLZiBYt9wQWB8X6/NQHNnavHt2zJDXD+fCawgoLABmncQG18TrQmegA2ehytGn+SkIlE5qXgGSLr2pXbaNmyRNu2cZyTtuv1fWSxlNWpkz6j06gR0TXXeLYfZwLSh8h+/pnFttL08vkogESMm6KdTvHF02iqJNxIZUgyFJF5vM0nA7/9RlSjRnCb97r1+2TId9NNgR/JF15IdOSI899WIW5IHyI7eJCOfTWa93tZEUm0JBZLR/H5iM49NzoZXIAiMo+3+WSgd29up507B7tGSRVbi8mQT4aOHYl27LDPq5AQpA+REbHDPIBo1y6+tiIS88jMCeI5dWGWwaMkRmTfOFTwUJtPJKSR60qViDZtCp4ij3Zq3enHYzLyxKusjz7Slcjq1GGXTQpJg927KvXU74FAFfwFC4ALL9RV78ePB4YODUzfsiU3w3Do0AEoLAQmTOD9JQAfJ0zg+A4dnOUjCpZh8GBnMigoJBN79gB33cXnY8YAJ50UvB/TaquKGUSR50tWnniWNX8+sHo1q+b//jtwySXAI4+oPWduw4rdrAJYU2s1gDna9TSEMa4KoIPhfgmAQwC6hyrH0ddp+/b8RfTxx4GenY2q92atxWRN7RlHcAMHBloU8ej0ItSIzNXg6ohMWre/+OLA/qPWyELnOXKEaMQI3bhyhw5qqjEJsHtXOe5sAIYAmGEispA26Uz5awDYA6BiqHSOOvV117Ho06eHNhhs1lo0+x0rLCQaNSp+hkiJgqcZ+/fXXxShpiddhCKyDCWy+fO5TZYvT7R+vTtagV7XWgyX55NP9K09deoQLVxoLYNCXGD3rirjZNQmhGgAoCuAxzRCiwa9AMwnooNR5tchpxa3b9en/gCeXgR46D9oEMcPH87mqPLyeDoA0KcMJk3iacdVq4AffwSaNgXGjuU08j4RW+bo0AG47LLwssnpxzFjWLYqVfj6yivZ9I3d9KSCQjLx119A//58/uijQLNmwW0X0PtX2bKhp9YjzZesPMZ8t93GU/z16wM9esSnrEsuAUpKgBtu4PdE585sIWjECCA721oehfjDit3MAcBMAHkA2iNwRBbWuKrhGYsBXGlzLzIDqtOm8RdQly56nN1CbrgpAqOKPsAaW1Jry3ge7Ujq3ns5/6RJ0eVPAqBGZK4GV0Zkd97J7bJFi8xQJd+3j6hChcB3RHY2W/Do2JFo5Ej2Gh2to82jR/kZxqlGs4cOhZhh964K28nAZnue086NRFYPbManHIDp0Cwj2DyjHoCdAHLCleeoU+/axY2wTBl2LxFOtTbcfZ8veDO1mdDMm6edQk4t/ve/kedNEhSRZRiRLV3KbTInh2jNmuSWnWz4fGxiSloskdqZxr5uDFlZRDVrEp1wAmtHlytHVLky0erVzspbuFBfl1dTjXFHLET2OIBtAH4B8BvYKOrrpjTHCM7mGfkAXghXFkXSqTt0YPFfe83ZQq55xCa/vELtFYnHqEwupr/8cuR5kwRFZBlEZAcOsAdogOiRR5JXbrLxwgtE55/PH7vG/jxuHN8vLWX3NAMG8EdstWr27wCA75eWOit7xw79/SQEj9SOHk1YVTMJURNZQGLTiEw72hpXNeT7AkAHJ2U47tTSeHCzZhR2UdaKrKQBYRkvF2ytiCwWjcM+fTj/G29EnjdJUESWQUQ2dCi3x7PPdv5iTjW0bh3Yh7OyiM48k+jFF0Pn8/mIlizh2ZchQ9hB7pQpgdOFTmE11ai0GmNGIogsrHFV7boRgO0AspyU4bhTb97M4pcvTzR4sL3W4bx5wVqNZieb8jovj1XmrcgsWrX5q6/mZ3jYlboisjQkMitDuV98wS91IazXbJNplDdReW66KbDvjh8fO2GPHas/b+JE+3RW8i1cqH8kV6+uphpjRFyILBkhok4tiWn2bPs0RjVas2KHMTRvrpOYmcxiITLpSv3DD6PLnwQoIkuhNu8UZvXxQ4eIzjhDb9Opot4eSZ4LLwzu11OmRPa72UFun8nKIlq71jqNnXy3367LIwTRww9nhoJNApCeRDZ6NFfh1lv52uqLyGzp3ufT95nJULs2UY8edIzEzLbmAKKePZ3LJfHPP0QnncT5P/008vxJgiKyFGrzTmFeK37oIb0t33OP9zccR5OnZk0+P+20wD2l8cA//+im8Y4/3pqIQsk3cCATmJxqbNtW9zCg4BjpSWTr1gU2rHBfbAUF9iMy6XuoeXP9GT6fft2kSeSjssce47xnnunpLzBFZCnU5iOB1dpw797eNsobS57rr+frypWJpk6lYyOgw4dDy+4UX32lE5Fx608kdVq0iKhuXTr23pozJz6yZQjSk8j8fl3ZY9my8F9s0qGmeURmNmdlVgSR9yPRWty2TVfz9bhhUUVkKdTmI8Xhw4Ft3etGeWPJ8/vvOtEYtRClpmI88Oij+nOffTa6Ov3+OxOhvD90aPoq3sQZ6UlkRPo04IgRfB3qi8joSdpMdHl53EitnHP6fJGbq7rhBs7fo0dk9XEBishSrM07hd9PdNFF1n0hVJ5UHZER8QZv82xL8+ahZY8ULVvyc7OzeTtDNHXy+YjGjOFnALxV4Oef4ytnGiJ9iezDD7karVvrcXZfRE4WmKP9GjVCerEuV45o48bI8ycZishSrM07gd8fqMH30UepZZQ32jyffcYKGZUq6Xs4K1aM6acMwoED+sjvs89iq9Nnn7GTYICoShWid96Jr6xphvQlsr/+4q+a7Gw2QxPqiyic+q5RVT+Sr1Ejjh4lOu88zvvww5HVxSUoIkuxNu8Ew4bpbVgqQyVCkzDafInM888/fNy5U/8N4r39RRLZDz/EVicioj179G06AG/SPngwvvKmCdKXyIj06ZPZs6P7oiSK/mvUjBdf5DwNGhDt3x95XVyAIrIUbPPhcMcd3A7r1WMzbhLx3tsVbb5k5TnzTP4dypfnD914QZLOzp2xyWdM8/TT7Pke4A3r330XP3nTBOlNZCNG0LF1rmi+iIii/5oyYu9efXP1m29GXg+XoIgsBdt8KPz8M0+nAWyGKZOxYwfblJTvh3jgyBE6RmTx1kb++muipk352RUqED3/fHTLG2mK9CayZcu4KiedFP0XUSxfUxKDB7McbdqkVOPLRCIDUB7AKgDfAFgH4FEtfhrCOIzV0o3T8n0P4CkAQotfqnmFkPlrh5MlrkS2fDnRySdzO+zTJ37PTWXUr68TTzz2cxqnLBOBv/4KXN/s3p0NpSukOZGVluqq7m65Tti2jQ2UZmU5t5TtEWQokQmDWbUcAF8CuBAOHMYCaA1gJdhrejaAzwG0J53IWkQiS1yIrLSUaPhwfe0mNzdw2itTMWOGTgjHH8+q77Hihx/4eULE/qxQePNNVgAB2Bq/Mm9l+67KQjqgbFmgXTs+X7Qo+P6CBcB993FzNoKI4xcsiF2GvXuBo0fZaV/z5rE/TyGh0PrFfu0yRwsUIktAdvCIrizYjVEOgN/jLqRTrFsHXHCB7gTygQeAL74AatZ0TSTP4IEH+Ni8ObBrF1C7duzP3IGDXhsAACAASURBVLOHj1kJfn326QN88w1w0UXAr7+yE8+hQ4HS0sSWm4JIDyID+E8GgIULg+8tWQKMG8denyWZEfH1uHF8P1acfjp7g966lRudguchhMgWQpQA+APAJ0T0pXbrMSHEGiHERCFEOXM+IvocwBIAO7TwERF9b0jyihCiRAjxsBDSvXBQ2f2FEEVCiKKdO3dGVwG/H5g4kb2fl5QAjRsDy5YBjz3GH3eZDr8f2LKFzx99NH7PvPdePs/Jic8zQ6FRI2DpUuA//2GP0xMmsFf7desSX3YqwWqY5maIepplzRo6pi1otc4VD43EcJC79VNsLwgycGrRGABUAxPTWXDgMBZAEwBzAVTWwucA2mj36mvH4wB8DOCmcOWHbPN2a7e//KLvPwKIbruN11bC5cskrcWPP6ZjU4DRen42orRUtyQE8FRurHWKJN8XX+i+5HJyiDp1CvZzFsmafgrC7l3l+kvEHKImMr9fd5dg3NthvB/rHrFwkOZr8vPj98wkINOJjH8CjAAwzBTXHhYOYwEUAHjYlLfQIt0tAJ4JV3bINm/WpvX7iV59VVfTrljR2vtDsvZ2hconLWDcdVdg+ief5LVkuWcqUfINH87x1atbyx0J9u0LVBoZNSp2+aLJ9/ffRP366XKcdBLR9u3Oy0pxpD+RERFddx1X6ZlnrO/Hw2pHKCxcSHFV800SMpHIANQCUE07rwBgOYAr4cBhLIBrASwEUAa8PrYIQDftuqaWJgfATAADwskSss0bX0533MFGf2X7bdyY6LffwudLpLUNu3z33KPL2a1bYPoGDQL7obmsW2+Nj3z336+XsXu3/W8cDjt2ENWooT8rlI3FZFlGmTmT98bJPXKzZiVmlsljyAwie+klOqauaobfz/FWHSjW4bicGjBaGfn775QZ5mcokZ0DYDWANWDnsCO0+LAOY8GailPAqvffAZigxVcCUKw9cx2AyQCyw8kSts1btd3OncNPl0UzCxHtzIVVPqOsRvz73/o96ebIGLKz+cUcq3ylpXqaVatCy2+HDRvYmr6conzrrdDp4/n7hcu3bVvg9LIc/aYpiRHZv6tcf6GYQ0xE9umnXKVTTgmMNzeS/Hzd2LDxPNrhuHFqQG7KXrgwZYb5mUhkXgph2/yoURT0snf6sopmFiLamQtzvvvuIzr9dLK0NzppEtHll/P+KCvyK1OGaMgQ3kJgtEwSqXxyO0KnTpFvR1i9mu2lSnJdsMBZvnj9fk7yHT0amKdJE6Ivv3RWXgoiM4jM6BzTCEk0dgRm/vqJdNHWSJTSf1mrVsHP9SgUkXm8zZ9/PgW96L1kkT6eZdmF55+PrhzjiC8ri71Iz5gRfjS7ZAmTKcDrkU5HdF74/bKz+ePHrAiSBkh/Ilu9mqtTvnzwpmgjMVn9+d27BzaYaBZtk6FMkiAoIvNwm/f7iTp04PZUr15i112SucZjleb22ynopSxDp07RlePzsSd5abLLOOJr1cqa1GbN0pVRKla0Vh6L1+8Qz9/PuC4JsA3aTZucyZ4iSH8iu+oqrs7gweHThhvCx9Igjc/1sFdoIxSRebjNy48qOcX12Wfx14SLJU+8y+rbN7APyVCrFk81RiMfEZPVxIlsRFhON8pQrZpOZlOn6verViXautX+mV78/WQeaemoShWiN95wXgePI72JrLiYq1Khgr0Wl4TTkVOkIyyr9Ndeq0ZkKsTW5uVsglQl79WL48MpEqXyPjLZf049lahtW/06FluqRpSWEo0fr1vGB3gd7rbb9Os6dSLXdPTK71dYyEop8uMeILrxRqI//4ysPh5EehNZt25clWHDQqeLdKTldPHV/Jz+/ckR+XkEishSoM1v386bYLOyrJUn0gmy79Sty5uA5XWVKvEt5/ffA/u3DCefHOz5ORXh9xO98II+rXrSSWxUOoWRvkS2ahUdm8sOZxA0kiF8JCMy83P/9z++lhsoldaiCvFo81JtfdAgZ+lTFVLJon59vp41i6/Llo1vOXJd3TzNmCJLAo7xww+6NnVWFtEDDxAdPuy2VFEhfYnsiivIMVk4HcJHOnIzP3fvXm4wOTmcXu0jUyEebf7rr7kNVq6cFtNElhg5ko6RSo0aHCcJJ97W5ufOpQASq1o1PqasvIjSUt4gLtf/WrQgWr/ebakiRnoSmZx2qFw5sj0i4QitVy9+bvfuRAUFwdqOBQXh5+QvuIDTzpnjXC6XoIgshdq81GB88knneVIJctsKwHuiiNi4gIyLJ9FMnao/Nycnvh6kvYply/QtCRUrEk2Z4vmlDyPs3lWpbf1+5Eg+3ntvZC4rwlnDB4DCQqBpU+CJJzgeYMvTBQXATz9xupdesncPI108fPxx9PVTUDC7IJJt8amnuI3auSCKxnVRtO6O4llW1ar6+W+/8bFyZT3u7rvjJ9/vBs87DzzA3isSUadk/n7hyjpwAOjeHbjhBuDgQeCOO/j6jz/i59LKDVixm5vB8dfpypX8VXHccZF7T3U6dRguXUFB6DwAWzbwOKBGZN5t8+b1V5+P6MQT9faVbur35r1k0qqHXDeLp3xyWaJSpfT5/SIp6403eDpVjs5SeD3f9U5sDo6J7JJLWPyHHgqf1jyVOH8+k5DRsgfA1wUFgVOGoZQ+QhHdvffqNtoi2YviAhSRebjNW7Wxs87i6zZt0mtDtPFaGsS97jq2kC/7XufO8ZWvRw+igwfT7/dzWtYvvwRa9b/rLk9rbMZMZGBDqauhubUAu4TfBKBEC81t8jUE+2WSBlYbhSrHEZEtX86iV6lCtGdP+PTmrxd5Lc1JGYnM6qsklBp+KKK78kq+fuWV8DK6CEVkHiYyIus2BgRblXeSL9x2kGjyJKKsChWC6wuwyxovyOf13y/Ssi6+mNcJ5SzS11+HfoZLiAeRDQEww0RkvRzkWwqgs3ZeGUDFUOkdEVnHjiz6iBHOam/+CvH5gknMSGZOicqYxoroJk/m6+uvdyanS1BE5nEiIwpuYwCrioezp+d0L2SseeJdVp8+fJ2Vxaa5LrmE6O23vSNfvPN4oaziYqLTTuPrnByisWM9Z68xJiID0ADsc6ljJEQG4AwAK5yUIUPYTj12LItdtWqgVWw50rLTJLT7qg1FZE6G76GI7rvv+LpWLU+r9Soi8ziRWbWx447j4+rVkeVLpRGF1WyLl+SLVx4vlXXgAE8vyvj27Yk2bw79vCQiViKbCSAPBo+5GpGtB/temgignEW+7gDmAJilTUs+AQv/TAD6AygCUNSwYcPQNZEL3RdeGNkip0xnJi2ra/mMcAuqBQWhic7n050IhnrhuAxFZB4mslAfUwDR009Hni9V13iUfMkra+5cotq16dig4c03rZ+XZERNZGCvuc9p50Yiqwf2olsOwHRojglNeXsB2AegMdh77v8A9AtVXshOvXu33oFj+QNlsFojMyp72O03mzeP3bjL/WbG0VlBge4EsbCQvd0CROPGRfJ/JRWKyDxMZHYfU3I/2dlnR5Yv1Aef17XulHzJLev333XzfwDRDTe4vhE/FiJ7HMA2AL8A+A3AQQCvm9IcIzhT/IUAPjVc/xvAs6HKCzvNcvnlegd2OqS2Ginl5vJ1bi5fOyFDCdkYnGyYnjGD48xecj0ERWQeJjK7j6mFC7ldnXFGZPlCTcF7yeitks8bZfn97AtOKt80bMgOjF1CXNTvzSMy7SgATAIwxiJ9NoBvANTSrl8BcHeoMsIS2XffseM46S9IhlDkY6e1KMlM/rFOpieJIhvWS8Ok5csT/fNP6Oe6BEVkHiYyO2zbxu3q+OMjz6ugECl++IHNWgFs5ur++9nsVZJh966KxbLHG0KIbwF8C6AmgNEAIIRoIYSYCgBE5AMwDMAiLa0A8GIMZQKnnw7ceSfg9wfGG610GLFgAbBxI1vkmDABEALo0IGvO3UCevXiayH4fmEhX4eCTDtoEDBpElvxmDSJr2UZcud9rVpAbi5w6BCwYgXLmMo76BW8gRNOACpWBHbvBvbudVsahXTHqacCn30GPPggv98efxxo3RpYv95tyRhW7OZmcKTBdccddGwk9sEHoacFo51LdgKz8oixbGO50gLIsGHxKTfOgBqRebvN20FOr3/1VXT5FRSigdFeY4UKRP/9b/jlmDjB7l3leic2h7CdWhJEu3Z8bNqU6NAhe4KIVrsnHMKpvhrvX301H+vWjb3cBEARmcfbvB1ku/KIRplCBuHPP3W3QgAbfwjnRisOSB8ik4uVpaXsQRZgb6+hFjmj3W9hB6fkaFXuPfd4isSI7BtHOgcA5QGsAq/hrgPwqBY/Dc4s1ozT8n0P4CkAwnR/NoC1TmSJmsjkSH/UqOjyKyjEirfe4o35AKvrJ9jbR/oQmRHSn1DVqkR//BE6bahpQKLItHkiddBpLPezz5zXL0nIUCITACpr5zkAvgRr2U5D+I3+rQGsBCszZQP4HEB7w/0eYCs4iSUy6YbkuuuC7yVCE27ePKJNm9hx7IMPEl10EdspVFqBqSFftAhX1quv8sZp+Y67806i995LiHzpSWR+P9Gll+o/Xqh04UZk4cipVy/9jzH+sfKPGTWKv5CN+9D8/uBN123bqhGZxwKAigC+BnCBQyJrBaAYQAUtbxGA07V7lQGsAFu1SSyRSYeTJ54Yn7XhcHmkIV9jOPnk8F4gMmWfltflixZOyvL5eK+stNdYo0ZC5EtPIiMiWrtWV8f/9tvg+9FMA1qlC9dZW7YMvG8mMUm4ao3MM0EbTZUA2A9grBY3DWEs1mjpngTwJ3jD/2OG+IkArgbQKOFE5vOxHUKAZyeMiGZtOFQeGWrWDGzL7drFvyy3rFmku3zRIpKyVq9mo8OAvkVq4MC4yZe+REZEdPfdXJVLLontyzTUyC3cn2neVC2JD2BCW7uWz6WNPKW16JkAoBqAJQDOgjOLNU0AzNVGX5W1qcU2AJoDmK2lCUlkiMQsWyg88QS3pxNOCLaJZ9WeBwxgrbPFi4n27w9+nlUe6bPq8cf5/pEjuuuPjz+2zxfuhZWsPEq+yHHwINHOndGVdfAgu7Eyj97jIF96E9nOnfqC44cfBt6LZie78Yc35gv3Z1rdl0aIfT79hXDnnfGdw44RmU5k/BNgBIBhprj2sLZYUwDgYVPeQgB3AvgVbAVnG4DDAJaGKzvqERkR0eHDPF0NEDVrxhtXv/6aaNYsogkT7F8oAM9knH8+t9mZM4l27OBnmvsAwBZEDh/m+//7n16e0Rh2qL5jh2TlUfIF4+OPiUaPZgW03r3Zv12zZuweSz6ve3c2DRhNWfPmEdWpE718FkhvIiMimjhR71zR7jh38tUR7s8Mdb9LF457993o5EsQMpHIANQCUE07rwBgOdiuqBOLNdcCWAi2H5oD9gzRzZQm8VOLEnv3Ep17bmC7swtlyvCWlaZNg63jAESnnKJPDRnD4sV6edKN0qRJepwa8aSGfERMTNJNjl3IySEqW5bPGzYk+uST4PV+J/L17x+Y59RTnfmQtEH6E1lpKZMYwKQWClajNOOaVvfu1lOJ0YzI8vL0r9aRIzluyJD4axbFgAwlsnPAHhnWAFgrpxABLAZbq1kL4HXomo0tAEzVzrMBTIHuLHaCxfOTR2RERL/9xs4RK1dmD9JXXqmTW9euRKtW6YYEZJv96y9+QT3yCE/LV6pEli81Yx7pmqhiRd2NklqDSg35iPidI9dVK1bktBMnsk3YxYv5/929m/Nv3Kiv/ctwwQVE338fmXz5+URTpvBHFMDLK0uWWMsXBulPZEQ8rQjwNKNxftcMq3Uz45pWQQHHGf8MJy5bzPfz8vhaGiaeP5+vZUPyyDpZJhKZl0JciMyMaLTahg3je+3a8VRTvXpsweG22/Q8cqry9ttjKysdtQK9LJ/fHzjN3Lo10U8/WctjRGkpb7Mwf9jk5elT2k7lW78+cKpx+PCIZ88yg8j8fn367u67g+/LkZiZdHw+oiZN6NjXg3mkVlgY7LJF3pPPMWstEgVa2c/LC3RDc++98VmUjQMUkaVwm7dDIvYmzZqlr5+UlCS2rFTcp+Vl+bZu5f8tOztyz8/z5/NHziefsFsq4xqaEER9+zqXr7SUqFUrfVr7vPN4hOcQmUFkRKwdmJXFf9jatYH3jF8JRjIzfmWEGpKHajijRlnf9/n0kZkxfPFFbPWMIxSRpXibTxaef57b7kUXuS2JQqTYvJn/u/r1Y3/WwYO8zt+9OxMZwNbwI/kwX7GC9yACbK/x2Wcd5c8cIiPSXXV36RI8ujKPxIzkYtTAihZ2629mIjMulLsMRWRp0OYTDb9fN1L8xhtuS6NgxJ9/Eq1cGZoIpNufevXiW/Ybb/CgAYjc/N6+fUQ33aS/EydPDpvF7l0VixsX7+LRR4GqVYGPPwamTdPjze5XsrMD8w0dyj+pU0hXLcY8S5YA48ax24PevdndzNVXB+edNYvzKZcuCqmAFSuAb78FatcGevZ0WxoFI4YOBS66CHjqKfs0Varw+2/XLuDAgfiVff31wP/+B5QtCzzzDAenqFIFmD4dGDaMrzdujF4OK3ZzM8Tt6/Tll5nly5UjKioKvGc1EnOi8WOG1YKmcV0MIGrePPBcalaWLatrSbqs9AE1IkuPNp9ISHXtBx90WxIFM848k/+b8uVDrzfJJY4FC+Ivw9tv87PLlIncnuzYsZx32LCwSe3eVa53YnOIa6eW6sYNG+pajH6/rtghg5xmlMTSq5ez5xunKrt3Z81G+QwjgQFEjRvr58b9O/HaiR8DFJGlUZuPBXbrwL/+ym02KyvYekiofJmkTJFs+SQOH9btG0qlsyNHrNPKD2+plR1vyHdh/fq8wd5pnUaPpmNajGGQmUR26BDvewCIOnXiP92o4JGfrxOP8dwpkRFZ7x3Lzw8e9VmFli1dJzEi+8ahQgq2+Vhgp9LdqhXHN20aWb5MUm9PtnwScl9fvXpEDRrw+ejRwemmTmU3K1J3IBE4fJioRQsuo3Nn53V65BGOe/jhsEVkJpERsdqp/AMlqdkRWLQjJLMyx8CBwbvgrYKDoXQyoIgszdp8tLDaVDtwoN5epV1FJ/kyccNxMuWTmDmT03XtyurxAI/QVq/W03z7ra5deN55rBiSKLzxBpfTqpXzOj34IMf/5z9hH5+5REZEtHSprllz1VX8I1qNpKRFDwknQ36r5xiDeYrRGE44QY3IVPAOkRHZt2ezXUUn+cK9hJOVJ13lI+KXv3G6UBpQb9iQPzyOHtVNig0YkPj3zcaNXFajRs7rJEekjz8e9vGZTWREui3GSpWI1q3jOPNIyu4HthseW1n7sNozNnCgPaF5YGO0IrI0bfPRItrtIuH6k5t50lU+qYTzyit8vX8/G4I2f0jXqsVmzBKNDRu4vMaN+dpJnQYP5ntPPBH28XbvqvRUv7dCfj6rih44AHTpwteDBwemGTwYKCzU1eHbtwfy8lhVf8gQ/iuGDOHrvDxg0yY+HzSI1foBoEGD4LKXLwdKSlh1WaJ6dT4+/TQwfHjcq6ugEBVkGzeiTBngppsizyf7jNt50lU+APjuOz6ecQYfK1UCVq4ExowBypfn905WFvD220CdOqGfFQ/4fHzMznZep6NH+ZiTE325VuzmZkjo1+n+/UTnnBP4hSBNUhnXtOQwXY7IjOr0xmuj12jjSC0/n++ZtSPlvYIComuu4etTTnHdeDDUiCx923wkMK9j3HOP3m4zbQ3K6/IRsXZiuXKc7q+/gu//9BOblJoxwzp/IiCVT0491XmdBgzg+GefDft4u3eV653YHBLeqTds0P/8cEQWav3L/GdYTUPOmxc81ShJ7+67dTnM+9ySDEVkad7mncLchrt25Wujp18n+YhSXyvQ6/IREf34I9+PxTFrvCEdCB9/vPM69evHcVOmhH283buqTPRjuRTFKacA/fvzlF7lysDkyRwAnm4sWxbo2JGvhQDGjwc+/RRYvVp/Rm4uxwuhx3XowMcxY/T4JUuA4uLA8uV05rPPAuefD3z1FTBqFPD++/Gvq4JCJDC34alTedr71Vc5/uWXgZNPBu66i/uJXT5At6JTtqx+P1RZicqTrvIBwLp1fJTTil6AnFo87jigXz9ndVJTi1Hi6NHg6UKrhchIRmRmmEd5ZjX//Hz2yFuhAl8b1WWTDKgRWfq3+Vjw9ddEHTrobbdJE/YS7QGN24zGY4/x/zFkiNuS6Fi9mmU691znea6/nvO89lrYpHbvqsxR9jAiO1sfhRnRowcre5C2GDl8OCtzNG8emC43l+NDKWkMH85l5OdzMI78AP4qqVsXGDCAr0eNiq1OCgqJQm4usGgRMHs22xDdsIHtLbZtC6xa5bZ0mQuzoocXIEdkWRFQy5EjfIxhRJaZREbERnslevVizcP33weeeIKn/4h0rcWSEk5XUMDpVq/m+Pbt7cvo0IFJceJEDkYUFOjTlwUFrF00axawZk08a6mgED8IAXTrxoaDn30WqFmTDQlfcAFrA2/Z4raEmQcvEpkkJbNB9lCQU4tlol/pykwikyOtfv2AihWBmTO5k+bn8/3Jk1lN1LjGlZ8PjB2rW88vLgaWLrUv47LLOD0QrIJ65Ahw6aV8Xq8er9kBwOjRcauigkJCkJPDa2QbNrDnhnLlgDffBE47DRg5Mr6W1RXs4fMB33/P56ef7q4sRkgL9iee6DxPHEZkrq8PmENS1guMFjukwcqzzmJbYQUFbOHDSkVfwmjZIxScqtVu28bW8IUIdgaaBECtkaV/m08UNm0iuvZava+ccALRq6/Gx7efgj1+/ln/vb2E++9nuR56yHmeSy/lPPPmhU1q965y3NkAZANYDWCOdj0NwCYAJVpobpPPZ0gzO1w5Se/U//yjeyp95hmOi3ZnvhmRqNVK0zLXXhtbfaKAIjJFZDFj+fLArSYtW4Z25+F16/Jel2/2bP6dO3e2fqZb6NaN5Ro+3HmdOnXiPHa2PA2IB5ENATDDRGS9HOTb77QMcqtTv/ce/xTVq7O7l2hsnlnB3EDltc8X+Gf6/bwpsEwZHpV991386uYAisgUkcUFPh/RtGlsiV32neuus3b94vV9Wl6Xb8wYOjZb5CVId1V9+zqvU9u2HLdkSdjHx0RkABoAWASgY1oSmd9PdMkl/HNIyx+RWqF2gnANVtpFu/762OsUATKRyACUB7AKwDcA1gF4lPR27WSmYZyW73sATwEQWvwCwzOfB5AdTpa0ITKJv/9mi+Zyw3/58uyiY/9+PY3XLWd4Xb6bbiKnm4iThv37WaacHHah5bRO0k3QihVhi4iVyGYCyAPQ3kRk6wGsATARQDmbvEcBFAH4AkB3mzT9tTRFDd3apb5une7q4IYbIvuicopwDXbTJm4EWVlEP/wQc5WcIkOJTACorJ3nAPgSwIVOPtAAtAawUptuzwbwOYD22r0qhuf/D0CfcLKkHZFJ/PJL6PUzq32aTvZnRjNjkqyykiWf9Pu1fHno5yYTX33FMp15Jl87rZOsyxdfhC0iaiIDcCWA57RzI5HV0zprOQDTAYywyV9fOzYG8AuAU0KV52qnvuoq/knato1OucMJwv25t9/Ocf/+d+xlOUQmEpkxAKgI4GsAFzgkslYAigFU0PIWATjdlCYHwIcArg1XftoSmcSKFfrLCmDr7NInVjTr0dGuYSerrETL5/OxFw+AaNcuZ89OBl55hWUyrvM7qZOciSouDltELET2OIBtGgn9BuAggNdNaY4RXJhnhX1JuNqp9+4lqlmTf5a3305cOaH+3I0bea0sK4uNfiYBmUpk2miqBMB+AGNJb6NOZhqeBPAngH0AHjPd+wjAXvCacuZNLVrBbv1MrqV4bcSTzDyR5vvlF75fp07oZyYb/fuzXKNG8bXTOp15Jt9bsyZsETFNLR5LbBqRaUcBYBKAMRbpq8sXAYCaAH4CcEaoMlzv1C+8wD/LiScSHTgQ/+c7+XNl577llviXb4FMJTIZAFQDsATAWU5mGgA0ATAXQGUtfA6gjSlNeW1qsbNNme5Pp7uBv/9m1ezy5fX237w50fr13lqD8vIa2bx5fK9Dh9C/dTLxzz+8hQggevrpyOrUrBnHf/992GISQWSLAXwLYC2A16GvN7QAMFU7b62l+UY79gtXhutEdvQouwMHiEaOjO+znf65GzawR+vsbN4vkmBkOpHxT4ARAIaZ4ixnGgAUAHjYlLfQIt1NAJ4JV7brbd4N3HEHBXzMAaxKfvnlfO62VqCXtRafeILj77nH+nluoWdPlqtKFdYzcFonqenoYAYqLkSWjOCJTj19Ov80HTvG97mRNNibb+a4fv3iK4MFMpHIANQCUE07rwBgOXg92MlMw7UAFgIoA14LWwSgmzY6k/nLAHgbwD3hZPFEm0825FaUkhL2mWUcoZUrx4T22WeBX+5qHxnjllv4d3ruOevnuYXSUvbDCLAx9Ouuc1anE0/kPL/8ErYIRWSRQE7tPfJIfJ8bSYP98UdeJytThrUZE4gMJbJzwBv812izCiO0eCczDdkApoBV778DMEGLrwPgK8MznwZQJpwsnmjzbmPPHp6SMiqFAERNm7L1Hau9aJkIn4+ofn3+bVatcluaYBw5oo/GKlQg+uST8Hnq1uX027eHTaqIzCkOHyaqUYN/mnXr3JXlxhtZjv79E1pMJhKZl4Lrbd5rWLuWTcUZFUOE4BmS6dMD96NlGlau5N+jYcPY97UmCkeP6s4yy5cPb3ZPKtj9/nvYR9u9qzLTaHAoLF4M7NnDFqXdtir90ENsdfyVV5R1cYXMwZlnAuPGcZufPx+49lp2e7R4MXDzzUCdOsAtt7BRb7/fbWmTi3fe4WPv3oGOfb2E7GzghReAG28EDh3i/0pauLeCcuOSALz7Lh9793ZXDoB9P/Xpw3/0mDFuS6OgkFyUKcNeJN56C9ixA3j+eaBVK7awP306u0Jq2pQ9WRw8LVc9pgAAGfdJREFU6La0iYffz546AG+8n0IhK4vd/Zx4IlBUBDz5pH1aRWRxxpEjwHvv8Xm0DWXBAnZvQRQYV1jIYcECPZ6I0xrjzJCjspdeArZti04mBYVUR/XqwB13AJ99Bqxfz/2iYUN2GzJ4MNC4Mfv9syI0qz4JhO5/ycoTSb5PPwW2b+d6t2xp/SwvoUoVYOpUPh85Eli3zjqdHK0pNy5xwkcf8Vzt6adH/wwrzcSCAn2uv6CA4yIxfXXNNZRIdVuoNbLMbfOpjKNH2eC30ep+3bpEEycSHTyop0sX9fsuXSghSmiJhrRWdMIJRJ9/HnxfmgY8ejTso+zeVa53YnNwtVPfdhv/JCNG2KcJp3k4b17w3rD8fL2jSd9mkRgj/vZbOqaW7ECzJ1IoIsvgNp8O8PuJPvzQntDSYUP0qlV8Xbky0e7d8f39Eo2//iJq04blL1uW369//833jh6lY8o8DqCILBwOHyY6/nj+Sb791j6dk68nK+sd+fmBhOaUxCTkZsMEuG1QRJahbT7d4Pezny5p0EAS2qRJbKUnlU1USWe/sRoudwuHDxPde69et1q12OOIfK+VLevoMYrIwuHjj/nnOO005w011FeXlT3FaI2JEvHGUanOumNHbHU1QRFZhrb5dIUVodWrx4QWaf/zgtHgtWvp2IxMnPt+0rF8OdGFFwbWU047OoDdu0ope0hIbcXq1YPvkWHRVQigSxcgL4+1pbKy+DhoEDB+PDB8OKsMDxkS+IzBgzkYMWRI8AKvhHkB+Nxzge7dWZ21W7fQCiIKCpkMIbiPFBUBH3wA5Oay1uOgQYHpQvU/gO+Z+3Ei8oTL9/jjfH3bbUDduqGf43VcfDEr7KxeDXz0EWukPv88MHdubM+1Yjc3gytfp0eO6JvynCy6yulFY/D59HRyrj6WNTKrKcziYv05d98dt+pDjcgyr81nEnw+om7dAvurdINyzz3eXiO79Vbdwo8DE07pDrt3leud2Bxc6dSffMI/xamn6qQTqiH6fES5uYEdQ16bSYwoOq3FUI0bIBo2LG7VV0SWgW0+kyA/CvPzWctR+r+SoVMntt5ulccLWotA0jxheB2KyEJB+tF56CFuQHJh1TxC8/uZiOR9KzKbOzdYq3H+fM5XUBBoTzGcMVGrBeDrruNjxYpEf/wRl+orIsvANp9JMGsa+/1MaOeeq/erunXZqvy+fdZ5JJJpNHjrVvaAIURSPcZ7GYrI7GCcVvzmG44zjqCMihpGUhk0iEdm5inGeGsVWS0Ad+3K58OHx6UIRWQZ1uYVGD4f0axZgR+kxx3H0/YOnDwmHIMHs0y9e7stiWdg965Syh6ffgrs2gU0awacfTY358OHA9NIRY1Jk1jhoqCAFTuGDg1Ml5cHtG+vXxOFt9wRCmSzAPzww3z+zDPA7t3RPVtBIdORlQVcfTVQXAzMmQO0awf8/TebVjrnHOD889layPbtyZdt1y5gyhQ+f+CB5JefarBiNzdD0r9OpYO/Bx8MHHVZ7fuSihrm9SvjCM4qTTSjtHALx5ddRvEalUGNyDKrzSvY45tviO68k51Dyj4tBFHbtkTPPuvIQntc8OCDXPYVVySnvBSB3bvK9U5sDknt1EeO8MY8gPdpmRddzdN6UlFDOo8zkpaR9Hr1isxyhxXCLQBLnz9ZWbxnJgYoIsugNq/gDAcPEr3zDlGPHrx/S/btrCz2ZP3SS+xDLRH480+iqlW5vJUrE1NGisLuXZXZU4vLlgE7d7IF7V9/ZQOkBQXAhAl83zyt99NP1vtB5N6yWrX4eubM4L1lkU4vdujARoYnTNDdNQjB14WF7CLhkUfYInafPsBXX0X2fAUFBXtUqMCGw//3P+CPP4DXXgO6duXpyE8+Afr1Y3cy3boBb7zBU5LxwnPPAfv28VRn69bxe246w4rd3AxJ/TodMIC/eh54IHAEZNwTZqVWb7wvR0xGO29We8sSYVrG79fdnteuTbRxY1SPgRqRZU6bV4gNu3cTvfgim1fKytL7evnyPBPz7ruBBosjxYED+izRxx/HT+40gd27yvVObA5J6dRS3VWaSlm8OHB6sGlTCiAxM3nZ2VOUDdC8tyza6UUnOHyYpzrkPrgoDIoqIsuANq8Qf/z2G9EzzxBdfHFgv69cmaf+P/yQqLQ0smdOnszPOP9873qAdhGKyIyQo6+TTuLj3LnB61zmPWREwftBzGtodnvLfL7E1mffPqJzzuHy2rQJ3twZBorIMqDNKyQWW7YQjR/PBGTs/9WqEfXty6OrI0dCP6O0lKhBA8733nvJkTvFoIjMCPNoavr0QG1F874tJ8+QpCXdEiRyb5kVtm4lql+fy+vTJyLyVESWAW1eIXnYsIHosceIzj478D1QqxZrGdsR2tSpnO7MMxP/8ZuiUERmht8fbKrGqasVsyr8vHn6Gpl5RJaXx/eTgZIS3tAJEN13n+NsisgypM0rJB/r1rH/rWbN9HfCDTcEO5E8coSoSRO+//rr7siaAlBEZoWRIymIyOz2bRnJzEo13mh/MS8vcueZ8cJHH7GBUYDoueccZVFElkFtXsEd+P1s07VyZe6bN98cSGYzZnB848bhpyAzGIrIzPD7idq3p6BpQLk3TKax0jq0so0myc04Aot1U3S0ePllOrbn5cMPwybPRCIDUB7AKgDfAFgH4FEtfhqATQBKtNDcJv84Ld/3AJ4CIABUBDAXwA/avTFOZFFElkFYtoztpAJs45WIP4LPOovjpkxxVz6Pw+5dVSaRqv2eBRHvEVu6lK+vuYbNUr3/fmA6uW+rbFne1yVx2WUcjJD3x4wJ3vdlzp9o3HorsHkz8OijwLXXshmuFi2SV35qoBRARyLaL4TIAbBCCDFfu1dARDPtMgohWgO4CMA5WtQKAO3AxPgkES0RQpQFsEgIcTkRzbd5lEKmoU0b9r11+eXACy8ArVoBNWoAa9cC9esDN9/stoSpCSt2czMkVWvx6qv52KmT9egpWmvWyYSUcd68QFmNX3mVKhFNm2b7CGTgiMwYwCOprwFcAB6R9QqTvhWAYgAVtLxFAE63SDcZwO3hylcjsgyEnDUpX563zQBEEye6LZXnYfeuykzLHtJqhjTGuXs3eyvNyWHLHnL0tGQJMG4cGw+dr31UEwE9enD84sWBzyWKzUhwNJAyPvwwH4cMYWsfQ4fyV95xxwEHDnD83r3JkysFIITIFkKUAPgDwCdE9KV26zEhxBohxEQhRDlzPiL6HMASADu08BERfW96djUA3QAssim7vxCiSAhRtHPnzjjWSiElcOutbB3k0CFg/XqgZk3g9tvdlip1YcVuVgFANoDVAOZo19PgYC1BS1sFwDYAz4QrJ6lfpxs38pfQiSdGrsDhdD0t0TCWK2U1H48/no9t2xIdOhT0CKgRWTUwMZ0FoB54vascgOkARlikbwJeC6ushc8BtDHcLwNgPoBBTspXI7IMxcGDuub06NFuS5MSsHtXRdLZhwCYYSKykFMwhryTtbzeIrJ9+/gnqFDB2uyUdKBpVtMfODDYE3Q8NBSjncq02tNm3D6weTPRCSeQ3R6zTCcy/gkwAsAwU1x72d5N8QUAHjblLTRcvwzgKadlKyLLYPzxB0/7Hz7stiQpAbt3laOpRSFEAwBdAUx1kt6UNw9AHQAfR5o34TjuODb0+88/wHffsWLGoEFs8DcrS1f+aNcuMJ8Q7MMoL09PK40EG438Rgo5TThkCNMQoCumjBvH960glUqsMGEC0LAhMG8e1/ett4CHHopOvjSCEKKWNv0HIUQFAJ0B/CCEqKfFCQDdAay1yL4FQDshRBlNUaQdWHsRQojRAKoCGJT4WiikPGrVYgWPnBy3JUltWLGbOQCYCSAPhi9U8IhsPYA1ACYCKGeRLwvAUgANANwCmxEZgP7gBfOihg0bJo/eidh8DEA0ahRfm81OWZmtkvFmD9Gx7hUL54MsEisjVhu633yTji0wG56FDByRgTUOV2vtdy20KUQAiwF8q8W9DqCyFt8CwFTtPBvAFDB5fQdgghbfAABp8XLK/bZwsqgRmYKCM9i9q5x0+CsBPKedG4nMyVrCPdCmXEIRGbnZqWfPJts1MKtpRaM5KivHmwUF+hRgNFOFVjI4JTG7NTKZXzrrM7lOz0Qi81JQRKag4AyxENnjYEWNXwD8BuAggNdNaY4RnCn+DfA0zC8AdgH4C2E2iSa9Ux88qG9Q7Ncv8MUv18HMZGY0QyUVPoykJh1whnOOaacUYh4VhhrpGTdiy7KM630yfvBgXeljxYqARygiU0SmoJAKsHtXhd0QTUT3A7gfAIQQ7cEL4jcKIeoR0Y5QawlEdIM8F0LcAqAFEQ0PV2ZSUaECb26eNQt46aXAda5Ro3g96euvA/NkOdy1MGYMb7SeNImvJ0zg9S65njZmTHAeomCHnj16sPPPsWMD19+I2Blor15A3768wVtuyJYbsdu35/jDh3mbQYsWylmfgoJCesGK3ewCAqcWw64lmPLeAi9OLRIRvfoqEcBuXYyjH+NIq3v3wNFUXh5rL4aaWpTPcDpVGGqNLBZ1f79f3xz92mtBt6FGZK4GNSJTUHAGu3eV4HveQYsWLaioqCi5he7ZA9Suzec7dwLVq/N5797AzJlAfj4wcSKPdIiAwYOByZN5JDTTYMnI77fWWCQKHMX5fIHXRMDw4Ty6mjkzcFRoLA/Q7/XowVqVV10FNGumj9bks+RIrEMHHpl16gTUrcumq8qWDRBPCFFMRMqGlUsQQuwEsDmCLDXBU/WpjFSvQ6rLD6RmHU4iolrmyMy0tWhGjRqsYr94MdtBu/FGZ/lKSgKvhwwJVr8ni6nCFi2AoiImM3l/0iQmxsLCYHuNEycy+fz0E6eTU5UAsG0b8MEHwJEjwPjxbNFj0iRg0SLeIgAA32tGJ+66K4jEFNyHVccMBSFEUap/eKR6HVJdfiA96nAMVsM0N4Nr0yxPPUUEEPXsqceFmuozawRaqcmb442WQqTnaHO+UJqOBQVkuTUglLbijz8SCUFUrhzR779bVh1qajGlQjr8X6leh1SXP13qcKwubgtgDq4R2ebN/HNUqsSajBJWa1xmqx7mdHLNKpzZK6s1MztNR6v9bFaOQM3PlOt4t95qW/V0atCZENLh/0r1OqS6/OlSh2N1cVsAc3B14fu88/gnmT07MN6sDm+2NG9MZ9wfZje6CrWR2mp0ZyQrqfAh09gRmd/PJrikx+iSEttqp1ODzoQAoL/bMmR6HVJd/nSpw7G6uC2AObhKZKNG8U/St68eF+kG5XBw8jw7Sx1mrUU7EgOIBgwgGj6czzt0CCmSIjIVVFAhlYPrApiDq0T27bf8k9SsyW7IozUZZYdInmceBRYU2G8NOOUUPjZqxMcKFQLzvv9+SLEUkamgggqpHFwXwBxcJTK/XyeFZcvs16ukVXxpwcOYP5SVeqeWPsKN2vx+3QqJtNZhDkLwUToNDQFFZCqooEIqh8x0rGkHIYDu3fn8/fd1B5xGlXoh2MoGwOrwRHxOFN5Kvd3zJkzg+A4dAtXxBw3ivWkDB/J1+/ZAz568H+yllzj/7t18PO444NJL2RrJddexpZLCQmDYsOit8Su4BiFEDSHEJ0KIn7RjdZt0C4QQfwoh5pjiTxZCfCmE2CCEeFsIkfR9FxHU4WYtzU9CiJsN8dcJIb7VnJwuEELUTJ70cZG/rBDiBSHEj0KIH4QQPZMn/TEZYqqD4f5sIYSVJwhvwG0mNQfXrRwsX04EEDVu7MxQbzz9kRHpo7aePYkefZSoc2eiypUpaMRVpw5Rr15EkycTFRcTHTkSdZFQIzLPBQDjAAzXzocDGGuTrhPYE/UcU/w7APpo588DuNOLdQBQA8BG7VhdO68O3uP6B4Cahmc9kirya/ceBTBaO8+SdUmlOmj3e4D9Sa5NtvyO6+m2AObgOpEdPUpUqxb/NGvW6PFmDcR4KoHs3s2akgUFRKedRpSVRUHE1aQJ0dln63vDYnUZY4AiMu8FsIuketp5PQDrQ6RtbyQysFeKXQDKaNetAHzkxToAuA7AFMP1FC0uB8BOACdp9XkeSdayi0V+7XwrgEpeb0dh6lAZwAoAZ3iZyJRlDzOys4F//Yun7t5/Hzj7bI6XTi8PH9anBsePD7Sy4dSp5tatwPLleli3LvB+VhaQmwtcfDHQpg0f69WLXx0VUgF1iGiHdv4b2DmtUxwP4E8iOqpdbwNQP57COYSTOtQHv/AltgGoT0RHhBB3gu25HgDwE4C7EymsBaKWXzptBTBKM7b+M4B7iOj3RAlrg6jroJ2PAjAe7PXEs1BEZoXu3XUie/hhjjNbsh8/nk1NGWFlosrvB374IZC4tmwJzFeuHNCyJZNWmzZAq1ZA1aqJq5+CJyCEWAigrsWtB40XRERCCEqOVJEhUXXQPG/fCSAXPNX1NNgLx+jopbUsJ1H/QRmwo9XPiGiIEGIIgCcB/DtqYW2QwP+gOYBTiGiwEKJRTEImGIrIrNCpE1CpErtv2bIFaNhQV8oAAu0d5uay3URp49DnA66/Hlixgklr5UpdIUOialXgoot04mrRgslMIaNARJfY3RNC/G5wlVQPvF7kFLsBVBNClNFGZQ0AbI9RXEvEoQ7bwVOjEg3AXuWba8//WXvWO+A1nrgigfLvBo9iZmnx7wLoFw+ZzUhgHVoBaCGE+AXMFbWFEEuJqD08BqW1aAXpowxgg7wSRjKTWLqUjQ0fdxzQoAHw9NM8oiooAGbPZhI74QTg2muBZ54BvvmG4+bOZSv1F12kSEzBCrMBSO2xmwF8ECJtAIgXN5YA6BVN/jjCSR0+AtBFCFFd06jrosVtB3CGEEIaVO4M4PsEy2tG1PJr/8GH0AmiE4DvEiuuJWKpw3+J6AQiagTgYgA/epHEAChlD1u89hoRQNSxox7n9xP170+W+7WMoUEDottuI5o2jejnn+OqmJEIQCl7eC6A17kWgdeGFgKoocUH+PsDsBysFPEPeG3jUi2+MYBVADaARwPlPFyHvpqcGwDcaogfACavNWBSOD7F5D8JwDJN/kUAGqbaf2C43wgeVvZQ/sjssHcvUEv7GHz+eeCLL3hv1t69wWnr1AH69AHatmXFDOnbLEWg/JEpKCikMtQamR2qV+cNyIsWAbffrseXKROoTfjBB8Bzz/H0YI8eromroKCgkKlQRBYKBQXsuLJZMx5llZYCU6cGrml17gxUrsxWORQUFBQUkg41taigphYVFBRSGkprUUFBQUEhpaGITEFBQUHBFkKI3kKIdUIIvxDCduZGCHGZEGK9Zqh6uCF+mhBikxCiRAvNw5TXWQhRrBmMLhZCdAwno1ojU1BQUFAIhbVgw8FT7BIIIbIBPAve77cNwFdCiNlEJPfOFRDRTIfl7QLQjYh+FUKcBd7nFtLEmhqRKSgoKCjYgoi+J6L1YZK1BLCBiDYS0WEAbwG4KlQGIUQlIcTLQohVQojVQoirtPJWE9GvWrJ1ACoIIUJajVBEpqCgoKAQK0IZHgaAxzS/chMNpPQggMVE1BJABwBPCCEqmZ7bE8DXRFQaqnBFZAoKCgoZDiHEQiHEWosQclTlEPcDOA3A+WCfZ/dp8V0ADBdClIBtO5YH0NAg05kAxgK4I1wBao1MQUFBIcNBIQwPO8R2ACcaro8ZqibdjUypEOIVAMO0awGgp9W0pRCiAYD3ANxEmuHoUFAjMgUFBQWFWPEVgKZCiJOFEGUB9AEbLIZmdR9CCAGgO1h5BGAljnu1eAghcrVjNQBzwZ6tVzop3HMbooUQOwFsjtPjaoI1YNIFiarPSURUK3wyBQWFTIMQ4mqwP7haAP4EUEJElwohTgAbHr5CS3cFgEkAsgG8TESPafGLtbwCQAmAAUS0XwhRQUvfGjyo2kREVwohHgJPR/5kEKMLEdm6MvIckcUTQoiidLJYkW71UVBQUIgH1NSigoKCgkJKQxGZgoKCgkJK4//bu3+XKqM4juPvD1hLZT8QxIaaKnIo8y8oImprCIKmsqkocA0aGhpqqCFwEgKlIEiQyiUIhxpqi7Ko6Cc4SBo5thh+G84xLqJJ6b3X8/R5gXju43Mezne5H57rfb6n6kHW3+wFrLCq1WNmtmyV/h+ZmZlVX9XvyMzMrOIcZGZmVrRigiw3l5yS9Lrm2BZJjyR9yL835+MbJY1Iepm3H+hZ5JprJfVLei/pnaRjhddzIm99MCbpoaS2RtVjZtYsxQQZMAAcmXfsAjAaETuA0fwa4BzwJiL2AvuB6/lp8/kuAlMRsRPoBB7XYd2LGWAF65HUAtwADkTEHmAMOF+31ZuZrRLFBFlEPAGm5x0+Cgzm8SCp/QlAABty65P1ed7PBS57GriSrz8bEQ3rAlKHepR/1uXzWoEJzMwqrpggW0R7TUPKr0B7HvcBu0lv5K+A3oiYrZ2Y+3kBXJb0XNKQpHaa65/riYgZ4Gz++wTpDvNmIxZtZtZMpQfZb5GeI5h7luAwqafXVqAL6JPUOm9KC6lD89OI6AaeAdcatNwl/W09ktaQgmxfPm+M1K/MzKzSSg+yyZrOyh3AXFPJHmA4ko/AF9J+OLW+Az+A4fx6COiu/5L/aDn1dAFExKccgndJzTjNzCqt9CB7AJzM45PA/TweBw4C5I8LdwGfayfmN/sR0pcnyOe/qe9yl/TP9ZD2/umUNNfF/hDwtq6rNTNbBYrp7CHpDil02oBJ4BJwj3TnsY209cvxiJjO2wsMAB2kL0BcjYjb+TovIqIrj7cDt4BNwDegJyLGC67nDNALzOT5pyLieyPqMTNrlmKCzMzMbCGlf7RoZmb/OQeZmZkVzUFmZmZFc5CZmVnRHGRmZlY0B5mZmRXNQWZmZkX7BXDABmtVKuLTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "####quick sanity check#######\n",
    "BUCKET = 'rio-mora'\n",
    "BUCKETS = ['ikh-nart','rio-mora']\n",
    "fig, ax = plt.subplots(1,2,figsize=(8,4))\n",
    "for i,BUCKET in enumerate(BUCKETS):\n",
    "    if BUCKET == 'rio-mora':\n",
    "        shp = './data/shp/RioMoraNWRBndry_2019/RioMoraNWRBndry_2019.shp'\n",
    "        pointsfile = f'/datadrive/{BUCKET}/wdpa-555609346_points_v1.csv'\n",
    "        \n",
    "    if BUCKET == 'ikh-nart':\n",
    "        shp = './data/shp/IN_boundary/IN_boundary.shp'\n",
    "        pointsfile = f'/datadrive/{BUCKET}/wdpa-99848_points_v1.csv'\n",
    "    df = gpd.read_file(shp).to_crs(\"EPSG:4326\")\n",
    "    pnts = pd.read_csv(pointsfile)\n",
    "    pnts = gpd.GeoDataFrame(pnts, geometry=gpd.points_from_xy(pnts.Longitude,pnts.Latitude))\n",
    "        \n",
    "    pnts.plot(ax=ax[i], marker='x', color='red', markersize=60)\n",
    "    df.plot(ax=ax[i], facecolor='none', edgecolor='red', lw=2)\n",
    "    fc,g = build_features(shp)\n",
    "    srtm = srtm_image_generator(fc)\n",
    "    REGION = g.getInfo()['coordinates']\n",
    "    PREFIX = 'processed/srtm'\n",
    "    DESC = f'{BUCKET}_srtm'\n",
    "    outbands = ['elevation', 'slope', 'aspect']\n",
    "    image_collection_export(srtm, outbands, BUCKET, PREFIX, DESC, REGION)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:15: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.\n",
      "\n",
      "  from ipykernel import kernelapp as app\n"
     ]
    }
   ],
   "source": [
    "#BUCKET = 'ikh-nart'\n",
    "\n",
    "BUCKET = 'rio-mora'\n",
    "BUCKETS = ['ikh-nart','rio-mora']\n",
    "\n",
    "for BUCKET in BUCKETS:\n",
    "    if BUCKET == 'rio-mora':\n",
    "        shp = './data/shp/RioMoraNWRBndry_2019/RioMoraNWRBndry_2019.shp'\n",
    "        target_months = (4,10) \n",
    "        #shp = './data/shp/MRW/Mora River Watershed July2011.shp'\n",
    "    elif BUCKET == 'ikh-nart':\n",
    "        shp = './data/shp/IN_boundary/IN_boundary.shp'\n",
    "        target_months = (6,10) \n",
    "\n",
    "    #need an automated way of determining growing season....\n",
    "    fc,g = build_features(shp)\n",
    "    REGION = g.getInfo()['coordinates']\n",
    "\n",
    "    years = np.arange(1986,2021) \n",
    "    # Define a FeatureCollection fc above\n",
    "    #Build landsat collection\n",
    "    cloud_thr = {'l5':25,'l7':20,'l8':10}\n",
    "    landsat = landsat_collection_generator(bounds=fc,cloud_thr=cloud_thr)\n",
    "    landsat = landsat.map(clipper)\n",
    "  \n",
    "    '''\n",
    "    ######## ALL BANDS #############\n",
    "    DESC = f'{BUCKET}_growingseason'\n",
    "    PREFIX = 'processed/allbands'\n",
    "    landsat_images = []\n",
    "    outbands = ['B5', 'B4', 'B3', 'B2', 'B1','NDVI','EVI']\n",
    "    if BUCKET == 'rio-mora':\n",
    "        target_months = (4,10) \n",
    "    elif BUCKET == 'ikh-nart':\n",
    "        target_months = (6,10)\n",
    "    out = landsat_aggregate(landsat, outbands, years, target_months)\n",
    "    #export these to bucket\n",
    "    image_collection_export(out, outbands, BUCKET, PREFIX, DESC, REGION)\n",
    "    del out\n",
    "   \n",
    "    ######## NDVI BAND #############\n",
    "    DESC = f'{BUCKET}_growingseason_ndvi'\n",
    "    PREFIX = 'processed/ndvi'\n",
    "    landsat_images = []\n",
    "    outbands = ['NDVI']\n",
    "    if BUCKET == 'rio-mora':\n",
    "        target_months = (4,10) \n",
    "    elif BUCKET == 'ikh-nart':\n",
    "        target_months = (6,10) \n",
    "    out = landsat_aggregate(landsat, outbands, years, target_months)\n",
    "    #export these to bucket\n",
    "    image_collection_export(out, outbands, BUCKET, PREFIX, DESC, REGION)\n",
    "    del out\n",
    "  \n",
    "    ######## PEAK BAND #############\n",
    "    DESC = f'{BUCKET}_growingseason_color'\n",
    "    PREFIX = 'processed/peakcolor'\n",
    "    landsat_images = []\n",
    "    outbands = ['B3','B2','B1']\n",
    "    if BUCKET == 'rio-mora':\n",
    "        target_months = (6,8) \n",
    "    elif BUCKET == 'ikh-nart':\n",
    "        target_months = (6,9) \n",
    "    out = landsat_aggregate(landsat, outbands, years, target_months).float() \n",
    "    #export these to bucket\n",
    "    image_collection_export(out, outbands, BUCKET, PREFIX, DESC, REGION)\n",
    "    del out\n",
    "    ''' \n",
    "    ######## EVI BAND #############\n",
    "    DESC = f'{BUCKET}_growingseason_evi'\n",
    "    PREFIX = 'processed/evi'\n",
    "    landsat_images = []\n",
    "    outbands = ['EVI']\n",
    "    if BUCKET == 'rio-mora':\n",
    "        target_months = (4,10) \n",
    "    elif BUCKET == 'ikh-nart':\n",
    "        target_months = (6,10) \n",
    "    out = landsat_aggregate(landsat, outbands, years, target_months)\n",
    "    #export these to bucket\n",
    "    image_collection_export(out, outbands, BUCKET, PREFIX, DESC, REGION)\n",
    "    del out\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}