Test_preprocessing_2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import rasterio\n",
    "from rasterio.mask import mask\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import mapping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the use of os in python "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Landsat image 2013/All_bands\\\\LC08_L1TP_141041_20130326_20170505_01_T1_B1.TIF',\n",
       " 'Landsat image 2013/All_bands\\\\LC08_L1TP_141041_20130326_20170505_01_T1_B2.TIF',\n",
       " 'Landsat image 2013/All_bands\\\\LC08_L1TP_141041_20130326_20170505_01_T1_B3.TIF',\n",
       " 'Landsat image 2013/All_bands\\\\LC08_L1TP_141041_20130326_20170505_01_T1_B4.TIF',\n",
       " 'Landsat image 2013/All_bands\\\\LC08_L1TP_141041_20130326_20170505_01_T1_B5.TIF',\n",
       " 'Landsat image 2013/All_bands\\\\LC08_L1TP_141041_20130326_20170505_01_T1_B6.TIF',\n",
       " 'Landsat image 2013/All_bands\\\\LC08_L1TP_141041_20130326_20170505_01_T1_B7.TIF',\n",
       " 'Landsat image 2013/All_bands\\\\LC08_L1TP_141041_20130326_20170505_01_T1_B8.TIF']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import os\n",
    "Landsat_folder_path=(\"Landsat image 2013/All_bands\")\n",
    "# print(Landsat_folder_path)\n",
    "\n",
    "# Landsat_band_path \n",
    "Landsat_bands_path=[os.path.join(Landsat_folder_path,i) \n",
    "               for i in os.listdir(Landsat_folder_path) if os.path.isfile(os.path.join(Landsat_folder_path,i))]\n",
    "\n",
    "Landsat_bands_path.sort()  \n",
    "Landsat_bands_path    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need a rasterio dataset object containing all bands in order to use the mask() function and extract pixel values using geospatial polygons.\n",
    "\n",
    "We'll do this by creating a new raster dataset and saving it for future uses.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data/products/Landsat_bands_composite.tif\n",
      "{'driver': 'GTiff', 'dtype': 'uint16', 'nodata': None, 'width': 7441, 'height': 7421, 'count': 1, 'crs': CRS.from_dict(init='epsg:32645'), 'transform': Affine(30.0, 0.0, 237585.0,\n",
      "       0.0, -30.0, 3146115.0)}\n",
      "{'driver': 'GTiff', 'dtype': 'uint16', 'nodata': None, 'width': 7441, 'height': 7421, 'count': 8, 'crs': CRS.from_dict(init='epsg:32645'), 'transform': Affine(30.0, 0.0, 237585.0,\n",
      "       0.0, -30.0, 3146115.0)}\n"
     ]
    }
   ],
   "source": [
    "# create a products directory within the data dir which won't be uploaded to Github\n",
    "image_dir= \"data/products/\"\n",
    "\n",
    "# check to see if the dir it exists, if not, create it\n",
    "if not os.path.exists(image_dir):\n",
    "    os.makedirs(image_dir)\n",
    "    \n",
    "# filepath for image we're writing out\n",
    "image_folder_path=image_dir+\"Landsat_bands_composite.tif\"\n",
    "print(image_folder_path)\n",
    "\n",
    "\n",
    "# Read metadata of first file and assume all other bands are the same\n",
    "with rasterio.open(Landsat_bands_path[0]) as Land_band_0:\n",
    "    meta_data=Land_band_0.meta\n",
    "    \n",
    "\n",
    "# Update metadata to reflect the number of layers\n",
    "print(meta_data) \n",
    "meta_data.update({\"count\":len(Landsat_bands_path)})\n",
    "print(meta_data)\n",
    "# OR\n",
    "# meta_data.update(count=len(Landsat_bands_path))\n",
    "# print(meta_data)\n",
    "\n",
    "# Read each layer and write it to stack\n",
    "with rasterio.open(image_folder_path, 'w', **meta_data) as composite_band:\n",
    "    for count, layer in enumerate(Landsat_bands_path, start=1):\n",
    "        with rasterio.open(layer) as Land_band_1:\n",
    "            composite_band.write_band(count, Land_band_1.read(1)) #read to convert into array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(7421, 7441)\n",
      "8\n"
     ]
    }
   ],
   "source": [
    "composite=rasterio.open(image_folder_path)\n",
    "print(composite.shape)\n",
    "print(composite.count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Let's clip the image and take a look where we know the training data was from the last lesson (on the NC Rachel Carson Reserve) just to make sure it looks normal:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3, 450, 1150)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x151008821d0>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAAD7CAYAAAA8XnP0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAc7UlEQVR4nO3df4xd5X3n8fenNgZalhgwEMsmtZM46lJUHDBepxvarAsJQdoFNnRDVQk3SWUtlLArRItRV23Z/JO06jbNthtKoCp0fwAlCbAEQrwEUikKBhOM+V2GFIoDC5ukprFSoeJ+94/zTLidzJ254xnsOeP3Szqac7/nec6958kZ88l5zpmbqkKSJEn98GMH+gNIkiRpdIY3SZKkHjG8SZIk9YjhTZIkqUcMb5IkST1ieJMkSeqRacNbksOSPJDkkSSPJ7mq1S9JMpakkiybpN9pSfYmOX+g9rYkX0nyZJInkqxq9dVJtiV5JslNSZa0+qHt9VjbvmpgX1e2+tNJPjDbgZAkSeqFqppyAQIc0dYPAbYBG4B3A6uA54BlE/osAr4K3AmcP1C/DzizrR8B/Hhbvxm4oK1fDVzU1i8Grm7rFwA3tfUTgUeAQ4HVwLPAommOo1xcXFxcXFxcerL8v2GZZtorb9XZ014e0paqqoer6rkh3T4OfB54ZbyQ5ERgcVVtbfvdU1U/SBJgI3BLa3o9cG5bP6e9pm3/hdb+HODGqnqtqv4aGAPWT3cskiRJPfH8sA0j3fOWZFGSHXRhbGtVbZui7QrgPLoraIPeBexO8oUkDyf5vSSLgGOA3VX1emu3C1jR1lcALwC07a+29j+sT9Jn8LNsTrI9yfZRjlOSJGm+Gym8VdXeqloLrATWJzlpiuafBq6oqr0T6ouB04HLgdOAtwO/Qjct+yNv2X4O2zZVn8HPfU1VrauqdVN8XkmSpN5YPJPGVbU7yX3AWcBjQ5qtA27sZjdZBpyd5HW6q2MPV9W3AJLcSnfv3J8CS5MsblfXVgIvtn3tAk4AdiVZDLwF+N5AfdxgH0mSpAVrlKdNj02ytK0fDpwBPDWsfVWtrqpVVbWK7j61i6vqVuBB4Kgkx7amG4Enqnua4F5g/KnUTcBtbf329pq2/aut/e3ABe1p1NXAGuCBEY9ZkiSpt0aZNl0O3JtkJ10A21pVdyS5NMkuuqteO5NcO9VO2jTq5cA9SR6lm/r8XNt8BXBZkjG6e9qua/XrgGNa/TJgS9vX43RPqD4BfBn4tUmmaSVJkhactD+jseAlOTgOVJIkLQQPDbtn329YkCRJ6hHDmyRJUo8Y3iRJknrE8CZJktQjhjdJkqQeMbxJkiT1iOFNkiSpRwxvkiRJPWJ4kyRJ6hHDmyRJUo8Y3iRJknrE8CZJktQjhjdJkqQeMbxJkiT1iOFNkiSpRwxvkiRJPWJ4kyRJ6hHDmyRJUo8Y3iRJknrE8CZJktQjhjdJkqQeMbxJkiT1iOFNkiSpRwxvkiRJPTJteEtyWJIHkjyS5PEkV7X6JUnGklSSZZP0Oy3J3iTnT6gfmeTbSf5ooHZqkkfb/j6TJK1+dJKtSZ5pP49q9bR2Y0l2JjlltgMhSZLUB6NceXsN2FhVJwNrgbOSbAC+DpwBPD+xQ5JFwKeAuyfZ3yeAr02ofRbYDKxpy1mtvgW4p6rWAPe01wAfHGi7ufWXJEla8KYNb9XZ014e0paqqoer6rkh3T4OfB54ZbCY5FTgeOArA7XlwJFV9Y2qKuAG4Ny2+Rzg+rZ+/YT6De2z3Q8sbfuRJEla0Ea65y3JoiQ76MLY1qraNkXbFcB5wNUT6j8G/D7w6xO6rAB2Dbze1WoAx1fVSwDt53EDfV4Y0mfwPTcn2Z5k+9RHKEmS1A8jhbeq2ltVa4GVwPokJ03R/NPAFVW1d0L9YuDOqnphQj2TveU0H2mkPlV1TVWtq6p10+xPkiSpFxbPpHFV7U5yH909aY8NabYOuLE9c7AMODvJ68B7gNOTXAwcASxJsgf4Q7pQOG4l8GJbfznJ8qp6qU2Ljk/D7gJOGNJHkiRpwRrladNjkyxt64fTPaTw1LD2VbW6qlZV1SrgFuDiqrq1qn65qt7W6pfT3bO2pU2Hfj/JhvaU6YXAbW13twOb2vqmCfUL21OnG4BXx6dXJUmSFrJRpk2XA/cm2Qk8SHfP2x1JLk2yi+6q184k187ic1wEXAuMAc8Cd7X6J4EzkzwDnNleA9wJfKu1/xzdlKwkSdKCl+4Bz4UvycFxoJIkaSF4aNg9+37DgiRJUo8Y3iRJknrE8CZJktQjhjdJkqQeMbxJkiT1iOFNkiSpRwxvkiRJPWJ4kyRJ6hHDmyRJUo8Y3iRJknrE8CZJktQjhjdJkqQeMbxJkiT1iOFNkiSpRwxvkiRJ80hVTbnd8CZJktQjhjdJkqQeMbxJkiTNE9NNmYLhTZIkqVcMb5IkST1ieJMkSZoHRpgxBQxvkiRJ88Ro6c3wJkmS1CPThrckhyV5IMkjSR5PclWrX5JkLEklWTZJv9OS7E1yfnu9Nsk32j52JvnwQNvVSbYleSbJTUmWtPqh7fVY275qoM+Vrf50kg/MfigkSZLmv1GuvL0GbKyqk4G1wFlJNgBfB84Anp/YIcki4FPA3QPlHwAXVtVPA2cBn06ytG37FPAHVbUG+FvgY63+MeBvq+qdwB+0diQ5EbgAGN/Xf2vvKUmS1Duj3u8GI4S36uxpLw9pS1XVw1X13JBuHwc+D7wysJ+/qqpn2vqLbduxSQJsBG5pTa8Hzm3r57TXtO2/0NqfA9xYVa9V1V8DY8D6EY5XkiRpHho9vY10z1uSRUl20AWurVW1bYq2K4DzgKunaLMeWAI8CxwD7K6q19vmXcCKtr4CeAGgbX+1tf9hfZI+g++zOcn2JNtHOU5JkqT5bqTwVlV7q2otsBJYn+SkKZp/GriiqvZOtjHJcuDPgY9U1T8Cmewtx5sP2TZVn8HPfU1VrauqdVN8XkmSpANmJlOmAItntvPaneQ+uvvMHhvSbB1wYze7yTLg7CSvV9WtSY4EvgT8p6q6v7X/DrA0yeJ2dW0l8GLbtgs4AdiVZDHwFuB7A/Vxg30kSZJ6ZGbpbZSnTY8df7AgyeF0Dyk8NfTtq1ZX1aqqWkV3n9rFLbgtAb4I3FBVfzHQvoB7gfNbaRNwW1u/vb2mbf9qa387cEF7GnU1sAZ4YMRjliRJ6q1Rpk2XA/cm2Qk8SHfP2x1JLk2yi+6q184k106zn38H/BzwK0l2tGVt23YFcFmSMbp72q5r9euAY1r9MmALQFU9DtwMPAF8Gfi1YdO0kiRJ89VMp0wBMsq31y8ESQ6OA5UkSb0xLIcleWjYPft+w4IkSVKPGN4kSZJ6xPAmSZLUI4Y3SZKkA2BfHzswvEmSJB0Q+5beDG+SJEk9YniTJEnqEcObJElSjxjeJEmS9rPZfEeC4U2SJGm/2/f0ZniTJEnqEcObJElSjxjeJEmS9qNhX0Y/KsObJElSjxjeJEmSesTwJkmStJ/MdsoUDG+SJEm9YniTJEnaD+bgohtgeJMkSdpP5ia9Gd4kSZJ6xPAmSZL0JpuLBxXGGd4kSZJ6xPAmSZLUI9OGtySHJXkgySNJHk9yVatfkmQsSSVZNkm/05LsTXL+QG1TkmfasmmgfmqSR9v+PpMkrX50kq2t/dYkR7V6WruxJDuTnDIXgyFJkjTX5nLKFEa78vYasLGqTgbWAmcl2QB8HTgDeH5ihySLgE8Bdw/UjgZ+G/gXwHrgt8fDGPBZYDOwpi1ntfoW4J6qWgPc014DfHCg7ebWX5IkacGbNrxVZ097eUhbqqoerqrnhnT7OPB54JWB2geArVX1var6W2ArXRBcDhxZVd+oLpreAJzb+pwDXN/Wr59Qv6F9tvuBpW0/kiRJC9pI97wlWZRkB10Y21pV26ZouwI4D7h6wqYVwAsDr3e12oq2PrEOcHxVvQTQfh43zb4mfpbNSbYn2T71EUqSJPXDSOGtqvZW1VpgJbA+yUlTNP80cEVV7Z1Qz2S7nqI+lZH6VNU1VbWuqtZNsz9JkqQ5N9f3uwEsnuEH2J3kPrp70h4b0mwdcGN75mAZcHaS1+mujr1voN1K4L5WXzmh/mJbfznJ8qp6qU2Ljk/D7gJOGNJHkiRpwRrladNjkyxt64fTPaTw1LD2VbW6qlZV1SrgFuDiqrqV7uGF9yc5qj2o8H7g7jYd+v0kG9pTphcCt7Xd3Q6MP5W6aUL9wvbU6Qbg1fHpVUmSpIVslGnT5cC9SXYCD9Ld83ZHkkuTjF8125nk2ql2UlXfAz7R9vEg8J9bDeAi4FpgDHgWuKvVPwmcmeQZ4Mz2GuBO4Fut/eeAi0c5WEmSpP3lTZgxBSBvxlzsfJTk4DhQSZI0L8wmYyV5aNg9+37DgiRJUo8Y3iRJknrE8CZJktQjhjdJkqQ59mY+UmB4kyRJmnNvXnozvEmSJPWI4U2SJKlHDG+SJEk9YniTJEmaQ2/2FyAY3iRJknrE8CZJktQjhjdJkqQ5sj++M97wJkmS1COGN0mSpB4xvEmSJM2B/TFlCoY3SZKkXjG8SZIk9YjhTZIkaZb215QpGN4kSZJ6xfAmSZLUI4Y3SZKkWdifU6ZgeJMkSeoVw5skSVKPGN4kSZJ6ZNrwluSwJA8keSTJ40muavVLkowlqSTLBtqfk2Rnkh1Jtid578C23237eDLJZ5Kk1U9N8mjb32D96CRbkzzTfh7V6mntxtp7nTLXAyNJkjSd/X2/G4x25e01YGNVnQysBc5KsgH4OnAG8PyE9vcAJ1fVWuCjwLUASX4W+JfAzwAnAacBP9/6fBbYDKxpy1mtvgW4p6rWtP1uafUPDrTd3PpLkiQteNOGt+rsaS8PaUtV1cNV9dwk7ffUGzH0J4Dx9QIOA5YAh7b9vJxkOXBkVX2j9bsBOLf1OQe4vq1fP6F+Q/ts9wNL234kSZIWtJHueUuyKMkO4BVga1Vtm6b9eUmeAr5Ed/WNqvoGcC/wUlvurqongRXAroHuu1oN4Piqeqn1fwk4rtVXAC8M6TP4OTa3qdvtoxynJEnSyPb/jCkwYnirqr1tGnQlsD7JSdO0/2JV/RTdlbJPACR5J/DP2z5WABuT/ByQyXYxzUcaqU9VXVNV66pq3TT7kyRJmpE6QOltRk+bVtVu4D7euCdtuvZ/CbyjPdBwHnB/m1bdA9wFbKC7arZyoNtK4MW2Pj6tSvv5SqvvAk4Y0keSJGnBGuVp02OTLG3rh9M9pPDUFO3fOfC06Cl097h9F/gb4OeTLE5yCN3DCk+26dDvJ9nQ+l0I3NZ2dzuwqa1vmlC/sD11ugF4dXx6VZIkaSFbPEKb5cD1SRbRhb2bq+qOJJcCvwG8FdiZ5M6q+lXgQ3TB6h+Avwc+XFWV5BZgI/Ao3RTnl6vqf7f3uAj4M+Bwuityd7X6J4Gbk3yMLvz9YqvfCZwNjAE/AD6yrwMgSZI0UwfgL4T8UA7E3yc5EJIcHAcqSZLedG92fkry0LB79v2GBUmSpB4xvEmSJM3AgZ60NLxJkiTNyIFNb4Y3SZKkHjG8SZIk9YjhTZIkaUTz4a90GN4kSZJ6xPAmSZLUI4Y3SZKkEcyHKVMwvEmSJPWK4U2SJKlHDG+SJEnTmC9TpmB4kyRJ6hXDmyRJUo8Y3iRJknrE8CZJkjSFeXS7G2B4kyRJmsb8Sm+GN0mSpB4xvEmSJPWI4U2SJGmI+fT33cYZ3iRJknrE8CZJktQjhjdJkqQemTa8JTksyQNJHknyeJKrWv2SJGNJKsmygfbnJNmZZEeS7UneO7DtbUm+kuTJJE8kWdXqq5NsS/JMkpuSLGn1Q9vrsbZ91cC+rmz1p5N8YK4GRJIkCebn/W5A98GmWoAAR7T1Q4BtwAbg3cAq4Dlg2UD7I4C09Z8BnhrYdh9w5kC7H2/rNwMXtPWrgYva+sXA1W39AuCmtn4i8AhwKLAaeBZYNM1xlIuLi4uLi4vLqMuBBGyfLM9U1fRX3to+9rSXh7SlqurhqnpukvZ72psC/EQbAJKcCCyuqq0D7X6QJMBG4JbW53rg3LZ+TntN2/4Lrf05wI1V9VpV/TUwBqyf7lgkSZJGUtM3OVBGuuctyaIkO4BXgK1VtW2a9ucleQr4EvDRVn4XsDvJF5I8nOT3kiwCjgF2V9Xrrd0uYEVbXwG8ANC2v9ra/7A+SZ/Bz7G5Td1uH+U4JUmSAGoep7eRwltV7a2qtcBKYH2Sk6Zp/8Wq+im6K2ifaOXFwOnA5cBpwNuBX6Gblv2RXbSfw7ZN1Wfwc1xTVeuqat1Un1eSJKkvZvS0aVXtprtv7awR2/8l8I72QMMu4OGq+la7inYrcArwHWBpksWt20rgxba+CzgBoG1/C/C9wfokfSRJkvbd/L3oBoz2tOmxSZa29cOBM4Cnpmj/znZfGklOAZYA3wUeBI5KcmxruhF4ot0fdy9wfqtvAm5r67e317TtX23tbwcuaE+jrgbWAA+MdsiSJEnDzecpUxjtytty4N4kO+kC2NaquiPJpUl20V312pnk2tb+Q8Bj7R65PwY+3B562Es3ZXpPkkfppj4/1/pcAVyWZIzunrbrWv064JhWvwzYAlBVj9M9ofoE8GXg19r+JUmSFrS88WDowpbk4DhQSZI0K/MhGyV5aNg9+37DgiRJUo8Y3iRJknrE8CZJktTMhynT6RjeJEmSesTwJkmS1COGN0mSJKAHM6aA4U2SJKnpR3ozvEmSJPWI4U2SJKlHDG+SJEn9mDEFDG+SJEnz/svoBxneJEmSesTwJkmS1COGN0mSpB4xvEmSpINaH77PdJDhTZIkqUcMb5IkST1ieJMkSQetns2YAoY3SZJ0UOtfejO8SZIk9YjhTZIkqUcMb5Ik6aDUx/vdwPAmSZIOWv1Mb9OGtySHJXkgySNJHk9yVatfkmQsSSVZNtD+nCQ7k+xIsj3Jeyfs78gk307yRwO1U5M82vb3mSRp9aOTbE3yTPt5VKuntRtr73XKXA2IJEnSfDbKlbfXgI1VdTKwFjgryQbg68AZwPMT2t8DnFxVa4GPAtdO2P4J4GsTap8FNgNr2nJWq28B7qmqNW2/W1r9gwNtN7f+kiRJC9604a06e9rLQ9pSVfVwVT03Sfs99cb3TPwEA9ckk5wKHA98ZaC2HDiyqr7R+t0AnNs2nwNc39avn1C/oX22+4GlbT+SJEnT6ttXYg0a6Z63JIuS7ABeAbZW1bZp2p+X5CngS3RX30jyY8DvA78+ofkKYNfA612tBnB8Vb0E0H4eN9DnhSF9Bj/H5jZ1u336o5QkSZr/RgpvVbW3TYOuBNYnOWma9l+sqp+iu1L2iVa+GLizql6Y0DyT7WKajzRSn6q6pqrWVdW6afYnSZIOEn2+6gaweCaNq2p3kvvo7kl7bIT2f5nkHe2BhvcApye5GDgCWJJkD/CHdKFw3Ergxbb+cpLlVfVSmxZ9pdV3AScM6SNJkrRgjfK06bFJlrb1w+keUnhqivbvHHha9BRgCfDdqvrlqnpbVa0CLqe7Z21Lmw79fpINrd+FwG1td7cDm9r6pgn1C9tTpxuAV8enVyVJkhayUa68LQeuT7KILuzdXFV3JLkU+A3grcDOJHdW1a8CH6ILVv8A/D3w4Zr++uRFwJ8BhwN3tQXgk8DNST4G/A3wi61+J3A2MAb8APjIKAcrSZIObn2fMgXIQjiIUSQ5OA5UkiQN1Zfck+ShYffs+w0LkiRJPWJ4kyRJB4W+XHWbjuFNkiSpRwxvkiRJPWJ4kyRJC95CmTIFw5skSVKvGN4kSdKCtpCuuoHhTZIkqVcMb5IkST1ieJMkSeoRw5skSVKPGN4kSZJ6xPAmSZIWrIX2pCkY3iRJknrF8CZJktQjhjdJkqQeMbxJkiT1yOID/QH2o+8Azx/oDzHEMrrPp33nGM6eYzh7juHsOYaz5xgOSDLTLvNl/H5y2IYsxKcw+ibJ9qpad6A/R585hrPnGM6eYzh7juHsOYaz04fxc9pUkiSpRwxvkiRJPWJ4mx+uOdAfYAFwDGfPMZw9x3D2HMPZcwxnZ96Pn/e8SZIk9YhX3iRJknrE8CZJktQjhrdZSHJCknuTPJnk8ST/YcL2y5NUkmXtdZJ8JslYkp1JThlouynJM23ZNFA/Ncmjrc9n0v5gTZKjk2xt7bcmOWp/Hfdc2ocxfF+SV5PsaMtvDbQ9K8nTbay2DNRXJ9nWxuqmJEta/dD2eqxtX7V/jnpuDRvDJL+T5NsDY3X2QJ8r23E/neQDA3XHcIQxTLIqyd8P1K8e2NeMfmen+nehT6b6XU7y8XZePZ7kdwfqnofNTMfPc/BHTfF7fNPAOD2XZMdAn36eg1Xlso8LsBw4pa3/M+CvgBPb6xOAu+n+MPCyVjsbuAsIsAHY1upHA99qP49q60e1bQ8A72l97gI+2Oq/C2xp61uATx3o8dhPY/g+4I5J9rMIeBZ4O7AEeGRgPzcDF7T1q4GL2vrFwNVt/QLgpgM9HnM5hsDvAJdP0v7ENj6HAqvbuC1yDGc0hquAx4bsa0a/s8P+XejbMsUY/ivg/wCHtm3HeR7Oyfh5Do44hhPa/D7wW30/Bw/4YC+kBbgNOLOt3wKcDDzHG8HjT4BfGmj/dDvZfgn4k4H6n7TacuCpgfoP2433HThhnz7Qx7+fxvB9TB7e3gPcPfD6yraE7i9lL57Yji4YvqetL27tcqDHYK7GkOHB40rgyoHXd7dxcQxHH8NVTPIfzn35nR3278KBHoM5HMObgTM8D+d8/DwHRxzDgdcBXgDW9P0cdNp0jrRLpO8GtiX5N8C3q+qRCc1W0J0443a12lT1XZPUAY6vqpcA2s/j5uRADqARxxDgPUkeSXJXkp9utWFjeAywu6pen1D/J33a9ldb+94aHMNWuqRNg/xp3phan+l56Bj+6BgCrE7ycJKvJTm91fbld3bYuPfWhDF8F3B6m0r6WpLTWjPPwyFGHD/wHBxqkt9jgNOBl6vqmfa6t+eg4W0OJDkC+DzwH4HXgd8EfmuyppPUah/qC84MxvCbwE9W1cnAfwVuHd/FJG2nG8MFNb6DY1hVfwd8FngHsBZ4iW66AOb2PDxYx/Al4G1V9W7gMuB/JjmSfRuPhT6Gi+luB9kA/Dpwc7sHy/NwEjMYP8/BISYZw3G/BPyvwaaTdO/FOWh4m6Ukh9CdJP+jqr5A9w/9auCRJM8BK4FvJnkrXUo/YaD7SuDFaeorJ6kDvJxkefsMy4FX5vbI9p+ZjGFV/V1V7QGoqjuBQ9I9zDBsDL8DLE2yeEKdwT5t+1uA771pB/ommmQMqaqXq2pvVf0j8DlgfWs+0/PQMZwwhlX1WlV9t60/RHd/zLvYt9/ZYePeO5ONId3xfaE6DwD/SPfF356HE8xk/DwHJzdkDMfPi38L3DTQvLfnoOFtFtr/+7kOeLKq/gtAVT1aVcdV1aqqWkX3P+gpVfV/gduBC9uTPRuAV9ul67uB9yc5qk3LvJ9uHv0l4PtJNrT3upBuDp+2r01tfdNAvVdmOoZJ3tr6kGQ93Tn8XeBBYE17EmgJ3Q2jt1d3A8K9wPntLQfHanAMzwe+2tr3ymRj2OrLB5qdBzzW1m8HLmhPR60G1tDd4OwYjjiGSY5Nsqitv51uDL+1j7+zw/5d6JVhY0h3dXxja/MuuhvAv4Pn4T8x0/HzHPxRU4whwBl09wIOTin39xx8s26mOxgW4L10l0V3AjvacvaENs/xxs32Af6Y7v8hPQqsG2j3UWCsLR8ZqK+j+w/Gs8Af8ca3YhwD3AM8034efaDHYz+N4SXA43RP/9wP/OxAu7Ppni56FvjNgfrb6X4hx4C/4I2ntg5rr8fa9rcf6PGYyzEE/rydZzvp/mFZPtDnN9s4PU17Es0xHH0MgQ8NnIffBP71wL5m9Ds71b8LfVqmGMMlwH9vY/JNYKPn4ezHz3Nw9DFs2/4M+PeT9OnlOejXY0mSJPWI06aSJEk9YniTJEnqEcObJElSjxjeJEmSesTwJkmS1COGN0mSpB4xvEmSJPXI/wcGVtgxjz8ksQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from rasterio.plot import show\n",
    "\n",
    "composite_in_array=composite.read([4,3,2])[:, 150:600, 250:1400]# Only read three of the bands and selective heightand width\n",
    "print(composite_in_array.shape)\n",
    "\n",
    "# Lets plot the raster now\n",
    "fig,ax=plt.subplots(figsize=(10,7))\n",
    "show(composite_in_array[:, :, :], ax=ax, transform=composite.transform) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Okay looks good! Our raster dataset is ready!\n",
    "Now our goal is to get the pixels from the raster as outlined in each shapefile.\n",
    "\n",
    "Our training data, the shapefile we've worked with, contains one main field we care about:\n",
    "\n",
    "   a Classname field (String datatype)\n",
    "\n",
    "Combined with the innate location information of polygons in a Shapefile, we have all that we need to use for pairing labels with the information in our raster.\n",
    "\n",
    "However, in order to pair up our vector data with our raster pixels, we will need a way of co-aligning the datasets in space.\n",
    "\n",
    "We'll do this using the rasterio mask function which takes in a dataset and a polygon and then outputs a numpy array with the pixels in the polygon.\n",
    "\n",
    "Let's run through an example:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CRS.from_dict(init='epsg:32645')"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "composite.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'init': 'epsg:32645'}"
      ]
     },
     "execution_count": 156,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kathmandu_shapefile=gpd.read_file(r\"C:\\Users\\user\\Downloads\\GIS_self_project\\kathmandu_classified_deep_learning.shp\")\n",
    "kathmandu_shapefile.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [],
   "source": [
    "kathmandu_shapefile=kathmandu_shapefile.to_crs({\"init\":\"epsg:32645\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17405"
      ]
     },
     "execution_count": 158,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(kathmandu_shapefile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['OBJECTID', 'Id', 'Classname', 'Shape_Leng', 'Shape_Area', 'geometry'], dtype='object')"
      ]
     },
     "execution_count": 201,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kathmandu_shapefile.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "26"
      ]
     },
     "execution_count": 198,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kathmandu_shapefile[\"Classname\"][14]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to extract the geometry of each feature in the shapefile in GeoJSON format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 1 1 180.000000002 1800.00000003\n",
      "  <shapely.geometry.polygon.Polygon object at 0x0000015101FE20F0>]\n",
      " [2 2 1 330.730401102 4658.2611777\n",
      "  <shapely.geometry.polygon.Polygon object at 0x0000015101FE2E80>]\n",
      " [3 3 1 119.999999996 899.999999944\n",
      "  <shapely.geometry.polygon.Polygon object at 0x0000015101FE2E10>]\n",
      " ...\n",
      " [17403 17403 1 120.0 900.0\n",
      "  <shapely.geometry.polygon.Polygon object at 0x00000151024DD748>]\n",
      " [17404 17404 1 807.749932817 14816.5819543\n",
      "  <shapely.geometry.polygon.Polygon object at 0x00000151024DD630>]\n",
      " [17405 17405 26 30004.513804 1309974.39518\n",
      "  <shapely.geometry.polygon.Polygon object at 0x00000151024DDF98>]]\n"
     ]
    }
   ],
   "source": [
    "print(kathmandu_shapefile.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<GeometryArray>\n",
      "[<shapely.geometry.polygon.Polygon object at 0x0000015101FE20F0>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x0000015101FE2E80>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x0000015101FE2E10>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x0000015101FE2F60>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x0000015101FE28D0>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x0000015101FE2EF0>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x0000015101FE2048>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x0000015101FE2A90>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x0000015101FE8160>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x0000015101FE8240>,\n",
      " ...\n",
      " <shapely.geometry.polygon.Polygon object at 0x00000151024DD4E0>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x00000151024DD518>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x00000151024DD550>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x00000151024DD588>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x00000151024DD5C0>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x00000151024DD128>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x00000151024DD668>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x00000151024DD748>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x00000151024DD630>,\n",
      " <shapely.geometry.polygon.Polygon object at 0x00000151024DDF98>]\n",
      "Length: 17405, dtype: geometry\n"
     ]
    }
   ],
   "source": [
    "# this generates a list of shapely geometries\n",
    "geoms=kathmandu_shapefile.geometry.values\n",
    "print(geoms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'shapely.geometry.polygon.Polygon'>\n",
      "<class 'list'>\n",
      "[{'type': 'Polygon', 'coordinates': (((347181.09350611316, 3077934.320210689), (347181.0935061127, 3077964.3202106897), (347241.0935061133, 3077964.3202106897), (347241.0935061136, 3077934.320210689), (347181.09350611316, 3077934.320210689)),)}]\n"
     ]
    }
   ],
   "source": [
    "# let's grab a single shapely geometry to check\n",
    "geometry=geoms[0]\n",
    "print(type(geometry))\n",
    "\n",
    "# transform to GeoJSON format\n",
    "from shapely.geometry import mapping\n",
    "features=[mapping(geometry)]          # can also do this using polygon.__geo_interface__\n",
    "print(type(features))\n",
    "print(features)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's extract the raster values values within the polygon using the rasterio mask() function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(8, 2, 3)\n",
      "<class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "source": [
    "out_raster,out_transform=mask(composite,features,crop=True)\n",
    "print(out_raster.shape)\n",
    "print(type(out_raster))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Okay those looks like the right dimensions for our training data. 8 bands and 6x8 pixels seems reasonable given our earlier explorations.\n",
    "\n",
    "We'll be doing a lot of memory intensive work so let's clean up and close this dataset.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "metadata": {},
   "outputs": [],
   "source": [
    "composite.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building the Training Data for scikit-learn\n",
    "\n",
    "Now let's do it for all features in the shapefile and create an array X that has all the pixels and an array y that has all the training labels.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8\n",
      "(8, 66, 76)\n",
      "(8, 1469)\n",
      "(8, 1469)\n",
      "(1469, 8)\n"
     ]
    }
   ],
   "source": [
    "X = np.array([], dtype=np.int8).reshape(0,8) # pixels for training\n",
    "y = np.array([], dtype=np.string_) # labels for training\n",
    "with rasterio.open(image_folder_path) as composite_1:\n",
    "    band_count = composite_1.count\n",
    "    print(band_count)\n",
    "    for index, geom in enumerate(geoms):\n",
    "        feature = [mapping(geom)]\n",
    "        \n",
    "    # the mask function returns an array of the raster pixels within this feature\n",
    "        out_image, out_transform = mask(composite_1, feature, crop=True) \n",
    "    # eliminate all the pixels with 0 values for all 8 bands - AKA not actually part of the shapefile\n",
    "        out_image_trimmed_1 = out_image[:,~np.all(out_image == 0, axis=0)]\n",
    "        \n",
    "    # eliminate all the pixels with 255 values for all 8 bands - AKA not actually part of the shapefile\n",
    "        out_image_trimmed_2 = out_image_trimmed_1[:,~np.all(out_image_trimmed_1 == 255, axis=0)]\n",
    "        \n",
    "    # reshape the array to [pixel count, bands]\n",
    "        out_image_reshaped = out_image_trimmed_2.reshape(-1, band_count)\n",
    "\n",
    "    #append the labels to the y array\n",
    "        y = np.append(y,[kathmandu_shapefile[\"Classname\"][index]] * out_image_reshaped.shape[0])  \n",
    "    # stack the pizels onto the pixel array\n",
    "        X = np.vstack((X,out_image_reshaped))     \n",
    "        \n",
    "print(out_image.shape)\n",
    "print(out_image_trimmed_1.shape)\n",
    "print(out_image_trimmed_2.shape)\n",
    "print(out_image_reshaped.shape)        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 220,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 9033  8993  8407 ...  8201  7709  7859]\n",
      " [13549 14019 12377 ... 10279  7432  7995]\n",
      " [ 9242  9131  9148 ...  9056  8464  8371]\n",
      " ...\n",
      " [ 8568  8342  8237 ...  8440  8286  8393]\n",
      " [ 8396  8493  8428 ...  8227  8501  8398]\n",
      " [ 8483  8170  8099 ...  8325  8315  8400]]\n",
      "[b'1' b'12' b'12' b'1' b'1' b'1' b'1' b'1' b'1' b'1' b'1' b'1' b'1' b'1'\n",
      " b'1' b'1' b'1' b'1' b'1' b'1' b'1' b'1' b'1' b'1' b'1' b'1' b'26' b'1'\n",
      " b'1' b'1' b'1' b'1' b'1' b'1' b'26' b'1' b'26' b'1' b'1' b'1']\n"
     ]
    }
   ],
   "source": [
    "print(X)\n",
    "print(y[10:50])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 223,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels=np.unique(kathmandu_shapefile[\"Classname\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 225,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n"
     ]
    }
   ],
   "source": [
    "print(labels.size\n",
    "     )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 226,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The training data include 3 classes: [ 1 12 26]\n",
      "\n",
      "Our X matrix is sized: (459360, 8)\n",
      "Our y array is sized: (459360,)\n"
     ]
    }
   ],
   "source": [
    "# What are our classification labels?\n",
    "labels = np.unique(kathmandu_shapefile[\"Classname\"])\n",
    "print('The training data include {n} classes: {classes}\\n'.format(n=labels.size, \n",
    "                                                                classes=labels))\n",
    "\n",
    "# We will need a \"X\" matrix containing our features, and a \"y\" array containing our labels\n",
    "print('Our X matrix is sized: {sz}'.format(sz=X.shape))\n",
    "print('Our y array is sized: {sz}'.format(sz=y.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "It all looks good! Let's explore the spectral signatures of each class now to make sure they're actually separable since all we're going by in this classification is pixel values.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 231,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[b'1' b'12' b'26']\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Band Intensities Higher Ref Subset')"
      ]
     },
     "execution_count": 231,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABJ8AAAHwCAYAAAAIF85ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeZwdVZn4/8/TS5bOviHZSFgCWRABA4JBDDCMgCyiIAQRFB2WrwzgiKwiMAOjoj9H0UFFg4BIIMqOoCwDIghiFMFsQGRJmkVCQtbO0sv5/VHVye1Od9KB3L6d7s/79bqvrnvqVN3n1u2u0/epc05FSglJkiRJkiSpGMpKHYAkSZIkSZI6L5NPkiRJkiRJKhqTT5IkSZIkSSoak0+SJEmSJEkqGpNPkiRJkiRJKhqTT5IkSZIkSSoak08quYiYHBHVpY5jc0TE/RFx8kbW/zgiLmnPmFrT/PhGxCsR8S+ljKmtOtJxlLTleN4XQERMiogXI2JFRHxiC+/70Yj44pbcp6Tis33YcjbndSPi+oi4otgxvVsRcXRELMjbiz228L63mu9GWzuTT2pR/ke4Kv8DfycifhMRI0sUS4qIndpYt11OHimlQ1NKN+Sv+bmIeLzZ+tNTSv+1JV8zbxTW5p9J4+O4Lfka+et8OCL+LyKWR8TSiLgnIsZv6ddpi2IcR0kt87y/cSU873fkLwOF7dLiiHgwIsZuxi7+E/hhSql3SunOFva/X0T8MW+LFkfEExGx15Z7B5unpc9d6gpsHzauo7QPETE6Pz4VxXrddyu/ILA6/x16OyJuj4ihm7GL7wBn5u3FMy3s/6iI+FtELMv3/3BEjN5S8W+uiLgsIm4q1et3VCaftDFHpJR6A0OBfwI/KHE8gqvyk27j49YtufOI2Bd4ALgLGAZsDzwLPBERO2zh14qI8BwkdSye99Wixi8zLbgq/50ZDrwGTN2M3Y4CZrXyen2Be8l+Bwfm+78cWLMZ+5e05dg+qE0ioryVVWfmv0M7Ab3JEkpttbH2YifgRuArQD+y7y/XAA2bsX+1A7/4aZNSSquBXwPrer9ExMcj4pk8u7wgIi4rWNeYdT85Iubn2eeLC9b3zLP170TEbKDNVzHzLPL0iLgx75kzKyIm5ut+AWwH3JNn1c/Ly/fJr5wuiYhnI2Jywf4ejYj/yq+mLo+IByJicL6uR0TcFBGL8m3/HBHvK9juixExDvgxsG/+mkvy9U2uRkTE4Xk2fkkey24F686PiNfy138+Ig5q6/Eo2EeTq0Dv4Wr5VcCNKaXvp5SWp5QWp5S+BjwFXJbve05EHF7wWhX5Z7xn/nxTx/vKiHgCqAEuiogZzd7LlyPi7pbeR2vHMSI+HxH3FNSbFxHTC54viIjd38XxkLokz/tbxXn/w3l8S/OfH87LD4iIvxfUeygini54/njkQ9wiYlhE3BYRCyPi5Yg4q9lx/3V+PJYBn9tYPCmlVcB0oMm5NiJOyduNdyLidxExKi//B7AD6z+77s12uXO+32kppfqU0qqU0gMppecK4rup4HWaXPHP7RgRT+fH6K6IGJjX3djn3C8ipkbEG/lndEVElLf2uUtdje1Dx28fCvbV/HXPy89tr+fxNu9FNiCyXm3LI+JPEbFjwbZjI+vdujiP69PNXudHEXFfRKwEDthYXCmlJcCdFLQXEVEWERdExD/yYzw9IgZGRPeIWAGUA8/mbUdzuwMvp5QeTpnlKaXbUkrzWzkOk2PDoZ17RcTs/Pfw5xHRI687OCLuzT+rxRHxh8gvnkcrbWhEHAJcBByX/x48u7Hj0ZWYfNImRUQVcBxZAqLRSuAkoD/wceCM2HC+hv2AXYCDgK/nJ2SAS4Ed88fHgFbHSLfiSOCW/LXvBn4IkFL6LDCf/MpMSumqiBgO/Aa4guzK6bnAbRExpGB/JwCfB7YBuuV1yOPqB4wEBgGnA6sKA0kpzcnLn8xfs3/zYCNLylwHnJbv5yfA3fnJdBfgTGCvlFKf/Hi8spnHY4vIP+cPA79qYfV04OB8eRowpWDdx4C3U0p/bePx/ixwKtCH7KrZLhExpmD9CcDNLcTX6nEEfg98JG+4hgKVwKR8ux3Irq4818ZDIXV5nvc79nk/siTKb4Cr8/1/F/hNRAwCngR2yv9hrgB2BUZERJ+I6Al8EGj85/kest6tw8k+s3Mi4mMFL3UU2ZfM/sAvNxFTL7K2YV5B2SfI/gH/JDAE+ANZG0JKaUeafnbNezS9ANRHxA0RcWhEDNicY5Q7CTiFrCdvHdnxgo1/zjfkdXcC9gD+FfhiWz53qSuwfejY7UNr8oTIfwD/QnZ++2gL1aaQ9TAdQHYuvzLfthfwINn/59vk9a6JiAkF256Q1+8DbHR4ct5WfZKC9gI4C/hEHtcw4B3gf1NKa/LeUgAfyNuO5v4KjI2I/4nsAkzvFupsymfIjveOZBc/vpaXfwWoJmvD3kfWpqWNtaEppd8C/w3cmv8efOBdxNMpmXzSxtyZZ+yXkSUevt24IqX0aErp7ymlhvwq5DQ2PIldnl+pfJbsD7PxD+/TwJV5r5oFrP9nsK0eTyndl1KqB35RsN+WnAjcl9dvSCk9CMwADiuo8/OU0gstXLWtJWsUdsqvuv4lpbRsM2MF+DfgJymlP+X7uYFs2MA+QD3QHRgfEZUppVdSSi1l9Budm2fel0TE2+8ilo0ZSHZOeKOFdW8Ag/Plm4Ej838+oGmyqC3H+/qU0qyUUl1KaSnZEL8pAHkSaizZPw/NtXocU0ovAcvJPruPAr8DXots7pGPAn9IKdn1Vto0z/sd77zfko8DL6aUfpGfS6cBc8m+ZK3O3+/+wESyxPvjZAn5ffLtFpH1LhiSUvrPlNLa/Dz6U+D4gtd5MqV0Z34cm3zJKnBu/juznOzL5WcL1p0GfCOlNCelVEf2z/jukfd+2pj8uO8HpDyuhRFxd+Q9DdroFymlmSmllcAlwKcjGw7S4uec7/tQ4JyU0sqU0lvA/9D0mEhdle1Dx2sfCr8XLGHjF1o/nb+3WSmlGrIkU3O3p5Sezs/Xv2T9ez8ceCWl9PO8zfkrcBtwTMG2d6WUnsiP6+pWYrg6IpYCb5N9r/j3gnWnARenlKpTdjHiMuCYaH3I9zp5+zWZLAk0HXg77+20OUmoH6aUFqSUFpMl0RovtNeSDTUdlVKqTSn9IaWUaFsbqmZMPmljPpGyjH13siz87yNiW4CI+FBEPJJ3M1xKluUf3Gz7NwuWa8h6n0CWzV5QsO7VzYyr+X57bOTENAo4ttmJeT+yk8im4vwFWRLjlsi6p14VEZWbGWtjDF9pFsNIYFhKaR5wDtkJ9q2IuCUihm1kX99JKfXPH82P93v1DtnY6JYm/xtK1lCQxzwHOCJPQB3J+uRTW4534WdPvm3jCf4E4M68UWyu1eOYr/89WcOzf778KNk/Ph/Nn0vaNM/7He+835JhbHgMXyX7xxvadj4cBQxrFuNFZFd2GzU/X7fkO/nvzGiyXgC7FKwbBXy/YP+LgSiIc6PypNXnUkojyHpwDQO+15ZtW4j/VbJesYNp/XMeldd5oyDmn5Bd6Ze6OtuHjtc+FH4v6A/stpG6zY9zS+f31t77KOBDzWL+DLDtJvbX3FkppX55nAOAEQXrRgF3FOx/Dlkyrk0XHFJKT6WUPp1SGgJ8hKz9u3gTmxVq/jvYeNy/TdZD64GIeCkiLiiId1NtqJox+aRNyrPyt5OdAPbLi28m650yMj+J/JjsH8q2eIPsJNtouy0VK9kV0kILyK589i949EopfXOTO8qy25enlMaTDUc7nKxL8aZes7kFZFd0CmOoyq9Uk1K6OaW0H9lJLAHf2lRsLagBqgqeb9taxdbkV4afBI5tYfWngYcLnjcOvTsKmJ03ltC24938eD0ADI5sTqYptDDkrmDfrR5H1n/Z+ki+/HtMPknviuf9Dn/efz3fttB2ZBN+w4bJp5bOhwvI5sgojLFPSqmwB8Cm3uf6itncGmeTJZt6FrzGac1eo2dK6Y9tfqfr9z8XuJ4sCQXZMJ9NtXvNf+dqyYaJt/Y5LyDrgTC4IN6+KaXGoSVtPh5SZ2X70OHbh9a8QdNkz+bcrXAB8PtmMfdOKZ1RUGdz2ou/kw19/N+IaPw9WQAc2uw1eqSUXmt9T63u/8/A7by39uL1fF/LU0pfSSntABwB/Edk83Btqg21vWiBySdtUmSOIstQz8mL+wCLU0qrI2Jvsh4rbTUduDAiBkTECJp2uXyv/kk2gWmjm8h66HwssglDe0Q2ydyIVrZfJ7Ixw+/Pu+gvI/untb6V1xwREd1a2dVPgdPzq0IREb0im5ixT0TsEhEHRjZv0Wqyq8Ytvcam/A04IX+Ph9DyOO62uAA4OSLOyuMbENkEffvStHvuLWTzYJxB02TRZh/vlHXt/TXZlYWBZGPKW9LqcczX/55sgsOeKaVqsrlFDiHrIr3BLVkltc7zfoc67ze+h8ZHN+A+YOeIOCGymz4cRzb57735Nn8k64G0N/B0SmkW+ZVr4LG8ztPAssgmt+2ZH6tdI6LNk/02l7IhLK+TzesH2RfQCyOfFySyybxbusCxgcgmt/1K4+cW2W3dp7B+npm/AftHxHYR0Q+4sIXdnBgR4yPrpfufwK9TSvWtfc4ppTfILoj8fxHRN7J5BHeMiMY2dVOfu9Tp2T50qPZhc0wHPh8R4/Jz4tc3Y9t7ydqcz0ZEZf7YK9bP2/Vu3EDWq/TI/PmPgStj/U0phuS/Z5sUEftFxL9FxDb587H5fgvbi8Mim8B8W7LeZc19KSJGRDan4kXArfm+Do+InSIiyD73+vyxqTb0n8Do8M7eTXgwtDH3RHZ3gWVkY19Pzv+BBfh/wH9GxHKyk9f0VvbRksvJujO+TPZP3i+2XMh8A/haZN0fz03Z2PGjyE4iC8my1F+lbb/725IlRZaRNa6/J2u0mvs/slt/vhktzMOUUppBNr77h2RD2+ax/q5B3YFvkg1pe5PsJHxRW95oM2eTZeMbu8He+S72QUrpcbLJ9j5JdoXkVbIJV/dLKb1YUO8Nsl5SHyY/Oefl7/Z430w2AeKv8mRUS7Ft7DiSUnoBWEGWdGqcL+Ql4ImUzQMgadM873e88/4FZF9AGh//l7I5mw4nmwh1EXAecHhKqXF49EqyCVhnpZTW5vt5Eng1ZfMYkZ8XjyC/S1Aez8/IJtR9L74NnBcR3VNKd5Bdtb8lsjvmzSSbU6ktlpMly/4U2d2Tnsq3/0oe/4Nk7c9zwF9Yn3gr9Auy3lJvAj3IJrSFjX/OJ5FNMjyb7LP7NeuH5Gz0c5c6OduHjtc+tFlK6X6y+bQeyV/zyXxV85s9tLTtcrKLzseTXWB4k+zc3vwupZsTz9o8nkvyou+T9Z57IP89eoqsDWiLJWTJpr/nv6O/Be4gu4s3ZL9Tz5JN3v4ABd9dCtycr3spfzTeHW8M8BDZd4wngWtSNsfZptrQxhs4LYqIv7bxfXR6kZI9wiRJkiRJ6gryXkszge6tXfiVtjR7PkmSJEmS1IlFxNER0S0iBpD1XLrHxJPak8knSdJWL5+34emIeDYiZkXE5Xn5QRHx14j4W0Q8HhE7lTpWSVLxRMR1EfFWRMxsYd25EZEiYnD+PCLi6oiYFxHPRcSe7R+x1G5OIxtu+A+yeYvO2Hh1acty2J0kaauXTwTZK6W0IrJbHz9ONhfajcBRKaU5EfH/gL1TSp8rYaiSpCKKiP3J5me5MaW0a0H5SLI5WcYCH0wpvR0Rh5FNcH0Y2fwy308ptXWeGUnSZrDnkyRpq5cyK/Knlfkj5Y++eXk/8lvnSpI6p5TSY8DiFlb9D9nk/IVX3o8iS1KllNJTQP+IGNrCtpKk96ii1AFIkrQl5Lc//guwE/C/KaU/RcQXgfsiYhXZHWr2KWWMkqT2FxFHAq+llJ7NOsquM5zsjmeNqvOyN9oxPEnqErpc8mnw4MFp9OjRpQ5Dkjqcv/zlL2+nlIaUOo53K7/t7e4R0R+4IyJ2Bb4MHJYnor4KfBf4YvNtI+JU4FSAXr16fXDs2LHtGLkkbR22xnYiIqqAi8luFb/B6hbKWpyTxHZCkjZtY+1El0s+jR49mhkzZpQ6DEnqcCLi1VLHsCWklJZExKPAocAHUkp/ylfdCvy2lW2uBa4FmDhxYrKdkKQNbaXtxI7A9kBjr6cRwF8jYm+ynk4jC+qOoJXh2bYTkrRpG2snnPNJkrTVi4gheY8nIqIn8C/AHKBfROycVzs4L5MkdREppb+nlLZJKY1OKY0mSzjtmVJ6E7gbOCm/690+wNKUkkPuJKkIulzPJ0lSpzQUuCGf96kMmJ5Sujci/g24LSIagHeAU0oZpCSpuCJiGjAZGBwR1cClKaWprVS/j+xOd/OAGuDz7RKkJHVBRUs+RcR1wOHAW423OY2IbwNHAGuBfwCfTyktydddCHwBqAfOSin9Li8/BPg+UA78LKX0zbx8e+AWYCDwV+CzKaW1xXo/kqSOK6X0HLBHC+V3AHe0f0SSpFJIKU3ZxPrRBcsJ+FKxY5IkFbfn0/XAD4EbC8oeBC5MKdVFxLeAC4HzI2I8cDwwARgGPFQwTOJ/yYZKVAN/joi7U0qzgW8B/5NSuiUifkyWuPpREd+PJEmSuqDa2lqqq6tZvXp1qUNpFz169GDEiBFUVlaWOhRJUidRtORTSumxiBjdrOyBgqdPAcfky0cBt6SU1gAvR8Q8YO983byU0ksAEXELcFREzAEOBE7I69wAXIbJJ0mSJG1h1dXV9OnTh9GjR5NPWt1ppZRYtGgR1dXVbL/99qUOR5LUSZRywvFTgPvz5eHAgoJ11XlZa+WDgCUppbpm5S2KiFMjYkZEzFi4cOEWCl+SJEldwerVqxk0aFCnTzwBRASDBg3qMr28JEntoyTJp4i4GKgDftlY1EK19C7KW5RSujalNDGlNHHIkCGbG64kSZK6uK6QeGrUld6rJKl9tHvyKSJOJpuI/DP5JH+Q9VwaWVBtBPD6RsrfBvpHREWzckmSJKnTeeWVV9h11103KJ88eTKvvPIKABdffDEjR46kd+/eTepcdtllXH/99e0QpSRJLWvX5FN+57rzgSNTSjUFq+4Gjo+I7vld7MYATwN/BsZExPYR0Y1sUvK786TVI6yfM+pk4K72eh+SJElSR3PEEUfw9NNPlzoMSZI2ULTkU0RMA54EdomI6oj4Atnd7/oAD0bE3/K71JFSmgVMB2YDvwW+lFKqz+d0OhP4HTAHmJ7XhSyJ9R/55OSDgKnFei+SJElSqdXV1XHyySez2267ccwxx1BTU8PAgQMpLy8HYJ999mHo0KEbbNe7d2969uzZ3uFKkrROMe92N6WF4lYTRCmlK4ErWyi/D7ivhfKXWH9HPEmSJKnoLr9nFrNfX7ZF9zl+WF8uPWLCJus9//zzTJ06lUmTJnHKKadwzTXXcPvtt29yu3PPPXdLhClJ0rtWyrvdSZIkSWqjkSNHMmnSJABOPPFEHn/88RJHJElS2xSt55MkSZLU2bSlh1KxNL8LnXelkyRtLez5JEmSJG0F5s+fz5NPPgnAtGnT2G+//UockSRJbWPySZIkSdoKjBs3jhtuuIHddtuNxYsXc8YZZzRZf9555zFixAhqamoYMWIEl112WWkClSSpGYfdSZIkSR3c6NGjmT179kbrXHXVVVx11VXtFJEkSW1n8knqZFbWriQIyqKM8rJyKqLCOSEkSZIkSSVj8knqJFbWruTrT3ydB159YIN1QaxLRBUmpRqXyyN7VJS1vL6luq0tr9tHYXkrdfcZtg8TBpVu4lZJkiRJUvGZfJI6gZeWvsQ5j5zD/GXzOXn8yQzuOZi6VEd9Qz0NqYG6VEdDaqC+oX7dcl1DHfWpft1y8/UtLa+tX0t9Qz31KXus265wuWB98+W6VNck7rJnyjhx3ImcuceZ9KzoWaKjJ0mSJEkqJpNP0lbu4Vcf5uInLqZ7eXeuPfha9h66d6lD2qjGZNXKtSv5wTM/4MbZN/LIgke4/MOXs9e2e5U6PEmSJEnSFubd7qStVH1DPd/7y/c459Fz2KHfDtx6+K0dPvEEUBZlVJZV0r9Hfy7Z9xKm/utUUkqc8rtTuOKpK1hZu7LUIUqSJEmStiCTT9JWaMnqJZzx0BlMnTmVY3Y+husPuZ5te21b6rDelb2H7s1tR97GieNOZPrz0zn6rqN54rUnSh2WJEmSJGkLMfkkbWVmL5rNcfcex4x/zuCyfS/j0n0vpVt5t1KH9Z5UVVZx/t7nc+OhN9KjogenP3Q6lzxxCUvXLC11aJIkdQivvPIKu+666wblkydP5pVXXgHg4osvZuTIkfTu3btJne9+97uMHz+e3XbbjYMOOohXX3113T4nT55c7NAlSTL5JG1N7pp3FyfdfxINNHDjoTfyqZ0/VeqQtqjdt9mdXx3xK774/i9yzz/u4ei7jub/5v9fqcOSJGmrcMQRR/D0009vUL7HHnswY8YMnnvuOY455hjOO++8EkQnSerKTD5JW4Ha+lqueOoKvvbE19h9yO7cevit7Dp4w6ufnUH38u6cvefZ3PzxmxnQYwBnP3I25/3+PBavXlzq0CRJKqm6ujpOPvlkdtttN4455hhqamoYOHAg5eXlAOyzzz4MHTp0g+0OOOAAqqqq1tWprq4GoLy8nIEDB7bfG5AkdVne7U7q4P658p985fdf4dmFz/L5CZ/nrD3PoqKs8//pjh80nls+fgtTZ07lJ8/9hKfeeIoLP3Qhh4w+hIgodXiSpK7q/gvgzb9v2X1u+3449JubrPb8888zdepUJk2axCmnnMI111zD7bffvlkvNXXqVA499FAARo4cudnbS5L0btjzSerAZrw5g+PuPY4X3nmB73z0O/zHxP/oEomnRpXllZz+gdOZfvh0hvceznmPncfZj5zNwpqFpQ5NkqR2N3LkSCZNmgTAiSeeyOOPP75Z2990003MmDGDr371q8UIT5KkVnWdb7HSViSlxC/n/JLvzPgOI/uMZOrHprJj/x1LHVbJjBkwhl8c9gtumn0TP/zbDznqrqM4b6/zOGrHo+wFJUlqX23ooVQszdu8zWkDH3roIa688kp+//vf07179y0dmiRJG2XPJ6mDqamt4YI/XMC3/vwt9h+xPzd//OYunXhqVFFWwed2/Ry/PuLXjOk/hkueuIQzHjqD11e8XurQJElqF/Pnz+fJJ58EYNq0aey3335t2u6ZZ57htNNO4+6772abbbYpZoiSJLXI5JPUgSxYtoAT7z+R+1++n3/f49/53gHfo0+3PqUOq0MZ3W80Pz/k51z0oYv461t/5ei7jubWubfSkBpKHZokSUU1btw4brjhBnbbbTcWL17MGWec0WT9eeedx4gRI6ipqWHEiBFcdtllAHz1q19lxYoVHHvssey+++4ceeSRJYhektSVOexO6iAeq36MC/5wAUHwo3/5EZOGTyp1SB1WWZQxZewU9h+xP5f98TKu+NMV/PaV33L5hy9nu77blTo8SZK2uNGjRzN79uyN1rnqqqu46qqrNih/6KGHihWWJEltYs8nqcQaUgM/+tuPOPPhMxneezi3Hn6riac2Gt57ONcefC2Xf/hynl/8PJ+6+1PcMOsG6hvqSx2aJEmSJCln8kkqoWVrl3HW/53FNc9ew+E7HM6Nh97IiD4jSh3WViUi+OSYT3LHUXewz9B9+M6M73DS/SfxjyX/KHVokiRJkiRMPkkl88I7L3D8vcfzxGtPcNGHLuLK/a6kZ0XPUoe11Xpfr/dx9YFX862PfIv5y+dz7D3Hcu1z11LbUFvq0CRJkiSpSzP5JJXAfS/dx4n3ncjqutVcd8h1TBk7ZbNul6yWRQSH7XAYdx51JwdudyA/eOYHnPCbE5i7eG6pQ5MkSZKkLsvkk9SOahtquerPV3H+H85n3MBx3Hr4reyxzR6lDqvTGdRzEN/56Hf43uTvsbBmIVPuncLVf72atfVrSx2aJEmSJHU5Jp+kdvL2qrf5twf+jV/M/gWfGfcZfvaxnzGkakipw+rUDhp1EHd94i4O2+Ewfvr3n/Lpez7NcwufK3VYkiRJktSlmHyS2sGzC5/luHuOY9bbs/jv/f6bC/a+gMqyylKH1SX0696PK/e7kmsOuoYVtSv47P2f5dt//jar6laVOjRJktrslVdeYdddd92gfPLkybzyyivU1NTw8Y9/nLFjxzJhwgQuuOCCJvWmT5/O+PHjmTBhAieccMK6fU6ePLk9wpckdXEmn6QiSikx/fnpfO63n6OyvJKbDruJI3Y8otRhdUkfGfER7jzqTj415lPcOPtGPnX3p/jzm38udViSJG0x5557LnPnzuWZZ57hiSee4P777wfgxRdf5Bvf+AZPPPEEs2bN4nvf+16JI5UkdTUVpQ5ga3H7i7fTs6InvSp7UVVRRVVlFb0qe6173rOipxNGq4k19Wu48qkruWPeHUwaPolvfeRb9Over9RhdWm9u/Xm6/t+nUNGH8Klf7yUU353Csftchxf/uCX6VXZq9ThSZK0UXV1dZx88sk888wz7Lzzztx4440MHDiQ8vJyqqqqOOCAAwDo1q0be+65J9XV1QD89Kc/5Utf+hIDBgwAYJtttgGgvLycgQMHlubNSJK6FJNPbVDbUMulf7x0o3WCyBJSFb2aJqYqq6iqaPp8gzoVG9Y3mbV1e33F63z50S8ze9FsTtvtNM74wBmUl5WXOizl9h66N7cdeRs//NsPuWn2TTxW/RiX7nspk4ZPKnVokqQO7ltPf2uL30V17MCxnL/3+Zus9/zzzzN16lQmTZrEKaecwjXXXMPtt9++Qb0lS5Zwzz33cPbZZwPwwgsvADBp0iTq6+u57LLLOOSQQxg5cmSL20uStKWZfGqDiqjgwWMeZGXtSmpqa1hZt3L9cu1KaupqmjxvLKupreGNFW80eb66fnWbXrMsytb1sNogeVXZa10CqzGZVbi+X/d+7NR/J3tylMiTrz/JeY+dR11DHT848AdMHjm51CGpBVWVVZy313n866h/5et//DqnP3Q6n9jpE5w78Vx7qEmSOqSRI0cyaVJ2oeTEE0/k6quv5txzz21Sp5mCEMUAACAASURBVK6ujilTpnDWWWexww47rCt78cUXefTRR6muruYjH/kIM2fOpH///u3+HiRJXZPJpzaICLbtte0W2VddQ926RNS6ZFVBMqt5cqsxsdX4/I0VbzRJdrWWzAqCUX1HMW7gOMYNGsfYgWMZN3Ac/Xv4T0axpJS4buZ1XP3M1ezQbwe+d8D3GNV3VKnD0ibsvs3u/OqIX/HjZ3/Mz2f+nMdfe5xL9rmEA7c7sNShSZI6oLb0UCqW5r3iW+olf+qppzJmzBjOOeecdWUjRoxgn332obKyku23355ddtmFF198kb322qvoMUuSBCaf2l1FWQV9u/Wlb7e+W2R/LSWzFq1exPOLn2fO4jk8u/BZ7n/l/nX1h/UaliWiBo1j/KDxjB04liE9hzjE7z1aWbuSS564hAdffZBDRh/C5R++nKrKqlKHpTbqXt6ds/c8m4NHHcwlT1zC2Y+czaGjD+WCD13AwB7OhSFJ6hjmz5/Pk08+yb777su0adPYb7/9mqz/2te+xtKlS/nZz37WpPwTn/gE06ZN43Of+xxvv/02L7zwwrpeUZIktQeTT1u51pJZhUO9lq5ZypzFc5izaM66n48seIREAmBQj0GMGzRuXS+pcQPHMbz3cBNSbfTS0pc455FzmL9sPudOPJeTxp/ksdtKjR80nls+fgtTZ07lJ8/9hKfeeIoLP3Qhh4w+xM9UklRy48aN44YbbuC0005jzJgxnHHGGevWVVdXc+WVVzJ27Fj23HNPAM4880y++MUv8rGPfYwHHniA8ePHU15ezre//W0GDRpUqrchSeqCIqVU6hja1cSJE9OMGTNKHUbJraxdua53VGNS6qUlL1GX6gDo061PlowaOI6xg8YyfuB4RvUd5aTZzTz86sNc/MTFdC/vznc++h322tbu653Fi++8yNef+DozF81k8sjJnLn7mewycJdSh1VUEfGXlNLEUsdRarYTkpqbM2cO48aNK3UY7aql92w7kbGdkKSWbaydsOdTF9Wrshd7vm9P9nzfnuvK1tSvYd4785i9eDZzF81lzuI5TJs7jbUNawHoWdGTXQbswtiBYxk/aDzjBo1jx347UlleWaq3UTL1DfX84JkfMHXmVN4/+P18d/J3t9i8YOoYxgwYwy8O+wU3zb6Ja569hkcXPMqe2+zJlHFTOGi7g6gs63q/95IkSZL0bph80jrdy7szYfAEJgyesK6stqGWl5e+zNzFc5mzaA6zF83mnpfu4ZbnbwGgsqySnfrvtG7+qHGDxrHzgJ3pWdGzVG+j6N5Z/Q7nP3Y+T77xJMfsfAwX7n0h3cq7lTosFUFFWQWf2/VzHD3maO6cdye3zL2Fr/7+qwzpOYRjdzmWY3c+lsE9B5c6TEmSJEnq0Ew+aaMqyyrZecDO7DxgZ47c8UgAGlIDC5YvaDKH1MPzH+a2F28DoCzK2L7v9k3mkRo7cCx9uvUp5VvZImYvms2XH/kyC1ct5LJ9L+NTO3+q1CGpHfTr3o+TJ5zMieNO5InXn+DmuTdzzd+u4drnruXgUQdzwtgT+MCQDzgvlCRJkiS1wOSTNltZlDGq7yhG9R3FIdsfAkBKiTdXvpklo/KE1NNvPs29L927bruRfUY2mdR87MCxDOq59Ux2eee8O/mvJ/+LgT0HcuOhN7Lr4F1LHZLaWXlZOfuP2J/9R+zPq8te5Za5t3DXvLu4/+X7GTdwHFPGTuHQ7Q+lR0WPUocqSdrCUkpd5iJDV5sTVpJUfE44rqJ6e9XbzF08l7mL5zJ70WzmLJpD9Yrqdeu7lXWjvKycsiijLMooj/KWf5a1Uh7rt23yvGwj+3oX+3x1+av85qXf8KFtP8RVH72KgT0GlvCoqiOpqa3h3pfuZdrcacxbMo9+3fvxyZ0+yad3+TQj+owodXibxYlkM7YTkpp7+eWX6dOnD4MGDer0CaiUEosWLWL58uVsv/32TdbZTmRsJySpZU44rpIZ3HMw+w3fj/2G77eubNnaZTy/+HlmL5rNotWLaGhooD7V05DW/yxcrk/1LdbZ4GdDPXUNda3u593ssyE1AFAe5Xx+wuc5a8+zqCjzz0brVVVW8eldPs2xOx/LjH/OYNrcadw4+0aun3U9Hx3xUaaMncI+w/ahLMpKHaok6V0aMWIE1dXVLFy4sNShtIsePXowYsTWdQFFktSx+S1a7a5vt77ste1e7LXtXqUOZZNSSjSkBhLJpJM2KiLW/V6/ufJNfvXCr/j1C7/m0YceZXTf0Rw/9niO3PHITjH3mSR1NZWVlRv0ApIkSW3npXhpIyKC8rJyE0/aLNv22pZ/3+PfefCYB/nGR75B3+59+ebT3+SgXx3EFU9dwbx35pU6REmSJElqN36jlqQi6VbejcN3OJzDdzicWW/PYtrcadzx4h3c+vyt7L3t3kwZO4XJIyeb3JQkSZLUqdnzSZLawYTBE7hivyt46NiHOGfPc1iwfAFffvTLHHLbIfz0uZ+yaNWiUocoSZIkSUVh8kmS2tGAHgP4wvu/wP2fvJ/vH/B9tu+3PVc/czUH//pgLvrDRfx94d9LHaIkSVutiLguIt6KiJkFZf8VEc9FxN8i4oGIGJaXR0RcHRHz8vV7li5ySercTD5JUgmUl5Vz4HYH8tN//Sl3HXUXx+x8DA/Pf5gT7juBKfdO4e5/3M2a+jWlDlOSpK3N9cAhzcq+nVLaLaW0O3Av8PW8/FBgTP44FfhRewUpSV2NySdJKrEd+u/ARR+6iIePfZiLPnQRK+tWcvHjF3Pwrw7m+3/9Pm+seKPUIUqStFVIKT0GLG5WtqzgaS8g5ctHATemzFNA/4gY2j6RSlLX4iy3ktRB9O7Wmyljp3D8Lsfzpzf/xLQ507hu5nVcN/M6Dhh5AFPGTmHvbfcmIkodqiRJW5WIuBI4CVgKHJAXDwcWFFSrzsu86iNJW5jJJ0nqYCKCfYbuwz5D9+H1Fa8z/fnp3PbibTw8/2F26LcDU8ZO4Ygdj6BXZa9ShypJ0lYhpXQxcHFEXAicCVwKtHQ1J7VQRkScSjY0j+22265YYUpSp+WwO0nqwIb1HsY5HzyHh459iCsmXUGPih5c+acrOehXB/GNP32Dl5e+XOoQJUnamtwMfCpfrgZGFqwbAbze0kYppWtTShNTShOHDBlS5BAlqfMx+SRJW4Hu5d05aqejuOXjt/DLw37JASMPYPoL0znyziM59YFTeWT+I9Q31Jc6TEmSOpyIGFPw9Ehgbr58N3BSfte7fYClKSWH3ElSETjsTpK2IhHBbkN2Y7chu/GViV/hthduY/oL0znrkbMY1msYx409js+O/yyVZZWlDlWSpHYXEdOAycDgiKgmG153WETsAjQArwKn59XvAw4D5gE1wOfbPWBJ6iJMPknSVmpwz8Gc9oHTOOX9p/DogkeZNnca9/zjHj4/wf+dJUldU0ppSgvFU1upm4AvFTciSRKYfJKkrV5lWSUHjzqYg0cdzMrald4NT5IkSVKH4pxPktSJeAc8SZIkSR2NPZ8kSVu9iOgBPAZ0J2vbfp1SujQi/gD0yattAzydUvpEicKUJEmSuiSTT5KkzmANcGBKaUVEVAKPR8T9KaWPNFaIiNuAu0oWoSRJktRFOexOkrTVS5kV+dPK/JEa10dEH+BA4M4ShCdJkiR1aSafJEmdQkSUR8TfgLeAB1NKfypYfTTwcEppWSvbnhoRMyJixsKFC9sjXEmSJKnLMPkkSeoUUkr1KaXdgRHA3hGxa8HqKcC0jWx7bUppYkpp4pAhQ4odqiRJktSlFC35FBHXRcRbETGzoOzYiJgVEQ0RMbGgfHRErIqIv+WPHxes+2BE/D0i5kXE1ZHfQzwiBkbEgxHxYv5zQLHeiyRp65FSWgI8ChwCEBGDgL2B35QwLEmSJKnLKmbPp+vJ//EvMBP4JNkdiZr7R0pp9/xxekH5j4BTgTH5o3GfF5ANoRgDPJw/lyR1QRExJCL658s9gX8B5uarjwXuTSmtLlV8kiRJUldWtORTSukxYHGzsjkppefbuo+IGAr0TSk9mVJKwI1A4y2yjwJuyJdvKCiXJHU9Q4FHIuI54M9kcz7dm687no0MuZMkSZJUXBWlDqDA9hHxDLAM+FpK6Q/AcKC6oE51XgbwvpTSGwAppTciYpvWdhwRp5L1nmK77bYrRuySpBJKKT0H7NHKusntG40kSZKkQh1lwvE3gO1SSnsA/wHcHBF9gWihbmqhbKOcSFaSJEmSJKk0OkTyKaW0JqW0KF/+C/APYGeynk4jCqqOAF7Pl/+ZD8trHJ73VvtFLEmSJEmSpLboEMmnfKLY8nx5B7KJxV/Kh9Utj4h98rvcnQTclW92N3ByvnxyQbkkSZIkSZI6iKIlnyJiGvAksEtEVEfEFyLi6IioBvYFfhMRv8ur7w88FxHPAr8GTk8pNU5WfgbwM2AeWY+o+/PybwIHR8SLwMH5c0mSJEmSJHUgRZtwPKU0pZVVd7RQ9zbgtlb2MwPYtYXyRcBB7yVGSZIkSZIkFVeHGHYnSZIkSZKkzsnkkyRJkiRJkorG5JMkSZIkSZKKxuSTJEmSJEmSisbkkyRJkiRJkorG5JMkSZIkSZKKxuSTJEmSJEmSisbkkyRJkiRJkorG5JMkSZIkSZKKxuSTJEmSJEmSisbkkyRJkiRJkorG5JMkSZIkSZKKxuSTJEmSJEmSisbkkyRJkiRJkorG5JMkSZIkSZKKxuSTJEmSJEmSisbkkyRJkiRJkorG5JMkSZIkSZKKxuSTJEmSJEmSisbkkyRJkiRJkorG5JMkSZIkSZKKxuSTJEmSJEmSisbkkyRJkiRJkorG5JMkSZIkSZKKxuSTJEmSJEmSiqai1AFI2nJWrKnjzaWr6F5RTs9u5VR1K6dHRTllZVHq0CRJkiRJXZTJJ2kr9c7Ktcx6fRkzX1/KrNeXMeu1pby8aCUpbVi3R2UZPSvLqepWQY/KMqq6VdCzMktQZeXl9Chczn82rVNBz25l9KysWJ/YqszWdauwE6UkSZIkqWUmn6QOLqXEW8vXMPO1pcx8bRmz8mTTa0tWraszvH9PJgzryyf2GM6oQVWsqWtg1dp6VtXWN/lZs7ae1bX11KytY1VtPUtq1vJGbUG9tfXU1Na3mMDamIqyWJ+oKkharU9gVdAzT3o1JqyqupUzuE83Rg6oYuTAKob07m4PLUmSJEnqhEw+SR1ISon5i2uyHk2v5T2aXl/K2yvWAhAB2w/uxQdHDeCkfUex6/B+jB/alwG9um3RGAqTV40Jq8blrLyOVWsbqFlblyezsvWrm9TJfr6zsrZJnVVr61lb37DB63avKGPEgJ6MHFjFdgOr8qRU9nzkwCr69qjcYu9RkiRJktR+TD610eKVa6ksD7pVlNGtvIwIe2jovamrb+Clt1cy6/WmPZqWr64Dst5EY97XhwN22YYJw/qy6/B+jBval17di/tnGxH0qMyG1A0o0mvU1Tewqraet5avYf7iGqoX1zB/cQ0LFq9iwTs1/OXVd9Ydh0b9elZmSanGhFTeY2q7gVUM79/ToX+SJEmS1EGZfGqDhobEnv/1YJOybuVldKsoW5+QypNS3SrK6bZBWVZeWR50b1JWRmX5+nrdK5qVF9TLtiunsiKalheUVZT75bujWlNXzwtvrsgSTXmSac4by1hdm/UA6l5RxrihfTlq92FMGNaPXYf1Y8z7etOjsrzEkRdHRXkZfcrL6NOjkh2H9G6xztKa2iwh9U4NCxqTU++sYu4by3lo9ltNek9FwLZ9exQkpXrmiars+TZ9HNInSZIkSaVi8qkNEnD5kRNYW9fA2voG1tQ1ZMt1DdTWN6wrX1uXr6tvoLaugdW1DSxbVddk/bqf+XJ9w2ZOrrMRZUGeqCqnd/cKenVv/FlBnx4V9OpWQe8eFfTuXrGuvHC5T4+mZT0q7eH1bqxcU8ecN5atGzo38/VlvPjP5dTln3Wf7hWMH9aXz3xo1LoeTTsM7mXysJl+VZW8v6of7x/Rb4N1DQ2Jfy5fzYLFq/IeU/njnRqemPc2by5b3aR+t8YhfQOq1veeyntOjRxYRb+eDumTJEmSpGIx+dQG5WXByR8eXZR91zckagsTWvVNE1utla+ta2BNYSKrsby+gdW19axcU8+KNbWsXFPPstV1vLF0NStW17FyTR0r1ta1aULp8rKgV7csgdW7WWKqecKqV/cK+jRLaGXblNOne2WnTWQtqVm7bl6mma9ld557+e31d5wb1KsbE4b344BdhmQ9mob3ZeSAKnvhvEdlZcHQfj0Z2q8ne28/cIP1q2vreW3JqoKk1CrmL8qSU3+dv+GQvr49KthuUFWThNTIAVnvqeEDetK9onP2QJMkSZKk9mDyqcTKy4LysvJ2HV6VUqJmbT0r19SxfE2ekFpdx4o12WPlmjpWFCSvljcmrdbUsTxPZDU+X7GmbYmssqBpgirvgdU9H45YUV5GRVlkj/JsOGN5WVDZrLyiPF8uK1u3XVYvKyusV1lWsK5xP03qFe4/33e+vryF5NBby1ZnQ+byJNPM15recW5Yvx5MGN6Poz4wfF2Ppvf17d4pk24dXY/KcnYc0nujQ/qaDufL5pt6/s3lPDxnwyF92/TpzsBe3RlQVUn/qkr6V3Wjf89KBlR1o19V9rN/VSUDqirp1zNbrrQnmyRJkiQBJp+6pIhY11tpm/e4r5QSq2rrmySvsgRWlrxasaZ+fY+rxsfqOlauzRJZb9etpa6+gbqGRF1DA3X1idr6RH3jcv6zbgsOT2yLCNYlryrKs+RRYW+Z7Qf3Yo/t+vPZfbOhcxOG9WPgFrzjnIqrX1Ul/ar6sevwlof0NU6E3picem3JKpbUrOWdmlqef3M5S2pqWbKqdqPDZnt3r8gTVXmSquf6JNW65FWvLFk1IC/r17OyxcSnJEmSJG3NTD7pPYkIqrpVUNXtvSeyNialRH1DypNUibr6hjxJlQ1brGvIEla19Wld0mrdumb1GrdvTGrVNTSsS3g1bl/X0NDkdVJKjB7ciwnD+jFuaB/69HCOoM6qrCzYtl8Ptu3Xo8UhfY1SSixfU8fSmlreqVnLkvzn0lW1vLOyliWrsrLGpFX1O6vWrd9Yb8G+PSoY0CtLTvVf16OqMXm1vqx/VZ606tmNPj0qHMopSZIkqcMy+aStQkQ+JM6pd9RBRAR9e1TSt0clIwdWtXm7hobEstW163pPvVOztkkCa0nN2rw8W3757ZUsqVnLsmbzVBUqC9b1rBo5sIobTtl7S7xFSZIkSdoiTD5JUjsqK4u899LmDdOsq29g2eq6pkmqwt5W+XMnR5ckdWURcR1wOPBWSmnXvOzbwBHAWuAfwOdTSkvydRcCXwDqgbNSSr8rSeCS1MmZfJKkrUBFeRkDe3VzbjFJkjbueuCHwI0FZQ8CF6aU6iLiW8CFwPkRMR44HpgADAMeioidU0r17RyzJHV63o5JkiRJUqeQUnoMWNys7IGUUuP49aeAEfnyUcAtKaU1KaWXgXmAY9clqQhMPkmSJEnqKk4B7s+XhwMLCtZV52WSpC3M5JMkSZKkTi8iLgbqgF82FrVQrcV70kbEqRExIyJmLFy4sFghSlKnZfJJkiRJUqcWESeTTUT+mZRSY4KpGhhZUG0E8HpL26eUrk0pTUwpTRwyZEhxg5WkTsjkkyRJkqROKyIOAc4Hjkwp1RSsuhs4PiK6R8T2wBjg6VLEKEmdnXe7kyRJktQpRMQ0YDIwOCKqgUvJ7m7XHXgwIgCeSimdnlKaFRHTgdlkw/G+5J3uJKk4TD5JkiRJ6hRSSlNaKJ66kfpXAlcWLyJJEjjsTpIkSZIkSUVk8kmSJEmSJElFY/JJkiRJkiRJRWPySZIkSZIkSUVj8kmSJEmSJElFY/JJkiRJkiRJRWPySZIkSZIkSUVj8kmSJEmSJElFY/JJkiRJkiRJRWPySZIkSZIkSUVj8kmSJEmSJElFY/JJkiRJkiRJRWPySZK01YuIHhHxdEQ8GxGzIuLyvDwi4sqIeCEi5kTEWaWOVZIkSepqKkodgCRJW8Aa4MCU0oqIqAQej4j7gXHASGBsSqkhIrYpaZSSJElSF1S0nk8RcV1EvBURMwvKjs2vSDdExMRm9S+MiHkR8XxEfKyg/JC8bF5EXFBQvn1E/CkiXoyIWyOiW7HeiySpY0uZFfnTyvyRgDOA/0wpNeT13ipRiJIkSVKXVcxhd9cDhzQrmwl8EnissDAixgPHAxPyba6JiPKIKAf+FzgUGA9MyesCfAv4n5TSGOAd4AtFeh+SpK1A3m78DXgLeDCl9CdgR+C4iJgREfdHxJjSRilJkiR1PUVLPqWUHgMWNyubk1J6voXqRwG3pJTWpJReBuYBe+ePeSmll1JKa4FbgKMiIoADgV/n298AfKJIb0WStBVIKdWnlHYHRgB7R8SuQHdgdUppIvBT4LqWto2IU/ME1YyFCxe2X9CSJElSF9BRJhwfDiwoeF6dl7VWPghYklKqa1YuSeriUkpLgEfJetJWA7flq+4Admtlm2tTShNTShOHDBnSLnFKkiRJXUVHST5FC2XpXZS3vHOvaEtSpxYRQyKif77cE/gXYC5wJ1lPWYCPAi+UJkJJkiSp6+ood7urJrsbUaMRwOv5ckvlbwP9I6Ii7/1UWH8DKaVrgWsBJk6c2GqSSpK01RoK3JDPFVgGTE8p3RsRjwO/jIgvAyuAL5YySEmSJKkr6ijJp7uBmyPiu8AwYAzwNFkPpzERsT3wGtmk5CeklFJEPAIcQzYP1MnAXSWJXJJUciml54A9WihfAny8/SOSJEmS1Khow+4iYhrwJLBLRFRHxBci4uiIqAb2BX4TEb8DSCnNAqYDs4HfAl/KJ46tA84EfgfMIbuSPSt/ifOB/4iIeWRzQE0t1nuRJEmSJEnSu1O0nk8ppSmtrLqjlfpXAle2UH4fcF8L5S+R3Q1PkiRJkiRJHVRHmXBckiRJkiRJnZDJJ0mSJEmSJBWNySdJkiRJkiQVjcknSZIkSZIkFY3JJ0mSJEmSJBWNySdJkiRJkiQVjcknSZIkSZIkFY3JJ0mSJEmSJBWNySdJkiRJkiQVjcknSZIkSZIkFY3JJ0mSJEmSJBWNySdJkiRJkiQVjcknSZIkSZIkFY3JJ0mSJEmSJBWNySdJkiRJkiQVjcknSZIkSZIkFY3JJ0mSJEmSJBWNySdJkiRJkiQVjcknSZIkSZIkFY3JJ0mSJEmSJBWNySdJkiRJkiQVTZuTTxHRq5iBSJLUyDZHkgS2B5LUWWwy+RQRH46I2cCc/PkHIuKaokcmSepybHMkSWB7IEmdTVt6Pv0P8DFgEUBK6Vlg/2IGJUnqsmxzJElgeyBJnUqbht2llBY0K6ovQiySJNnmSJIA2wNJ6kwq2lBnQUR8GEgR0Q04i7z7qyRJW5htjiQJbA8kqVNpS8+n04EvAcOBamD3/LkkSVuabY4kCWwPJKlT2WTPp5TS28Bn2iEWSVIXZ5sjSQLbA0nqbDaZfIqInwOpeXlK6ZSiRCRJ6rJscyRJYHsgSZ1NW+Z8urdguQdwNPB6ccKRJHVxtjmSJLA9kKROpS3D7m4rfB4R04CHihaRJKnLss2RJIHtgSR1Nm2ZcLy5McB2WzoQSZJaYJsjSYI2tgcRcV1EvBURMwvKjo2IWRHREBETm9W/MCLmRcTzEfGxIsQtSaJtcz4tJxtvHfnPN4HzixyXJKkLss2RJMF7ag+uB34I3FhQNhP4JPCTZq8xHjgemAAMAx6KiJ1TSvXvNX5JUlNtGXbXpz0CkSTJNkeSBO++PUgpPRYRo5uVzQGI+P/bu/coS+v6zvfvD32hu7k1lwabS3ORHtTgCNiDt4QooALHDJjEDJ5MJI7amRw0mKzMBM85K4zO0tEVJ84YjXNQiZg4KKNxYJSoLd7PERAQEERDh2vbSLfSXBqw6W6+54/nKamurqouqvdTu2rv92utvfbev/3svb+7Keq76rN/v9+TsYefDXy6qrYAdyVZC5wMfHc67y1JmtiE4VOSkyZ7YlXd2PtyJEnDyJ4jSYIZ7weHAdeMur+uHRuvrtXAaoAVK1wNLknP1GQzn/7zJI8VcGqPa5EkDS97jiQJZrYf7DQVqn2PnQerLgYuBli1atW4x0iSJjZh+FRVr5jJQiRJw8ueI0mCGe8H64AjRt0/HFg/g+8vSUNjl3s+ASQ5HngesGhkrKo+OfEzJEmaHnuOJAlmpB9cCfz3JH9Js+H4SuC6Hr6+JKk1lbPdXQS8nOYX/1XAmcB32PEMEpIk7TZ7jiQJpt8PklzWPu+gJOuAi4AHgb8ClgFfTHJTVb26qm5LcjnwQ2AbcL5nupOkbkxl5tNvAy8Avl9Vb0xyCPCxbsuSJA0pe44kCabZD6rq9RM89PkJjn838O5pVylJmpI9pnDME1X1FLAtyb7ABuCYbsuSJA0pe44kCewHkjRQpjLz6fokS4GPAjcAm3EttCSpG/YcSRLYDyRpoOwyfKqq/6O9+d+SfAnYt6pu6bYsSdIwsudIksB+IEmDZsJld0l+mOT/SvLskbGquttf+pKkXrPnSJLAfiBJg2qyPZ9eD+wNfCXJtUnenuTQGapLkjRc7DmSJLAfSNJAmjB8qqqbq+odVfVs4ALgSOCaJF9L8pYZq1CSNPDsOZIksB9I0qCaytnuqKprquqPgTcA+wMf6rQqSdLQsudIksB+IEmDZJcbjif5FzTTX38LuBu4GPgf3ZYlSRpG9hxJEtgPJGnQTBg+JXkP8K+ATcCngZdV1bqZKkySNDzsOZIksB9I0qCabObTFuDMqvrHmSpGkjS07DmSJLAfSNJAmjB8qqp3zmQhkqThZc+RJIH9QJIG1ZQ2HJckSZIkSZKmw/BJkiRJkiRJndll+JTGv07y5+39FUlO7r40SdKwsedIksB+IEmDZiozn/4aeAnN4JSssAAAIABJREFUqU4BHgU+3FlFkqRhZs+RJIH9QJIGymRnuxvxoqo6Kcn3AapqU5KFHdclSRpO9hxJEtgPJGmgTGXm09Yk84ACSLIMeKrTqiRJw8qeI0kC+4EkDZSphE8fBD4PHJzk3cB3gPd0WpUkaVjZcyRJYD+QpIGyy2V3VfWpJDcApwEBzqmq2zuvTJI0dOw5kiSwH0jSoNll+JTkxcBtVfXh9v4+SV5UVdd2Xp0kaajYcyRJYD+QpEEzlWV3HwE2j7r/WDsmSVKv2XMkSWA/kKSBMpXwKVVVI3eq6immNmPqkiQbktw6auyAJGuS3NFe79+OvzzJw0luai9/Puo5ZyT5cZK1SS4cNX50kmvb1/qMZ7+QpIEwrZ4jSRo49gNJGiBTCZ/uTPJHSRa0lwuAO6fwvE8AZ4wZuxC4uqpWAle390d8u6pOaC/vAmjPcPFh4EzgecDrkzyvPf59wAfa19oEvGkKNUmSZrdp9Zwki5Jcl+TmJLcleWc7/okkd436cuOEzj+BJKkXpvs3iCRpFppK+PRvgZcCPwHWAS8CVu/qSVX1LeDBMcNnA5e2ty8FztnFy5wMrK2qO6vqSeDTwNlJApwKfPYZvJYkafabVs8BtgCnVtULgBOAM9r9QgD+3agvN27qomhJUs9Ntx9IkmahqZztbgNwbo/e75Cqur993fuTHDzqsZckuRlYD/xpVd0GHAbcN+qYkcZzIPBQVW0bNX5Yj2qUJPXJdHtOuzRjZG+QBe2lJn6GJGk26/HfIJKkPpvK3k3LgLcAR40+vqr+TQ/ruBE4sqo2JzkL+J/ASprTqo5Vk4yPK8lq2m9KVqxYsfvVSpI6sTs9p12qfQNwLPDhqro2yR8C7273ErwauLCqtozzXPuEJM0iM/Q3iCRphkxl074rgG8DXwW27+b7PZBkeTvraTmwAaCqHhk5oKquSvLXSQ6imdF0xKjnH04zM+pnwNIk89vZTyPj46qqi4GLAVatWuU34ZI0e02751TVduCEJEuBzyc5HngH8FNgIU0f+DPgXeM81z4hSbNLL/8GkST12VTCpyVV9Wc9er8rgfOA97bXVwAkeRbwQFVVkpNp9qL6OfAQsDLJ0TTrvc8F/vf2uK8Dv02zD9QvX0uSNKftds+pqoeSfAM4o6re3w5vSfI3wJ/uboGSpBnRy79BJEl9NpUNx7/QLoV7RpJcBnwXOC7JuiRvogmdXpnkDuCV7X1oQqRb2z2fPgicW41twFuBLwO3A5e3e0FB8+31nyRZS7MH1MefaY2SpFlnuj1nWTvjiSSLgdOBH7WzbGlPVHEOcGsvi5UkdWZa/UCSNDtNZebTBcD/mWQLsJVmv6Wqqn0ne1JVvX6Ch04b59gPAR+a4HWuAq4aZ/xOmrPhSZIGx7R6DrAcuLTd92kPmi8rvpDka+2+IQFuojl7kiRp9ptuP5AkzUJTOdvdPjNRiCRJ0+05VXULcOI446fudlGSpBnn3yCSNFimMvOJJPvTnH1u0chYVX2rq6IkScPLniNJAvuBJA2SXYZPSd5MM+31cJolCy+m2cvJb5MlST1lz5Ekgf1AkgbNVDYcvwD4F8A9VfUKmmUNGzutSpI0rOw5kiSwH0jSQJlK+PSLqvoFQJI9q+pHwHHdliVJGlL2HEkS2A8kaaBMZc+nde3pq/8nsCbJJmB9t2VJkoaUPUeSBPYDSRooUznb3Wvbm/8hydeB/YB/6LQqSdJQsudIksB+IEmDZpfL7pL87cjtqvpmVV0JXNJpVZKkoWTPkSSB/UCSBs1U9nz6ldF3kswDXthNOZKkIWfPkSSB/UCSBsqE4VOSdyR5FPjnSR5pL48CG4ArZqxCSdLAs+dIksB+IEmDasLwqar+U1XtA/xFVe3bXvapqgOr6h0zWKMkacDZcyRJYD+QpEE1lWV31yXZb+ROkqVJzumwJknS8LLnSJLAfiBJA2Uq4dNFVfXwyJ2qegi4qLuSJElDzJ4jSQL7gSQNlKmET+MdM7/XhUiShD1HktSwH0jSAJlK+HR9kr9M8uwkxyT5AHBD14VJkoaSPUeSBPYDSRooUwmf3gY8CXwGuBx4Aji/y6IkSUPLniNJAvuBJA2UXU5drarHgAuT7F1Vm2egJknSkLLnSJLAfiBJg2aXM5+SvDTJD4EftvdfkOSvO69MkjR07DmSJLAfSNKgmcqyuw8ArwZ+DlBVNwOndFmUJGlo2XMkSWA/kKSBMpXwiaq6b8zQ9g5qkSTJniNJAuwHkjRIpnK60vuSvBSoJAuBPwJu77YsSdKQsudIksB+IEkDZSozn/4tzZklDgPWASfgmSYkSd2w50iSwH4gSQNlKme7+xnwuzNQiyRpyNlzJElgP5CkQTNh+JTkr4Ca6PGq+qNOKpIkDR17jiQJdr8fJLkEeA2woaqOb8cOAD4DHAXcDfxOVW1KEuC/AmcBjwO/X1U39uBjSJLGmGzm0/UzVoUkadjZcyRJsPv94BPAh4BPjhq7ELi6qt6b5ML2/p8BZwIr28uLgI+015KkHpssfDq9qn4vyQVV9V9nrCJJ0jCy50iSYDf7QVV9K8lRY4bPBl7e3r4U+AZN+HQ28MmqKuCaJEuTLK+q+6dbvCRpfJNtOP7CJEcC/ybJ/kkOGH2ZqQIlSUPBniNJgm76wSEjgVJ7fXA7fhhw36jj1rVjO0myOsn1Sa7fuHHjNMuQpOE12cyn/wZ8CTgGuAHIqMeqHZckqRfsOZIkmNl+kHHGxt1vqqouBi4GWLVq1YR7UkmSxjfhzKeq+mBVPRe4pKqOqaqjR138I0CS1DP2HEkSdNYPHkiyHKC93tCOrwOOGHXc4cD6aRcvSZrQZMvuAKiqP0zyq0neCJDkoCRHd1+aJGnY2HMkSdDzfnAlcF57+zzgilHjb0jjxcDD7vckSd2YbNkdAEkuAlYBxwF/AywE/g54WbelSZKGjT1HkgTT7wdJLqPZXPygJOuAi4D3ApcneRNwL/C69vCrgLOAtcDjwBt7/kEkScAUwifgtcCJwI0AVbU+yT6dViVJGlb2HEkSTLMfVNXrJ3jotHGOLeD83SlSkjQ1u1x2BzzZ/mIugCR7dVuSJGmI2XMkSWA/kKSBMpXw6fIk/w+wNMlbgK8CH+u2LEnSkLLnSJLAfiBJA2WXy+6q6v1JXgk8QrPm+s+rak3nlUmSho49R5IE9gNJGjRT2fOJ9hf9GoAk85L8blV9qtPKJElDyZ4jSQL7gSQNkgmX3SXZN8k7knwoyavaU5C+FbgT+J2ZK1GSNOjsOZIksB9I0qCabObT3wKbgO8Cbwb+Hc0pTs+uqptmoDZJ0vCw50iSwH4gSQNpsvDpmKp6PkCSjwE/A1ZU1aMzUpkkaZjYcyRJYD+QpIE02dnuto7cqKrtwF3+0pckdcSeI0kC+4EkDaTJZj69IMkj7e0Ai9v7Aaqq9u28OknSsLDnSJLAfiBJA2nC8Kmq5s1kIZKk4WXPkSSB/UCSBtVky+4kSZIkSZKk3WL4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkznYZPSS5JsiHJraPGDkiyJskd7fX+7XiSfDDJ2iS3JDlp1HPOa4+/I8l5o8ZfmOQH7XM+mCRdfh5JkiRJkiQ9M13PfPoEcMaYsQuBq6tqJXB1ex/gTGBle1kNfASasAq4CHgRcDJw0Uhg1R6zetTzxr6XJEmSJEmS+qjT8KmqvgU8OGb4bODS9valwDmjxj9ZjWuApUmWA68G1lTVg1W1CVgDnNE+tm9VfbeqCvjkqNeSJEmSJEnSLNCPPZ8Oqar7Adrrg9vxw4D7Rh23rh2bbHzdOOOSpCGTZFGS65LcnOS2JO8c8/hfJdncr/okSZKkYTabNhwfb7+mmsb4zi+crE5yfZLrN27cuBslSpJmqS3AqVX1AuAEmhmyLwZIsgpY2s/iJEmSpGHWj/DpgXbJHO31hnZ8HXDEqOMOB9bvYvzwccZ3UlUXV9Wqqlq1bNmynnwISdLs0S7ZHpnZtKC9VJJ5wF8A/75vxUmSJElDrh/h05XAyBnrzgOuGDX+hvasdy8GHm6X5X0ZeFWS/duNxl8FfLl97NEkL27PcveGUa8lSRoySeYluYnmS401VXUt8FbgypHl3pIkSZJm3vwuXzzJZcDLgYOSrKM5a917gcuTvAm4F3hde/hVwFnAWuBx4I0AVfVgkv8IfK897l1VNbKJ+R/SnFFvMfAP7UWSNISqajtwQpKlwOeTnELTY16+q+cmWU1z9lRWrFjRZZmSJEnS0Ok0fKqq10/w0GnjHFvA+RO8ziXAJeOMXw8cvzs1SpIGS1U9lOQbwCuAY4G1zQRZliRZW1XHjvOci4GLAVatWjXu/oGSJEmSpmc2bTguSdK0JFnWzngiyWLgdOCGqnpWVR1VVUcBj48XPEmSJEnqVqcznyRJmiHLgUvbDcb3AC6vqi/0uSZJkiRJGD5JkgZAVd0CnLiLY/aeoXIkSZIkjeKyO0mSJEmSJHXG8EmSJEmSJEmdMXySJEmSJElSZwyfJEmSJEmS1BnDJ0mSJEmSJHXG8EmSJEmSJEmdMXySJEmSJElSZwyfJEmSJEmS1BnDJ0mSJEmSJHXG8EmSJEmSJEmdMXySJEmSJElSZwyfJEmSJEmS1BnDJ0mSJEkDL8kFSW5NcluSt7djByRZk+SO9nr/ftcpSYPI8EmSJEnSQEtyPPAW4GTgBcBrkqwELgSurqqVwNXtfUlSjxk+SZIkSRp0zwWuqarHq2ob8E3gtcDZwKXtMZcC5/SpPkkaaIZPkiRJkgbdrcApSQ5MsgQ4CzgCOKSq7gdorw/uY42SNLDm97sASZIkSepSVd2e5H3AGmAzcDOwbarPT7IaWA2wYsWKTmqUpEHmzCdJkiRJA6+qPl5VJ1XVKcCDwB3AA0mWA7TXGyZ47sVVtaqqVi1btmzmipakAWH4JEmSJGngJTm4vV4B/CZwGXAlcF57yHnAFf2pTpIGm8vuJEmSJA2DzyU5ENgKnF9Vm5K8F7g8yZuAe4HX9bVCSRpQhk+SJEmSBl5V/do4Yz8HTutDOZI0VFx2J0mSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM70JXxKckGSW5PcluTt7dh/SPKTJDe1l7NGHf+OJGuT/DjJq0eNn9GOrU1yYT8+iyRJkiRJkiY2f6bfMMnxwFuAk4EngS8l+WL78Aeq6v1jjn8ecC7wK8ChwFeT/LP24Q8DrwTWAd9LcmVV/XAGPoYkSZIkSZKmYMbDJ+C5wDVV9ThAkm8Cr53k+LOBT1fVFuCuJGtpgiuAtVV1Z/s6n26PNXySJEmSJEmaJfqx7O5W4JQkByZZApwFHNE+9tYktyS5JMn+7dhhwH2jnr+uHZtoXJIkSZIkSbPEjIdPVXU78D5gDfAl4GZgG/AR4NnACcD9wH9un5LxXmaS8Z0kWZ3k+iTXb9y4cfc+gCRJkiRJkqasLxuOV9XHq+qkqjoFeBC4o6oeqKrtVfUU8FGeXlq3jqdnRgEcDqyfZHy897u4qlZV1aply5b1+uNIkiRJkiRpAv06293B7fUK4DeBy5IsH3XIa2mW5wFcCZybZM8kRwMrgeuA7wErkxydZCHNpuRXztRnkCTNHkkWJbkuyc3tmVTf2Y5/vB27Jclnk+zd71olSZKkYdOPDccBPpfkQGArcH5VbUryt0lOoFk6dzfwBwBVdVuSy2k2Et/WHr8dIMlbgS8D84BLquq2mf8okqRZYAtwalVtTrIA+E6SfwD+uKoeAUjyl8Bbgff2sU5JkiRp6PQlfKqqXxtn7PcmOf7dwLvHGb8KuKq31UmS5pqqKmBze3dBe6lRwVOAxUywN6AkSZKk7vRl2Z0kSb2WZF6Sm4ANwJqqurYd/xvgp8BzgL+a4LmemEKSJEnqiOGTJGkgtCetOIHmBBQnJzm+HX8jcChwO/CvJniuJ6aQJEmSOmL4JEkaKFX1EPAN4IxRY9uBzwC/1aeyJEmSpKFl+CRJmvOSLEuytL29GDgd+HGSY9uxAL8B/Kh/VUqSJEnDqV9nu5MkqZeWA5cmmUfzxcrlwBeBbyfZFwhwM/CH/StRkiRJGk6GT5KkOa+qbgFOHOehl810LZIkSZJ25LI7SZIkSZIkdcbwSZIkSdLAS/LHSW5LcmuSy5IsSnJ0kmuT3JHkM0kW9rtOSRpEhk+SJEmSBlqSw4A/AlZV1fHAPOBc4H3AB6pqJbAJeFP/qpSkwWX4JEmSJGkYzAcWJ5kPLAHuB04FPts+filwTp9qk6SBZvgkSZIkaaBV1U+A9wP30oRODwM3AA9V1bb2sHXAYf2pUJIGm+GTJEmSpIGWZH/gbOBo4FBgL+DMcQ6tCZ6/Osn1Sa7fuHFjd4VK0oAyfJIkSZI06E4H7qqqjVW1Ffh74KXA0nYZHsDhwPrxnlxVF1fVqqpatWzZspmpWJIGiOGTJEmSpEF3L/DiJEuSBDgN+CHwdeC322POA67oU32SNNAMnyRJkiQNtKq6lmZj8RuBH9D8HXQx8GfAnyRZCxwIfLxvRUrSAJu/60MkSZIkaW6rqouAi8YM3wmc3IdyJGmoOPNJkiRJkiRJnTF8kiRJkiRJUmcMnyRJkiRJktQZwydJkiRJkiR1xvBJkiRJkiRJnTF8kiRJkiRJUmcMnyRJkiRJktQZwydJkiRJkiR1xvBJkiRJkiRJnTF8kiRJkiRJUmcMnyRJkiRJktQZwydJkiRJkiR1xvBJkiRJkiRJnTF8kiRJkiRJUmcMnyRJkiRJktQZwydJkiRJkiR1xvBJkiRJkiRJnTF8kiRJkiRJUmcMnyRJkiRJktQZwydJkiRJkiR1xvBJkiRJkiRJnTF8kiRJkiRJUmfm97sAST3w5GNw17dh7VfhwTth6Qo44JinL/sfBQuX9LtKSZIkSdIQMnyS5qIq2PhjWLumCZzu+f9g+5OwYAkceCys/z488eCOz9nn0DaMOnrHYOqAo2HPffrzOSRJkiRJA8/wSZorfvEI3PVNuGMNrL0aHlnXjC97Dpy8Go49HVa8BBYsasaf2AQP3tXMhPrl9Z3wj1+Gxzbs+Np7HbxzIDVyvXj/mf2ckiRJkqSBYvg0VXd9Cxbs1SxdWrgXLNy7uZ63EJJ+V6dBVAU//UEzs2ntV+G+a+GpbbBwHzjm1+GUP20Cp6VHjP/8xfvDYfvDYSft/NiWR3cMpEYCqju/ATf/951fZ4dgatRlyYH+/EuSJEmSJmX4NBVPbYdLf2P8x/aY34ZSk132bpZDjQ6tFi4ZdXvsMXvBvAUz+xk1Ozz+INz59WZm09qvwuYHmvFnPR9e+rYmbDriRbv/87HnPrD8nzeXsZ58HDbd/XQotakNqe67Fn7wWaBGvc6+4yzjay97H2IwJUmSJEkyfJqy3/9is6nzDpfNzfXWx5++/eRjzR/vmzc8fX/rY7BlM9T2qb/fvD13DqgWjBNYjZ2JtWBJEyzssxz2XQ6LlhoAzGZPPdXszzQyu+kn10M9BYv2g2efCse+Eo49DfZ51szVtHAJHPK85jLWti3w0L1jZkzdCfffDD+8csef8QVLnt7sfGwwte9hsIcn25QkSZKkYWD4NBV7zIOjfnX3XqOq2RB6ygHWmGO2Pt7cfmTd0wHXyGOjZ6KMNX9xE0Ltc2h7vRz2PXTU9bNg72fB/IW79/k0dZs3wj99rdks/J++Bo//HAgceiL82p/CylfCoSfBvFn4v+f8PeGglc1lrO1b4eH7dt5j6mf/CHd8pfn5HzFvzyaU+uXlSFh65NO33QBdkiRJkgbGLPzrdkAlzR/u8/eEJQf07nWrYOsTO4ZXv3gYHr0fHrm/vV7fXN93HTz6U9i+ZWxxsNdB4wRTy3cMrpxFNT3btzUzmtZ+tdks/P6bmvElBzXL6I49vZnltNdB/a1zd81b8PTMprGe2g6P/GTnfaY23d2cqe/JR3c8fvEB44dSS4+E/Y4wLJUkSZKkOcTwaa5L2qV3S6YWXlQ1+wo9ur4Np9Y3gdRIQPXwT2Dd99rZOGPschbV8mYmlftVNf+eI/s23fn1JhDMHnD4yfCK/7tZSrf8hOFZerbHPFi6orkc8+s7PlbVnJlv012w6R546J4mlNp0T7Oc7/b/1Wy0PiJ7NMv2lh7ZBFL7H7XjbfeakiRJkqRZxfBp2CSw14HN5VnPn/i4bVtGzZ5a7yyqXdn2JNx3Tbt309XwwK3N+D7L4bm/0cxuOublzZnjtKOkmQ245AA47IU7P/7U9ubnbnQo9dA9zfXaq2HzT3c8fv6iJuQaG0qN3F603wx8KEmSJEnSCMMnjW/+qD15JjLeLKrR11OZRbXkoGaj9D33bjdN33uC+/vsuLH6yP0FS/oXYm265+mw6a5vNsse95gPK14Cp7+zCZwO+ZXBCtn6YY95sPSI5jLe3mtbn4CH7muCqV8GVO3te6+BLY/sePyipRMs6TuqeY/5e3b9iSRJkiRpqBg+afqmOotq6y+a2SnjzaJ64sEmtHn0p+0G6o8216M3p568iCaQ2nPvUeHUM7m/z6iwq70/0UbfW38B9/y/T5+Z7mf/2IzvdwQ8/3XNRuFHn+Jm2TNtwWJY9s+ay1gjS/pGZkqNnj31wG3w438Y87OWZqbe2H2mlq5oZq0t2hf23Lf5GRqWJZOSJEmStJsMn9S9BYt2PYtqrG1Pthuot5uob9n89P3Rt8c+NnL/kfU73t/62NTfe96eO4dTe8yH9d+HbU80jx/1Mnjh78Oxr2zO/Obsptlp9JK+Q0/c+fGnnmqC0J2W9N0Nd36jeWzcs0mmCaFGwqjR14v2GzM29r4BliRJkqThYvik2Wn+Qph/QO/ODPjUU00ANW5Y1c622uH+mGBr2y/gpDc0S+mOelkTSGnu22MP2O+w5nLkS3d+fNuWZknfw/fCEw81S/h+8cj414/eDz/78dP3R2+SPq6JAqz9xh8zwJIkSZI0Rxk+aTjssUezHG7PfcBVcZqq+XvCQcc2l2eiqtmL6hcPjwqoHp44uOplgLXfEfBbH532R5YkSZKkXjN8kqReS2DhkubC8um9xkiAtUNI9fDT1xMFWFse7elHkSRJkqTdZfgkSbPR6ABrn2f1uxpJkiRJmjY3C5EkSZIkSVJnDJ8kSZIkSZLUGcMnSZIkSZIkdcbwSZIkSZIkSZ0xfJIkSZIkSVJnDJ8kSZIkSZLUmb6ET0kuSHJrktuSvL0dOyDJmiR3tNf7t+NJ8sEka5PckuSkUa9zXnv8HUnO68dnkSRJkiRJ0sRmPHxKcjzwFuBk4AXAa5KsBC4Erq6qlcDV7X2AM4GV7WU18JH2dQ4ALgJe1L7WRSOBlSRJkiRJkmaHfsx8ei5wTVU9XlXbgG8CrwXOBi5tj7kUOKe9fTbwyWpcAyxNshx4NbCmqh6sqk3AGuCMmfwgkiRJkiRJmlw/wqdbgVOSHJhkCXAWcARwSFXdD9BeH9wefxhw36jnr2vHJhqXJEmSJEnSLDF/pt+wqm5P8j6amUqbgZuBbZM8JeO9zCTjO79AsppmyR4rVqx4RvVKkiRJkiRp+vqy4XhVfbyqTqqqU4AHgTuAB9rldLTXG9rD19HMjBpxOLB+kvHx3u/iqlpVVauWLVvW2w8jSZIkSZKkCfXrbHcHt9crgN8ELgOuBEbOWHcecEV7+0rgDe1Z714MPNwuy/sy8Kok+7cbjb+qHZMkSZIkSdIsMePL7lqfS3IgsBU4v6o2JXkvcHmSNwH3Aq9rj72KZl+otcDjwBsBqurBJP8R+F573Luq6sGZ/BCSJEmSJEmaXF/Cp6r6tXHGfg6cNs54AedP8DqXAJf0vEBJ0pySZBHwLWBPmt722aq6KMmngFU0X3ZcB/xBVW3tX6WSpH5IchzwmVFDxwB/DnyyHT8KuBv4nfZM2pKkHurLsjtJknpsC3BqVb0AOAE4o12q/SngOcDzgcXAm/tXoiSpX6rqx1V1QlWdALyQZkXF54ELgauraiVwdXtfktRjhk+SpDmvGpvbuwvaS1XVVe1jRTPz6fC+FSlJmi1OA/6pqu4BzgYubccvBc7pW1WSNMAMnyRJAyHJvCQ30ZwtdU1VXTvqsQXA7wFf6ld9kqRZ41yaEx4BHNKezIj2+uC+VSVJAyzNl8HDI8lG4J5pPv0g4Gc9LKcrc6VOsNYuzJU6wVq7sDt1HllVy3pZTD8kWUqzlOJtVXVrO/ZR4LGqevsEz1kNrG7vHgf8eJpvPww/JzPNWntvrtQJc6fWuVInDHmfSLIQWA/8SlU9kOShqlo66vFNVbX/OM+zT8xec6XWuVInWGsX5kqd0FHTOiWpAAAIVUlEQVSfGLrwaXckub6qVvW7jl2ZK3WCtXZhrtQJ1tqFuVJn15JcRBM2vb+9fSLwm1X1VMfvOyf+/edKnWCtXZgrdcLcqXWu1Alzq9YuJDmb5mzbr2rv/xh4eVXdn2Q58I2qOq7D958T//5zpU6YO7XOlTrBWrswV+qE7mp12Z0kac5Lsqyd8USSxcDpwI+SvBl4NfD6roMnSdKc8HqeXnIHcCVwXnv7POCKGa9IkobA/H4XIElSDywHLk0yj+aLlcur6gtJttEstf5uEoC/r6p39bFOSVKfJFkCvBL4g1HD7wUuT/Im4F7gdf2oTZIGneHTM3NxvwuYorlSJ1hrF+ZKnWCtXZgrdfZUVd1Cs7Ru7PhM97m58u8/V+oEa+3CXKkT5k6tc6VOmFu19lRVPQ4cOGbs5zRnv5spc+Xff67UCXOn1rlSJ1hrF+ZKndBRre75JEmSJEmSpM6455MkSZIkSZI6Y/g0BUkuSbIhya39rmUySY5I8vUktye5LckF/a5pIkkWJbkuyc1tre/sd02TSTIvyfeTfKHftUwmyd1JfpDkpiTX97ueiSRZmuSzSX7U/ry+pN81jSfJce2/5cjlkSRv73ddE0nyx+3/T7cmuSzJon7XNCzsE70113oE2Cd6zT7RDftEf8yVHgH2iS7ZJ3rLPtF7XfcIl91NQZJTgM3AJ6vq+H7XM5H29LDLq+rGJPsANwDnVNUP+1zaTtLs/LtXVW1OsgD4DnBBVV3T59LGleRPgFXAvlX1mn7XM5EkdwOrqupn/a5lMkkuBb5dVR9LshBYUlUP9buuybQbWf8EeFFV3dPvesZKchjN/0fPq6onklwOXFVVn+hvZcPBPtFbc61HgH2i1+wTvWef6J+50iPAPtEl+0Rv2Sd6ayZ6hDOfpqCqvgU82O86dqWq7q+qG9vbjwK3A4f1t6rxVWNze3dBe5mVSWiSw4H/DfhYv2sZBEn2BU4BPg5QVU/O9kbROg34p9nWKMaYDyxOMh9YAqzvcz1Dwz7RW3OpR4B9otfsE52yT/TBXOkRYJ/oin2it+wTnem0Rxg+DagkR9Gc+ena/lYysXbq6U3ABmBNVc3WWv8L8O+Bp/pdyBQU8JUkNyRZ3e9iJnAMsBH4m3bq8ceS7NXvoqbgXOCyfhcxkar6CfB+mtNE3w88XFVf6W9Vms1me5+YQz0C7BO9Zp/ogH1Cz5R9oqfsE71ln+ixmegRhk8DKMnewOeAt1fVI/2uZyJVtb2qTgAOB05OMuumISd5DbChqm7ody1T9LKqOgk4Ezi/neY928wHTgI+UlUnAo8BF/a3pMm1U3n/JfA/+l3LRJLsD5wNHA0cCuyV5F/3tyrNVnOhT8yFHgH2iY7YJzpgn9AzYZ/oHftEJ+wTPTYTPcLwacC0a54/B3yqqv6+3/VMRTtF8hvAGX0uZTwvA/5lu/b508CpSf6uvyVNrKrWt9cbgM8DJ/e3onGtA9aN+nbqszTNYzY7E7ixqh7odyGTOB24q6o2VtVW4O+Bl/a5Js1Cc61PzPIeAfaJLtgnumGf0JTYJ3rOPtF79one67xHGD4NkHbjvY8Dt1fVX/a7nskkWZZkaXt7Mc0P+4/6W9XOquodVXV4VR1FM03ya1U1K78lTLJXuzEk7bTTVwGz7qwqVfVT4L4kx7VDpwGzahPLcbyeWTpFdpR7gRcnWdL+LjiNZp8G6ZfmSp+YKz0C7BNdsE90xj6hXbJP9J59ovfsE53ovEcYPk1BksuA7wLHJVmX5E39rmkCLwN+jyZNHzmV41n9LmoCy4GvJ7kF+B7NOu1ZfdrROeAQ4DtJbgauA75YVV/qc00TeRvwqfa//wnAe/pcz4SSLAFeSZP+z1rtNz+fBW4EfkDz+/3ivhY1ROwTPWeP6IZ9ogP2Ce3KHOoRYJ8YdvaJDsyFPjETPSJVs/akAJIkSZIkSZrjnPkkSZIkSZKkzhg+SZIkSZIkqTOGT5IkSZIkSeqM4ZMkSZIkSZI6Y/gkSZIkSZKkzhg+SdOQZHt76tmbk9yY5KU9et2jkty6i2Mua497e5Jze/G+kqTesk9IkiZjn9CwMXySpueJqjqhql4AvAP4TzP43kdX1d3ArwPfnsH3lSRNnX1CkjQZ+4SGiuGTtPv2BTYBJNk7ydXttxc/SHJ2O35UktuTfDTJbUm+kmRx+9gL2288vgucP9GbJPlUkh8CxyW5CXgV8MUkb+78E0qSdod9QpI0GfuEBl6qqt81SHNOku3AD4BFwHLg1Kq6Icl8YElVPZLkIOAaYCVwJLAWWFVVNyW5HLiyqv4uyS3A26rqm0n+Ajizqo6f4H1/BzgC+BzwF1X1uq4/qyTpmbNPSJImY5/QsHHmkzQ9I9NknwOcAXwySYAA72kbwFeBw4BD2ufcVVU3tbdvAI5Ksh+wtKq+2Y7/7S7e90TgJuD57bUkaXayT0iSJmOf0FCZ3+8CpLmuqr7bfiuxDDirvX5hVW1NcjfNtxkAW0Y9bTuwmKa57HL6YZKzgPcARwOvad/jsSSnV9UrevVZJEm9Z5+QJE3GPqFh4MwnaTcleQ4wD/g5sB+woW0Ur6CZHjuhqnoIeDjJr7ZDvzvBcVcBLwRurarnA7cBJ9ooJGn2s09IkiZjn9AwcOaTND2L2036oPm24byq2p7kU8D/SnI9zTTWH03htd4IXJLkceDLkxx3InBzkoXAgqp6ZDfqlyR1yz4hSZqMfUJDxQ3HJUmSJEmS1BmX3UmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM4YPkmSJEmSJKkzhk+SJEmSJEnqjOGTJEmSJEmSOmP4JEmSJEmSpM78/yM6KiRqwEK6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x576 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,3, figsize=[20,8])\n",
    "\n",
    "# numbers 1-8\n",
    "band_count = np.arange(1,9)\n",
    "\n",
    "classes = np.unique(y)\n",
    "print(classes)\n",
    "\n",
    "for class_type in classes:\n",
    "    band_intensity = np.mean(X[y==class_type, :], axis=0)\n",
    "    ax[0].plot(band_count, band_intensity, label=class_type)\n",
    "    ax[1].plot(band_count, band_intensity, label=class_type)\n",
    "    ax[2].plot(band_count, band_intensity, label=class_type)\n",
    "# plot them as lines\n",
    "\n",
    "# Add some axis labels\n",
    "ax[0].set_xlabel('Band #')\n",
    "ax[0].set_ylabel('Reflectance Value')\n",
    "ax[1].set_ylabel('Reflectance Value')\n",
    "ax[1].set_xlabel('Band #')\n",
    "ax[2].set_ylabel('Reflectance Value')\n",
    "ax[2].set_xlabel('Band #')\n",
    "#ax[0].set_ylim(32,38)\n",
    "ax[1].set_ylim(32,38)\n",
    "ax[2].set_ylim(70,140)\n",
    "#ax.set\n",
    "ax[1].legend(loc=\"upper right\")\n",
    "# Add a title\n",
    "ax[0].set_title('Band Intensities Full Overview')\n",
    "ax[1].set_title('Band Intensities Lower Ref Subset')\n",
    "ax[2].set_title('Band Intensities Higher Ref Subset')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [],
   "source": [
    "def str_class_to_int(class_array):\n",
    "    class_array[class_array == b'1'] = 0\n",
    "    class_array[class_array == b'12'] = 1\n",
    "    class_array[class_array == b'26'] = 2\n",
    "    return(class_array.astype(int))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GaussianNB(priors=None, var_smoothing=1e-09)"
      ]
     },
     "execution_count": 259,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.naive_bayes import GaussianNB\n",
    "\n",
    "gnb = GaussianNB()\n",
    "gnb.fit(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 260,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rasterio.plot import show\n",
    "from rasterio.plot import show_hist\n",
    "from rasterio.windows import Window\n",
    "from rasterio.plot import reshape_as_raster, reshape_as_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 261,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(8, 450, 1150)\n",
      "(450, 1150, 8)\n"
     ]
    }
   ],
   "source": [
    "with rasterio.open(image_folder_path) as src:\n",
    "    # may need to reduce this image size if your kernel crashes, takes a lot of memory\n",
    "    img = src.read()[:, 150:600, 250:1400]\n",
    "\n",
    "# Take our full image and reshape into long 2d array (nrow * ncol, nband) for classification\n",
    "print(img.shape)\n",
    "reshaped_img = reshape_as_image(img)\n",
    "print(reshaped_img.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {},
   "outputs": [],
   "source": [
    "class_prediction = gnb.predict(reshaped_img.reshape(-1, 8))\n",
    "\n",
    "# Reshape our classification map back into a 2D matrix so we can visualize it\n",
    "class_prediction = class_prediction.reshape(reshaped_img[:, :, 0].shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\user\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:2: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  \n",
      "C:\\Users\\user\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:3: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n",
      "C:\\Users\\user\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:4: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  after removing the cwd from sys.path.\n"
     ]
    }
   ],
   "source": [
    "class_prediction = str_class_to_int(class_prediction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 265,
   "metadata": {},
   "outputs": [],
   "source": [
    "def color_stretch(image, index):\n",
    "    colors = image[:, :, index].astype(np.float64)\n",
    "    for b in range(colors.shape[2]):\n",
    "        colors[:, :, b] = rasterio.plot.adjust_band(colors[:, :, b])\n",
    "    return colors\n",
    "    \n",
    "# find the highest pixel value in the prediction image\n",
    "n = int(np.max(class_prediction))\n",
    "\n",
    "# next setup a colormap for our map\n",
    "colors = dict((\n",
    "    (0, (48, 156, 214, 255)),   # Blue - Water\n",
    "    (1, (139,69,19, 255)),      # Brown - WetSand\n",
    "    (2, (96, 19, 134, 255)),    # Purple - Emergent Wetland\n",
    "\n",
    "))\n",
    "\n",
    "# Put 0 - 255 as float 0 - 1\n",
    "for k in colors:\n",
    "    v = colors[k]\n",
    "    _v = [_v / 255.0 for _v in v]\n",
    "    colors[k] = _v\n",
    "    \n",
    "index_colors = [colors[key] if key in colors else \n",
    "                (255, 255, 255, 0) for key in range(0, n+1)]\n",
    "\n",
    "cmap = plt.matplotlib.colors.ListedColormap(index_colors, 'Classification', n+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\user\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:8: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.\n",
      "  \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASQAAAD8CAYAAADe49kaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAcRElEQVR4nO3da6xcZ33v8e9/rTWzL96+bMeXbF/ANnECTguJMeAUaFHTUAiFIJ2iJIrA0CBLLZVCQS0JfdWjSodULamQKsAi9KQohdAQkSgHFFkhSKHnYHIh5OZcnPtOnDiOL3t732bWWv/zYp7tTHb2bM94z+xZdn4fydprPeuZmWcejX961vUxd0dEpAiibjdARGSaAklECkOBJCKFoUASkcJQIIlIYSiQRKQwOhJIZvYxM3vczPaZ2dWd+AwROf1Yu69DMrMYeAK4CBgG7gEud/dH2/pBInLa6cQI6f3APnd/2t0rwI+ASzrwOSJymkk68J5rgRfq1oeBD8ysZGY7gZ1h9b0daIeIFNdBd185s7ATgWSzlL1pv9DddwG7AMxM96+IvLU8N1thJ3bZhoH1devrgJc68DkicprpRCDdA2w2s41mVgYuA27rwOeIyGmm7bts7p6a2V8DdwAx8H13f6TdnyMip5+2n/Y/qUboGJLIW8197r5tZqGu1BaRwlAgiUhhKJBEpDAUSCJSGAokESkMBZKIFIYCSUQKQ4EkIoWhQBKRwlAgiUhhKJBEpDAUSCJSGAokESkMBZKIFIYCSUQKQ4EkIoWhQBKRwlAgiUhhKJBEpDAUSCKyoF554ZmG2xRIIrKg/vdNP224TYEkIgtqYLa5rQMFkogsmO99+7vkcdxwe9snihQRaSSdGCWeYxZGBZKILJinnn+RLG+8z6ZdNhFZMD3lHqYmxhpu1whJRBbEdf/rWoZWrmD4pf0N62iEJCILIjJ47dBrTFUmG9dZwPaIyFvYyNERJqYqVKppwzoKJBFZEP19vaRpRl//QMM6CiQRWRDDLx/g2ReGKZfLDevooLaIdNz5522jt1zi7HdsYt3aNQ3raYQkIh23aNEAURyzZvUqItOV2iLSRZd98mJWrljB+qEzGZ2YalhPgSQiHXXd//xHkjhmbHyMQ0dGGDk22rDuCXfZzGy9md1lZnvN7BEzuyqULzez3Wb2ZPg7GMrNzL5lZvvM7EEz29q2byYip5znX3yJcilh2ZKl9Pf1smb1qoZ1mxkhpcBX3f1+M1sM3Gdmu4HPA3e6+zfM7GrgauBrwMeBzeHfB4Bvh78i8hZ04NARfvJ/fs7ZmzZy1qZN9JUbx84JA8nd9wP7w/Kome0F1gKXAB8J1W4AfkktkC4B/sPdHfi1mS0zs6HwPiLyFrLj8s9yaGSEszdt5H3v3cqqM84gbtfjR8xsA3A+sAdYPR0y7r7fzKbHYWuBF+peNhzK3hBIZrYT2NnK54vIqWXTujW8Z8k5rB86kyWLl1CtVLn/0Qca1m86kMxsAPgJ8GV3HzFr+AiB2Ta86Qko7r4L2BXee44npIjIqWqyUuXM3j4OHnyNOE548ulnWD442LB+U4FkZiVqYXSju98Sil+Z3hUzsyHgQCgfBtbXvXwd8FLL30RETnlTE+P0lBKGDx2ikqZkWUYUN46dZs6yGXA9sNfdv1m36TZgR1jeAdxaV/65cLZtO3BUx49E3nou/x+XsnZoDb+67wEOvnaIl145QG9vH/bmHabjmhkhfRD4LPCQmU3v/H0d+AbwYzO7Enge+EzY9jPgYmAfMA584eS+joicykZHj7H7//4/znn7Ot65eTNpmrLmzNVMVuZxYaS7/4rZjwsBXDhLfQe+1GyjReT09Oqh1/ij928jLpV47sVhNr3tbbgZpTnOsuleNhHpiIEly3joqWepViqsWbmSVw8dwnDWDA01fI1uHRGRtrvik59iaO0a1q5Zw+EjIxw6fJiDI8c45x3v4Oho41tHFEgi0naDS5eQ5s6xY2P09paBiA1vexsO8zvLJiLSqmpUojI5yeixMZKkTKncQ55WWb1yBeeetbHh6xRIItJWZ5+9hf++517++/7fkVarpNUqkcHYyBHGxic5cODVhq/VLpuItFeekwJuxhPPPscj+55m1eBStr773eSVSZ57/rWGL1UgiUhblUpl8jwnq1Z56rnnKZXLeJ4zOnqUX+35DcuXLW34WgWSiLTNti2/z1iWkWcZF37gvZy5ejW/2/sE7zvvPUxUKixeupT169c3fL2OIYlI25R7+8CdUhyz8ozlTFYqvO/891BJU+IkYeP69cQlzToiIgugL3LKy5YwMjHFwcMjZA5T1ZTzf+/3yD1nfGKCt69c0fD1CiQRaZukVKJsEdbr5HnO75+7BdwZHR1hYNEi3ja0GkurDV+vXTYRaYtt52zh2MQkk2mVxX19vHb4KI899gQTk1MsO2MFa9etZ2RiilJPT8P30AhJRNpiSRSR5045z4g8Z8WyJRwbG2P02DEGFw/g7rz7rI1MVNKG76FAEpG2yKoppXJMlPTyzOERvPdVxitVSqUEz53eJGZ0bIwjR3Uvm4h00CfeeS7j1QwyZ3xygqXliFde3s8Z/b08se9JHnp0L+e8YyOVNCOb40iRAklE5m2kmpEYVN2pVlL6IyNJYqbSFHLYc89vuPvXv+b9572H1Wee2fB9dFBbROZt0p0xoGeyyrJKxjhQnUqppE6UZRweHWf1qlUM9Pfz4O/aMOuIiMjs3kkvzhF3BuOYKYMykGU5ZBk9AwO8dPQIyctlDr7yMhs2ndXwnTRCEpF52b7JySOjDyCO6OstU0pzyg6xQWVijE1L+2BslKX9i+iZ43lIGiGJyLyUIuMozooowvOcSXc8jkhj6DGDKGI8KhNPTXB05BBjk2MN30sjJBGZl9RhpUWUSzGexCTVlFLu9JVK5HFENXeiNCXu72OympFY49hRIInIvJQi4xgwledM5U6eQ5YYXk0pT1YpGVQqk4wfG2Ogtx/0CFsR6YSLzt5CnxlrB/qwKGJRNaOc50RTtb9ZZFTTnCwHB46OHSPp72/4fjqGJCInLU+i2sWQE1Pkec54bGQGSRJD7qRZDmZkSZlyHFGZmmRs5EjD99MISUROWgnAnWrupO5gUDLDspw0iRg3qBhQqUCeU+7ppTI22fD9FEgiclIuOOudjFczYq+FUWIRPdRCKnOnUslI0wx3J7eINE3J3Dh8TGfZRKTNeqOIJI7IzegHeiOjN4dKmoPXdtf6o4gcsDwn7unF8pRe84bvqWNIInJSUgPPnRLO0iX9TKQZIxNT9AFHcifNnapFkOYMLuqjL3aII6qmCyNFpM3cIc1z0iTm6NgkcSkhdhiLjMnMSYFlGD09CZ7lTIxMYD0JZHnD91QgiUjLPrT5XWTulMwomzGeZnhWIY1qd/xHBuswIndSDPKMKYPFZnhZ1yGJSBvFDrk75Tgiz3MsNjIgzh1yZw3GFJCZYZUqiUNfHIFDpidGikg7TYVR0JFqSmxGYobnOaU4Zk1slL1WJ3OHUoy5E01WiS2l8SP+NUISkRZddNa7MGrHj3rNyIHIndUDfbx9UR+WxIy4YziWRHg1Y6qaUYoixpKIqjV+b42QRKQlxyIoW0Q5gxRYhLFiUbl2a8jkFJY5ce64QSV3YjN6oojXPGdyKmdxpJtrRaRNktxxh0mHEsaKJX1ESQLVlKicUIoML8e4Gf1ee1jbSJ7XHvQfRYznjc+yNR1IZhab2W/N7PawvtHM9pjZk2Z2k5mVQ3lPWN8Xtm+Y39cXkSIpY5g75cgY95xj1ZRKNYVyicSMqBQxNZVSAQ66sz/LyNyJgAxnMGnPCOkqYG/d+rXAde6+GTgMXBnKrwQOu/tZwHWhnoicBj646RwmqN0q4u6s7uuhH8OiiL5KRmWiwnjmpAZH0gxzJ3EYSGKmgL7coa/xRJFNBZKZrQM+AXwvrBvwx8DNocoNwKfD8iVhnbD9wlBfRE5xpTiiP47oS2IG+3qYiIw4iSjnzkQ1xXPn5akKGZCYYUAV51g1pT+uPcQtmWx8nq3Zg9r/CvwdsDisnwEccffpCwqGgbVheS3wAoC7p2Z2NNQ/WP+GZrYT2Nnk54tIAfQnMQlGAlSznJ4kIptMSc14Natd/NhnxtIoqj1P253cjd7IqBossYhjjW9lO/EIycz+DDjg7vfVF89S1ZvY9nqB+y533+bu207UBhEphjSHNHfKZoxUUqKJKpPAi5NT9EQR68plBsolAA5lOdXc6YljSGKWl0ukSUQWNd5hamaE9EHgU2Z2MdALLKE2YlpmZkkYJa0DXgr1h4H1wLCZJcBS4NDJfX0RKYpPnrUFz3IcOJRm9MRGta8MY1MMxjFLywnjacZIlhGb0WvGhDuRO31VZ0lfzKHRydrzkRo44QjJ3a9x93XuvgG4DPiFu18B3AX8eai2A7g1LN8W1gnbf+HucwzSRORUUAXcnSzPyQzickIpzcmA/jiiUs2YzHLKccTyckJK7T63UhKR9pd4tpqyZOVi+uO44WfM5zqkrwFfMbN91I4RXR/KrwfOCOVfAa6ex2eISAF89qxzca9d5OhRxADGosRYXMnIcqcvzfE8p+zOQA4Hp1LcnYo7scOiDCbGKzz9ylEOVCoNP6elK7Xd/ZfAL8Py08D7Z6kzCXymlfcVkWJ7GaeP2g20uTtVM9KxFPecRWa86DkRkJsR4fQmEVQzVvSVoZrzYqXKQGQkODbrYeYa3ToiIicUhRDKATeYApaaMRXFJLmzHAczUiCL4LA7G/p7SaaqvJBluNVuM+mNozef4aqjQBKREyonMb3ARFaLk8SMgyWjzyIqWUZazSn3lOirpFRyZzlGWs2ohmuRViUxU5lTDZcMNKJ72UTkhEpmvOw5SRjfTGQZFtXu3C9VnX6M8lTKEXc8jlmxcgl5ZDyfpnhkjFYzjuYZE1nG0TRr+DkaIYnInD517ntISjEbj2U8T21Wkf4zBpisZvhklRLOAZx+d/LeEgPlhGdeOYIBQ+WEOHNKiZFGxkieY3OcdFcgicicDJiaqFCOIlLP6ekvU61m9KQ5WZpzFGcpRtoTs6qS88zkJMsxYnemplKiyDjsDuGCyGSOO8msCJcImdko8Hi329GCFcy4Fabg1N7OUntb93Z3XzmzsCgjpMdPpVtIzOxetbdz1N7OKnJ7dVBbRApDgSQihVGUQNrV7Qa0SO3tLLW3swrb3kIc1BYRgeKMkEREFEgiUhxdDyQz+5iZPR5mKen6o0rMbL2Z3WVme83sETO7KpQvN7PdYZaV3WY2GMrNzL4V2v+gmW3tUrtPmVlhzGyZmd1sZo+Ffr6gyP1rZn8TfgsPm9kPzay3SP1rZt83swNm9nBdWcv9aWY7Qv0nzWzHbJ/VcR5mD+jGPyAGngI2UZu+6XfAli63aQjYGpYXA08AW4B/Aq4O5VcD14bli4GfU7ugdTuwp0vt/grwn8DtYf3HwGVh+TvAX4blvwK+E5YvA27qQltvAL4YlsvAsqL2L7VnxD8D9NX16+eL1L/AHwJbgYfrylrqT2A58HT4OxiWBxf8t7HQHzijIy8A7qhbvwa4ppttmqWNtwIXUbuSfCiUDVG7mBPgu8DldfWP11vANq4D7qQ2E8zt4cd2EEhm9jNwB3BBWE5CPVvAti4J/8FtRnkh+5fXJ61YHvrrduBPi9a/wIYZgdRSfwKXA9+tK39DvYX61+1dtuMzlAT1s5d0XRhunw/sAVa7+36A8HdVqFaE7zA9K8z0lKBNzwoDTM8Ks1A2Aa8C/x52Mb9nZosoaP+6+4vAPwPPA/up9dd9FLd/p7Xan0X4HXc9kJqaoaQbzGwA+AnwZXcfmavqLGUL9h06NStMByXUdi++7e7nA2PM/ZjjbvfvILW5BjcCa4BFwMfnaFO3+/dEGrWvEO3udiBNz1AyrX72kq4xsxK1MLrR3W8Jxa+Y2VDYPgQcCOXd/g7Ts8I8C/yI2m7b8VlhZmnT8fZ2aVaYYWDY3feE9ZupBVRR+/dPgGfc/VV3rwK3AH9Acft3Wqv92e1+BrofSPcAm8MZizK1g4C3dbNBYZbd64G97v7Nuk31s6nMnGXlc+HsxXbg6PRQeSH4KTYrjLu/DLxgZueEoguBRylo/1LbVdtuZv3htzHd3kL2b51W+/MO4KNmNhhGhR8NZQtroQ9azXIw7mJqZ7KeAv6+AO35ELWh6oPAA+HfxdSOA9wJPBn+Lg/1Dfi30P6HgG1dbPtHeP0s2ybgN8A+4L+AnlDeG9b3he2butDO84B7Qx//lNpZncL2L/APwGPAw8APgJ4i9S/wQ2rHt6rURjpXnkx/An8R2r0P+EI3fsO6dURECqMju2xWsIsdReTU0PYRkpnF1HbBLqI2fLyH2vUMj7b1g0TktNOJEdL7gX3u/rS7V6id+bmkA58jIqeZTjzCdrYLrD4ws5KZ7QR2ApRj3ruqv/F83yJyehkezQ76Aj1Tu6kLrNx9F+FBUeuXJH7V9mUdaIqIFNHf7n7tudnKO7HLVogLrETk1NOJQCrcxY4icmpo+y6bu6dm9tfUrvKMge+7+yPt/hwROf10ZF42d/8Z8LNOvLeInL66fS+biMhxCiQRKQwFkogUhgJJRApDgSQihaFAEpHCUCCJSGEokESkMBRIIlIYCiQRKQwFkogUhgJJRApDgSQihaFAEpHCUCCJSGEokESkMBRIIlIYCiQRKQwFkogUhgJJRBbUI/d/teE2BZKIFIYCSUQWzFyjI1AgiUiBKJBEpDAUSCJSGAokESkMBZKILIgTHdAGBZKIFIgCSUQKQ4EkIoWhQBKRwlAgiUhhKJBEpONuvPTupuopkESk48772k+bqqdAEpGOaub6o2knDCQzW29md5nZXjN7xMyuCuXLzWy3mT0Z/g6GcjOzb5nZPjN70My2nvQ3EZG3lGZGSCnwVXd/F7Ad+JKZbQGuBu50983AnWEd4OPA5vBvJ/DttrdaRE4pD1z76abqnTCQ3H2/u98flkeBvcBa4BLghlDtBmD6Ey8B/sNrfg0sM7Oh1povIqeD6d21jhxDMrMNwPnAHmC1u++HWmgBq0K1tcALdS8bDmUz32unmd1rZvceq3orzRCR01TTgWRmA8BPgC+7+8hcVWcpe1PiuPsud9/m7tsGSrO9RETeapoKJDMrUQujG939llD8yvSuWPh7IJQPA+vrXr4OeKk9zRWRU0UrZ9emNXOWzYDrgb3u/s26TbcBO8LyDuDWuvLPhbNt24Gj07t2IiJzSZqo80Hgs8BDZvZAKPs68A3gx2Z2JfA88Jmw7WfAxcA+YBz4QltbLCKnhOkza80e0IYmAsndf8Xsx4UALpylvgNfaroFInJauuKmD7e826YrtUWkI2689O6mrz+a1swum4hISx65/6ucd3/zu2rTNEISkcJQIIlIYSiQRKTtzt36L5y79V9afp0CSUQ65oFrP93SgW0Fkoh0zBU3fZgrbvpw06Gks2wi0hH11yA1e8ZNIyQRKQwFkoi0TbMP829Eu2wi0lbToaQLI0Wkq6646cPzer0CSUQKQ4EkIm0zvbs2faq/1ZtrFUgi0nY3Xnr38d23VkJJB7VFpG3qjyHVn3Fr9mFtCiQR6YiZ4dTMkyO1yyYibVc/OprefWvmmJICSUTaqv740XQwHb82SbtsIlIUx0dIX/z6rNs1QhKRtpo5OmrlYkmNkESk7WbbbWuGRkgiMm8zQ6c+jKafiTRtruNICiQRmbeZu2Uzd9dmuyZpNgokEWm7K2768BtCqNnjSAokEemImSE0vT5XOCmQRGRBzLwmaTYKJBGZtxsvvXvOoGn2TJsCSUTmba7dsFauR1IgiUhbzPdpkaALI0Wkw1q5QFKBJCILon4E9UCDOgokEWmLk7nuaCYdQxKReWtHGIECSUTmab6TQ9ZrOpDMLDaz35rZ7WF9o5ntMbMnzewmMyuH8p6wvi9s39C21opI4c0noFoZIV0F7K1bvxa4zt03A4eBK0P5lcBhdz8LuC7UE5HTUKO7/Ft5Tb2mAsnM1gGfAL4X1g34Y+DmUOUGYPoW3kvCOmH7haG+iJxGGoXRiUKpHfey/Svwd0Ae1s8Ajrh7GtaHgbVheS3wAkDYfjTUfwMz22lm95rZvceq3mQzRKQI6h/A1urr5nLC0/5m9mfAAXe/z8w+Ml08S1VvYtvrBe67gF0A65ckSiSRU9BsodQorN5wc+3uLbO+XzMjpA8CnzKzZ4EfUdtV+1dgmZlNB9o64KWwPAysBwjblwKHmvgcETlFzHb19YnuWTvvaz9909MjZzphILn7Ne6+zt03AJcBv3D3K4C7gD8P1XYAt4bl28I6Yfsv3F0jIJHTzFwHtGduu/HSu3ng2k+f8KkA87lS+2vAj8zsH4HfAteH8uuBH5jZPmojo8vm8RkiUmCNRjuzlTdz64gVYfCyfkniV21f1u1miEibNTqe9Le7X7vP3bfNLNeV2iLSESdzJk6BJCInZebxoJO5SHImBZKItMXMg9oncwuJHj8iInOa7XR+o1P8JzN9dj0Fkoic0Fyn9OvL5/sY20KcZTOzUeDxbrejBSuAg91uRAvU3s5Se1v3dndfObOwKCOkx2c7BVhUZnav2ts5am9nFbm9OqgtIoWhQBKRwihKIO3qdgNapPZ2ltrbWYVtbyEOaouIQHFGSCIiCiQRKY6uB5KZfczMHg+zlFxdgPasN7O7zGyvmT1iZleF8uVmtjvMsrLbzAZDuZnZt0L7HzSzrV1q9ykzK4yZLTOzm83ssdDPFxS5f83sb8Jv4WEz+6GZ9Rapf83s+2Z2wMweritruT/NbEeo/6SZ7ZjtszrO3bv2D4iBp4BNQBn4HbCly20aAraG5cXAE8AW4J+Aq0P51cC1Yfli4OfUHt27HdjTpXZ/BfhP4Paw/mPgsrD8HeAvw/JfAd8Jy5cBN3WhrTcAXwzLZWBZUfuX2jPinwH66vr180XqX+APga3Aw3VlLfUnsBx4OvwdDMuDC/7bWOgPnNGRFwB31K1fA1zTzTbN0sZbgYuoXUk+FMqGqF3MCfBd4PK6+sfrLWAb1wF3Unu88O3hx3YQSGb2M3AHcEFYTkI9W8C2Lgn/wW1GeSH7l9cnrVge+ut24E+L1r/AhhmB1FJ/ApcD360rf0O9hfrX7V224zOUBPWzl3RdGG6fD+wBVrv7foDwd1WoVoTv0PZZYTpoE/Aq8O9hF/N7ZraIgvavu78I/DPwPLCfWn/dR3H7d1qr/VmE33HXA6mpGUq6wcwGgJ8AX3b3kbmqzlK2YN+hflaY+uJZqrY0K0wHJdR2L77t7ucDY9R2KRrpdv8OUptrcCOwBlgEfHyONnW7f0+kUfsK0e5uB9LxGUqC+tlLusbMStTC6EZ3vyUUv2JmQ2H7EHAglHf7O5xqs8IMA8Puvies30wtoIrav38CPOPur7p7FbgF+AOK27/TWu3Pbvcz0P1AugfYHM5YlKkdBLytmw0yM6M2UcFed/9m3ab62VRmzrLyuXD2YjtwdHqovBD8FJsVxt1fBl4ws3NC0YXAoxS0f6ntqm03s/7w25hubyH7t06r/XkH8FEzGwyjwo+GsoW10AetZjkYdzG1M1lPAX9fgPZ8iNpQ9UFqkyM8ENp4BrUDx0+Gv8tDfQP+LbT/IWBbF9v+EV4/y7YJ+A2wD/gvoCeU94b1fWH7pi608zzg3tDHP6V2Vqew/Qv8A/AY8DDwA6CnSP0L/JDa8a0qtZHOlSfTn8BfhHbvA77Qjd+wbh0RkcLo9i6biMhxCiQRKQwFkogUhgJJRApDgSQihaFAEpHCUCCJSGH8f0zH6k9ozdwqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(2,1,figsize=(12,4))\n",
    "\n",
    "img_stretched = color_stretch(reshaped_img, [4, 3, 2])\n",
    "axs[0].imshow(img_stretched)\n",
    "\n",
    "axs[1].imshow(class_prediction, cmap=cmap, interpolation='none')\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}