3. Vectorization and Postprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import rasterio\n",
    "from rasterio.windows import Window\n",
    "\n",
    "import cv2 \n",
    "import geopandas as gpd\n",
    "\n",
    "from shapely.geometry import MultiPolygon, Polygon\n",
    "from shapely.ops import cascaded_union\n",
    "from collections import defaultdict\n",
    "\n",
    "from pathlib import Path\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creating polygons from pixel masks\n",
    "# from https://www.kaggle.com/lopuhin/full-pipeline-demo-poly-pixels-ml-poly\n",
    "\n",
    "def mask_to_polygons(mask, epsilon=10, min_area=50.):\n",
    "    # first, find contours with cv2: it's much faster than shapely\n",
    "    image, contours, hierarchy = cv2.findContours(\n",
    "        ((mask == 1) * 255).astype(np.uint8),\n",
    "        cv2.RETR_CCOMP, cv2.CHAIN_APPROX_TC89_KCOS)\n",
    "    # create approximate contours to have reasonable submission size\n",
    "#     epsi = lambda x: epsilon*cv2.arcLength(x,True)\n",
    "    approx_contours = [cv2.approxPolyDP(cnt, epsilon, True)\n",
    "                       for cnt in contours]\n",
    "    if not contours:\n",
    "        return MultiPolygon()\n",
    "    # now messy stuff to associate parent and child contours\n",
    "    cnt_children = defaultdict(list)\n",
    "    child_contours = set()\n",
    "    assert hierarchy.shape[0] == 1\n",
    "    # http://docs.opencv.org/3.1.0/d9/d8b/tutorial_py_contours_hierarchy.html\n",
    "    for idx, (_, _, _, parent_idx) in enumerate(hierarchy[0]):\n",
    "        if parent_idx != -1:\n",
    "            child_contours.add(idx)\n",
    "            cnt_children[parent_idx].append(approx_contours[idx])\n",
    "    # create actual polygons filtering by area (removes artifacts)\n",
    "    all_polygons = []\n",
    "    for idx, cnt in enumerate(approx_contours):\n",
    "        if idx not in child_contours and cv2.contourArea(cnt) >= min_area:\n",
    "            assert cnt.shape[1] == 1\n",
    "            poly = Polygon(\n",
    "                shell=cnt[:, 0, :],\n",
    "                holes=[c[:, 0, :] for c in cnt_children.get(idx, [])\n",
    "                       if cv2.contourArea(c) >= min_area])\n",
    "            all_polygons.append(poly)\n",
    "    # approximating polygons might have created invalid ones, fix them\n",
    "    all_polygons = MultiPolygon(all_polygons)\n",
    "    if not all_polygons.is_valid:\n",
    "        all_polygons = all_polygons.buffer(0)\n",
    "        # Sometimes buffer() converts a simple Multipolygon to just a Polygon,\n",
    "        # need to keep it a Multi throughout\n",
    "        if all_polygons.type == 'Polygon':\n",
    "            all_polygons = MultiPolygon([all_polygons])\n",
    "    return all_polygons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "INPUT = Path('/media/splash/00F2D259F2D25310/Data/final/outputs')\n",
    "COG_URL = Path('/media/splash/00F2D259F2D25310/Data/final/arid.tif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## If needed, reproject to EPSG:4326"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "src_tif = 'arid1x_merged.tif'\n",
    "dst_tif = 'arid1x_merged_reproj.tif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'driver': 'GTiff',\n",
       " 'dtype': 'uint8',\n",
       " 'nodata': None,\n",
       " 'width': 6431,\n",
       " 'height': 6503,\n",
       " 'count': 3,\n",
       " 'crs': CRS.from_epsg(3785),\n",
       " 'transform': Affine(0.14931428515517545, 0.0, 8135024.072190155,\n",
       "        0.0, -0.14931402813080957, 3982066.4113659593)}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raster = rasterio.open(INPUT/src_tif,'r')\n",
    "raster.meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://rasterio.readthedocs.io/en/latest/topics/reproject.html\n",
    "\n",
    "from rasterio.warp import calculate_default_transform, reproject, Resampling\n",
    "\n",
    "dst_crs = 'EPSG:4326'\n",
    "\n",
    "with rasterio.open(INPUT/src_tif) as src:\n",
    "    transform, width, height = calculate_default_transform(\n",
    "        src.crs, dst_crs, src.width, src.height, *src.bounds)\n",
    "    kwargs = src.meta.copy()\n",
    "    kwargs.update({\n",
    "        'crs': dst_crs,\n",
    "        'transform': transform,\n",
    "        'width': width,\n",
    "        'height': height,\n",
    "        'compress': 'JPEG',\n",
    "        'tiled': True, \n",
    "        'blockxsize': 256, \n",
    "        'blockysize': 256\n",
    "    })\n",
    "\n",
    "    with rasterio.open(INPUT/dst_tif, 'w', **kwargs) as dst:\n",
    "        for i in range(1, src.count + 1):\n",
    "            reproject(\n",
    "                source=rasterio.band(src, i),\n",
    "                destination=rasterio.band(dst, i),\n",
    "                src_transform=src.transform,\n",
    "                src_crs=src.crs,\n",
    "                dst_transform=transform,\n",
    "                dst_crs=dst_crs,\n",
    "                resampling=Resampling.nearest,\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or with gdalwarp on CLI:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gdalwarp -co compress=JPEG -co PHOTOMETRIC=YCBCR -co TILED=YES -co \"BLOCKXSIZE=256\" -co \"BLOCKYSIZE=256\" -s_srs EPSG:32737 -t_srs EPSG:4326 {INPUT}/{src_tif} {INPUT}/{dst_tif}   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Open and Read Windows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "dst_tif = 'arid1x_merged_reproj.tif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster = rasterio.open(INPUT/dst_tif,'r')\n",
    "data = raster.read()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'driver': 'GTiff',\n",
       " 'dtype': 'uint8',\n",
       " 'nodata': None,\n",
       " 'width': 6997,\n",
       " 'height': 5890,\n",
       " 'count': 3,\n",
       " 'crs': CRS.from_epsg(4326),\n",
       " 'transform': Affine(1.2328306174800615e-06, 0.0, 73.07816460728642,\n",
       "        0.0, -1.2328306174800615e-06, 33.651230626615664)}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raster.meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(256, 256), (256, 256), (256, 256)]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raster.block_shapes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pad_window(window, pad):\n",
    "    col_off, row_off, width, height = window.flatten()\n",
    "    return Window(col_off-pad//2, row_off-pad//2,width+pad,height+pad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# padded windowed reads with blocks and rasterio: https://rasterio.readthedocs.io/en/latest/topics/windowed-rw.html\n",
    "\n",
    "pad = 64\n",
    "windows = []\n",
    "for ji, window in raster.block_windows(1):\n",
    "    assert len(set(raster.block_shapes)) == 1\n",
    "    window_padded = pad_window(window, pad)\n",
    "    windows.append((ji, window_padded))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "672"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(windows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[((5, 10), Window(col_off=2528, row_off=1248, width=320, height=320)),\n",
       " ((5, 11), Window(col_off=2784, row_off=1248, width=320, height=320)),\n",
       " ((5, 12), Window(col_off=3040, row_off=1248, width=320, height=320)),\n",
       " ((5, 13), Window(col_off=3296, row_off=1248, width=320, height=320)),\n",
       " ((5, 14), Window(col_off=3552, row_off=1248, width=320, height=320)),\n",
       " ((5, 15), Window(col_off=3808, row_off=1248, width=320, height=320)),\n",
       " ((5, 16), Window(col_off=4064, row_off=1248, width=320, height=320)),\n",
       " ((5, 17), Window(col_off=4320, row_off=1248, width=320, height=320)),\n",
       " ((5, 18), Window(col_off=4576, row_off=1248, width=320, height=320)),\n",
       " ((5, 19), Window(col_off=4832, row_off=1248, width=320, height=320))]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "windows[150:160]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def poly2coords(poly_mask, transform):\n",
    "    geo_coords = [(coord[0],coord[1])*transform for coord in poly[0].exterior.coords for poly in poly_mask]\n",
    "    return Polygon(geo_coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 672/672 [00:01<00:00, 662.94it/s]\n"
     ]
    }
   ],
   "source": [
    "mask_thres = 0.7\n",
    "epsilon = 10\n",
    "min_area = 100.\n",
    "erode = 5\n",
    "dilate = 7\n",
    "\n",
    "polys = []\n",
    "\n",
    "for window in tqdm(windows):\n",
    "    win_tnfm = rasterio.windows.transform(window[1], raster.meta['transform'])\n",
    "    win_img = raster.read(window=window[1])[0,:,:]\n",
    "\n",
    "    mask = win_img > (mask_thres*255)\n",
    "    # erode and dilate\n",
    "    mask = cv2.erode(mask.astype('uint8'), np.ones((erode,erode),np.uint8), iterations=1)\n",
    "    mask = cv2.dilate(mask.astype('uint8'), np.ones((dilate,dilate),np.uint8), iterations=1)\n",
    "\n",
    "    # label via connected components\n",
    "    _, instances = cv2.connectedComponents(mask.astype('uint8'))\n",
    "\n",
    "    # make polys from instances\n",
    "    uniques = list(np.unique(instances))\n",
    "    for b in uniques[1:]:\n",
    "        poly = mask_to_polygons(instances==b,epsilon, min_area)\n",
    "        try: \n",
    "            if poly.type == 'MultiPolygon': \n",
    "                geo_coords = poly2coords(poly, win_tnfm)\n",
    "                polys.append(geo_coords)\n",
    "            #else: print('not a MultiPolygon')\n",
    "        except Exception as exc: print(f\"{exc}: {window}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"73.08036654139569 33.650764715268714 5.281446365756892e-05 0.00019705564589145297\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,67.30172648618333)\"><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"3.941112917829059e-06\" opacity=\"0.6\" d=\"M 73.08041205750209,33.65095447255735 L 73.08039479787344,33.650943377081795 L 73.08037383975295,33.65083858647931 L 73.08038740088973,33.65077201362597 L 73.08041205750209,33.65078434193214 L 73.08041205750209,33.65095447255735 z\" /></g></svg>"
      ],
      "text/plain": [
       "<shapely.geometry.polygon.Polygon at 0x7f40d4799208>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polys[5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "357"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(polys)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dedupe windowed polys\n",
    "merged_poly = gpd.GeoDataFrame()\n",
    "merged_poly['geometry'] = gpd.GeoSeries(cascaded_union(polys))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f40d4835860>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAIICAYAAADE0iKgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3zb1b3/8deR5D3jbcdJnO0sZwdCCCOkEGbLKlBKGeXSlvaWtlxK2193772lu6WTfVsopUBL2WUTAmTvvZfjON57Sjq/P6yYJJZtOZYt2X4/Hw8/sM53fQRY/viMzzHWWkREREQk9ByhDkBERERE2igxExEREQkTSsxEREREwoQSMxEREZEwocRMREREJEy4Qh1AMKSlpdm8vLxQhyEiIiLSrbVr15ZZa9P9HRsUiVleXh5r1qwJdRgiIiIi3TLGHOzsmIYyRURERMKEEjMRERGRMKHETERERCRMKDETERERCRNKzERERETChBIzERERkTChxExEREQkTCgxExEREQkTSsxEREREwoQSMxEREZEwocRMREREJEwoMRMREREJE0rMRERERMKEEjMRERGRMKHETERERCRMKDETERERCRNKzERERETChBIzERERkTChxExEREQkTCgxExEREQkTSsxEgqCktondx2pDHYaIiAxwrlAHIDKQVda3cNOjK9lypIaZI5N57s4FoQ5JREQGMPWYifTCd57fwpYjNQCsP1RFUVVjiCMSEZGBTImZyGl6cWMRL2062v7aYSDCqR8pERE5ffotInIaSmqa+M7zW05qmzc6hfSEqBBFJCIig4ESM5HT8M1/bqaqofWktkumZYcoGhERGSy6TcyMMdHGmFXGmI3GmK3GmB/42h/xtW0yxjxrjInv5PoCY8xy37WbjTHRvvZ3jTE7jTEbfF8ZvvZzjDHrjDFuY8w1wXyzIsFQWd/CWztKTmozBi6akhWiiEREZLAIZFVmM7DIWltnjIkA3jfGvAp81VpbA2CM+SXwJeC+Ey80xriAJ4CbrLUbjTGpwIndDDdaa9ec8rxDwC3Af53OGxLpa5uPVHdomzNqGJmJ0SGIRkREBpNuEzNrrQXqfC8jfF/2hKTMADGA9XP5hcAma+1G373KA3jeAd99vQHEL9LvthR1TMwunqphTBER6b2A5pgZY5zGmA1ACfCGtXalr/0xoBjIB37r59IJgDXGvOYbnvz6Kccf8w1jfseX4ImEvS1+eswunqZhTBER6b2AEjNrrcdaOwPIBeYZY6b62m8FcoDtwHV+LnUBZwM3+v55pTHmAt+xG62104CFvq+behK4MeYOY8waY8ya0tLSnlwq0iunDmXOGplMdlJMiKJps6O4hte3FrN0V6lqqYmIDGA9qvxvra0yxrwLLAG2+No8xpi/A/cAj51ySSGw1FpbBmCMeQWYBbxlrT3iu77WGPMkMA/4Sw9ieRB4EGDOnDn+hlFFgq6qoYXDFScnPqFcjdnY4uEXr+/k0Q/24/X9FCyelMnDN88JWUwiInL6AlmVmW6MSfZ9HwMsBnYaY8b52gxwObDDz+WvAQXGmFjfQoBzgW3GGJcxJs13fQRwGb5ETyScHa/yf6Klu0ppavWEIBr4yb938PD7HyVlAG/vOKZeMxGRASqQocxs4B1jzCZgNfAG8DLwZ2PMZmCz75wfAhhjrjDG/BDAWlsJ/NJ33QZgnbX2ZSAKeM13zw3AEeAh3/VzjTGFwLXAA8aYrcF6syK95W/i/77SeqIjnNQ3u6mob+m3WPaW1vHkykMd2r0W/r76cL/FISIiwRPIqsxNwEw/h/zu1mytfQF44YTXT9BWMuPEc+qB2Z1cv5q2uWwiYcdf4jV9RBK1Ta08+v4Bqhtb+fy5Y8joh9IZr24+ivWzGNoYKKxUj5mIyEDUozlmIkNddlLHhKsgN5mE6AjuWjy+X2O5aEoWP399V4d2a+G2s/P6NRYREQkOJWYiPbC5sONQZkyEk61F1UQ4HTgdhmGxkaTERfZ5LF0NV76wsYgpOUl9HoOIiASXEjORHvBX9f97L5w8DfKOc8bwrUsm9WkczW4P/1hX2OnxlzYe5d6L8nE4VB5QRGQg0SbmIgFqaHGzt7SuQ/uk7MSTXrv6IRl6Y9sxKk/ZRP1ER6oaWX+4ss/jEBGR4FJiJhKgvSX1J5WlyEmK5t9fWcjL/3k2V80c3t7ucvb9j1VXw5hp8VH86BNTKchN7vM4REQkuDSUKRKg+OiTf1xuO3s0+VltvWU/uaaAs8encaiigen9kBDtK63v0BYX6eTWBaP5zFmjyEjQhuoiIgOREjORAGUmRrV/Hx/l4rq5I9pfRzgdXDWr/6q85GclcMRXRDbCabjxjFF8adE40uKjurlSRETCmRIzkQDFRrpIiHZR2+TmtrNHkxAdEbJYLpmWzd7SOubkpfDlReMZmRobslhERCR4lJiJ9MB3L5vMqNQ45o1OCWkcV8/O5erZqsMsIjLYKDET6YFr54zo/iQREZHTpFWZIiIiImFCiZmIiIhImFBiJiIiIhImlJiJiIiIhAklZiIiIiJhQomZiIiISJhQYiYiIiISJpSYiYiIiIQJJWYiIiIiYUKJmYiIiEiYUGImIiIiEiaUmImIiIiECSVmIiIiImFCiZmIiIhImFBiJiIiIhImlJiJiIiIhAlXqAMQEQnUg+/tpbKhlXl5KZyfnxHqcEREgk6JmYgMGP9cd4QdxbUcmd6oxExEBiUNZYrIgPH9K6YwJj2O0trmUIciItIn1GMmIgPGmWNSefWuhby08WioQxER6RPqMRORASXK5eTq2bmhDkNEpE8oMRMREREJExrKFJF+VVzdxMr95eQOi2Xq8ESiXM5QhyQiEjbUYyYi/crt9XLXUxu4+o8fsqekLtThiIiEFSVmItKvhifHEBfZ1kvW4vaGOBoRkfCioUwZ0t7YdoyC3CQyE6NDHcqQYYzh4ZvnsqO4hrT4qFCHIyISVpSYyZBVUtPEPc9upLbJzbcvncStC0aHOqQhY/7YVOaPTQ11GCIiYUeJmQxZZXUteDwWr7WcPS4toGueXnMYgDNHpzIyNbYvwxMRkSFIiZkMWZNzEln5/y5g/aEqxmcmBHTN61uP8eb2Y3zu3DF88+JJfRyhiIgMNUrMZEiLjXSxIMDeMoCHPjOb+haPJq2LiEifUGIm0gPGGOKjXKA56yIi0geUmIkE0dJdpby0sYjs5BhykqLJSoomJzmG7KRoEqIjQh2eiIiEOSVmIkHU4vbyzNpCv8fio1xkJ0Xz/Sum9Gj4VEREhg4VmBUJouykzuuh1TW72V1Sx5j0uH6MSEREBhIlZiJBtGx3WZfHx2fEk50U00/RiIjIQKPETCRINhVW8YvXd3Z63GHgmtm5/RiRiIgMNJpjJhIE9c1u7npqA26v7XAsOsLBJ+eM4Pazx6go7RBgreVIVSNHq5soyE0iyuUMdUgiMoAoMRMJgh++uI39ZfUd2q+elcu3LsknVXtCDlr1zW42FVaz/nAl6w9VseFwFaW1zQAkxURwaUE2V84czuyRw3A4TIijFZFwp8RMpJd2FNfwd99WTae6+axRSsoGqbd3HONnr+1iZ3ENfjpKAahubOXJlYd4cuUhcofF8K1LJnHJtOz+DVREBhTNMRPppWM1zX7bI5yGiVmBbfUkA09pbTPbj3aelJ2qsLKRh5ft69ugRGTAU2Im0ktVDS1+2ydmJWh+0SA2LDayx9d8bHJWH0QiIoOJhjJFeqm6sdVv+7ThSf0ciQTiJ//ewfPrj9DisUS5HHz27NHcdvboHt9nck4inzt3DMXVTewtrWPLkZr2YwlRLhZOSGN4cgw5vq8Rw2KZnJMYzLciIoOQEjORXqpq8J+YTVViFlbcHi+HKhr495Ziiqqb2ttf2lR0WolZ7rBYvnnxJAAaWzwcq2kiKymaumY39c1uRqWqkLCI9Fy3iZkxJhp4j7Ztm13As9ba7xljHgHmAAbYBdxira3zc30B8ACQCHiBudbaJmPMu0A20Og79UJrbYkxJgr4CzAbKAeus9Ye6NW7FOlD/hKzacOTWDwpMwTRyKmstXzn+S08vbqQFo+3w/GNhdXUNLWS2Iu9TGMineSltSVi0RFO0rTgQ0ROUyA9Zs3AImttnTEmAnjfGPMq8FVrbQ2AMeaXwJeA+0680BjjAp4AbrLWbjTGpAIn/ha70Vq75pTnfRaotNaOM8ZcD/wEuO503pxIf6hqPHmO2cTMBP5551lEODWFMxzc9+oOnlhxqEP7sNgIcofFEuE07C+tZ/qI5BBEJyJysm4TM2utBY73hEX4vuwJSZkBYgB/a5MuBDZZazf67lUeQEwfB77v+/5Z4HfGGOOLQyTsVJ/SY1Ze36ykLEyU1DTxwHv+V0LOHpXCwzfP6eeIRES6FtBvD2OM0xizASgB3rDWrvS1PwYUA/nAb/1cOgGwxpjXjDHrjDFfP+X4Y8aYDcaY7/gSPIDhwGEAa60bqAZSe/rGRPqDtZZjtU0ntaUnRNPs9oQooqHHWsuuY7W8uvkotU0nJ8npCVFkJHQcVoxwGs4ep48VEQk/AU3+t9Z6gBnGmGTgOWPMVGvtFmvtrcYYJ21J2XXAY37ufzYwF2gA3jLGrLXWvkXbMOYRY0wC8A/gJtrmlvkrjd2ht8wYcwdwB8DIkSMDeRsiQWeM4dGb57KlqJo1ByqZnJPIJVOzVeG9HxyuaODxFQd5fWsxB8obAIh0Ojh7fBq3nJXHORPSMcZw64LR/PnDAxTXtCXQY9Pj+M31M7U4Q0TCkunpCKEx5ntAvbX25ye0nQvcY6297JRzrweWWGtv8b3+DtBkrf3ZKefdAsyx1n7JGPMa8H1r7XLfHLViIL2rocw5c+bYNWtOnaomIoPZT/+9gz+8u9fvsTmjhvHsF85qf13X7GZTYRUzRwzDmLYJ+iIioeLrpPI7l6LboUxjTLqvpwxjTAywGNhpjBnnazPA5cAOP5e/BhQYY2J9Sda5wDZjjMsYk+a7PgK4DNjiu+YF4Gbf99cAb2t+mUjf+3BPGf9cV8iaAxWU1TUT7j92DtN5r2TyKcVf46NcnDU2jZhIp5IyEQlrgQxlZgN/9g1ZOoCngZeBZcaYRNqGHjcCXwAwxlxBW+/Xd621lb4Vm6tpG458xVr7sjEmDnjNl5Q5gTeBh3zPewR43BizB6gArg/SexWRTny4p4xPP7LypO2FEqNd3HPRRD595ihMF0lQqFR0suMCwNc+NqEfIxERCZ5AVmVuAmb6ObSgk/NfoK3X6/jrJ2grmXHiOfW01Snzd30TcG13cYlI8DS2ejrs+VjT5OY7z29l6a4yfnpNASlxPd+CqC8VVjb6bb+0IFsV9kVkwNKafhHh0Q/2d3rsze3HeGCp/7lcoWKtZVNhVYd2h4GvLh4fgohERIJDWzKJDHHbj9bwwZ6uSwxmJEb3UzRtrLV4LXitxWst1oLH2/Z9Y4uHtQcr/e648ImZwxmXkdCvsYqIBJMSM5Eh7rEuesuOy07qPjHbV1qHx2uJi3IRH+0iLtKFM4CyIU+tOsTBigbuXZIPwBeeWMurW4q7D/wULofhrgvUWyYiA5sSM5EhrKyumX9tKOr2vKwAErMP9pTxnee3ntQWE+HkoimZ/PeV04iP8v9xs+FwFftK62lq9RDlcrBqf0VgwZ/i2jkjtHG4iAx4mmMmMoQ9ufIQLe6OG3ufKiuAocxr54zoUGW/sdXDvzYUcflv32drUbXf6+67uoCnPz+f6Agndc1urpiRw7zRKYG9gRO0ery8u7OE+mZ3j68VEQkXPS4wG45UYFak55rdHhbc9w5ldc1dnhfpdLD1hxcFtP/nw8v28d8vb/d/H5eDb12czwWTMskdFtNlCY5XNh/lzr+u6/Z5/rgchoLcJOaPTWX+mDRmjxpGTKRql4lI+OiqwKyGMkWGqNe2Hus2KQP4+pKJAW/KfuMZo/jju3spr+9YY6zF7eX7L27j+y9uIykmgm9enM/18/xvp9ab0hxur2XdoSrWHari9+/sJdLpYMaIZM4cm8r8ManMzRuGawBsMv/ge3tZvrecRfkZ3DQ/L9ThiEg/Cf9PJxHpE4vyM5g/puuNvM+ZkM5tC0YHfM+YSCe3LxzT7XnVja3MGjWs0+OTcxL54cencP3cEUwbnkRkLxKpFo+XVQcquP+t3dzw0AqeXHXotO/VX/aU1PLz13fxzs5S3txeEupwRKQfqcdMZIiKj3Lx5H+cQW2zm+Lqpo++apo4Wt1EWV0z//2JKT3ekP2m+aP4/Tt7qOtirteMEclMyOy8rEVidASfOaGXqNXjZWdxLTuLa9l2tIatRdVsK6rhW5dMYsPhKp5afTjg+H7z5m6unDmchOiIgK/pT26Pl7uf3tg+929cRnyIIxKR/qTETGQIM8aQGB1BYnREh0TpcEUDdz+9ifL6Fj4+I4ePz8ghOymm23vGR3VfJuO6uSPav39g6V6W7irl3iX5TB+R7Pf8CKeDqcOTmDo8iat9bcfnx2443LHQbFfy0uKorG8N28TsrysPsbGwbaHEdXNG8GWVABEZUpSYiUgHNU2tXPG796n0FXHdfrSGn/x7B2eOTuUPN85iWBdzwJrdHqobOxZ/PS4mwsllBdlAW1L15vZj/Piqgh73DB1fPPDB3rKAzh+dFse9S/K5aEpmWO79eVyMb5P1/7lyKjeeMSrE0YhIf9McM5EhZu3BCt7YdoymVk+n5yRGR/DNiyed1GYtrDpQ0W1B2vK6zjcXB7hkWjYJ0RF4vZbyumaeumP+aQ/XHa5o4HCF/z0zT+R0GP56+xksmZoV1kkZwFxfqZCcAHonu1NS20RJbROV9S00tKiMiMhAoB4zkSHmiRWHeG79EeIinVwwKZOLp2Zxfn4G0REnl5T45NwR5CTHEOlykBofSWpcJInREd3OOatpaiUtPoqfXN1WVPZgRQMHy+s5WN5Ai9vLLWflAeBwGC6YlNmr97J8b9dbSbkcBo+1XDw1i5zk3ic6/SEvNZYpOYkEI3+884l1rDlY2f76ypnDueeiiQPm34XIUKQ6ZiJDSH2zm3N/9m6HMhn3XTWt09IVfe29XaXUNbs5d0I6cZ3sDuCPtZY7Hl/LG9uOdXrO7z81i+kjkvB6YWRqbDDCHTBaPV4+9dAKVh+oPKk9yuXgGxfnc2sPVtuKSHCpjpmIAPDs2kK/tcsefn8/6w5VkpkYTUZiNBkJUbR6vEzITOhy9WQwvLqlmL+tOsTjn53HwvHpJx1ravXwh3f3suZABfffMJO0+I92FnhixcFOk7KUuEh+dk0Bi/Izwn7o8kQNLW5aPZakmN4vTIhwOkiK6TgXsNnt5QcvbqPZ7eXz547t9XNEJLiUmIkMEdZa/tZJDa89JXXsKanr0G4MnDE6hbsvnMjcvJ5vk9Sdf6wtpKqhhfgoF02tXlbtr6Cx1UNDs5vKhlYeeG8vB8sbAPj47z7gwc/MZkpOEgDjO0kYhyfH8Phn5zEmfeCVmXh1czF3P7ORCZnxzB41jLsvnHhSMtpT9y6ZSGZiFIWVjRyraSuFUuVb0HHfqzuIiXBys29oWUTCgxIzkUGqor6FYbER7T1Gq/ZXsKO4tkf3sBZW7Kvg689u4p3/Oi/oMSbHRvDqlmIA/uMvXU9HOFLVyDV/XM6vrpvOkqnZnDkmlWtm5/Ls2kIcBubkpbBwXBqfnDuCzAD29gxHK/e3zZnbdayOXcfquPmsvF4lZnXNbtwey/kT07llwWiW7y3nhodWtB8vqW3qdcwiElxKzEQGmaqGFv7w7l7+78MDXDw1i0+fOYo/f3igPQHqCWMgPyuRCZl90/t0waRMLsjP4K0dgVW3b2z18Pkn1rHs6+czIiWWn1xdwJcXtdX5GgxzyFburwjq/Z5ceYhn1hZy5czheL2Wh5ftaz8W6XJwWUFOUJ8nIr2nxExkEHn0/f38+s1d1DS1lUZ4fkMRz28o6vF98lJjuWl+HgvGpZKflRjsME/y8ZnDA07MjnM523oBnQ4zKBIyaJtPlxgdQYTT0OqxZCZGkd6L3jKAablJPLO2kCOVjfxx6d6T/j1/Y0k+k7L79r+tiPScEjORQeTDvWXtSdnpumHeSL5xcT4NLe6AKv331oZDPavcD4Rt1f7eiI5w8uJ/nk2z28Mj7+8nNS6S1F4mZosnZfLd57eytaiaNQc/6o07d0I6ty7I62XEItIXlJiJDCLnTEjv1abXF07O5H+vnMo1f1rOgbJ61n7nY0GMzr/0hCiMaZvPFojzJqYTF+ns/sQBKsrl5M7zxgXlXtlJ0aQnRFFa27YSNyMhijvPG8v180YOqNWqIkOJKv+LDBKHKxp46IQ5RKdj7cFKnllTyNqDlWQn988E+i+cN5ZnPz+fMWlxzMtLIdLV+cfS+Ix47r9hppKKABljmJ7btor1jNEpvPf187llwegOxYRFJHyox0xkEDhQVs+nHlpBUXX3q+yGJ8cwfUQSk7MTGZkax7qDlby14xiHKxopr2/hW89tBiArsf+qw88elcKbXzuX/3pmI6sO+J8Af+XM4Xx18QQSB+EwZl9we7y8taOEQxVt5UYuLchWQiYyACgxExkECisbSY2P6jYxG5MWxyt3LTzpF/QV03P43uWT2V1Sx1V/+JC65rY5ajn91GN2omvnjCAm0smrW4qpqP9oz82vL5kYtOG9oeBvqw5x/1u7OXrC/w/7SuvZcLiKacOTcHazrVZ/qW928+B7+8hMjGbx5AwyEgZmmRORYNKWTCKDhLWWZbvL+NuqQ7yx7Rhu78k/28bA05+b32Wh2Bk/fL29AGkokyG3x8vyfeV88a/rqGly8+bXzj3tjc6HGmst8/73rfZ5ZaeaNTKZf965oJ+jgn9vOUphZSNea/F4weP1crS6ib+ubCt6bAzMGjmMs8elERvpxOkwOIxhxshkZo0c1u/xivQlbckkMgQYYzhnQjrnTEjnpU1FfOnJ9Scdv3l+XpdJWX2zm/ysBBKiI4h0OrhoSlZfh9wpl9PB+IyE9hWmCdEuvvjkOoYnx5CXGkdeaix5aXFkJ0Vrvtkp9pbWd5qUAewvq+/HaD6SlxbHPc9uoraTVcPWts1xXHvw5L09U+Ii+fdXFqo3TYYMJWYig9BlBTlsK6rhD+/uBSB3WAz3XDSxy2violw8dcf8/givS+sPVfKL13extai6ve1TD61gb2nHhGLOqGFMzklkUnbb18TMBGIG8YrNU1lreeT9/UzKTmT2qGFERzhZvq+8y2sqG1ppaHETG9m/H//5WYk8/Jk53PToKlrc3oCvq6hv4Zev7+K+qwv6MDqR8KHETGSQuvvCieRnJ5KZEMX4zATiogbGj/tjHxzg/T1lJ7X5S8oA1hysZM0JPSwOA6PT4rhh3khuXzimT+MMB+/uLOW/X94OQJTLwdy8FL+b1B+XHBvBx6fnBFyaJNjOGJPKb2+YydOrDxMd4WRLUXX7XqhdmTo8qR+iEwkPA+OTWkR6zOkwXDF9YG25U9vUymtbe7511HFe25bE3ffqDhaOT2dilv+NzgcDr9fys9d2tr9udns7JLTQ9v/B+RPTuWZ2LufnZxDlCm2P4kVTstqHyasbWrn2gQ/Zdayu0/MjnQ4uK8jur/BEQk6JmYiEjVc3F9Pcg2Guzri9lh+9tI0nbj8jCFGFp6LqRrYdrWl/bQzcc9FE5o9JZcPhKjYVVjMlJ5GPzxhOekLvdhDoK0mxEXzz4knc+n+rOz3ngkkZJMdG9mNUIqGlxExEwsY/1hUG7V71Lb3bmirc5Q6L5d4l+by4sYjCygZ+c/1Mzs/PAGDmAFrFGB/d9a+hq2bl9lMkIuFBiZmIhIXDFQ2s3O+/uOzpGJ0aF7R79VaL20uE0wR9BekXzhvLF84bS6vHS4RzYG7kEtXFTg8pcZGcOyG9H6MRCb2B+ZMsIoNOXJSLl/7zbC4P0ry4oupGqhpauj+xj7W4vXzzn5v7tKzHQE3KABxd/Hu5YnpOl1t0iQxG+j9eRMJCSlwkU4cncduCvKDcb8W+Chb/cikvbCzC7en9vLXT0eL2cudf17GpsCokzx8ISrtYRXrVrOH9GIlIeFBiJiJhJSUueBO9y+pa+PLf1nPRr9+j2e0J2n0D9du3d/Pm9mNkJak4ameO+dlG7PhChmkqkyFDkOaYiUhYSYsP/grC2iY3kf083Get5fkNRQBkJg7cxGxTYRXHapo5PuCYmxJDflZi0O5/6sbqCdEu7j9hIYPIUKPETETCSlyUi7hIJ/UtwevhOnNMar9v3bS1qIZLpmVz0ZRMcofF9uuzg2FPSR0vbSrit2/vwXPKvqu3nJXHNy7O75BUnY5LpmUTG+kkMzGarKRo0uKjwmaTdZFQUGImImGlxe2lJchzws4ckxrU+wVif1k9/3XhBFwDdGL+vf/Y1GHfyuP+78MDvLXjGE/dMZ/hyTG9ek6ky8GFIdyXVSTcKDETkbBQ1+xm6c5Slu0uJcLpoNUTnB6zCKfh6tn9P4n8nPHpAzYpa2hxs+Fw1wsWDlc0sr+0vteJmYicbGB+aojIoPPA0r188cl1rDpQwZy8lKDd9/aFY0KyDVFSbES/PzNYNhyq6jB86c/R6ka8AZzXmXd2lHDlHz7gxodXcLDc/36oIkONesxEJCysPVjJ/1w5levmjMBjLZ95ZFWvC85++9JJQ2Iz82BbEeC/93ue3cS3nttMVlI0OUkxDE+OIcf3lZ0c3f46Psr/r5oRKTHsKamjtsnNtX9azvv3LlLdMhnylJiJSMh5vJaHPjOHON8vcBfwvcuncMn9y3p8r+HJMSzKz+DSguyQzC0bDFbtL+/2HIeBCZkJZCdFs6mwussk+tfXzeATMzsOJ4/LSGDNtxfzwZ4yVu2vVFImghIzEQkDTodpT8qOC3Roy+kwzB41jEX5GSzKz2B8Rny/r8AcTJrdHtYf6np+mcPARVOyOD8/g/ljUskdFsOB8gbWHqxk7cEKKupbSI2PIjUuktS4SMakd749VpTLyaL8TBblZwb7rYgMSErMRCQsbS+u7fac+2+Yybnj0wf0fK5w09ji4csXjKe4uoljNU0cq22mpKaJktrm9nlnXpg/ub0AACAASURBVAuvbinm1S3FRDodbP7BhYxOi2N0WhzXzPa/6fhLm4q479UdNLV6aWr1MCY9jiduP4PEaP23EzmREjMR6RMHy+sZ1YuNxN0eL1EuB81u/6UznA7DFUHaV1M+khwbyRfPH9eh3eO1lNc3c7SqiaPVjZTUNlNS00yUyxHQ4opLp2Wzp6SOd3aWtrf94Z293Ltkono4RU5grD39FTXhYs6cOXbNmjWhDkNEfHYW1/LCxiPcc1F+r+7T4vaypaiadQcr+evKQ+wv+2h4MzrCwY4fXdzbUEVE+p0xZq21do6/Y+oxE5GgG58R3+lqyMMVDaTGRxIb2f3HT6TLwayRw5g1chi3LRjNkapGdpfUsutYHYWVDcEOW0Qk5JSYiUifuOQ3y4iOcJKflUBCtIvr5o5keHIM331+C3eeP465PaxV5nAYRqTEMiIlVhPFRWTQUmImIkG39lAlR6ubANqHH7MSo/nahRN57NZ5oQxNRCSsqWiMiATdy5uOdmgrrmkKQSQiIgOLEjMRCbrcYR33TzzegyYiIp1TYiYiQffZs0dz05mjTmo7XNEQ0P6LIiJDmRIzEQk6YwzzRp88uf9AeQOff2ItjS2ek9o3Hq6isLKBwVC6R0Skt7pNzIwx0caYVcaYjcaYrcaYH/jaH/G1bTLGPGuMie/k+gJjzHLftZuNMdGnHH/BGLPlhNfTfedvNsa8aIxJ7O2bFJH+56937I1tx/jUwyuoqG9pb4uJdHL7n9ew4XDX2wCJiAwFgfSYNQOLrLXTgRnAEmPMmcBXrbXTrbUFwCHgS6deaIxxAU8An7fWTgHOA1pPOH4VUHfKZQ8D37DWTgOeA+7p8bsSkZBzdzJsuf5QFVf/8UNW7itny5FqWtxefn5tARFOB0VVjf0cpYhIeOm2XIZtG184njxF+L6stbYGwLTtpRED+PsUvhDYZK3d6LtX+fEDvh62rwF3AE+fcM1E4D3f928ArwHfCfwtiUgoNbV6iI5w4vH630oJ2kpoXPfgig7tucNiePvu84h0aZaFiAxNAX36GWOcxpgNQAnwhrV2pa/9MaAYyAd+6+fSCYA1xrxmjFlnjPn6Ccd+BPwCOLV89xbgCt/31wIjAn0zIhJa9c1uzvzxW1x6/zL+se5Ij68vrGzkqdWH+iAyEZGBIaDEzFrrsdbOAHKBecaYqb72W4EcYDtwnZ9LXcDZwI2+f15pjLnAGDMDGGetfc7PNbcBXzTGrAUSgBY/52CMucMYs8YYs6a0tNTfKSLSj7xey/1v7aaqoZWtRTWs2l9xWve5/609tHo6720TERnMejReYK2tAt4FlpzQ5gH+Dlzt55JCYKm1tsxa2wC8AswC5gOzjTEHgPeBCcaYd33322GtvdBaOxv4G7C3k1getNbOsdbOSU9P78nbEJE+YIF/by3u9X3K6pp5c9ux3gckIjIABbIqM90Yk+z7PgZYDOw0xozztRngcmCHn8tfAwqMMbG+hQDnAtustX+01uZYa/No60nbZa09z3e/DN8/HcC3gT/17i2KSF+qrG/hpU1FfPXvGzhYHpyNxQ9WaINyERmaAtkrMxv4szHGSVsi9zTwMrDMV8rCABuBLwAYY64A5lhrv2utrTTG/BJYTdsf1K9Ya1/u5nk3GGO+6Pv+n8BjPX1TItI/yuqamfs/bxLsEmSltc3BvaGIyAARyKrMTcBMP4cWdHL+C8ALJ7x+graSGZ3d/wAw9YTXvwF+011cIhJ6lfUtQU/KQImZiAxdWpMuIqetsqG1+5NOgxIzERmqlJiJyGmravC7aLrXnA7TJ/cVEQl3SsxE5LRV9VGP2c1n5fXJfUVEwp0SMxE5bZNzErm0IJu4SGev7hMT8dH1Y9PjuCA/o7ehiYgMSErMROS0TR2exO8/NYt13/0Yf/r07NMagvzyBeNZ/e3FfHnROKJcDj53zlgcGsoUkSHK2L5YUtXP5syZY9esWRPqMESGvP99ZTsPvrcvoHONgV9fN4OPzxje3nakqpG0+EiiXL3rgRuMPF6Lw0Bb6UgRGciMMWuttXP8HVOPmYgEzV0XjMcVYG+XtbBkatZJbcOTY5SUdWLprhI+9dBK9pXWhToUEelDSsxEJGgaWz24vYH3wq/eX8lg6LXvD8+tL2L5vnKW/HoZz64tDHU4MsBpP9rwpcRMRIJmX2l9j87/9CMrmf/jt/nGPzax61htH0U18JXUNPG6bx/SFo+XLUeqQxyRDGTWWm77v9Xc/ufVrDtUGepw5BSBbMkkIhKQQIbZhsVG8PlzxzJ/bCprD1ay+kAFb24/RnpCFHdfOLEfohx4fvraTprdH/VwREdouFdO39JdpSzbXQbAm9tLuGRaFn+4cXaIo5LjlJiJSND8c/2Rbs8ZmRJLRX0LWYnR3LpgNLcuGI21lpomdz9EOPB4vZa3d5Sc1BYdocEOOX2TcxIpyE1iU2Fbz+vbO0qobmwlKSYixJEJaChTRIIokCG2jYXVPPDePvaXfTTsaYzRL4VObDtaQ0X9yTssDIuNDFE0MhhkJETz+GfPYEpOIgBxka6AF+1I31OPmYgETWykk4YWT6fHF45P445zxjApO5G0+Kh+jGzg2nxKsjs8OYZr5+SGKBoZLJJiInj8s2fws9d28rlzxhAXpXQgXOi/hIgEzdWzc8lNjmFsejxVja3c+dd1Jx2vaXKzcHx6iKIbmD4xYzipcZEs3VXK0l2lfPeyycRG6qNbei8lLpIfXzUt1GHIKfTTLSJB882LJ7V/b61l8aRMqhtbyEyM9n1F0erxEuHULIpAxUQ6uXBKFhdOycJaqwKzIoOcEjMR6RPGGB6+2W9hazlNSspEBj/92SoiIiISJpSYiYiIiIQJJWYiIiIiYUKJmYiIiEiYUGImIiIiEiaUmImIiIiECSVmIqfh9a3F7CyuDXUYIiIyyKiOmUgPVdS3cP/bu9lypIbZo4Zxw7yRXFaQTXSEM9ShiYjIAKceM5EuNLV6OFBWj7W2vS0lLpLn7lzAiJQYRqfF8e7OEp5ceSiEUYqIyGChHjORTlz3wHLWHKzE47WMz4jnkmnZXDwti4mZCRypbORwRSM/uXo4kU4HK/aV8/t39vDF88eFOmwRkUGjxe2ltqmV1PioUIfSb5SYiXSisLIRj7etp2x3SR2/eWs3v3lrN6PT4kj3fUjc9Miq9nPio1zccc4Y7QMpIhIEX3t6Ay9tPMq3LsnnlgWjQx1Ov1FiJtIJl9P/voT7y+rZX1YP0J6UnTkmhW9dMklJmYhID7W4vby65SjrD1Wxu6SWn10znazEaFo9loumZjFvdGqoQ+xXSsxEOhFokhXhNDz0mTkkREf0cUQiIoOP02G49x+baGr1AvCbN3fzk2sK+O0NMzu9ptntIco1OBdcKTET6USgiVmrx/LOzlKumJ7TxxGJiAw+TodhYlYiGw9Xtb32M1pR1+zmRy9uaxuxKK+ntLaZKTmJXDw1iyVTsxmXEd/fYfcZjbuI+LG3tI6dxTUBn//06sMnrdyUoanF7eX93WWU1jaHOhSRgG05Us0j7+8PaQzjT0isov30hL2zo4S/rznMqgMV7T9fW4tq+Pnru7j0/mX9Fmd/UGIm4seYtLhuVwFFOh38x8LRvPils3ns1rkY439Omgwdm49U85lHVxIdoY9WGThqm9y8veNYSGO4fu4IUuIiiY100uT2dDj+763FIYgqNMxg+Ct/zpw5ds2aNaEOQwaZzYXV/Pz1naw/VEl8lIvU+ChS4yNJiYskIyGaG88YyYiU2FCHKWGkscVDWV2z/r8Q6QWP1+J0fPSHbk1TK/P+5832OWinmpeXwtOfn99f4QWFMWattXaOv2OaYybSiWm5Sfz5tnmhDkMGiBc3FlHT1ApAZkI0iydnhjgike55vZZXthylurGV+CgXlxfk4HCEtvffecrzX9hQ1GlSFhPh5GfXFvRHWP1GiZmISC+5PV7ufnojLZ62Xx6zRw1TYiYDgsNheHt7Cf9cf4QbzxjJ5QXht4jp76sPd3rsW5dOYlRqXD9G0/eUmImI9NK+svr2pAwgNS4yhNGI9MyNZ45kbl4K188bEXZzZasbW9l8pNrvsYXj0/j0GSP7OaK+p8RMRKSXNhee/ItDG9pLuHt8xUGW7y1j/aEqjlY3UZCbRFSEg0sLssOrPlgn0+ATol389JqCsEskg0FLh0REemlL0UeJmcPAHeeMCWE0IoF5ZXMxR6ubANhUWM3Xnt7Igvve5hev72RncW2Io2vT0Or22/6DK6aQnRTTz9H0DyVmIiK9tOWEoZZ7Lspn6vCkEEYj0jVrLe/vLvV7rKyuhd++vYdL71/Gf/5tPWsPVoS0RmN5XQuj0+JIiPpogO/CyZlcOXN4yGLqayqXISLSS1uOVLNiXzm1TW7uumB8yFe1iXTlkff386OXtgV8/rThSdx8Vh6XTw/tMGdTa1s5muTYSOKjBvZMLJXLEJE+d7C8nmv/tJyPTc7kf66cFupw+tXU4UnqJZMBYf2hSn78yvYeXbP5SDUr95XziRmhXbEZHeEkd9jgrxGoxExEguLlzUcpqW1mbPrg2bNOZLD5++rDeHo4UnbneWO556KJg3KifTjSHDMRCYopOUn8+KppXDwtK9ShiEgn7ru6gOXfuIBvXZJPflZCt+cvys/g60vy+ywp83otb24L7XZQ4UY9ZiISFOdOSA91CCISgKykaO44Zyx3nDOW7UdreOT9/Ty7ttDvuUum9O0fWp96eAXrDlWx678v7tPnDCRKzERERIaoSdmJ/PTqAv5j4RjWHqxk5f5y3th2jIYWDw4DF0zK6NPnN7Z6GafpDyfRqkwRERFp19DiZvvRWprdHs4am9anz9p9rJZxGfFDbv6aVmWKiIhIQGIjXcweNaxfnjU+s/t5bkONJv+LiIiIhAklZiIiIiJhQomZiIiISJhQYiYiIt1qaHHz4Z4yvN6Bv2BMJJxp8r+IiHSppqmV2x5bTU1TK81uL9+7fDKL8jNDHZbIoKQeMxER6VKr28tlBdm88KWzOXN0KoahVdpApD8pMRMRkS6lxkdxy4LRuByGmEgnI1MH/0bSIqHSbWJmjIk2xqwyxmw0xmw1xvzA1/6Ir22TMeZZY4zf0r3GmAJjzHLftZuNMdGnHH/BGLPlhNczjDErjDEbjDFrjDHzevsmRUSk91xOB19dPIHGFk+oQxEZtAKZY9YMLLLW1hljIoD3jTGvAl+11tYAGGN+CXwJuO/EC40xLuAJ4CZr7UZjTCrQesLxq4C6U573U+AH1tpXjTGX+F6fd1rvTkSGvA/3lPHP9UeYnpvE9BHJ5GclEunSYMHpSoqNICk2KdRhiAxa3SZmtm3PpuPJU4Tvy56QlBkgBvC3VOdCYJO1dqPvXuXHD/h62L4G3AE8feIjgUTf90lAUQ/ej4hIuxa3l28/v4V9pfXtmzRHOh1Myk6gIDeZAl+yNjY9HqdD86ZEJPQCWpVpjHECa4FxwO+ttSt97Y8BlwDbgLv9XDoBsMaY14B04Clr7U99x34E/AJoOOWarwCvGWN+TttQ61mdxHQHbUkdI0eODORtiMgQ8+gH+9lXWn9SW4vHy8bCajYWVre3xUU6GZkax9njUrlh3kjGaFNlEQmRgBIza60HmGGMSQaeM8ZMtdZusdbe6kvafgtcBzzm5/5nA3NpS8DeMsasBcqBcdbarxpj8k655gu0DZP+wxjzSeARYLGfmB4EHoS2TcwDercig8Cv3tiF11pcDgcupyHCaXA6HEQ4DS6Hg4Xj0xiRosnZ1loe+2B/QOfWt3jYfrSG7UdreGjZfhaMS+XGM0bxscmZRDg17Cki/adHdcystVXGmHeBJcAWX5vHGPN34B46JmaFwFJrbRmAMeYVYBZtQ6OzjTEHfDFkGGPetdaeB9wM3OW7/hng4Z6/LZHB69H391Pb7O70uDHw1H+cyRljUvsxqvBjjOFT80bxqzd39fjaD/aU88GectITorhuzgiunzeC3GFKdkWk7wWyKjPd11OGMSaGtt6rncaYcb42A1wO7PBz+WtAgTEm1rcQ4Fxgm7X2j9baHGttHm09art8SRm0zSk71/f9ImD36b45kcHI0c1cqNgIJzNGJvdTNOHtswtHkxIXedrXl9Y287t39nDOT9/h8RUHgxjZ6Wt2e9hXWofb4w11KCLSBwLpMcsG/uwbsnTQNlH/ZWCZMSYRMMBG2oYgMcZcAcyx1n7XWlvpW7G5mrZJ/a9Ya1/u5nn/AfzGl8g14ZtHJiJtupujvmBcGlEuZ/8EE+bio1xckJ/BM76J/6fLa+E7/9pCbVMrd543LkjRnZ79ZfUs+fUyIp0O8tJi+f7lUzhrXFpIYxKR4AlkVeYmYKafQws6Of8F4IUTXj9BW8mMzu5/AJh6wuv3gdndxSUyVHW3evC8iRn9FEn4a/V4eWdnSdDu99N/7yQnKYZPzBwetHv21JHKRqBtEcNZY9OUlIkMMprVKjLAtM0e6FxclHrLjlu6s5Syupag3rOr+X19qdntobK+hYPlHy1kf3P7MVrcGtIUGUyUmIkMMM5uErP1h6r6KZLw9491vRvC9CcxukdrpoLmLx8eZOaP3uCHL21rbyusbGTzkeourhKRgUaJmcgA091Q5pqDFf0USXirrG/hze3Hgn7f+KjQJGadJZmjtG+lyKASmk8YkSFoT0kdb24/xmUF2X5LL3i8NqDq8w/fPIdn1xbyr/VHKK/vOEw3PiMhKPEOdC9uKqLVc/olDkemxPKVxeNpdntpOf7l8TIhs//+/dY3u3l7RwmvbD7KjuLaDseTYiJI7cWqU5GB7qlVhxiZGstZYwfPXEslZiJ9bE9JHb97ezcvbCzCa+GlTUW89J8LO5z3+3f24PFavnj+uC73cpyUnch3LpvMNy7O599binlxYxFv7yjB7bUkRLv44vmhXTUYLp7tZiWmy2FweztP3C6emsVVs3KDHVZAnt9whJc3HWXprlKau5hDNi4jvts5hyKD2TWzcwfddmpKzET60Naiam5+dNVJE9C3FtVQ3dhKUkzESecermjgmbWFvLa1mKfuOJPk2K57QiKcDi6fnsPl03Mor2vmYEUD+VkJxEbqx3rXsVo2FXY+9+qiKZn86ONTSY6NpKK+hU2FVTy/sYh9pfVsP1rDxMwELp+e048Rt6luaOWZtYf531e200XO2G6cto6SIc41CHfm0Ce4SB+akpPELz85g888uopIp4P7b5jBzJHDiI3suHJyb2kdADuP1ZIQHdHheFdS46NIjY8KSsyDwWMfHOjy+Gtbj/Ha1mNMyUnk5S8v5E9Ly3l509H24zuP1XLLY6tY/f8W96pHqrCygQt/9R5uj+Vz547h7gsndjhnZ3EtNz2ykqykaHYdq6WpNbBVlvFRLiZmJdDi9nbZwyoiA4t+mkX62J8/PAC01Z26++mNrD1Y6Xf/xXMmpPOJGTn88pPTuy0iK50rrm7i2bWHAzr3eK/lzWflcWr+VVbXwtHqpl7FkjsslhvmjcRiGdnJ/qUWS0ltM5sKqwNOygDqmt388KVtnPPTd3hg6V5qmlp7FauIhAf1mIn0oWa3hxEpscweNYz4KBe3LMhj1shhfs/9yuIJ/Rzd4PTwsn0BT/of5ps4PzotjounZvHK5uKTjm8rqiEnOaZX8dy1eDzXzM5lUnZir+7TmeKaJn786g5+9/YefnXdDBZPzuyT54hI/1CPmQwJ1lre3VnCS5uKaGr19Ntzo1xOvn/FFL50/jiaWj2cPzGjw9wyCZ7K+haeXHUo4PNTTpjH97WPTezQU/mvDUc4XNFAbyRGR/RZUnai2mY33/7XFu2hKTLAqcdMhoRHPzjAj3yFOReOT+OBm2Zz11MbcBj4zfUziY7o22r5O4prqW8JTcX4oeTPyw/Q0BJ44j3shFIT4zLiuawghxc2FrW3vbTpKC9tOsr03CSunzeSG+aNDGa4QRXlcvDTawo6nQxdUd/C5iPVtLi9nDNB+6mKhCslZtJvvvb3DWwtqiEtIZL0+CjS4qNIS4giPT6KlPhIEqJcjM9MCHqP0qr9FfzvK9vbXy/bXcb8H79NdWMrY9Li/M73Craaplai9YuwT9U3u/k/33y+QEWdMGm+sr6F8vpmv+dtLKxm85HNnDkmldFpcfzgxa08vfowG753Yb/8/xOIe5fkc86EdL/HqhtbueeZjby1o23f0MRoF5cWZHPF9OGcMToFhyY1ioQNJWbSL6y1LN1VSnl9Czu7KMb+vcsnc+uC0UF7bklNE198ch2eU2oPVDe2TZS+YFJGv9TAuXdJfp8/Y6j726pDVDX0bAL8r9/cRX5WAo99cIAdxTVd7qvptfDA0r3ctXg8f11xiNT4yKAkZZFBSuw2H6mm1ePtENO6Q5Xc/uc1VJxQjLimyc3fVh3mb6sOk5UYzeXTs/n4jOFMyUlUXTSREFNiJv2isLLRb5X6U3W2cu10tHq83PnXdZTW+u8FATg/PyNoz5PQaXZ7eHjZ/h5f1+qxfP6JtQEvFvjHukIq6lto8XjJHda7RQHHjUmP50+fnsXXnt7Yo2HYUz23/gjF1U38/sZZpMRF0tDi5q8rDvHAe3tPSspOVVzTxEPL9vPQsv2MSY/jlS8v7POhfRHpnBIz6RfrD3e/sfZNZ47igknBW1H2l+UHWXOwstPjCVEu5ualBO15Ejqvbi6muOb0Slv0ZNumVo/l9W1tXb4j/GyrdbqWTM3G44UvPrmuV/dZvq+cK373PtfMzuWJFQe77AH0p7HFo6RMhpTDFQ1EuRxkJEaHOpR2SsykX6w/1HmCdFxTq4fGFg8xfoqvno7abuo6nTMxPWzmB0nvhKKG17pDlfz41e1Mzk5kcnYiY9LjT3tY3O3x8uB7e4MSV2FlI79+c/dpXTslp+9Xj4qEisdr2VZUw/biGjYVVrGpsJqtRTV4vJa0+Cim5CRyaUE2n5wzIqRxKjGTfrEhgB6zZ9YWsv5wFb+9YWZQygt0l3QtmqhhzMEiIbr/P8oOlDfwwNJ97a/vu2oa15/mqs0nVhxkYxdbSPWXyTlJoQ5BpE+U1TVzzR8/5EC5//I3ZXXNLN1Vyui0uH6OrCN1F0ifa3Z72FpUE9C5e0rqeH1rF6sDeiDC2XnvhTFw3kT/K9hk4Ens4RZW/kzJSeSWs/JIiOqY5EVHOJiYmdDl9X94dy8r9pVT19zzsihrD3X/h0t/UI+ZDEaNLR4+++c1nSZlJxqVGrwpCqdLPWbS52qb3Jw/MZ2V+ysCWjV34ZTgzDNblJ/JH9/dS6WfZ16Qn9mne0v+a/0RVu6vYHRaLHmpcYxOi2NkaqxqR/WRnu4teqoIp+HHV02jIDeZey6ayHPrj/D48oMcq23iM/PzuHn+KIprmrj0/vc7vcehigauf3AF18zO5efXTu/R84/1cuunYFFiJoNRWV0zGwMYtQHISw19j5kSM+lzafFRPHDTHLxey/biGpbvLWfFvnJW7q+gtunk3oXcYTHkZ3XdMxGocRnx/OW2M7jhoRUn9WLcuiCPb10yKSjP6MzSXaU8t/7ISW3GwMLx6fzltnl9+uyhKDEmsI8yY2BeXgrzRqfwyuaj7C2tJ8rl4E83zaYgNxmAuCgXnz5zFDeeMZJmt7d9MnxKXCQTMxPYeay2y2c8u7aQa2fncsaY1IDjP92FC8GUFBPB8F5uPyUSjnKHxXDbgtFsLapmW1ENtV30aqvHTIYUh8MwJSeJKTlJ3L5wTPtEzOX7ynjwvf2U1TXzscmZQa2jNC03iUdvmcsvXt/Z/gP51KrDrDlQyfiMeC6Zlt0newvuLa3r0GYtDE8On5U/g0l3PWazRiZz+fQcLpmWTaZv9dVXFk9g9YEKIpyG2aM6rs41xpy0QtEYw8M3z+H5DUd4adNRdhR3nqB9+19beO0r5wRUuNVay7EwSMyumjVcNcxkUDLG8N3LJwOwubCaDYfbFqPVNLnZdayWrKRocn1/lOQGcbX16VJiJiHjdBim5SYxLTeJ288eQ0VDS5+skpw3OoW/f24+Xq/lcGUDW4tq2HKkbTVOX2yTZK1lb0nHxAzodANz6Z3ETib/Lxyfxj0XTWzvDTuR02E4swe9WgAjUmL50qLxfGnRePaU1PLMmkIeeG9fh/NKapsJNMdpaPEQH+Wi2d2z0hbBNDk7UUWQO9Hq8fKv9UcYmRLL6PQ4MhL0x9VAdvx3TjhTYiZhweEwpPXhnK/jzxiVGseo1DgumZbdZ88pqW2mvpNCobNHKTHrC3GRLoxp65U80bLdZSzbXcY9F03ki+ePC+ozx2Uk8M1LJjEsLpL7Xt1x0rEx6XEB9z7FRblY8+3FFFU3sbmwio2F1by2pZh9ZfVBjbczsZFOfvupvt8vdqB6f08Z9zy7CYBIl4Ozx6UxY0QyM0YkM31EctC3kJO+5fFavNaGdakkJWYiQdZZb9mw2IiwWIo9GDkchoQoFzVN/ntAT2elZKA+f+5YKupbeNDXczYyJZbPnTOmR/cwxjA8OYbhyTEsmZrNvUvyOVzRwId7y/hwbzkf7i3vcgeL3vjhx6cyNj2+T+49GLy44aNN7VvcXt7eUcLbvj1HAX581bSw3txeTtbq8TL/x2/R4vaSFBNBYkwE80an8IMrpoTNUL4SM5Ega/VaFoxLpanVS2OLh6bWtq85eSlh84M/GM3NS2F/eT31zW7qmz3Ut7jbe9Ba3d4+ffY3L87nsoK2XtjJ2Ym4gvDX+IiUWK5LGcl1c0dirWVPSR0f7GlL1FbsK+80Ce2JxZMyuHrW8F7fZzC7Z8lEzhybyvK95Xy4t4xjNScnyD97bSeXTMtWz9kAER3h5PaFY9hRXMvMEcnMGJnMo+/vD6vPZmNP7fsfgObMmWPXrFkT6jDkNNQ1u4mJcPbLRuIytHi9loZWD/XNblwO06flUfqbx2vZWlTNtqdSAgAAIABJREFUh3vL+WBPGasPVNDU2vPk8/efmsWlBX03rD/YWGvZV1bPh3vLWb63jFX7K7msIJu7L5zQ65It0r8aWty4HA4iXQ4OlNWT18+jGcaYtdbaOX6PKTGTUDpc0cDVf/yQzMRoJmYlMDdvGNfN1bCASE80uz1sOFTF+3vKeGZNYUDlN4yBD+5dRI5KZMgQVN/s5vNPrOUri8f7XZXd17pKzMJ39psMCSNSYvnhx6eyp6SO5XvLeXHj0ZOOt3q8eL0D/48Hkb4U5XJyxphU7r5wIu99/Xw+e/bobq85d0K6kjIZsuKiXDz+2TNCkpR1R3PMJOSWTM1iydQlJ7U9vvwAb+0oYdX+CpJjIrh8Rg6fmDE8KHtoigxmkS4H9y7J5+nVh7sspHn93NBu1Cwi/qnHTMLSjuJa3t1ZSkOLh6LqJh5Yuo+Lf7OMi371Hn94d09YFOQUCVeRLgd3LR5PbKT/EhipcZEsyg9+YWUR6T31mElYun3hGJ5cdYj4SBeLJmWQlxrHw8v2sfNYLX9ffZirZ+WGOkSRsHb7wjFcNSuXx5cfpLCygWa3l2a3h2a3l/ljUol06e9ykXCkyf8StjYermJSdmL7L5CiqkZ+9cYu7lo8Piy2zRDpS40tHjYcrmJUaqzmgokMMl1N/lePmYSt6SNO3kYnJzmGn107PUTRiPSPoqpGvvy39WwsrKLVYzEGpucmc+uCPC6akqUK/SKDnBIzEZEw8syaQtYcrGx/bS1sOFzFXU9tICHaxeXTc/j2pZOIjdTHt8hgpEkGIvL/27vP8Liqc//73zWj3rvV3OSCu2UsF4zpYAihBEggnZBCICHkQEh/DuknhIRAThohIQkHEkr4QyAxnWCwKbbl3nuTbRVLtqrVZtbzYsZCskZWsTR7RvP7XNdcnll7r617SdbonrVXkRDy9s4jPR6rb27n6VVlIb3Pn4icHv12i4iEkF9/fBa/+mgxiT3MqDxjRLISM5FhTH3hIiIOa27z8PKmcp5YcYAZhal8+/LJTCtI5cmVB4h2G1bsqWHlXt/tzSlay09kWFNiJmFtd1UDr22pYOmOI5w9PotbzhvndEgi/WKtZVdVA6nx0Xx24VhmjfJNehmXncR3Lp8M+PbG/NTDy3lnVzVT8pWYiQxnSswkrH37mQ0s31MDQHby8NmkWiKHMYap+amnPMftMvzqo7P44P8uVWImMsxpoIKEtTaPt+P5gZomByMRGVrZybH89hNnMik32elQRGQIKTGTsPbgJ2eTHOfr+N1XrcRMhrc5YzJIjot2OgwRGUK6lSlhy1rL4g2HuWZWAW6XIcbtwuu1uFzG6dBERGQAqupb2Hiwlgsm5TgdimOUmEnYMsZw09ljnQ5DREQGwZGGFs659z80t3n59Fmj+e8rpkTk0jCR12IREREJOVlJsVw+LQ+A/3t3H596eDn1zW0ORxV86jETERGRkHDHJRMxxnDWuEzOGpcZkWMqlZiJiIhISBiZkcB91890OgxH6VamiIiISIhQYiYS5g4dO87W8jqstU6HIiIip0mJmUgYs9bypb+t5rIHlvLwsj1OhyMiIqdJiZlIGHttSyVrDxwD4H9e2MKza8ocjkhERE6HEjORMNXm8fLTF7d0vPZauOsf61m+u5rWdu8paoqISKjSrEyRMPXEiv3srmrsUubxWm546D2iXIaxWYlMzE3mKxeMZ1KeNr4WEQkH6jETCUN1zW3c/9qObuXjshMBaPdadlQ2sHj9YV7YcJhP/Ok93eYUEQkDvfaYGWPigLeAWP/5T1trv2eMeRgoAQywHfiMtbYhQP0ZwB+AFMALzLHWNnc6/jxQZK2d5n/9JHCG/3AacMxaWzzwJooMPw8u2UVNY2u38vtvKGbl3qOs3FPDpLxkJuWmMGdMOi9sLOeeF7dyzaxCB6IVEZG+6sutzBbgQmttgzEmGlhmjHkRuMNaWwdgjPklcBtwT+eKxpgo4DHgU9badcaYTKCt0/FrgS7JnLX2hk7H7wNqB9QykQiTGh/N1PxUZhSm8bmFXfcQfe3O87SchohIGOj1Vqb1OZE8RfsftlNSZoB4INC7/iJgvbV2nf9a1dZaj79eEnAn8ONAX9d/3euBx/vVIpEI8KULxpOdHNul7KyiTNwu02Md36+UiIiEsj6NMTPGuI0xa4FK4FVr7XJ/+V+AcmAS8OsAVScC1hjzsjFmtTHmG52O/Qi4D2jq4cueA1RYa7sPpPF97ZuNMaXGmNKqqqq+NENk2EiKjeKZWxfwP9dM55wJWcwsTGXR1BFOhyUiIqepT7My/b1cxcaYNOBZY8w0a+1Ga+1Nxhg3vqTsBuAvAa6/EJiDLwF73RizCqgGxltr7zDGjOnhy36MU/SWWWsfAh4CKCkp0T0aiTgjMxL4+LxRfHzeKKdDERGRQdKvWZnW2mPAEuCyTmUe4EngugBVyoA3rbVHrLVNwAvAmcBZwGxjzF5gGTDRGLPkRCX/2LRr/dcVERERiQi9JmbGmGx/TxnGmHjgYmCbMWa8v8wAVwJbA1R/GZhhjEnwJ1vnAZuttb+31uZba8fg61Hbbq09v1O9i4Gt1lrN7xcREZGI0ZdbmXnAI/5bli7gKWAxsNQYk4JvuYx1wK0AxpirgBJr7d3W2qP+GZsr8U0OeMFau7gPX/OjaNC/iIiIRBgzHKbQl5SU2NLSUqfDEBEREemVMWaVtbYk0DGt/C8iIiISIpSYiYiIyIA0t3nYUVFPRV0zbR6v0+EMC9rEXERERAbksff28ePFWwCYmp/Cnz8zhxEpcQ5HFd7UYyYiIiL9Vtfcxq9ee38N+E2H6rj2d++wvaLewajCn3rMREREpN/+vnw/9S3tXcoOHjvO5b9aytSCVOaOSWf26HRmjkwjNyVO28L1kRIzERER6RdrLX99e2/AY+1ey7oDx1h34Bh/XLoHgHHZidz74ZnMHp0exCjDk25lioiISL8YY7D0fbmt2uPtTM1PGcKIhg8lZiIiItJvozMT+3zusaZWvvPsBnZVNQxhRMODEjMRERHpF6/X9mt5jHav5ZnVB/nec5uGMKrhQYmZiIiI9Msf3trNmv3H+l1v+Z5q6pvbhiCi4UOJmYiIiPTZyr01/OKVbQOq2+axvLX9yCBHNLwoMRMREZE+qWtu4yt/X4PHO/B9tls9nkGMaPhRYiYiIiJ9svFgLeV1zad1jQk5yYMUzfCkxExERET65EBN02lfoyi777M5I5ESM3GM12s5dOw41g68S1xERILnQM3x06pfkBZPQozWtj8VfXfEMZsO1XHlb5aREOPmp9dO5+riAqdDEhGRU9jfzx6zhBg3507IJj8tnuNt7RSmJwxRZMOHEjNxxDOry3hm9UEAmlo9qNNMRCT09Scxcxn4040lLBiXNYQRDT9KzCTorLX86vUd7Kt+/xc8NT7awYhERKQv+jPG7OZzxykpGwCNMZOg23CwtktSBpASr88IIqHOWsv+6tMf/C3hqbGlnerG1j6dW5gez1cvmjDEEQ1PSswk6P617lC3MvWYiYQ+Ywx/W7GPytNcLkHC04GjvSfluSlxnDMhi198ZCbxMe4gRDX8qJtCgsrrtfx7/eFu5SlxSsxEwsG3PzDZ6RDEIb31lqYnRLPk6+cTF62E7HSox0yCavX+oxyu7fppe0JOEplJsQ5FJCIifXG87dQr9n983iglZYNAiZkETbvHy8PL9nQrv/X8cbhdxoGIRESkr66amd/juDG3y/DJ+aN7rNvm8bJybw0PvLadlXtrhirEYUG3MiVolu08wosby7uUjcpI4IoZ+Q5FJCIifWWM4Y5LJjK9IJVtFfVU1bdwpMH3mJSbQl5qfI91XcaQnRTL1PxUphekBjHq8KPETILmubXdB/3fdekZxESp41ZEJFxcPGUEF08Z0a86bpdhTFYiY7K0HVNvlJhJUDS1tvPypq69ZTMKU7liep5DEYmIOOtATRPff34TUW5DlNtFdlIs379qatDjaPd42VZRT7vH0u71Eh8dRVJsFH9bvo87LpnYbdzYxoO1pMZHk5UUq5mXQ0CJmQyKv769h0O1zUwvSGVGYSqjMhIw5v1xY69urqCptevA0W9dNgmXxpaJSIRKjoti7tgMSvcd5d1d1eQkOzMJygJX/eZtPF7fFiw3n1tEYkwUf3hrN7eeP65LYra/uokrf7OsY7eWxBg3aQkxfOOyM7St3iBRYianraXdw32vbqe+ub2jLCUuioUTsvjpNTNITYjmgkk5/OFTsznS0MKGslqMgQXjtSK0iESutIQYvnjeOL6Ib/HeqoYWR+KIdrvIT4vr2KD8+pKRPFV6gIzEmG5DTVLio4h2uWj1eAFobPXQ2HqcZTuOcPn0PKobfAvQjkiJ7fLhXPpOiZmctn+uOdglKQOoa27nhQ3l7Kxs4K83zSU/LZ5Lp+b6Ds5zIEgRkRBmjCEnOc6xr//Y5+bxyqYK9lQ3Mj4niZmFaTy+fD/xJ93GTEuIweUCTlo5o6nVQ3ltM+fc+wYA//zy2czwD/LXnZH+UWImp+2PS7svgXHC9ooGrv3dOzz5xfmMztSgTxGRUDQ6M5EvnFvU8fqDM/I4a1xmwHNP3PIE3xIa47KTiIlyMSLl/cTyZy9uJTrKxdcumcjMkWkd5VvL65iUmzIELRg+lJjJaXv0c3P509I97KtuYn9NI/trmmhu83Ycn1Zw6mnUIiLSu6bWdr76xFoaW9qJdrv44rlF/OSFLSTEuHngo7MoSBvc99mMxJiA5Z0TsxEpsXz14vfXNstMjKG6sZXle6rJSIxlf01Tl8Ts5y9t42cfnkGWFhXvkRIzOW15qfH89xVTOl57vZbK+hb2VTdyuLaZy6fnaUkMEZHT1NTq4dXNFQD8+mOz+M0bO9l0qA7wbTAeLL/66KyO5/knJYN3XzmFaLeLKXkpjM5M6DbO7OHPzMFai/RMiZkMOpfLkJsaR26qc+MlRESGmzb/gPvikWmMz0nknV3VgG9m5LjspKDFceXMnhcF78vMTE0KODUlZiIiImEgNsrNzz88gytn5vODf23qKJ9WkKpt7YYR3V+KYJV1zepSFhEJExmJMXykZCSxUS6W7TzSUa7xWsOLeswi1LGmVs7/xRKS46I4Z0I2507MZuH4rB4He4qISGh4b3dNx5pjAC9uPMyKPTXMHZvhYFQyWJSYRajH3ttHU6uHplYPT68q4+lVZRgD0/JTuee66UzN1yazIiKh6MmV+7u89lr46hNruO7MQmKiXL6H2/fvgnGZFAVx/JmcPiVmEarxpO2RAKyF2uNtTMhJdiAiERHpTW1TGy9sLO9Wfri2md+8sbNb+W0XjOeuS88IRmgySDTGLEIt7GE7pK8tmqilLUREQtQ/1x6ktd3b+4l+Gw7WDmE0MhTUYxaBKuuaeXHj4W7lk3KTuXJGz9OgRUTEWZdMGcGU/BRe31LJsabWjvJlO49QdvR4t/NX7KnhK4+v4b6PzNSH7jChxCzC1De3cc69b9By0ieu60sK+c7lk7WnmYhICMtPiyc/LZ45Y7oO9L/0/rcCnn+8zcMVM7TIdzjRTyrCuF2GQAtkPFVaxn8/t4mq+pagx9Qbay0t7d3HxImICDS0tLO9sr5Lmcv4ltH41PzRXDo1t6N81b4ajja2nnwJCSFmOKxjVVJSYktLS50OI2y8t7uajz70XsBjqfHRvHHX+Y4tm9Hc5mF7RT1bDtex5XA9mw/VsaW8jp9cM52rTrHatIhIpDrS0MKLG8sZnZFAVlIs2cmxZCTGdFt0dtOhWq76zdt4vJbclDiyk2MxBsZmJfL5hUVMzksmyq3+mmAwxqyy1pYEOqZbmRFoflEmJaPTKd13tNux9IRo0hOiHYgKDtceZ+HP3uiyQS5AtNtw3oRsR2ISEQl1J3rGepOdFNvx/lpe10x5XTMA68tqeW7tIWKjXNwwZyQ/vHrakMYrp6bUOEJ9duHYgOUfmJ7n2D5meanxTM7rvlTHuOwkUh1KFkWcdv+r23nknb1OhyHDQHZyLIkx7h6PT81P4csXjA9iRBKIErMItWjKiIC/oJdPy3Mgmk5ff7rv64/KSOCea6dzzoQsdlQ2sPbAMUfjEnHKkm2V/PDfm/n8I6X8c81BGlranQ5JwpQxhvhTJGZXFxcwIiUuiBFJILqVGaGi3C5mjkzjnV3VHWWF6fFMK0hxMCrfG0NhegKXTB5BfIybj84dRVV9C+W1zY7GJeKUwvQE1pXV8tqWCl7bUkFMlIsLzsjmgzPyuWhSDomxehuX0/ffV0zhxgVjnA5DUGIW0X7z8TOpbmihoaWdhpZ2UuKiHbuNeUJBWjwFafFdyrKTfYNZRSLRx+aO4q0dVdQ3txMT5WJ6QSqF6Qm4DDS1epSYSb98+wOTeXlTObHRbmKjXMRGuZg1Kp0Pzy50OjTx06xMEZEgW7bjCPe9uo3x2UlcM6uABT3sxHFCeW0zlfXNTMpN0XpUIsPAqWZl6jdcRCTIPNayZv8x/rGqrE+LOuemxjGjME1JmUgEUB+4iEiQRbt9yVhBWjxzT1rBXWSoHG/18K1n1pOeEENaQjTpCTGkJ8b4l0nyleWmxGktM4cpMRMRCbLZo9N58ub5eKzVNmgSNMbAc2sPnfKckRnxPH3LAs3OdJDSYhGRIIuNcjOvKJMF4049tkxkMEX3oSfsQM1xfvX6jiBEIz3p9adkjIkzxqwwxqwzxmwyxvzAX/6wv2y9MeZpY0xSD/VnGGPe9dfdYIyJO+n488aYjSeVfcUYs81f597TaaCIiPRffXMb/9lawXCYICY+bpfptk1TINvK63s9R4ZOX25ltgAXWmsbjDHRwDJjzIvAHdbaOgBjzC+B24B7Olc0xkQBjwGfstauM8ZkAm2djl8LNJxU5wLgamCGtbbFGJMz8OaJDK3K+mbuf3UHVfXNFKYnMDIjgZHp8YzMSGBcdpIGa0tY2V3VwH+2VvLGtkpW7Knh7PFZXDhphNNhySCKdptu296dbH3ZMWoaWx3bMznS9ZqYWd/HpRPJU7T/YTslZQaIBwL9pBcB66216/zX6ljN1N/DdidwM/BUpzq3AvdYa1v8dSr72SaRIWWtpanVwyPv7uW3/9lJY6sn4HlnjEjmgY8WMznP2UV7RfpiybZKPvOXlV3KLpqspGy4iXa7aG7znvIca6HNc+pzZOj06eO8McZtjFkLVAKvWmuX+8v/ApQDk4BfB6g6EbDGmJeNMauNMd/odOxHwH1AU4A65xhjlhtj3jTGzOkhppuNMaXGmNKqqqq+NEPktNU1t3HLY6s45943uPelbT0mZQDbKuq5+jdvc/DY8SBGKDIw68tqu5VdNEk3LIab2D704mtrJmf1KTGz1nqstcVAITDXGDPNX34TkA9sAW4IUDUKWAh8wv/vNcaYi4wxxcB4a+2zPdRJB+YDXweeMgGWo7fWPmStLbHWlmRnZ/elGSKnpby2mat+vYyXN1VQ09japzoxUS5y9QYnYeBEYmYMjM9J4jMLxpB/0i4cEv76MgHglvOKghCJ9KRfy2VYa48ZY5YAlwEb/WUeY8yT+JKov5xUpQx401p7BMAY8wJwJr5bo7ONMXv9MeQYY5ZYa8/313nGfwt1hTHGC2QB6hYTR22rqGdv9ckdvKdWMia9T4NtRZx25cw8bjmviMl5KdrmaRjrKTFzGbh48ghGZyYwYURykKOSznr97TPGZANt/qQsHrgYuNcYM95au9Pfm3UlsDVA9ZeBbxhjEoBW4DzgfmvtYuD3/uuPAf7tT8oA/glcCCwxxkwEYoAjA2+iyOCYOCKJe6+bQVNrO01tHh57dx+Hetlcfe5YLR4q4eHq4gKnQwDA47W0ebzERbudDmVYunjyCP6ztaLbh8wPFRfwyxuKHYpKOuvLx6I84BFjjBvfrc+ngMXAUmNMCmCAdfgG7WOMuQoosdbeba096p+xuRLf5IAX/EnZqfwZ+LN/CY1W4Ear+doSAvJS47l+zsiO100tHl7bUkF1YytV9S0B68wbmxms8ETCTnObh98v2cXOqgYOHj1O2dHjVDe28PgX5jO/SL87p9La7uVnL21lan4Ks0alMyYzgQCjfrq5+8op3H3lFI42trL+YC3rDhxj7YFj3Hbh+CBELX2hTcxFBsFb26v4yuNrqD3esRoMcdEu1n/vUi2ZIRJAS7uH7z67kadXlXUpz0iM4d1vX0hslHrMTuWZ1WXc+dS6jtcZiTF878opXDUzv08JmjhLm5iLDLFzJ2bzr9sWEu+//TIiJZY/frpESZlIAPuqG7n2d+90S8oAPjAtlxjt1XhK1lr+uHRPl7Kaxla++sTaPk9MktClEZ4ig2RUZgIvfvUcnlt7iOvnFJKXqhltIid7Y2sltz+xhvrm9oDHP33WGPX49OKdXdVsOVwX8Nje6kYyk2IDHvvzsj00trTT6vHS2u6lpd1Lq8dLm//f1nb/w+PFay0lozO4/aIJmsAUZErMRAbRmKxEvnrxBKfDEAlZu6oaekzKJuUmc0auZgT2Zs3+oz0e23Okidmju086OtrYyg//vblfX+ftndXERLn48gUafxZMSswkpFhreXFjOa9vqSQ22sW8sRnML8rUYociw8RNZ4+l9ngbD765izZP1zHOH5oVGjNDQ92MwrQej+090hiwfOsA9798Z9cRvnT+OPViBpESs2HsWFMrf357L3HRLm4+p4ioEB+30ebx8sk/LWf5npqOsr8v3w9AUVYiC8Zn8qHiAqYXpmpgsEiYcrsMX1t0BndeMpGaxla2VzTQ5vHS5vEye3S60+GFhbljM4iPdnO8rfvOI2kJ0QHr7KpqCFjem8YWD42tHpK0tl3Q6Ds9jK3Zf4xnVpfxt8/PC/mkDHwLH7p6+FS2+0gju4808th7+4mJcjGzMJWvXzpJ64SJhCljDJlJsZzVw3go6VlctJvikWm8u7u627ELethGa3dV4J603iTEuJWUBZm+28NYS7uXf9xyVlgNQr9+TmHAN5sTPrNgDBmJMbR5vGQlxQQxMpHwVl7bjDGQkxyr21LDwMVTRnR7rxyVkUBRVmK3c1/aWM7Gg7VMyk2mvrmd+uY2Glra8fZhtayjTW29nySDSonZMHbZtFynQ+gzr9fyVOkB7n152ynPWzRlBAvGZ/Xpmu0eL01tHpJionBpVpFEsHd2HeGL/7eK+pZ20hOiOSM3mUm5KUzOS+aM3BQmjkgiIUZ/DsLJoikj+JF/MH9SbBTpidF8ZkHgGa2XTcvt9vfAWktTq4f65nZqj7fR0NJGY4vvdUNLG/ExURSmx2vpEgfoN1FCwi9e2cbvluzq9byqhsAr7AfyVGkZ33l2A2MyE/jk/NHcdPZYTfuWiPTLV7ZT3+KbCXm0qY33dtfw3u73x3I+efN85mml/bAyMiOBFd+5iNSE6AGNuTXGkBgbRWJsFLmpmlwVSpQKiyOW767mwTd3UevvJs/r4xtDT1sfBTJhRBLf+sAkzhqXSXZybI/r/ogMZ29ur6J0X8/LK3xmwRglZWEqJyVOE6GGIfWYiSPy0+LJSo4l1T+D6OIpI6hrbmdHRT07KhvYc6SRESlxlNc2d5l5lBIfeMZRIHPGZDBnjCYHSGRbu/9Yl9dul+H7V02lKCuRwvR4CtLCZwxqpGv3eFm28wjnnxF4gL8MD0rMxBEjMxK6vM5Lje+yiKHXa3G5DNZa6lvaqahtps1jmZKfEuxQRcLa8+sOdnn9oeICPjV/tEPRyEBZa/ne85v4+4r93PeRmVx7ZqHTIckQUWImIenEYH1jDClx0aTE9b2nTETe97VFZ/Cjf29mZHoC43ISueW8cU6HJP3w9KoyNpQdo9VjeXyFb13Hu/6xjiXbqshIjCElPpqUuCiyk2O1gfkwocRMRGQYu3x6HpdPz3M6DBmgmsYWHnl3X5cyr4Xn1x3qdu6Fk3JI1ofYsKfB/xISvv/8Jr78t9VsLdcA/VB2tLHV6RBEhq2axtaOXrETcpL7PmOyoq7vk6MkdKnHTIacx2sp3VtDU6unY1Vqay07Kxt4ZXMFL28qZ31ZLQCfOXuMg5HKqeyrbuSCXyxh7tgMPjg9j+vnjNSMMJFBUFHXzP2vbufZNQdpafcyKTeZWaN821PlpPR9Z4TKumbG5yQNVZgSJErMZMi5XYZxOUlsPOhLvv60dDd/X76f3QE2242LoD/0G8pqmTAiibjo8Gjz82sP4bV0rIE1vyiTCSOSnQ5LQlRLuwdrIcbt0gLPvUiMjeLN7VW0tHsB+Os7e99PzPrTY1bfPCTxSXApMZOgyEqK7Zji/ci7ezlQczzgefExkXF3/UhDC598eDnzxmbw0KdLnA6nV9Za/rm26+y+qvoWJWbSo5c3VXD742sAeO/bF2kR01NIio3i3g/P4PUtlUwrSKV4ZGrHsRGdesxchlNuo1SpW5nDghIzGTLHWz3Ex3TvDRqTmYi1UHa0e3IWKbfGfrJ4C7XH22hq9fR+cgjYfLiOXSdtgtyfXRgk8swsTKUgLZ6vXjxB+9r2wTkTsjlnQna38qTYKB761GzG5yRRe7yNa373To/X0Biz4SEyuifEEceOt/Lzl7eyZFslFXXvd7E/+rl5vPX1C4gKcHsjutO+bDWNrf1a6T9cvLSxnGfXvN/7dKShhf99fYeDEZ3a3iON/PSFrd3K91U3ORCNhItRGQn8567zuL5kJFHab3HAjDEsmppLUXYSMVGn/j6+vrWCGk3QCXv6bZEhk5cazyVTcrn/tR38/KTNyY3xjavobEZhKinx75c98Np2rD1Fv30Yam33cs+LW7qUeb2WX766PWST0OfXHWLZziPdyv+0dLcD0Ui4MMZETA94sPS21+++6ibavd4gRSNDRbcyZUgVj0zjuS+f3a3cGMPFk0ew6VAtV87M54oZeYzOTOw43u7xMjozkfTE4XULpKXdw6HargN0T/SsRnJ/AAAgAElEQVQSrtp3lMum5ToR1im9vrUyYLl6QUSCK9Bdhs6ykmJJitWf9XCnn6A45r7rZ/Z4LMrt4nMLxwYxmuBIio2iuDCNjMQYLpkyggsn5RAf4+bq4nwOHgs8IcJJ1lpqmwLfGuntj4SIDC6Xf1X/GLcLi2VMZiKzR6dz5uh0Zo9OpygrUSv/DwNmONwqKikpsaWlpU6HITIsbTxYy21/X82h2mZa29+/TZKfGsc7377IwcjEKdbaiEoAPF6Lxz8d0u0yvd5SHCrHWz3sr2miKDuRNo+XhBj1rYQrY8wqa23AKfn6qYrIKU0rSGXJ1y/AWsvRpjbKa5upqGumoaXd6dDEAburGvjIg++SFBfFmMxExmYlcu7ELC6cNMLp0IbMo+/u5fv/2gzAL693bgPx+Bg3Z+T6lqiJ1lCCYUuJmYj0iTGGjMQYMhJjmJKf4nQ44pD/fm4j1Y2tVDe2sq+6iTe3V/Hoe/t49LNzWTA+y+nwhkTnpWLCZYkbCV9KzAZofdkxbn98DWkJMfz4Q9MYmZHAX9/eS2Ksm5goF9bCzJFpTM5L1sykIGpsaWfpjioq6lqYkp/CpNxkbeorMogamrv3lHq8ltX7jw67xGz57moOHjvO5sPv7+Hb1KqeYhlaSswG4LXNFXzl8TUcb/NAdRP3vLiV9WXHqAvwhpWeEM3/u3UBRdnav2yoVNY189Kmcl7bUsl7u6pp9XSdLj4mM4G5YzOYX5TJwvFZ5KRoBXKRgahrbmNbRX3AY6M6zaoeLhZvOMz/vbuvS9krmyqYNzaTmSPTHIpKhjvdpO4nr9dy+xP+pMxv2c4jAZMygKNNbXz2ryu16N8Q+sG/NnP3c5t4a3tVt6QMYG91E0+VlnHnU+t48E2tvSUyUM+vPURzW+B1skZnJAQ5mt7trKxn6Y4qdlTUU9fc1u91EWuPt3UrK913lM/8ZcWwW2NRQod6zPrJ5TKMy05ig39D7r7YW93Ezf9XymOfnxc2G1aHi3aPl6U7qno9Ly7axdO3LGBKnsZGiQzUkysP9HhsTAj2mP3gX5tZuuP9xZETYtzkpsSRkxJLbkocI1LjGJEcx3VnFpKa0H3Iw48+NI3tFQ1s6XQrE3zDVCJpVqoElxKzfqo93kZ9c/dPUb0p3XeUrz+9nl/dUIxL6z8NmnVltT32VnZ25yUTmVaQ2ut5IhLYxoO1PX4gTY2PDpjYOOmx9/Z1ScrAN3B/95FGdh/puu9r8ag0zhyV3u0aKXHRPPLZOXz49++yv8a3BZkxMDvAuSKDRYlZPzS3efjCI6XsHeAegf9ad4hRGfF8/dJJgxxZ5BqdmUCM2xXwFuYJ84sy+NzCoiBGJTL8PFXac2/Z6MzQu4353u5qAGaPTufq4nyeXHmAheOz+PSCMXzz6fUd24xdX1LIjFN8aMtJjuO5L59NZX0LCTFuCtLi9eFahpQSsz7yeC23P76GFXtrTus6v31jF6MzErl+zshBiiyyZSXFcsXMPJ5ZfTDg8aKsRB64YZZjC0KKDAfNbR7+uSbw7xjArBAbCP/ypnLW7D8GwHc/OJkzR6Xz6bPGdBz/040lfOlvq5man8Kdl0zs9bZkemLMsNseTkKXBv/30aPv7uWVzRWDcq3fLdk5KNcRn5sWdN+6KSMxhh9ePZWX7ziX3FTNwhQ5Hcea2rhyZj7jsruPI1s4PovvfHBy0GK586m1fPB/l57ynN+9sbNji7POKVdlnW+fWrfL8IdPzuZri87QWDEJOeox66PBXO5izpiMQbuWwPTCVP7ztfPYUdnAzsoGvF7LjWePIUXrl0kY21nZwKceXs63PjCJq4sLHI0lNzWOn1wzHfCNNSvdW8P6slrqW9p54IbioK3VaK3lQE1Tl310/7Z8H1cXF3Rs3u3x2i5LerS0ezlQ08RPFm9he0U9r3/tPG77+2oe/OTsoMQs0l9KzPpo7tgM4qJdPU4V76uirERuXDBmcIKSDkXZSRRlJ3HpVKcjERkcGw/WMjkvhStn5DsdShfTClIdnUjz1BfP6tLL9dzaQ/z85W187uyxfOHcIg4dO97xPp0Y46a+uZ2fvbSVlzaVA7D7SCM3n1uknjIJWUrM+igu2s38okyWbOt9aYZAYtwuvnTBOG49f5x2AhCRXm04WEtlfbMGmncSKJn67uWTMQbGZScRF+0mKTaKv39hHuOyk2hu9fDgW7s7Vu6/ZlYBRVmJGKMFvyV0KTHrh4XjswacmN1yXhH/dfHEQY5IRIarL55bxIzCVI63eoiP0Ye5npy8An9OSlzH7h61TW08uXJ/x3n3XDddPWUS8jT4vx92VDQMuO6LG8u1UrSI9FlOShxXFxcoKTsNqQnRnD0+izNyUyguTCXGrT95EvrUY9ZHOyrqefIU6/j0Wr+ygaNNbWRoyrWISNA8+rl5Tocg0i/6+NBHsVG+hQUHKtptSA+xlbFFREQktCgx66NRmQm8+F/ncN2ZhV3KE2LcTCvoff/FrKRYjW0QERGRU9KtzH5IiYvmvutncvHkHJ5be4gPzsjjosk5xEe7+e0bO/ndkl0d57qNIcptiHK7iHIZfvyhaQ5GHnk8Xstf3t7D58+J3K2YdlTU43IZxmYmamafiEiYMMNhQHpJSYktLS11OgwJIesOHOP6P7zLK3ecy+jM7quVR4Lr//AuK/bUkBwXxYzCVGYWpjGjMI3ikWnaDUFExEHGmFXW2pJAx9RjJsPSsp1HMIaITcoq6ppZ6d/Xtb65nbd3VvP2zuqO4znJsZw5Kp3ff/JM3WIXEQkhGmMWAprbPFpKYxDVNLby13f2ctPZYyP2+/rChsOcqumV9S28tKmc+pb24AUlIiK9UmLmsEPHjvPDf29Wr8Ugqjvexr+/spBvXjYpYr+v/15/uE/nndjUWUREQoMSM4dlJsVw9xVTnA5jWBmTlciIlMgdQ3Xo2HFW7Tvap3Mr61qGOBoREekPjTFzmPbNDJ7Wdi8xUcP/s8irmyv6fG5lvRIzEZFQMvz/Son4bTxUy+I+3uILZyfvHXgqlfW6lSkiEkqUmEnEmJafSrvXS31zm9OhDKnikWnMHp3e63nRbkNKnHajEBEJJbqVKREjJsrF1cUFTocRFD+4air/KD3Aqv1H2XyoDm+nGZpZSTF8Yt5oPjF/FDnJkTsWT0QkFCkxExmGphWkMq0gFYDGlnbWHjhG6d6jFKbHc8XMPI1tFBEJUb0mZsaYOOAtINZ//tPW2u8ZYx4GSgADbAc+Y61tCFB/BvAHIAXwAnOstc2djj8PFFlrp/lffx/4AlDlP+U71toXBtxCkQiXGBvF2eOzOHt8ltOhiIhIL/rSY9YCXGitbTDGRAPLjDEvAndYa+sAjDG/BG4D7ulc0RgTBTwGfMpau84Ykwm0dTp+LdAtmQPut9b+YkAtEhEREQlTvQ7+tz4nkqdo/8N2SsoMEA8EWmd8EbDeWrvOf61qa63HXy8JuBP48Wm3QkRERGQY6NOsTGOM2xizFqgEXrXWLveX/wUoByYBvw5QdSJgjTEvG2NWG2O+0enYj4D7gKYA9W4zxqw3xvzZGBNwepkx5mZjTKkxprSqqirQKSIiIiJhpU+JmbXWY60tBgqBucaYaf7ym4B8YAtwQ4CqUcBC4BP+f68xxlxkjCkGxltrnw1Q5/fAOKAYOIwveQsU00PW2hJrbUl2dnZfmiEiIiIS0vq1jpm19hiwBLisU5kHeBK4LkCVMuBNa+0Ra20T8AJwJnAWMNsYsxdYBkw0xizxX6/Cnwh6gT8Cc/vZJhEREZGw1GtiZozJNsak+Z/HAxcD24wx4/1lBrgS2Bqg+svADGNMgn8iwHnAZmvt7621+dbaMfh60rZba8/3Xy+vU/1rgI0DbZyIiIhIOOnLrMw84BFjjBtfIvcUsBhYaoxJwbdcxjrgVgBjzFVAibX2bmvtUf+MzZX4Jge8YK1d3MvXu9d/q9MCe4Ev9r9ZIuFp1b4aMhJjGZuV6HQoIiLiAGNtoMmU4aWkpMSWlpY6HYbIabv8V0v5SEkhN5091ulQRERkiBhjVllrSwId016ZIiFkWkEKBWnxTochIiIO0ZZMIiHk3g/PdDoEERFxkHrMREREREKEEjMRERGREKHETERERCREKDETERERCRFKzERERERChBIzERERkRChxExEJIi8Xkt9c9ugX3c4LBYuIkrMQkabx0vZ0aYuj9Z2r9Nhicgg21nVQMmPX+PWx1bx0sbDNLd5Tvuadc1tPLnyAABPrNjPCxsOU93QctrXFZHg0wKzIcLjtVx035u0+JOxURkJvHbneQ5HJSKDbfmeGlravby4sZwXN5aTHBvFpdNyufe6Gbhchkfe2UtCjJsPTM8jKbZvb9Hbyuv5yQtbuGBSDkeb2vjWMxsAmDgiice/MJ/MpNihbJKIDCL1mIWIuGg3c8dmdLy+85KJxETpxyMy3KQnRHd5Xd/Szo7KBlwuw87Ker73/Ca+/vR65v7kNZ5Ysb9P14x2u/jw7ELWHTjGJVNyAIhyGeqOt5MaH91LbREJJeoxCyFnj89i6Y4jTMpN5qqZ+U6HIyJD4NKpueSmxFFe18zZ4zO5/4ZicpLjADpuRwI0tXrYWdnQp2sWj0yjeGQa4Btrtu57i0iJi8IYM/gNEJEhpS6ZEHL2uCwA7lp0Bi6X3lBFhqNot4svXTAOt8tQXtvckZS1tnv5f6sPdjn3quL+f0AzxpAaH62kTCRMqccshEzJT+GiSTlcNDnH6VBEZAh9+qwxXDR5BLs69Yi9urmCmsbWjtdnj89kRmGaE+GJiIOUmIUQt8vwyxuK9UlXJAIUpMVTkBbf8fqJlV3Hk9163vhghyQiIUCJWYjRQF2RyLNqXw3Ldh7peJ2WEM38oowu59z3yjb+s7USa2HBuEwWjM9k7tjMPs/cFJHwoN9oERGHPbHiAJ3Xh71hzkii3O8PAV6z/ygPL9tDU6tvzbPNh+v407I9uF2G33xsFh+YnhfskEVkiGjw/yDyeG3A58PdsabW3k8SkYAaW9pZvOFwx2tj4JPzRnc5Z3dVIwkxbq6dVcDnF45lSl4K4HufKd13NKjxisjQUo/ZIGpobucL/1fKHZdM5K0dVXit5ZuXThrWMyxrm9oo/uGrfPfyyXzh3CKnw5Egs9ZqTORpWrzhcEdPGMBFk0YwMiOhyzkfnJHHNbMKuryXVDe08O7uaiblpgQtVhEZeuoxG0SpCdHExbh5ceNhzirK5KG3dtPqGd7bKrW0e7hyZj45KVpZPBK9urmCxesP936i9Ojp0rIur29cMLrbOXHR7m4f8DKTYrliRj7jc5KGND4RCS71mA2yb102ieZ2D2eOSucnH5pOXLS7y/Gt5XWs2neUHRUNXDEjj5IxGT1cKTzkpMTx64/NcjoMcciiqblOhxDWqhtaWLG3puN1UXYiC8dnORiRiDhNidlp8notXmtp91riot2MyUogLsqXjH10zshu5y/bcYQfL94C+KbH//hD0/nw7MKgxiwioSEjMaZjFwCAG88ao1vDIhFOidlp2FZez/+8sIXy2mY+Pm8UHykp5Ka/rOTJL54FEHBs2fVzRpKXGs8ZucmMyUzoMvNKRCKLMYavLZrIhoO1AFx7ZoHDEYmI04y14T97sKSkxJaWljrytRtb2nl2zUEyEmO4fHoeG8pqmV6Y6kgsIiIiEvqMMaustSWBjqnH7DQlxkbxyfnvD9ZVUiYiIiIDpftoIiIiIiFCiZmIiEgA9c1t3PCHd2lu8/R+ssggUWImIiISwMFjx6lvbqehpd3pUCSCaIyZQ+qa23h+7SH2HGkkPy2e/NQ48tPimVGYqunyIiIhYGxWIotvX6j3ZAkqJWYO+OmLW3jknb00t3XfFeCl/zpHW6yIiISA2Ch37yeJDDLdygwyr9fyxIoDAZOyb142iTNGJDsQlYhIZNtZ2cALGw7TNsy30ZPQpx6zIHO5DKX/38UcqGliz5FG9hxpZPeRRs4el8UHZ+Q5HZ6ISER6fUsFZ4/PIlqLfovDlJg5INrtoig7iaJsbT4sIhIKbj63SGPJJCToo4GIiEQ8JWUSKpSYiQyQx2tp13iUiPDW9iqtZSUiQaHETGSAFm84TPEPX+XGP6/gt2/sZOXeGlra9cd7uKlvbuNzj6zkxY2HnQ5FRCKAxpiJDNATK/bT0NLOm9ureHN7FQAxUS6KR6Yxb2wGc8dmcOaodBJj9WsWzt7eWU2bx7JgXJbToYhIBNBfDJEB2FfdyDu7qruVt7Z7WbGnhhV7agBwuwzT8lOYOzaDuWMzmTsmg9SE6GCHK6fhja2VZCfHkpkY43Qog2JHRT1/fWcvizf4egBj3C5io13ERrk7nse4XfziIzMZmZHgcLQikUeJmcgAPLnyQJ/O83gt68pqWVdWyx+X7uF/rpnOx+eNGuLoZLAcPHacf649SEu7l9ufWMMDN8wiJip8R4Dc+dRanll9sE/n/n3Ffm6/cAI7KxuYmp+Cy6XB8SLBoMRMZAC2ltf3u05MlEtr1YWZn7+0lZZ23wSPFzaU09hSyoOfnE18THiuCF/b1Nbnc//41m4efHMX1sKHZxdy73UzlJyJBIESM5EBON7a/0H+i6aMIDVetzHDxdoDx/jn2kNdyt7cXsVrWyq4cma+Q1GdntGZiX0+t91rO54/vaqM5jYPV83MZ39NEzcuGKOFWEWGiBIzkX6qb25j5d6afte7bnbhEEQjQ8Fay08Wb+5WHhPl4vwzsh2IaHDccckEPjF/FPXN7dQ3t/Gr13ZQuu9on+r+e/1h/r3eNy5tZ2UD91w3Y8jiPHTsOHmpcVpbTCKSEjORftpe0cBZ4zJJiY8mNT6alLhoDhxtYvH695dTcBmYOCKZ6QWpzBiZxoyCVKYXpDoYtZzKo+/u5d3d1fz242dijOHlTeWs3Ns9YTl/YjbJceHb65kcF90Rf9nRJtaVHevx3OzkWGaNTOOVzRXdjg32mm4er2V3VQPrympZsaeap0rLyE2JY+k3LyDa7cJay2PL95OeEM3EEcmMzUpUj50MW0rMQkBlfTNZibEavxEmZo9O59HPzetS5vVasL61zf72+XnMGpVGQszp/Xq1tnvZfaSBljYvE0ckh+24pnDwxrYq/rO1kk2H6piSl8JPX9wa8LzhNEbwN//ZSZvH9nj8KxeOx+O1AROzwZwA8evXd/Dgm7toPGl4QEyUqyP5evS9fdz93KaOY9Fuw9isRBZNyeWuS88YtFhEQoESM4cdrj3OWT/9D0/cPJ/5RZlOhyMD5HIZ7rt+JsbA2eMHZ72r460eLntgKQB/+/y8QbuudHfLeeP46kUTGJedhMtlqD3efZB8bJSLiyePcCC6wePxWn76whY2HqrtWNIlkML0eD46ZxSbD9d1OxYb5aJkTMagxWSBnJQ49hxpZOH4LNISonlvdw0Tcnx7Ca8vO8aP/t31tnKbx7K9ooEJIxoHLQ6RUKHEzGEpcdE8efN8Zo9OdzoUOU1x0W4euKF40K6XmuC7VVp7vI3tFfVKzIbQ3LFdE42irERW73//Nl9RdiJfOKco7BcLdhmYMzaDqQUpfPos3wB+X8+UIdbfQxUT5SIjMYaYKBeT85KJcbto9XjJT43j8+cUcf2ckSQN4vfh9osm8JULx7OjsoHkuChykuNobfeyq6oBgMq6FmaPTuf2Cyew+XAdb26vYvmeGlrbvRSmxQ9aHCKhwljbc1d2uCgpKbGlpaVOhyEy6GqPt2GtJTbKrVuZQfS1p9bx0sbDXDkzn4+UFHLmqPSIHYj+zafXM68ogytn5jsyruvE36jO3//mNg/L99SQEhfFrFH6UCvhxxizylpbEvCYEjMRka52VzWQmxp32uMERUQCOVVipncdEZGTFGUnOR2CiEQozTcWERERCRFKzERERERChBIzERERkRDRa2JmjIkzxqwwxqwzxmwyxvzAX/6wv2y9MeZpY0zAQRnGmBnGmHf9dTcYY+JOOv68MWZjgHp3GWOsMUZrBIiIiEhE6EuPWQtwobV2JlAMXGaMmQ/cYa2daa2dAewHbju5ojEmCngMuMVaOxU4H2jrdPxaoCFAvZHAJf7rioiIiESEXhMz63MieYr2P6y1tg7A+BaXice3gPPJFgHrrbXr/NeqttZ6/PWSgDuBHweodz/wjR6uKSIiIjIs9WmMmTHGbYxZC1QCr1prl/vL/wKUA5OAXweoOhGwxpiXjTGrjTHf6HTsR8B9QNNJX+sq4OCJZO4UMd1sjCk1xpRWVVX1pRkiImFtW3k9G8pqnQ5DRIZQnxIza63HWlsMFAJzjTHT/OU3AfnAFuCGAFWjgIXAJ/z/XmOMucgYUwyMt9Y+2/lkY0wC8F3g7j7E9JC1tsRaW5Kdnd2XZoiIhLXV+4+SFKflJ0WGs37NyrTWHgOWAJd1KvMATwLXBahSBrxprT1irW0CXgDOBM4CZhtj9gLLgInGmCXAOGAssM5/rBBYbYzJ7VerRESGoY/NHcXYrESnwxCRIdSXWZnZxpg0//N44GJgmzFmvL/MAFcCWwNUfxmYYYxJ8E8EOA/YbK39vbU231o7Bl9P2nZr7fnW2g3W2hxr7Rj/sTLgTGtt+ek3VURERCS09aVPPA94xBjjxpfIPQUsBpYaY1IAA6wDboWOMWIl1tq7rbVHjTG/BFbiG8j/grV28RC0Q0RERCTsaRNzERERkSA61SbmWvlfREREJEQoMRMREREJEUrMREREREKEEjMRERGREKHETERERCREKDETERERCRFKzERERERChBIzERERkRChxExEREQkRCgxExEREQkRSsxEREREQoQSMxEREZEQocRMREREJEQoMRMREREJEUrMREREREKEEjMRERGREKHETERERCREKDETERERCRHGWut0DKfNGFMF7OvhcBZwJIjhhIpIbTdEbtvV7siidkcWtXt4GW2tzQ50YFgkZqdijCm11pY4HUewRWq7IXLbrnZHFrU7sqjdkUO3MkVERERChBIzERERkRARCYnZQ04H4JBIbTdEbtvV7siidkcWtTtCDPsxZiIiIiLhIhJ6zERERETCghIzERERkRAR8omZMeYMY8zaTo86Y8x/GWN+ZIxZ7y97xRiT30P9G40xO/yPG/1lySdd84gx5gH/sVHGmDeMMWv81788mO3tFHew2z3aGPO6/9pLjDGFwWxvp7gHvd3+8o8ZYzb4r/GSMSbLX55hjHnVf/6rxpj0YLX1pLiD3e6PGGM2GWO8xhjHpqI70O6fG2O2+sufNcakBautJ8Ud7Hb36bpDLdjt7nT8LmOMPbk8WBz4eX/fGHOw09cbNn/H/OU9/ryNMV8xxmzzv7/dG4x2Djprbdg8ADdQDowGUjqV3w48GOD8DGC3/990//P0AOetAs71P38IuNX/fAqwN0La/Q/gRv/zC4FHh0u7gSigEsjyn3cv8P1Oz7/lf/4t4GcR0u7JwBnAEqDE6TYHsd2LgCj/859F0M+71+sOx3b7X48EXsa3CHlWJLQb+D5wl9NtdaDdFwCvAbH+1zlOt3sgj5DvMTvJRcAua+0+a21dp/JEINAshkuBV621Ndbao8CrwGWdTzDGTABygKX+Iguk+J+nAocGMf6BCka7pwCv+5+/AVw9iPEP1GC12/gficYYg+/ne+LnejXwiP/5I8CHBr8Z/Tbk7bbWbrHWbhvKRgxAMNr9irW23V//PcCRnuGTBKPdfblusAXj9xvgfuAbPVzTCcFqd6gJRrtvBe6x1rYAWGsrh6YpQyvK6QD66aPA4ydeGGN+AnwaqMWXKZ+sADjQ6XWZv6yzjwFPWn96je+TxivGmK/g+w9z8aBEfnqC0e51wHXAr4BrgGRjTKa1tnpQWjAwg9Jua22bMeZWYAPQCOwAvuw/Z4S19jCAtfawMSZn0FvRf8FodygKdrs/Czw5OKGflqC0uw/XDbYhb7cx5irgoLV2ne9veEgI1v/z24wxnwZKga/5kxsnBaPdE4Fz/NduxtdruHKwGzLUwqbHzBgTA1yF75YbANba71prRwJ/A24LVC1A2cmZeZf/LPgSlr9aawuBy4FHjTGOfZ+C2O67gPOMMWuA84CDQDsOGcx2G2Oi8X2SmgXkA+uBbw960INA7Q5Ou40x38X3//tvg9KAAQpmu/tw3aAJRruNMQnAd4G7Bzn8AQviz/v3wDigGDgM3DdYbRiIILY7Ct/tzvnA14GnTAhl5H0VNokZ8AFgtbW2IsCxv+Pr7TlZGb7xBScU0qmr1xgzE994k1Wdzvkc8BSAtfZdIA7fJqpOCUq7rbWHrLXXWmtn4Xszw1pbOwjxD9RgtrsYwFq7y99D+BSwwH9OhTEmD8D/r9Nd38Fqd6gJWrv9g4ivAD7RqcfYKU78vHu6bjAFo93jgLHAOmPMXv/5q40xuYPViAEIys/bWlthrfVYa73AH4G5g9eEAQnW//My4BnrswLw4uzf7wEJp8TsY3TtBp3Q6dhVwNYAdV4GFhlj0o1vtt0if1nAa/rtx3cvHGPMZHyJWdVpRz9wQWm3MSarU8/gt4E/D0Lsp2Mw230QmGKMyfafdwmwxf/8eeDEbJ8bgecGrQUDE6x2h5qgtNsYcxnwTeAqa23ToLei/4LV7r5cN5iGvN3W2g3W2hxr7Rhr7Rh8f7TPtNaWD35z+ixYP++8TvWvATYOWgsGJljva//EN3kNY8xEIAY4MojtCA4bAjMQensACUA1kNqp7P/h+8+2HvgXvnvPACXAnzqd91lgp/9x00nX3Q1MOqlsCvA2vjFXa4FFEdLuD+O7V78d+BP+WS3Dpd3ALfh+eU/Uz/SXZ+Kb9LDD/29GhLT7Gnx/qFqACuDlCGn3TnzjVtb6H47NTgxyuwNed7i3+6SvuxcHZ2UG+ef9KL4xWOvxffjMi5B2xwCP+a+9GrjQqXafzkNbMomIiIiEiHC6lSkiIvC5eS4AAAA4SURBVCIyrCkxExEREQkRSsxEREREQoQSMxEREZEQocRMREREJEQoMRMREREJEUrMRERERELE/w9K84/AZBi4NQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "merged_poly['geometry'].plot(figsize=(10,10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://gis.stackexchange.com/questions/271733/geopandas-dissolve-overlapping-polygons\n",
    "# https://nbviewer.jupyter.org/gist/rutgerhofste/6e7c6569616c2550568b9ce9cb4716a3\n",
    "\n",
    "def explode(gdf):\n",
    "    \"\"\"    \n",
    "    Will explode the geodataframe's muti-part geometries into single \n",
    "    geometries. Each row containing a multi-part geometry will be split into\n",
    "    multiple rows with single geometries, thereby increasing the vertical size\n",
    "    of the geodataframe. The index of the input geodataframe is no longer\n",
    "    unique and is replaced with a multi-index. \n",
    "\n",
    "    The output geodataframe has an index based on two columns (multi-index) \n",
    "    i.e. 'level_0' (index of input geodataframe) and 'level_1' which is a new\n",
    "    zero-based index for each single part geometry per multi-part geometry\n",
    "    \n",
    "    Args:\n",
    "        gdf (gpd.GeoDataFrame) : input geodataframe with multi-geometries\n",
    "        \n",
    "    Returns:\n",
    "        gdf (gpd.GeoDataFrame) : exploded geodataframe with each single \n",
    "                                 geometry as a separate entry in the \n",
    "                                 geodataframe. The GeoDataFrame has a multi-\n",
    "                                 index set to columns level_0 and level_1\n",
    "        \n",
    "    \"\"\"\n",
    "    gs = gdf.explode()\n",
    "    gdf2 = gs.reset_index().rename(columns={0: 'geometry'})\n",
    "    gdf_out = gdf2.merge(gdf.drop('geometry', axis=1), left_on='level_0', right_index=True)\n",
    "    gdf_out = gdf_out.set_index(['level_0', 'level_1']).set_geometry('geometry')\n",
    "    gdf_out.crs = gdf.crs\n",
    "    return gdf_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_out = explode(merged_poly)\n",
    "gdf_out = gdf_out.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>level_0</th>\n",
       "      <th>level_1</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>137</th>\n",
       "      <td>0</td>\n",
       "      <td>137</td>\n",
       "      <td>POLYGON ((73.08167443854619 33.64971041125181,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>138</th>\n",
       "      <td>0</td>\n",
       "      <td>138</td>\n",
       "      <td>POLYGON ((73.08048286196323 33.64842963545275,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>139</th>\n",
       "      <td>0</td>\n",
       "      <td>139</td>\n",
       "      <td>POLYGON ((73.08590979531148 33.64842425733259,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>140</th>\n",
       "      <td>0</td>\n",
       "      <td>140</td>\n",
       "      <td>POLYGON ((73.08602513530349 33.64849867396733,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>141</th>\n",
       "      <td>0</td>\n",
       "      <td>141</td>\n",
       "      <td>POLYGON ((73.08637402636823 33.64746556190988,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>142</th>\n",
       "      <td>0</td>\n",
       "      <td>142</td>\n",
       "      <td>POLYGON ((73.07970317989704 33.64744830228123,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>143</th>\n",
       "      <td>0</td>\n",
       "      <td>143</td>\n",
       "      <td>POLYGON ((73.08325249924476 33.64740392037901,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>144</th>\n",
       "      <td>0</td>\n",
       "      <td>144</td>\n",
       "      <td>POLYGON ((73.07990289845706 33.6474174815158, ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>145</th>\n",
       "      <td>0</td>\n",
       "      <td>145</td>\n",
       "      <td>POLYGON ((73.08633580861908 33.64752966910199,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>146</th>\n",
       "      <td>0</td>\n",
       "      <td>146</td>\n",
       "      <td>POLYGON ((73.08451491779707 33.647542315558, 7...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     level_0  level_1                                           geometry\n",
       "137        0      137  POLYGON ((73.08167443854619 33.64971041125181,...\n",
       "138        0      138  POLYGON ((73.08048286196323 33.64842963545275,...\n",
       "139        0      139  POLYGON ((73.08590979531148 33.64842425733259,...\n",
       "140        0      140  POLYGON ((73.08602513530349 33.64849867396733,...\n",
       "141        0      141  POLYGON ((73.08637402636823 33.64746556190988,...\n",
       "142        0      142  POLYGON ((73.07970317989704 33.64744830228123,...\n",
       "143        0      143  POLYGON ((73.08325249924476 33.64740392037901,...\n",
       "144        0      144  POLYGON ((73.07990289845706 33.6474174815158, ...\n",
       "145        0      145  POLYGON ((73.08633580861908 33.64752966910199,...\n",
       "146        0      146  POLYGON ((73.08451491779707 33.647542315558, 7..."
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gdf_out.tail(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fill holes\n",
    "gdf_out.geometry = gdf_out.geometry.apply(lambda x: Polygon([coords for coords in x.exterior.coords]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_out.drop(columns=['level_0','level_1'], inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>POLYGON ((73.08295538706595 33.6439754184318, ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>POLYGON ((73.07865527387217 33.6441800683143, ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>POLYGON ((73.08394535005179 33.64416157585504,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>POLYGON ((73.08381836849819 33.64418993095924,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>POLYGON ((73.08082628858956 33.6441800683143, ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>POLYGON ((73.08228349437942 33.64417143849997,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>POLYGON ((73.07858870101883 33.64438841668865,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>POLYGON ((73.08445327626619 33.64446361935632,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>POLYGON ((73.08436204680049 33.64436992422939,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>POLYGON ((73.0832330775162 33.64415745609443, ...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                            geometry\n",
       "0  POLYGON ((73.08295538706595 33.6439754184318, ...\n",
       "1  POLYGON ((73.07865527387217 33.6441800683143, ...\n",
       "2  POLYGON ((73.08394535005179 33.64416157585504,...\n",
       "3  POLYGON ((73.08381836849819 33.64418993095924,...\n",
       "4  POLYGON ((73.08082628858956 33.6441800683143, ...\n",
       "5  POLYGON ((73.08228349437942 33.64417143849997,...\n",
       "6  POLYGON ((73.07858870101883 33.64438841668865,...\n",
       "7  POLYGON ((73.08445327626619 33.64446361935632,...\n",
       "8  POLYGON ((73.08436204680049 33.64436992422939,...\n",
       "9  POLYGON ((73.0832330775162 33.64415745609443, ..."
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gdf_out.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"73.08228211360914 33.64417005772968 3.6047967242325285e-05 3.728079787634897e-05\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,67.28837739625723)\"><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"7.456159575269794e-07\" opacity=\"0.6\" d=\"M 73.08228349437942,33.644171438499974 L 73.08228349437942,33.64420225926541 L 73.08230568533054,33.644205957757265 L 73.08230568533054,33.64420558790808 L 73.08230815099178,33.644205957757265 L 73.0823167808061,33.644176369822446 L 73.08228349437942,33.644171438499974 z\" /></g></svg>"
      ],
      "text/plain": [
       "<shapely.geometry.polygon.Polygon at 0x7f40d29eeba8>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gdf_out.iloc[5].geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"73.08227995372462 33.64416789996328 4.036742893731571e-05 4.1597679654614694e-05\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,67.28837739760621)\"><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"8.319535930922939e-07\" opacity=\"0.6\" d=\"M 73.08228378747674,33.644169460093046 L 73.08228358640072,33.64416944061808 L 73.08228338438585,33.64416944152691 L 73.0822831834932,33.64416946281028 L 73.08228298577244,33.644169504251046 L 73.08228279324084,33.64416956542638 L 73.08228260786275,33.64416964571215 L 73.08228243152953,33.64416974428921 L 73.08228226604028,33.644169860151806 L 73.08228211308341,33.64416999211783 L 73.08228197421948,33.64417013884088 L 73.08228185086533,33.644170298823965 L 73.08228174427946,33.64417047043484 L 73.08228165554935,33.64417065192259 L 73.0822815855803,33.64417084143557 L 73.08228153508617,33.64417103704022 L 73.08228150458213,33.64417123674085 L 73.08228149437943,33.644171438499974 L 73.08228149437943,33.64420225926541 L 73.0822815044518,33.644202459735006 L 73.08228153456749,33.64420265818539 L 73.08228158442313,33.64420285261771 L 73.08228165351657,33.64420304107355 L 73.08228174115189,33.644203221654735 L 73.08228184644638,33.644203392542366 L 73.08228196833946,33.644203552015206 L 73.0822821056034,33.644203698466974 L 73.0822822568556,33.64420383042257 L 73.08228242057262,33.64420394655287 L 73.08228259510544,33.64420404568818 L 73.08228277869608,33.64420412682996 L 73.08228296949534,33.64420418916093 L 73.08228316558144,33.64420423205326 L 73.08230535653256,33.644207930545114 L 73.08230554802405,33.64420795303843 L 73.08230574079163,33.64420795698813 L 73.08230593304376,33.64420794235752 L 73.08230612299371,33.64420790928255 L 73.08230630887616,33.64420785807063 L 73.08230648896351,33.644207789197694 L 73.08230657909698,33.64420774434789 L 73.08230785431087,33.64420793562997 L 73.08230804471425,33.64420795493154 L 73.08230823609075,33.64420795594599 L 73.08230842668806,33.64420793866402 L 73.08230861476095,33.644207903243874 L 73.08230879858733,33.64420785000989 L 73.08230897648401,33.644207779449495 L 73.08230914682207,33.644207692208774 L 73.08230930804181,33.64420758908655 L 73.08230945866703,33.64420747102705 L 73.08230959731851,33.64420733911131 L 73.0823097227267,33.6442071945472 L 73.08230983374331,33.644207038658415 L 73.0823099293518,33.64420687287237 L 73.08231000867673,33.644206698707066 L 73.08231007099178,33.64420651775726 L 73.0823187008061,33.64417692982244 L 73.08231874728213,33.64417673447569 L 73.08231877393591,33.6441765354532 L 73.08231878049875,33.644176334761134 L 73.08231876690449,33.64417613442249 L 73.0823187332902,33.6441759364567 L 73.08231867999466,33.644175742859275 L 73.08231860755512,33.644175555581704 L 73.08231851670179,33.644175376511754 L 73.08231840835046,33.64417520745447 L 73.08231828359334,33.64417505011396 L 73.08231814368798,33.64417490607623 L 73.08231799004464,33.644174776793214 L 73.08231782421205,33.64417466356807 L 73.08231764786186,33.64417456754214 L 73.08231746277164,33.64417448968336 L 73.08231727080715,33.64417443077656 L 73.08231707390341,33.64417439141552 L 73.08228378747674,33.644169460093046 z\" /></g></svg>"
      ],
      "text/plain": [
       "<shapely.geometry.polygon.Polygon at 0x7f40d2ad67f0>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gdf_out.iloc[5].geometry.buffer(0.000002, cap_style=1, join_style=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# buffer polys\n",
    "gdf_out['geometry_buffered'] = gdf_out.geometry.buffer(0.000002, cap_style=1, join_style=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Affine(811141.4381028487, 0.0, -59276727.53347101,\n",
       "       0.0, -811141.4381028487, 27295907.604403656)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "~raster.meta['transform']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_out['coords_pixel'] = gdf_out.geometry_buffered.apply(lambda x: Polygon([np.round(coords*~raster.meta['transform']) for coords in x.exterior.coords]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_out['coords_geo_rounded'] = gdf_out.geometry_buffered.apply(lambda x: Polygon([np.round(coords,6) for coords in x.exterior.coords]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geometry</th>\n",
       "      <th>geometry_buffered</th>\n",
       "      <th>coords_pixel</th>\n",
       "      <th>coords_geo_rounded</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>POLYGON ((73.08295538706595 33.6439754184318, ...</td>\n",
       "      <td>POLYGON ((73.08295550451 33.64397342188305, 73...</td>\n",
       "      <td>POLYGON ((3886 5887, 3886 5887, 3886 5887, 388...</td>\n",
       "      <td>POLYGON ((73.082956 33.643973, 73.082955 33.64...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>POLYGON ((73.07865527387217 33.6441800683143, ...</td>\n",
       "      <td>POLYGON ((73.07865335514117 33.64418063264694,...</td>\n",
       "      <td>POLYGON ((396 5719, 396 5718, 397 5718, 397 57...</td>\n",
       "      <td>POLYGON ((73.078653 33.644181, 73.078653 33.64...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>POLYGON ((73.08394535005179 33.64416157585504,...</td>\n",
       "      <td>POLYGON ((73.08394356119742 33.64416068142785,...</td>\n",
       "      <td>POLYGON ((4688 5735, 4687 5735, 4687 5734, 468...</td>\n",
       "      <td>POLYGON ((73.083944 33.644161, 73.083943 33.64...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>POLYGON ((73.08381836849819 33.64418993095924,...</td>\n",
       "      <td>POLYGON ((73.08381655613191 33.64419077673016,...</td>\n",
       "      <td>POLYGON ((4585 5710, 4585 5710, 4585 5710, 458...</td>\n",
       "      <td>POLYGON ((73.083817 33.644191, 73.083817 33.64...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>POLYGON ((73.08082628858956 33.6441800683143, ...</td>\n",
       "      <td>POLYGON ((73.0808243073473 33.64417979503951, ...</td>\n",
       "      <td>POLYGON ((2157 5719, 2157 5719, 2157 5719, 215...</td>\n",
       "      <td>POLYGON ((73.08082400000001 33.64418, 73.08082...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                            geometry  \\\n",
       "0  POLYGON ((73.08295538706595 33.6439754184318, ...   \n",
       "1  POLYGON ((73.07865527387217 33.6441800683143, ...   \n",
       "2  POLYGON ((73.08394535005179 33.64416157585504,...   \n",
       "3  POLYGON ((73.08381836849819 33.64418993095924,...   \n",
       "4  POLYGON ((73.08082628858956 33.6441800683143, ...   \n",
       "\n",
       "                                   geometry_buffered  \\\n",
       "0  POLYGON ((73.08295550451 33.64397342188305, 73...   \n",
       "1  POLYGON ((73.07865335514117 33.64418063264694,...   \n",
       "2  POLYGON ((73.08394356119742 33.64416068142785,...   \n",
       "3  POLYGON ((73.08381655613191 33.64419077673016,...   \n",
       "4  POLYGON ((73.0808243073473 33.64417979503951, ...   \n",
       "\n",
       "                                        coords_pixel  \\\n",
       "0  POLYGON ((3886 5887, 3886 5887, 3886 5887, 388...   \n",
       "1  POLYGON ((396 5719, 396 5718, 397 5718, 397 57...   \n",
       "2  POLYGON ((4688 5735, 4687 5735, 4687 5734, 468...   \n",
       "3  POLYGON ((4585 5710, 4585 5710, 4585 5710, 458...   \n",
       "4  POLYGON ((2157 5719, 2157 5719, 2157 5719, 215...   \n",
       "\n",
       "                                  coords_geo_rounded  \n",
       "0  POLYGON ((73.082956 33.643973, 73.082955 33.64...  \n",
       "1  POLYGON ((73.078653 33.644181, 73.078653 33.64...  \n",
       "2  POLYGON ((73.083944 33.644161, 73.083943 33.64...  \n",
       "3  POLYGON ((73.083817 33.644191, 73.083817 33.64...  \n",
       "4  POLYGON ((73.08082400000001 33.64418, 73.08082...  "
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gdf_out.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_out.geometry = gdf_out.geometry_buffered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_out.geometry.to_file(f'/media/splash/00F2D259F2D25310/Data/final/outputs/arid.geojson', driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Crop detected building images for further classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cog_rst = rasterio.open(COG_URL,'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cog_rst.meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = gpd.read_file(f'/media/splash/00F2D259F2D25310/Data/final/outputs/arid.geojson')\n",
    "df.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cog_rst.meta['crs']['init']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyproj import Proj, transform\n",
    "\n",
    "condition = 'test'\n",
    "grid_num = '042'\n",
    "CLASSIFY = Path('/media/splash/00F2D259F2D25310/Data/final')\n",
    "CLASSIFY.mkdir(exist_ok=True)\n",
    "(CLASSIFY/f'{condition}_{grid_num}').mkdir(exist_ok=True)\n",
    "\n",
    "for i,row in tqdm(df[(df['geometry'].type=='Polygon')].iterrows()):\n",
    "    \n",
    "    poly = row['geometry'].buffer(0.00001) # padding around detection to crop\n",
    "#     print(poly.bounds)\n",
    "\n",
    "    inProj = Proj(init='epsg:4326') \n",
    "    outProj = Proj(init=cog_rst.meta['crs']['init']) # convert to cog crs\n",
    "    \n",
    "    # convert from geocoords to display window\n",
    "    minx, miny = transform(inProj,outProj,*poly.bounds[:2])\n",
    "    maxx, maxy = transform(inProj,outProj,*poly.bounds[2:])\n",
    "    ul = cog_rst.index(minx, miny)\n",
    "    lr = cog_rst.index(maxx, maxy)\n",
    "    disp_minx, disp_maxx, disp_miny, disp_maxy = lr[0], (max(ul[0],0)+1), max(ul[1],0), (lr[1]+1)\n",
    "\n",
    "    if disp_maxx-disp_minx <= 150: disp_maxx += 25; disp_minx-=25; \n",
    "    if disp_maxy-disp_miny <= 150: disp_maxy += 25; disp_miny-=25;\n",
    "\n",
    "    window = (max(disp_minx,0), disp_maxx), (max(disp_miny,0), disp_maxy)\n",
    "    data = cog_rst.read(window=window)\n",
    "    \n",
    "    pk = str(row.id).zfill(5)\n",
    "    tile_bgr = cv2.cvtColor(np.rollaxis(data,0,3), cv2.COLOR_RGB2BGR)\n",
    "    cv2.imwrite(f\"{str(CLASSIFY)}/{condition}_{grid_num}/{grid_num}_{pk}_{condition}.jpg\", tile_bgr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:fastai1] *",
   "language": "python",
   "name": "conda-env-fastai1-py"
  },
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}