notebooks/validation_with_rivgraph/1 - Orinoco Products.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create Orinoco products. This creation of the network products mirrors the GRWL approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from pathlib import Path\n",
    "import networkx as nx\n",
    "from orinoco import (get_distance_in_channel,\n",
    "                     get_distance_segments,\n",
    "                     get_undirected_channel_network,\n",
    "                     direct_channel_network,\n",
    "                     export_edges_to_geodataframe,\n",
    "                     export_nodes_to_geodataframe,\n",
    "                     get_map_centroid_from_binary_mask,\n",
    "                     add_flow_attributes,\n",
    "                     split_tuple_pairs,\n",
    "                     get_segment_df,\n",
    "                     get_geo_width_df,\n",
    "                     update_graph_with_widths,\n",
    "                     get_width_features_from_segments,\n",
    "                     add_width_features_to_graph,\n",
    "                     get_array_from_features\n",
    "                     )\n",
    "from skimage.color import label2rgb\n",
    "import random\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Initialize the Tile Directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "site_name = 'mackenzie'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "site_dir = Path(f'out/{site_name}')\n",
    "site_dir.exists()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read Data from the Tile Directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(site_dir/f'water_mask_{site_name}.tif') as ds:\n",
    "    water_mask = ds.read(1)\n",
    "with rasterio.open(site_dir/f'ocean_mask_{site_name}.tif') as ds:\n",
    "    ocean_mask = ds.read(1)\n",
    "    profile = ds.profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fc841250940>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALcAAAD8CAYAAAA18TUwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAWa0lEQVR4nO2dbcwm1VnHf1d3eW9JAdHwpoCuVCCywGZLxTSG7cpLCe2XGppUm0hC0tAW1KQu4ZMfMIjGYGNKQmgFUwoixdgQCq1YY0wQWGBLge3CLkVYWXnpFkFNoGwvP8zcZTo7L+fMy5lzZq5f8uS5n7ln7jn37n+u+5rrnOt/i6piGHPkPVMPwDDGwsRtzBYTtzFbTNzGbDFxG7PFxG3MluDiFpELRGSHiOwUkS2hz28sBwlZ5xaRNcAzwGZgN/AI8ElVfTrYIIzFEDpybwR2qupzqvo2cAfwscBjMBbC2sDnOw54sfD3buCD5Z1E5HLgcoDDDpWzP/ArB4YZnZEkz7/4Y17bu0/K20OLe78BAPvlRap6E3ATwIYzDtaH7z9h7HEZCbPx/Bcrt4dOS3YDRaUeD7wUeAzGQggt7keAdSJykogcCFwKfCPwGIwKzj92/dRD6MwzTxxauT1oWqKq74jIZ4H7gTXAV1T1qZBjMKq5/6VtUw9hcELn3KjqvcC9oc9rzJPsE+eByudshtKYLSZuI1na7hOCpyVGXBQFMre82yK3AWTCTqli4jJWE/eM8BXn+ceu5/6Xtv00Ypcjd0pir8LSkhnQR4TlY1NITVzfr0XuGbCKvr7CdNk/5eht4l44TalIlwtmbHwuNktLjP2YSwXFIreRjIB9UyQTdwTEnNfGPLY2LC2ZkBSEk0pUr8Ii90SshB3DTVvTReZ6AY59oXZ5fYvcRi3Fi65KXKvnV8+tJoViwSK3UUtVWbA8m1kW/RgRvOtrmrgDU0xHUqVO2LHRmpaIyFeAi4FXVPX0fNuRwN8BJwLPA7+jqj/Kn7sauAzYB3xeVe/Pt58N3AIcQtascKUu0By8uEApBoE3CdRnfGO9lz4XkEvkvgW4oLRtC/CAqq4ja4PYAiAip5L1RZ6WH/Ol3IgH4EYyu4Z1+U/5NRdDDDeR0C6c849d731D6XPM2LSKW1X/Fdhb2vwx4Nb88a3Axwvb71DVt1T1B8BOYKOIHAMcrqoP5tH6bwvHGBPR5wIrrih0ycWLx/mcow9dc+5fUNU9APnvn8+3V5nuHJf/7K7YXomIXC4iW0Vk66s/3NdxiEYVXQTje9O4EntThSUEQ5cC60x3nMx4fvpEyZRnmKEtl3KJbgiBt+3TN/UaIrXpGrlfzlMN8t+v5NvrTHd254/L240SY+SsoW72XJfQhsrJu4r7G8Cn88efBv6xsP1SETlIRE4iu3F8OE9d3hSRc0REgN8rHDM7iv+Bvv+RocpsQ93UxlTXLtMqbhG5HXgQOEVEdovIZcB1wGYReZbMjvg6gNxg507gaeA+4ApVXSXNnwFuJrvJ3AV8c5B3EDG+Aq/KbceO4kMJvDzWqnGHLoEG9efuQupGmL7/oUOspS7m2OXz140nllSoyzge0gd4Q/fud19nM5QjU6wcuJTHiqW1ushfFSWL25qEFUN9PRQm7kBUiapLzbdK1OX9XNOELuNwIWTUbsLEHZCqclzTSrquIqk6LrZp/xCYuCOjLeK6ltuqcuq+nx4uxBK1wcQ9KVXluOK2PstJm6oVXcuUqWHinoimCOdSAqxLPdqOGVPYMUVtMHFHg8vS0zZxuqzyWxLWZjYB5dKf6+RHVf7tsm1sfCN2qPHZJM4ErG7u6lbNpRhh2wQ+5nuySZyIaBJwSsKuc4eNBRP3RLjc/MWOa7f7VBes5dwT0rQCMPYI3lYzj+HitMgdAUMtPw0tqLY1LVP3U5q4A1O39gP6d5tPHe1d17KEwsQ9Ej5rnLv0Gk4pnhgvrCpcmhVOEJHviMh2EXlKRK7Mtx8pIt8WkWfz30cUjrlaRHaKyA4ROb+w/WwR+V7+3BfzrpzZ0aVnsW7ddYxM3fjrikvkfgf4I1X9NeAc4Ircn8S8Sypom2lsm3Z32S8UrmOItTbv4luyR1Ufyx+/CWwns2Uw75ISVSmIb8eLa8UhhPhd1reshB3DxVjGK+cWkROBM4GHGNm7ZC50iWguVYapptirOouSjdwrROS9wNeBq1T1jaZdK7Z5eZekbsoz1A1XLM0FVTORLu1zU+M0iSMiB5AJ+zZVvTvf/LKIHKOqe4b2LknVlGesiBaDeGIYgy8u1RIBvgxsV9W/LDxl3iU5TetE+jo9TZ13p4xLWnIu8LvAeSKyLf+5CPMuAdrbwPpGvKkjZqz5tAutaYmq/hvV+TLApppjrgWurdi+FTjdZ4AxU66OzC2Spv5+bIayA23LVcdKJVIXW2hM3I64+oC0CbwujenS1Z5quhAKE/cI1Am8XCsu7ts2eTPVtHzKF5CJu4Gu65PL9eDy67XVhuueK3+LgdGMNSu0UEwFfMwsu0bacnUihYVUsWKRu4E+wlpF566d4cVzTyHslEuAK0zcjpRnHrukBqmsspv6/ENh4q6hLvJ2nWX0EUyXVGYo5pTPm7hrqFok1LRPHcWOnLFuBoeOsnOI2mDibqUq5+36n1/1adBUAQmxrmROkbqMOU610NZZM0aUa3vdoc7rOtEUO3WOU1YKbKEopGIFY8yI11TnbqqeuIp+KeVFE7cDdTON5cdTi8VlHflShA2Wcw+KT6mvy7GunxixtahNhYnbEdfJFJeKSJ1/ies45uAzGAITtyN1H+dtH/9dz1OV+pTHULcsYCmRuQ2XNrODReRhEflubsrzJ/l2M+XBLb/1oWqVYNW2uvO3lQ+XJHyXyP0WcJ6qngGsBy4QkXNYqClPec3H1B04KTlVhcbFlEdV9X/yPw/If5QFmvL4piYur9UH1/MuVfhOObeIrBGRbWT2Dd9W1VFNeWL1LWm6qXR1ZyrvP1TUdzn/0nASt6ruU9X1ZF4jG0Wkqcm3tymPqt6kqhtUdcPRR62p2mVyXI13YjCvWVrEXuE1iaOqr4vIv5DlyqOZ8sRK8YasriFhaiFV3WBOPaapcKmWHC0i788fHwJ8BPg+CzXlKacWbb4lba81djrRtVVuDrhE7mOAW/OKx3uAO1X1HhF5ELhTRC4DXgA+AZkpj4isTHneYX9TnluAQ8gMeZI35Snjc5Pnsq9P1K2b+WwqJc6Z6FcFHi5H6ut7jgLiMYZcUTUe1/UdQwtt9Zpta13mmKYk/z2UMeW1ZYo17+K2pv3HHEsdsQSFUCQh7j557djUrRB0PabL81VUTSYtfVo+ibTkg7Kp9uPUx27BZ39X2tKQLuPu0qhQtd13TKlSl5YkI+46+lonDIWvIENF1SUIfLbiBr+10GOnOK6lt65T9lXvI1RbWqwkf0PZRNvsX7kUVje7OMQsYl3Db/G1fct7Q9lLuGybE0lE7tf3HOW9SAjaV8y5TqL0TWXqPjVCMufFU7OO3GWqInBd94tP10zXyD5E1Kzb3/eGuu84UiKJyN2Wc/vQJYJWRb0+ovD5tBiS8qeQy4RTCiwqco9xk1iV4nQVRF1eXjxf19d1JWUxuzJLcY85OzhEo8LquHJa5NJcbLgzS3GHpk8Ur8vp+/RfGhmzF3ddh/gY+Ih8qgi9pE+G2YvbpcY9FUNP8rTRtUUuVWYv7jIh/iNjFktd+TOmi34oFifuUPikJ11WFja9XhtVAnep+aeGs7jzDvjHReSe/G8z5WnAd0awzk2q+LyrcF1YgsB9IveVwPbC34sx5elauSjm+74VlaqSY4h+yDkJ3NW35Hjgo8DNhc2LM+Vpw6UJYYjFWa7na6NOyCEuohC4Ru4bgC8APylsC2LK82PechzieHSdrm96nS6NF239kUMR6jxj42LtcDHwiqo+6viag5ryHMBBjqedni5LWsu4lAfb8vO+xNzW54NL5D4XuEREngfuAM4Tka+Sm/IALMWUZ2i6rtQbagnA3HExwrxaVY9X1RPJbhT/WVU/xUJNeYambiluF4ZqLB5qPFPTp859HbBZRJ4FNud/o6pPAStTnvvY35TnZrKbzF3MwJSnrePc9bnyfuVlA0OWAdvOORcWt557DIZYG123Zrwu/y1G3LGqL6k0Fs+6QTg2fGwWmo51rVr4Ngz7EqOgiyyqWWFqXIVdV2OuE7Vrc4NrFaVrGpJK+mKReyT6RruyjcPUbWExR2+L3B3p0/Heh3LNvGnNiVGNiXsC+jYXl2c6h6jQuJDahWXibiFk00Cf16tbNjvUuYasx4fCxB0AlwhXJ06XOncxL4fhRTfkevOQ2A3lBPSphbvUwX2qNT7rzWON1HZDGRE+a0qKkdtlTUldE4LPOFxIIZqbuEdgyApL15WGQ99IprhYy8Q9Aq7dLD6TLL5uUiHEN1aOPxQm7pFw+Q/v03bmypgpQ6yiXmHijoChPQNdX8O1bh27iOswcUfAVOJxWWw11vqUEJi4A+CzjrvLcUMwhNdhkRhE79r9/nzuN7JNRLbm28y3ZGDGXLbah1UjQ6jGiaFwmsTJ+yc3qOprhW3XA3tV9ToR2QIcoap/nPuW3A5sBI4F/gn4VVXdJyIPk/mf/DtwL/BFVW3sxpnjJM4ciOniG2MSx3xLPIjhY9on8lYdV/d3rLiKW4FvicijInJ5vm0035I5EoPLbNcJoaqGiViidhOu4j5XVc8CLgSuEJEPN+zb27ckNlOeIeliVzbEBdH3NVY5d3GCqOsnQSicxK2qL+W/XwH+gSyfHs23JFVTHhdS7KKpi9axR28Xx6nDROR9q8fAbwNPYr4lnWhae910jAt1C69c69ddJ3FiFXlrtURETiaL1gBrga+p6rUichSZP8kvAi8An1DVvfkx1wC/D7wDXLWqiIjIBuAW4BAyz5LPacsA5lgt6dIJ3/U1VvtCfxFWNUHEUDUxa4eF0XZhQP8lr7EI3NZzJ8oYN2tDrVmJvXJi4o6coafF+77uilgrJEVM3BEQk1BcPAPLlm4xjb+I5dwzZQxLtRWxpSF1OffaKQZj7M/QYlySV2AdlpZEwlDT87GmCFNg4o4MF+PLpud9ImzM+fIQmLgjZuyKR6hG4iIhLyYTd4J0XXgVQ5QOeTHZDWXkVE2UdF2yOqX3yRSYuCMn5BR32XOweP4UMXEnQghhr36nLOgilnMnhKuLVRenqrkIuoiJOyFc1oJXmdP7vH4ddRdJzHm5iTsxius52sTYp45dd6yrg2zd/iExcSeIazRu6nVsE2752CGai0PjasrzfhG5S0S+LyLbReRDZsozDeUVeXUUxVglzCqD+vLrDuXfPRWukfuvgPtU9QPAGcB2YAvwgKquAx7I/yY35bkUOA24APiSiKzJX+dG4HKyvsp1+fOGB2NaRJSn/vucI4YbVJcG4cOBDwNfBlDVt1X1dcyUJzjltMG1D9OXqmgfQyT2xSVynwy8CvyNiDwuIjfnXfCjmfLM2bekD12iYR9x9pkVjQEXca8FzgJuVNUzgf8lT0Fq6G3KM2ffkjFIMaqGwEXcu4HdqvpQ/vddZGIfzZTH8GPob3EoE9PCKx9axa2q/wW8KCKn5Js2AU9jpjzBCCWqNlOesXL8sXCtlnwOuE1EngDWA38KXAdsFpFngc3536jqU2RmPU8D9wFXqOq+/HU+A9xMdpO5i8yYx2ghRKf6UDl1TLm5q1fgtjwH/nVV/biq/khVf6iqm1R1Xf57b2H/a1X1l1X1lKL/tqpuVdXT8+c+2+Y2tXSGioIu5psus5FdxjNlSmPd7wtjKIu2vucaEnOcSowhFir59FjG6mzVBxN3pBRTia7NA3VT7K77p46JO2KKVYohImvdhE6X9d8pYOJOBB/Lhy6vNXTUjuEiMXHPnKbcfUyX1hhSHBN35PQ13qnbt+7bFupyfdexxYSJO3LGrG7UfceN64xkyMpLF0zciTJke1fV2m2f78zxHVsobBJnwczFn8QmcRaIT107llRiSEzcEdHFl6SJIaKx60xp3277MTBxR8QUqUHbLGbdBFJVd/wUrrFNmLgTo0pAQ0ZL3wmemMRcxsQ9A4YSmK/hTuyYuI3kRVyHi7XDKSKyrfDzhohcZaY88RDbjVwsuPRQ7lDV9aq6Hjgb+D+y74I3U55IqLuR69M547tvjBeXb1qyCdilqv+BmfJERd81KF3OF7uvia+4LwVuzx+bKc9INEXBtg71sc8/xvnGwlncInIgcAnw9227VmwzUx4Ppii9xVajHgKfyH0h8Jiqvpz/baY8M6SL/VqM+Tb4ifuTvJuSgJnyzJ7UVgGWcfXnPpTMeOfuwmYz5UkE36aDtgVVsUbqMrbkdcHU+Yq0VV5i+8YzW/JqOOHyXTsxCbsJE/dMGaKcV6ygjNlMPBYm7pniI8K+X7Aaaw5u4k6AscXT9hUhXRuFp8bEnQBjT6P3OV+sURtM3EkwpoB87Y1TWvNt4k6Asc3nXRuFU6qUgIk7WXyiuY95ziqSF7+GO1VM3AlR9hkZy+4s9qWsrpi4E6LK/szlK7J9ST1irzBxJ0zfHLjJJDPliL3CxJ0wbdPkS8fEnTC+KUmXNdopXyQm7gVRztF9Lo4U0xQTd8J0EVxR4F3Xn6SCa7PCH4jIUyLypIjcLiIHm29J+viaaqYWvV1MeY4DPg9sUNXTgTVkXfDmW5IQdRM2TeXE2L85oQ3XtGQtcIiIrAUOJWvsNd+SRPGtl7cdHysujlP/CfwF8AKwB/hvVf0W5luSFC7LVstpStXjob4TMwQuackRZNH4JOBY4DAR+VTTIRXbzLckAK5mPk1pyBzWlKxY67DPR4AfqOqrACJyN/Ab5L4lqrrHfEvioM3Mx7VKUrdKcI6rAl8AzhGRQ/PqxiZgO+ZbkhxdKh7lY1IqCbrk3A8BdwGPAd/Lj7kJ8y2JljEEmGIp0HxLjE7E1A1vviULIVS6kMLNp4nb6EUX48xQmLhnxpRf9xebwE3cM2MqgcUocBO3MRhta1VCY+I2BicWgZu4Z0YMpbkVUwvcxG2MQgwXmYnbGI2pBW7iNmaLiduYLSZuY7aYuGfE1NWJ2DBxG7PFxD0jpq5OxEb067lF5E1gx9Tj6MnPAa9NPYiexPwefklVjy5vdOmhnJodqrph6kH0QUS22nsIj6UlxmwxcRuzJQVx3zT1AAbA3sMERH9DaRhdSSFyG0YnTNzGbIlW3CJyQe7vvVNEtkw9niIicoKIfEdEtue+5Vfm25PzLBeRNSLyuIjck+p7qEVVo/sh8wDfBZwMHAh8Fzh16nEVxncMcFb++H3AM8CpwPXAlnz7FuDP8sen5u/hIDJD0V3Amvy5h4EPkRmFfhO4MPB7+UPga8A9+d/JvYe6n1gj90Zgp6o+p6pvA3eQOc1GgaruUdXH8sdvknknHkdinuUicjzwUTKLuxVJvYcmYhV3ncd3dIjIicCZwEOM6Fk+EjcAXwB+UtiW2nuoJVZxO3t5T4mIvBf4OnCVqr7RtGvFNi/P8qERkYuBV1T1UddDKrZN+h7aiHVtSZ3HdzSIyAFkwr5NVe/ON6fkWX4ucImIXAQcDBwuIl8lrffQzNRJf81NzlrgObIbl9UN5WlTj6swPiHLLW8obf9zfvZm7Pr88Wn87M3Yc7x7M/YIcA7v3oxdNMH7+S3evaFM8j1Uvq+pB9DwD34RWRViF3DN1OMpje03yT56nwC25T8XAUeRfbPbs/nvIwvHXJO/lx0UqgnABuDJ/Lm/Jp81Dvx+iuJO8j1U/dj0uzFbYr2hNIzemLiN2WLiNmaLiduYLSZuY7aYuI3ZYuI2Zsv/AxKpLW1ZWRSLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(water_mask, interpolation='none')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate Orinoco Products\n",
    "\n",
    "## Obtaining the Distance Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fc8801507b8>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALcAAAD8CAYAAAA18TUwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deZQURbb/vzer0XGdEYeZx7hAg8i+t4DsA91d1bijILgAimwiKiDQiGyisjeKOiouA24wuIz61K7uFkVcQVbZbNkd3/BG37wzM/jGoz+67u+PyKrKqsqt9szs+JyTpzIjIzOj4FvRN27cuEnMDInEiyj5boBEki2kuCWeRYpb4lmkuCWeRYpb4lmkuCWeJefiJqIAEdUS0UEiKs/18yX1B8qln5uIfAC+BlAC4FsAXwAYzsz7ctYISb0h1z13NwAHmfkwM/8MYB2Aq3LcBkk9oSDHzzsPwF80x98C6B5fiYjGAhgLAGeccUbXVq1a5aZ1Eleybdu2/2HmRvHluRY36ZQl2EXMvArAKgAoKirirVu3ZrtdEhdDRMf0ynNtlnwL4ALN8fkA/prjNkjqCbkW9xcAWhBRIRGdAmAYgLdy3AZJPSGnZgkznySiOwBUAfABeI6Z9+ayDZL6Q65tbjDzuwDezfVzJfUPOUMp8SxS3BLPIsUt8SxS3BLPIsUt8SxS3BLPIsUt8SxS3BLPIsUt8SxS3BLPIsUt8SxS3BLPIsUt8SxS3BLPIsUt8SxS3BLPYiluInqOiL4joj2asoZEVENEB9TPczTnZqoJd2qJyK8p70pEu9VzK4lIb7GwRJIx7PTcqwEE4srKAWxg5hYANqjHIKI2EOsi26rX/EFNxAMAT0Cka2ihbvH3lEgyiqW4mXkTgP+NK74KwBp1fw2AqzXl65j5J2Y+AuAggG5E1BjA2cz8GYsUV89rrpFIskKqNvdvmfk4AKifv1HL9ZLunKdu3+qU60JEY4loKxFt/f7771NsoqS+k+kBpVHSHVvJeCInmFcxcxEzFzVqlJBISCKxRari/ptqakD9/E4tN0q68626H18ukWSNVMX9FoCR6v5IAG9qyocR0alEVAgxcNyimi4niKiH6iUZoblGIskKlnlLiGgtgP4Afk1E3wKYC2ARgPVENBrANwCGAAAz7yWi9QD2ATgJYCIz16m3mgDheTkNQKW6SSRZI6f5uVNBJsKUWEFE25i5KL5czlB6jBLfUJQoQ/LdDEeQ83RqkuxSU7c+301wDLLndhF2e+TS027KckvcgRS3m1B8KD1luGW16h9fhP/0m3PQIGcjxe0iak6uA8jef1nVv1/IcmucjxS321DEZK//rFEJp/y/vDXHjXE2UtwuICxk/1mjUP3jiwCAqhOrE+pV/fO5HLbK+Uhx5xntINFowBgWcrygA+eOTagbaHhbxtrmduQkjsMoKRgmbGuJbeQkjkuQws4cUtx5wspn7T9zpOl5iTVyhjKPlBQMA/l8qP7ppYRzVT+s0bnCOZQUDAM4BDCjJvRKvpuji+y580CJMgQgAoXdeqffnJGeuuzCu9O+h11IIdTUrXessAEp7pwSNkVqQq+gpm49qn9eG/FbgygpP7WekCu/eTgj7bSFzcmkfOL8FnqAEmUISpQhur1c9Y8vitlERYE220Wg0TjTe+ZUyHH4zxwJ8jlfOnbyllxARB8Q0X4i2ktEd6nlMneJCSW+oZH9mtArgOIzqa1OwCgU8V0Hv38qq+1Lh6of1gC+2O9TUjAsT60xxs7P7ySAqczcGkAPABPV/CQyd4kJFPefTwqh9NQbTSP2gv/7jBD4b29HoPHEbDfRkrLCKcYnlah0Sk+90ZEuTDt5S44z83Z1/wSA/RBpGWTuEhOqf14LANHoPJ8PUCgyfW5E8PunAIVARCg7b5JhvUCr8oy11YjKIxWm5/1njYL/9JtR/dNLjgyzTcpwIqKmADoD2Iws5y7xCuHoPPL5EnpzI4LHH0flXx+L6R0T6ny1KCPts0v4L0mg0TgEGt6Gqn8+J8IBfD74zxxp+aPNB7bFTURnAngNwN3M/C+zqjplSeUu8WJSHj071YrKvzwCAChrNjUbTbJN2XmTQCTMJRABPh8CDW9D4FejUXVitWN98rYmcYioAYSwX2Lm19XivxFRY2Y+nuncJcy8CsAqQMSW2PwueSXsuzZb5pXq+Lny8PJUm5UZFCXyQ3MTdrwlBOBZAPuZWWuEydwlWogsfb/Bfzyb9G0DLWcYnvN3npP0/VLCpU4tO2ZJLwA3AxhARDvVbRBE7pISIjoAoEQ9BjPvBRDOXRJEYu6SZyAGmYfgodwlYZs60wOrYO1iw3NsYpNnirLCKXn1qaeDpVnCzB9D314GgIEG1zwI4EGd8q0A2iXTQLdQ/fNalJ52U8qmR0r4svussoummQ5qnY57W+5Aqn98EVAU3SVgYcp+d0fK9/d3mB1bkOUfEivkWpMEkOJOGqtQ1aof1oBUb4IelX99LOY4cPF08dl6puWzq75cAH+XuZHj6s1ZtrkVBeyCaXYj3NvyHFGiDImZSreaRgfUgSMpCfEhZU0mRz7Ds3/Br5eIz/0LI/UC7WYl3LP0knkAgKrt85Nqf1oQmdr8TkfGc1uhekFKTxmO6p/XRsJUzQg0niii/YgQaDwRweOPAwAqj62I+TQiuCdhuAIQobT7/WAC2KfgvU/uS/67JIuLe21AituSmrr1IhEOKcITYmMiJnj8cTHhoVBE2Hbxd5iNqi8XRI5Lu98P9hFqtsw1uSpLuNjeBgAws6O3rl27cr4oVoZE9ktPu4lLzxjBpWeMSOoegd/dkdKzSzvPSem6TFHa8b68Pj8ZAGxlHe24++9ODqn69wtimplIBAzZXVig6enLmt9j65LSrvNya1vrPX/nAuuKDkeK24SISaKh6sRqMWFj9082EcoKp6Cs2VRUHlpm65LqbfOSamf/0swFUZV2m+8ZVUib2wIOMUpPvRHkU6L599TeOHDuWAT/vsr0eu3gseyiaWCfknEPxMbqDIa/EqE6H/Z9FvDIbzR71Jxch+qfXgJrkhcF/74KqKsTWxJUHlyKYO3iGJ+2v/Mc+LvMjbj6tBT3esDwXv0DmfmBDBiwMOaYszzrmUtkxqk08P/y1qzk5yvu/QDe+9jY1devbDE+rDQOqLLL74sXgUIMqmOAGRQSWjB7thMxyjiVd2+I1ZZPb0k8xb7rueSUG5iZhefkzJE5fX7/gQtNz3cfvszWffr5F3H/YvN7uQlIb0l6lChD1CQ0IeHvDoUAiN470PA2y9XqGcHCYvj8ZbGo4ZKRxsvD+g1aAibgg5rsL1PLN1LcNolJyxDSN+UCvxlv617+TrMTykp63G953Qc15bZs7S/WGC/sZUJGTBo3IMWdBJFVNnpT8KrgA40nWkb+6fmQaz63FwS1MTgD/QYtMa1j1nNb9f5eQoo7ScKr2gEAzIYDyrIL7krp/gP7irgSrRej72VCzL2vWQoA+PBdEUnY50px3OvaWP+5ac9tIzbGK9hZZvYLItpCRLvUpDzz1fL6m5QnxBHXoP+sUSLfCAAwR02WUCil3H0bNomIwPffj7oLN70jxPzxn6fFViYh+E9ei858dr3NuNfufc1SfPTWNMPzXsNOz/0TgAHM3BFAJwABIuqBepqUp8Q3FOCQSIP2wxpUnViNwK9G56Ut7IsKvucQ0Xtve0b22mHsJOVhZv5BPWygbox6nJQnxjSBGr8dXrUSFpDF8iyzhb92+fj1aC/MJo/rMVysntf28PUBWzY3EfmIaCdE+oYaZs5qUh4n5y0xSt0Q/P6p6Or3sLWlTWx58fQYQYen4O0GU1nx2Z+M70MM9BxqL67FS9gSNzPXMXMniFwj3YjIbJFv2kl5mHkVMxcxc1GjRo3sNDHn6L3ENPi3P0QDqtTPsmZTUXbRNAS/XqIfU5KDGeLP1k2tdyYJkGTgFDP/g4g2QtjKWUvK4waq/v2CWAisKDEek/g1kpbUhTLart7XLI0ZePYcugyfrr/H1GzxKna8JY2I6Ffq/mkAigF8BZmUR6QSi3MFJpOdteyCuyyXnPW5amlSbaK4WK5P1wtzhQnofmOeM1flGr05ee0GoAOAHQC+BLAHwBy1/FwIL8kB9bOh5ppZEEl3agGUacqL1HscAvAY1MAts81JsSWZxP8ft2f8nn0uX8x9Ll8cU9brmiXMzHzJTfbiTtwIDGJL8h4YZbU5UdzFdB2XNBgWOS49cyT7f3kr+88ZHSnT7sfjbzQu423qWyZE3WfQYu5zWaLIwxSNWJ7xZ+cbI3HXQ0ssfWpCr4A18SVVJ1Yj+I9no5M5QMx+PMHvnsx4m8LhqpvemS6G6Trj1KJbKuT0u8Sa8JsEtLkB7a6rtHqFdSozmxuDMxKWm4Wn58Ns/eOUyI+gXqDXnTtpc6JZoodRbHfJL27Ur3/2LcKUaTiG/b8ey/7fTOBA44mWq+WNbGe9WO/eVyzRrdvtRm/Z3zAwS+RKnBwQTuijxX/WKPGmBUUkc8/WC556D14aM5MZpscNyyPx327HaCVOvRJ3iW+oaXL4TBOThg3qS6BIiV1snAT+LnNRtX0+SrvNj1nEO7DfQ9jw4b2G1xkJ3CsYibte2dxhYZeeemOkLLySJlvP027hl6omI2x/h9kJ2V0TVqdbdVCcGBYLeH9Kvl713PH4zxolXlpUD+g9WAwuvdiDy55bh3yGq2r/etjBaHlZeHGDFR+/Pi1h9hIALh3m4VlLvVGmkza3eEsyib/hmKzcVzux0+vqqCfF7qp5pwI5iZMZrJLPJ3WvuAFnGKssVkZYLh7WWKDa4Cr26MSOFHeSxKyCV0lV8DV162Peme4/c6T9BJs6NPjXz6bnySAA0asRgx79WrlFT/C2r9W8M73qhzWRKMNU8qDUfJqYMkKL4eyk7LklOSUba6dDrGu6eHUhgxS3Q8lGcNXG6nJd+/qL1caLit2MFHc948PKGZE8KFrMUkK4FSluhxGf7D7p67vOs6yz6Z3p6HtFrMC3PTMFXcZ6TOB6/kG9DYAPYkXO2+pxQwA1ECtxagCco6k7EyKlQy0Av6a8K4Dd6rmV8PBKnGJlSMw7dZKl9LSbYo7DUYOZROvrDtN5rPsWMyDdlTgApgB4WSPuJQDK1f1yAIvV/TYAdgE4FUAhxJIyn3puC4BLIcbnldAsQTPanCBuf9t7892EBAJNJ6d9j/AStHg63+YugRuJ227ekvMBXAZAu7yk/iTlsRl/YzWVn67JEYOiRN4+nCofvz5NdwHy9qejA0w3myp2be6HAUwHoJ0GqDdJeYJ7H7JXMfyuHAMfdXxMd7JoUyRXHloG+n8n07ofoO/7Lro1Kujtq9zrSbGT2uFyAN8x8zab96wXSXn0CE+bZ2vhQdg9WHbeJABA5REhwrR6cJ1Zy63PxQq6yzh39t52eu5eAK4koqMA1gEYQEQvQk3KAwD1MSlPNjGKOQlT+V+PxhyH3x+fCmTH4sps3qDcoWeIG20A+iM6oFyK2AHlEnW/LWIHlIcRHVB+AaAHogPKQVbPdMKAMpcU03UpX5vqwLdvYFHKz3QCyEJU4CIAJUR0AECJegxm3gtgPYB9AIIAJjJzOJJ4AsSg9CCEF6Uyjec7gkDb2OVdZivXA7+93dY9S3xDI4PPcE5CO3HnwT32YrsTIELfy5Yk+L5dj57inbS5oef2t5sV2Q9ceHdK9yim6xJ6bb2V8/5G49j/mwniWeffyYEL7krpeXr0vmIJ975yCfe6Zgn3HLxUt44Tk/pAZpzKHYHCKYllFikbwhT7ro/sl552k3gl4Nm3JPU8f6tyW8+yS/dhzl7MYCRuOf2eBSoPxy7dKmsyOWEQCIj47XhqTq6LLkHz+VD1w5pIWmQjF2P884L7xft04hcWxzPw9/ZcnPHx3kW3uMR7oqd4J21u7LmZmf0XT0/r+tIzRoj7/PJW8alZepaNJJpWODmRJmTPnRp9L09tkKWbaD4Jqn4Qk7/BfzwrPlUfeqDxRFBBNK16WWGOJllcGPItxW1BNgL5k8nhHU/w+OOo/MsjkePKIxWmrx6xEyVoB1YowRxx/OSOXnfupM2tZkmu8beYFtkv7TQ7K8/oeovzPCXM0izJK4FW1u9Z1/rAw9Pr2n1T/3nrmQh+vQSB1uLdlVU7rF+1nQzhNxIzuWxRg57inbR5sefWcxXaui4uzNXfZmZCHa3PPczAngt079evJDEzrBFdbnVmr81s3HPnXbxWmxfFbUTg/DsTy1RBB5rfY/s+pR3vy1ibwsTHeHcaHz3uMKki489LBiNxS7MkCxT3eiCl67QDxUiZGvlXeTCJFz8Rwd9lrnU9m3QZVxET4w0AO56IHjs174lDm+Vu3vvkPpT0XGBZLz4uJeG8JpTVjt0epmrH/ajaPt92fSuMMlJ1nCTexPblw5Mz9qxMIsWdJawS5AA2FkFocpcEv1pkUtGYkh4ZGFwaiHvXo84UdRgpbgfg76T/Q0h3IggAaj6fY3guPsWDUb5urQniJqS4HUDVTmsTJhtseidq9vS6dlnkhaxarNZQtp9i/pLYfCLFnQN+X2zPpCjtFmsnB9rfl43m6PLJawaznBYrdXZX6JsmrWflX/R2V78fJaLdRLSTiLaqZQ2JqIaIDqif52jqzySig0RUS0R+TXlX9T4HiWglUTYS4jkQm99S+zoQf4fZCO5OzeuSURSg8/gKdJ5QgU4TRS/e8U5r4e5/0AH2uJ5/MH4DcBTAr+PK6k3eEkki7e/Mr29bC7Lg564/eUsywIABC/PdhEiGV71cgWZ0nhBrd3ectMLewuI8Y1fcDKCaiLYR0Vi1LGt5S7zI+++LuA+777DJBhuDMwDEDiTtoPWWdJpYATCwa6UDzA4L7Iq7FzN3AVAGYCIR9TWpm3beEqcl5ckkGzbNSroXL+6dvu2dalx6mC7jKtB5fAV2Pj4Fux4Twm43TdjebWbmf/Cohy1xM/Nf1c/vAPwZQDdkMW8JuzQpjx0GDFwY6cXtkm5MeZ8rl2LT26kn7ukyRvTWO56M9XfvWSpE7lgTRc8Q59iB4xkAztLsfwogAJm3JGX6F0ej8fpcttikZnJcOiRxxbpRssue1ybW7XG98VKyzmM8GBUIoJkq1l0A9gKYpZafC2ADRArjDQAaaq6ZBeElqYXGIwKgCMAe9dxj8HAKYzPshpr2H2hcT/sDseLS65bqCj9ZOo1bzp3GL+eOt0eF3m5y/r0mKYs735sXxZ0Legw1FnOPoUtNe2k7dLgjKupW98qQV0kKWL5b0gD2Gdvpn/3pnrTDVHc9Ohkd7hIDSUdM2OhQr9/9Lkmd9pOFsHevyL+w5bvfHUzvwUksRMgyRbdUoOtoi2CpqSsAEsJuN20F2pY70xUoe26Pcsmoioy+gq9t+QogJNx+YRegUzDquQv0KktyT88hy/DpK8b5R5Ilk8Juc+8K7FvkLEHbQZolDiEs7B7Dl1vUNMfKpEgFo2VmTkeK22F8vnZqZL/bCHOht74v0dbd9qz9HrvT7RUJQVG6uFTc0ixxMFuen2p6fv8D6ZkKO/9g74eQyZ67yaqlODZ2WuZuaILsuV1Im3uT806EFxl0uDtFr0YGVZIrYQPSW+J4uo6uAHHiG8ZyQau5K0SvTUDtXOcOKKW3xKWEbeiiWyqw9Y/ZFXibmSuwb+FktJ61QsxgKkDtPOeK2grZc0sARE0dJudOpxshZyg9gJ0Mq+2nrhAziCrhBQVW7HtosmtdfoboRVM5aZNRgbFo46o7TMx8NF7L+/Tv2XyRfjy3UXkugUFUoLS5Xcb2VVPQeUIFmIBdJq68tuUrEPKJ6fJ9DyVhZqi9d8t5K4TNPWdyTHmziuU4PCXqojw4w3wcUPjSQzhyo3lOxGwhzRIXsuOJKZYTK6wI25kpmiCn5byoidJyfqK50nLeikgobO08ce3FDwhT6OB0IWKtsO2QL2EDsJ235FcAXgXwFYD9ELlHGgKogViJUwPgHE39mRApHWoB+DXlXQHsVs+tRD1diZMO7aZEzQbtgoF0uXh+BV88X9yvxf3R+170QGpmx4VP6i9vywZIc7HCIwCCzNwKQEdV4OUANjBzC4hlZuUAQERtAAyDWEsZAPAHIvKp93kCwFgALdQtkOyPsb6ze7kwEzrcvSKjebGZoiYIK4wWD1agxUMVKf9tPzYud5M1hugpnmN77bMBHEFcLwvRKzdW9xsDqOVorz1TU68KoqdvDOArTflwAE9ZPV/23FHaTEu+p7bb8+rVa74wWlZY4c33UDYD8D2APxLRDiJ6hojOQBaT8ng5b0k67F2SvP/5wCxhKzdfYuFG1FHCwfLoYPHw5ORsbSdgR9wFALoAeIKZOwP4P6gmiAFpJ+VhD+ctyQa2Yk2sBqBe83HDnri/BfAtM29Wj1+FEHvWkvJIksOOq+/QNIupexMlNKsQobdNH00v1jzXWIqbmf8bwF+IqKVaNBDAPgBvARiplo0E8Ka6/xaAYUR0KhEVQgwct6imywki6qGmLh6huUZigl7cdqZhApov1Tddwu6/o5OsTZMmz6aXti2j6Bni8RuATgC2AvgSwBsAzoFMyuMamj5sbzBYuCz/s42pAJmUx31cPC9zfuwmK80FXrh8GRcuj63T5LHYxD4XPpG877pw7QPMzNzytXlJX2sXKW4JM5tPrjRdkVl3X/N1+m8uzjRG4pbT7w6lxUP69m/pRuN3wMejNwA0mlwpfGR5xoMxDl6fu3f66CHF7VAO3DsFFy0WAg97KwCguv/Dtu8RPwBs8ozxYI8V4Mid7vNlmyHF7WDCEXeHp0zNiBvu2G0iR3e8R6Pp48sAhaP7HkGK2yVoe+HuVfZfla3HsdHRRPRNnlwKKMDRiWpCoAwp4uJXM/Dm4jSR8dwu5GSdz7qSSpPVi3Fs1IzE8qeXCOc2AUfHRzNdcYZek/D1dcZvLs4Vsud2OE1WJSbJ3DbI/kuj9IQNAMfGTBcrGTTT7k1WLcWx8dMiz2zynHn65Jav5b93NkXPheKkTboC9Wn35uyM3OfCVZl7bUmYthlqm10gXYHeYveV5r1m4csP2buRIgaYCYPM5xNf6d35nVm2brnHom25QqZ2qMc0+ePiSFzmsVv1zRc3IFM71EMuWr/A9PyxW2YIBShi4Ok1pLgdhJ4pEE+zdQ/ioj/Ze+nqwaGzLeuQwpFNDyOXXuDDu2KOe9dMx6XVM9J2U2YSKW4HcXSEtTAKCuoy+lbToyPKxf1UcTdbG2urf33dHN2/AMF+jwAALv9oEgDg45Il+Kx0MTb7rX+guUKK22V8fd2chB653Vvp+ZRJYZD6gzk8PJqKoc0b8wCY/wV4u8+jaT07m0hxe4B0vROkiC2+19539by07ptvpLgloremaO/tFSzFTUQtiWinZvsXEd1NRA2JqIaIDqif52iumUlEB4moloj8mvKuRLRbPbdSXW4mSZNuwZm4tDo9V55SH8XNzLXM3ImZO0FkjPo3gD9DJuVxDFsCC/FZaaIrr1twpu17hD0mLV4xdx9qGfzJhMj+0M/G2b4uVyRrlgwEcIiZjwG4CsAatXwNgKvV/asArGPmn5j5CETqtG7qCvmzmfkzdcr0ec01kjQZ+EHi6vYtgYW2r1cUjmx2uOHzMXi91xOR4/WXPmX7WbkiWXEPA7BW3ZdJebKEma/YaCXOht+n94o+Us0SRQmh9Z/nWdZ/ucfTaT0vF9gWNxGdAuBKAK9YVdUpk0l5ksDMV5zMSpxk+GrwXHw1eC6IYLv3djrJ9NxlALYz89/UY5mUx4Psu3oeiBjtk/Cdj/5iVPYalAbJiHs4oiYJIJPyeBYfMXxKCJ3etrfA99lLVme3QSliS9xEdDqAEgCva4oXASghogPquUUAwMx7AayHyEoVBDCRmevUayYAeAZikHkIQGUGvoPEgvHbbjY9Hy/iXVcsABHDpzCKKhOTx4/ZOiKj7csWMuS1HtO9qlzXvu/67qyIzzvs/9bWu+HzMY4aUMqQV0kCis6kTbfgTGwb9CAU1TQh4ph6gz+Z4ChhmyHF7VGsTBEAuhM/PiUU+dzsX4TN/kXwEaPveyKZj9a37XSkuD3Kk11fsF233wax+r13zfSI4H2a3tpqat4L3hJJnhi55das3v/DgSIRz8clYh1l3/emxZgim4qX6powYVztLZHklzXdnsvavfVmPDcVL42YJ2EKKJRQDwBu2jw6K+3KBFLcLmDqzqFZu3d1/4cxaNOdCeVaMZduvDtillz50R0x9V7s/mzW2pYuUtwuYHmn9Wldf/+eK0zPv9t3ZWQ/vDZSa2MXKCEUKCFc/tEkvNXnsbTakkukuF3KR0eb2647p91/6pZf/fHtCWXBfo/gyo/uiPTol380CQoxCigERT8UyLFIcbuIiv0lkf0+TQ9h+7ELTes/sr/Y9HyBom9Hh3vnAgrh7T6PooDqoFAIBUqdbn2nIhNhuogprWtijrs0+QafH2uKHk2O6ta/q/V7pvcz8oBc9+l4nAwpUNQ4TqMfgdOR4nYxG4+2QP+mR1O+3miBwas9nwQgRA5EfwTSLJHkjP5NDxieCx5uk/b9wx6TAgqJzWVmiRS3i6k50srwXKDZvoSyhXsH2b73DZ+PgUKMmzaPhkIhNFDqYmYt3YAUt4tpQMn1pDPbvou5u6+KHJuFrr7c4+mIoBtQCGu6PWc4keNUpLhdjJlZYsT89m+ifNe1AICni543ras1R8ZsHQGFQrYCspyC3cUKk4loLxHtIaK1RPQLmbfEvdSp/+137hiOO3cMN6ynNUeeLnoeTxc976pBpZ2kPOcBuBNAETO3A+CDWAUv85a4CK3P2wdhXqzsvBYrO681nN4Pi7lAY/4UKHWYvOP67DY2Q9g1SwoAnEZEBQBOh1jYK/OWuAhFYy8v6vhazLnlndZHTJV4CpS6GH94A6rDis5/yk4jM4ydjFP/BWAZgG8AHAfwT2auhsxb4iomtXrf9Pyijq9h7u6rMHt3tL+ZvON6NKA6NKC6SO9u9kNwGnbMknMgeuNCAL8DcAYR3WR2iU6ZzFuSA8ym45+u7RPZX7xP3xqc3/5NLGj/Bu79cjDKd12LU5WTUIgjAg9jFtvtJOwlKjQAAAPvSURBVOzMUBYDOMLM3wMAEb0OoCfUvCXMfFzmLXEGXZp8Y3huTMuP8GRtP/wUaoAZbcyn5R/qEE1ycO+XgwEAPgph9u6rsaD9G0m7IPOF5ep3IuoO4DkAlwD4EcBqAFsBXAjg78y8iIjKATRk5ulE1BbAywC6QfT0GwC0YOY6IvoCwCQAmwG8C+BRZn7X7Ply9bszmLv7KvhUuz3sMbmv3dv5bFKElFe/M/NmAK8C2A5gt3rNKsi8JY4lPHP5xqGOGbvn/PZvYk67/0QDqoOPQmignMzYvbOFzFsiSYll+/xQKJQQqZgPjHpuGRXoMT462hx9mh7K+nPuaVMFAHj0qwHwIYTbW23M+jOTRU6/ewxfjmcQJ7V6HwpxjDfGKUhxe4yeTQ/n/JnjW36IBnQSLxzokfNnmyHF7TF2fnOBdaUsMOriz+BDCOsOJpi+eUPa3B4j12aJlhtabAEgvDRXN9+Vt3aEkT23x2jggJjrq5vvwruH2+W7GVLcXqPVBc6Y9FUQMl0plJs2SCRZINBsX15NJECKW5JFBhTW5vX5UtwSzyLFLfEsUtwSzyLF7SF+PN40301wFFLcHqIAPutK9Qgpbg/RoHH2owHdhOPjuYnoBID8+pTS59cA/iffjUgTJ3+HJsycsNjWDbEltXqB6G6CiLbK75B7pFki8SxS3BLP4gZxr8p3AzKA/A55wPEDSokkVdzQc0skKSHFLfEsjhU3EQXU/N4H1YxWjoGILiCiD4hov5q3/C613HU5y4nIR0Q7iOhtt34HQ5jZcRtEDvBDAJoBOAXALgBt8t0uTfsaA+ii7p8F4GsAbQAsAVCulpcDWKzut1G/w6kQCUUPAfCp57YAuBQiUWglgLIcf5cpEOnv3laPXfcdjDan9tzdABxk5sPM/DOAdRCZZh0BMx9n5u3q/gkA+yHSMbsqZzkRnQ/gMogUd2Fc9R3McKq4jXJ8Ow4iagqgM0Ryz6zlLM8SDwOYDkC7qtht38EQp4rbdi7vfEJEZwJ4DcDdzPwvs6o6ZUnlLM80RHQ5gO+YeZvdS3TK8vodrHBqbIlRjm/HQEQNIIT9EjOHE1q7KWd5LwBXEtEgAL8AcDYRvQh3fQdz8m30GwxyCgAchhi4hAeUbfPdLk37CMK2fDiufCliB2NL1P22iB2MHUZ0MPYFgB6IDsYG5eH79Ed0QOnK76D7vfLdAJN/8EEQXohDAGbluz1xbesN8af3SwA71W0QgHMhku0fUD8baq6ZpX6XWmi8CQCKAOxRzz0GddY4x99HK25Xfge9TU6/SzyLUweUEknaSHFLPIsUt8SzSHFLPIsUt8SzSHFLPIsUt8Sz/H9TVJDcyrBpZQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "transform = profile['transform']\n",
    "dx, dy = transform.a, -transform.e\n",
    "dist = get_distance_in_channel(water_mask,\n",
    "                               ocean_mask,\n",
    "                               dx=dx,\n",
    "                               dy=dy,\n",
    "                               # removes areas with less than percentage of total size\n",
    "                               min_rel_area=0.025,\n",
    "                               # Applys a 1 pixel buffer for distance computation\n",
    "                               # for artifical 8-connectivity.\n",
    "                               # See: https://github.com/scikit-fmm/scikit-fmm/issues/32\n",
    "                               # Removes buffer after computation complete.\n",
    "                               # Only recommended for validation here\n",
    "                               apply_mask_buffer=True\n",
    "                              )\n",
    "plt.imshow(dist, interpolation='none')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the distance product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = profile.copy()\n",
    "p['dtype'] = 'float32'\n",
    "with rasterio.open(site_dir/f'distance_{site_name}.tif', 'w', **p) as ds:\n",
    "    ds.write(dist.astype(np.float32), 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Segmentation of the Channel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "rag neighbors: 100%|██████████| 840/840 [00:00<00:00, 23876.66it/s]\n"
     ]
    }
   ],
   "source": [
    "pixel_step = 5\n",
    "segments, interface_adj_segments = get_distance_segments(dist,\n",
    "                                                         pixel_step,\n",
    "                                                         dx=dx,\n",
    "                                                         dy=dy,\n",
    "                                                         connectivity=8,\n",
    "                                                         min_size=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fc8d0702198>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALcAAAD8CAYAAAA18TUwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAZsklEQVR4nO2da4wkV3XHfye72LzseB0C2qwRGAkhQRQxs/ZiwmvlhNiMEbYUPhgJ4UiJJx9CBEFZZAthaodPYaIIWUQRs4BiQsAiPIJtnMASMSRKkGd3vAbbmMWLjfDGxpsHBIcPQMTJh3uru7q7qrq6nreqz1866urb9bj31r9On3vuuadEVTEYhohf6roCBkNTMHIbBgsjt2GwMHIbBgsjt2GwMHIbBovWyS0iV4vIaRE5IyI3tX19w/JA2vRzi8ge4DvA64GzwAngLar6rdYqYVgatK25DwFnVPURVf0ZcDtwbct1MCwJ9rZ8vQPAY4nvZ4FXTO8kIuvAuv96sIV6GXoOVZXpsrbJPVMBYMYuUtUtYAtARCw+wFAKbZslZ4HnJ75fAjzech0MS4K2yX0CeLGIXCoi5wHXA3e0XAfDkqBVs0RV/09E3g58CdgDfExVH2yzDoblQauuwDIwm9tQBGkDSpuhNAwWRm7DYGHkNgwWRm7DYGHkNgwWRm7DYGHkNgwWRm7DYGHkNgwWRm7DYGHkNgwWRm7DYGHkNgwWRm7DYGHkNgwWRm7DYDGX3CLyMRE5JyIPJMouFpHjIvKw/9yX+O1mn3DntIhclSg/KCL3+99uFZG0xcIGQ31Q1VwBXgusAg8kyj4A3OS3bwL+zG+/FPgGcD5wKfBdYI//bQd4JW4F/D8Ab5h3bX+cmpjMkzTuzNXcqvrPwH9PFV8L3Oa3bwOuS5Tfrqo/VdVHgTPAIRHZD1yoql9Xx9iPJ44xGBpBWZv7ear6BID/fK4vT0u6c8DL2ZTyVIjIuoicFJGTJetnMNS++j0r6U6hZDyjHywpj6EGlNXcT3pTA/95zpdnJd0567enyw2GxlCW3HcAN/jtG4AvJMqvF5HzReRS4MXAjjddnhKRK7yX5G2JYwyGZlDAW/Ep4Ang5zgN/PvArwD/BDzsPy9O7P8enJfkNAmPCHAZ8ID/7UP4nCnmLTGpQ9K4Y0l5DIOAJeVZAuimU2UGS6dmGAhMc/cca2zUut/QYZq7R4g7okhQjm6CHGmyNmEhTXMbuXsGRZFC9F4umFkyAIwe9Wj2N93eabMqwcM0dw+geFNEFSxSOBWmuQOEJhx3WQPB0V2bIramOP0C11WtwjR3YFi2gWBdsAGlYbAwsyQg6ObcPdqoxqBh5O4Kl++wtbuS6vVwCH3gqGztrhR4SDtEkci8LoUAIs7qljU2lMh/+oZqXedtKwpvs/t+nKhPGne6Ju9Skds1aEK2dld0a3fF3YkFCNM1uXR7p/v+TNbHyN2NbO2uFCCjuz+j7ykPQiiiin8gu6/LuE4lyI1bNvZV4CHgQeAdvvxi4DhuwcJxYF/imJtxK99PA1clyg8C9/vfbqXAgoWuO610Z0+RWYscpxocaTIlmvO97f4uSe79wKrfvgD4Di4/SSu5Szq/iXXdfHxva/5x8S0JQXPn/dskzZI2bf3M+tRhluDWPr4ep5X3Jx6A0wmtfXNi/y/hCL0f+Hai/C3AhwdL7qkbv8ZGYTtZt3d0jY05Wrxb8uv2jqOPr0fXNnhlcgMvBL4PXAj8aOq3H/rPDwFvTZR/FHgzbg3lVxLlrwHuGjq5R53vu3uxY7qv96gu/sFcY2Pi30c36dwkgXRyF85bIiLPBj4LvFNVf5yT6q9y7hIRWQfWi9atDxDE+bUPnlrgGIc1NribW5qpWAHo9g5y+HIUge2ruUYAfS+KhB3HVVBjPw1nXrwrUWZmSVKi+X/NvRksTmvFAOowt44p3CmS5VVwpsVDqvoXiZ8sd8mCWF9Aa8fIjdGOytdlsUq0dJ26UUBzvhrXvG8C93lZo6XcJQSgFYqIs0lbduVFA7lGHf2fwh2LCqwRa2zwRW5pLypEtdHFC/G/hhw+1Ng16oKmRAUu7ApsWwhAKywiuplveyc9DQvLtBZtWKsqYfiwC9U1jTtdk7d35HaVypVxvMj88yX94IWuH7Xd1vntDUGM3GU6bRPV5A0uQG4cI2YIHn9XJX/2L6086qbtXfe/kbvRTsPPxmlhkjktrAmppx5xiKzWdM5i7eieuEbuJiVy2luVwpob372LXmta249ivjtpe1fXNXI332FTpNOpsiYJototWfo24WTkXrTDUmxOxfuyo4LnyNjOlYLnbqzdgS1EKFRnI/fikmZ3xrUr1ule6zdImDptY1Xiju+VpHHHJnHmQlEVXDXcPEE8uXHN4X9cKKBJUdiUoPOSjFvZL6hN4lSRSW3mC8qdK6EZFWf+pNm4ef8OdbnpZs4Tdd3PJdthZknNHVrQNFlY5hBMa7pOvLZT0ZG7s65zt34vjNxVyYziCb21u9K6H3ge8YrWR7d3ejlozG2TkbsKsXyFIhTV0URNHEvSzmxe0X+K7P1UiTt2UGLkrtyBjhiqTkuuseG0qSd7UW2Yal8X1aQViakB9KORO0Byx8QcmyROc8caPS5vnkD5BM+7fvN1C4fclitwQSRX06h6t97RsRfqi7wXKP/SJfc8M7HKJj7XOC/fpEtSUZLooyuvERTQnE/H5Rv5Bi4pz1FfvrRJeZyHwZsX3kzQzUS1iafrq5kQcyVixrMScjqIRu9JGbMEpwie7befBtwDXMGyJuXxJsh02TS54+5u9IYmzl9kQNv4w9ahVLa5gWcC9wKvYElXv2tGeZxER70dHg842yJarsbW4ZJ61J9lbW4R2SMi9wHngOOqeg/wPHUr2vGfz/W7HwAeSxx+1pcd8NvT5WnXWxeRkyJyskj92kSWPRtPwx/bXYUjcN3u57lu9/PjHSImclnHr9vb2l2ppV65K+uPCtN2+VJgQc19ES4p5q+TnXHqL5nNOPW7wOXMZpy6s2+aeyQZ2lBJ2NuR2y9Xg7c0mTK0SZuZ9lU1SzzZ3gf8KUtqlkx06Gb1iDzVZkk8NleGbZqUIjfwq8BFfvsZwL8AbwQ2mRxQfsBvv4zJAeUjjAeUJ3CD0XhAudZncqfJIkH+RbTporZ55jkjlAETvCy5fwM4hUvK8wBwiy+3pDxVJGrivOolcdM3m7xeOFKLWdK2dN1pqR25yYTrTcFVNkGgXHNF69eg6j/H/vX0awzV9k7jji1WKImt3ZVSuf+agzL25aj/NvXGYe+tCXmxRFmoLVZoQKLxdnF7O19z66J1GGmv8fFxD9Z17tDFzJImJcoiagaRI8ZvJ/Dx1Uq+2xDIMWmKEzmrvM9i5O6q49Ps3GgcYZj5ANQi6efuUzapQn2cwp0ls7mVNmPmkrOP63eeQiNgewded6hkclZX/+lWzLP/2211N1CzueOnfGo7alHDLKilNXFMlk0/z5SZ9u6UrUvIYmbJlPQtq1IVmX5R09DEyJ0q3WivRf3NmQ+iFq9/mqbuU7LL3LYZufshTf2jJM+rSWIMYGInjTtLNqCsDvfaunpeo6GbdU+oKHlDR0VnJnbAMT3oV+4VgKYMKG0N5YJII7YzZUuc6wgTR2qsU8siymfosd3V9B++lvPGtB7DyF0D0rRh8WMT28JIhSYXNhRGVLoag4SRe4mwfuep9AfgdeG/rawMzOZeMnT9qu2mkGZzG7mXEGnDzjoHyl3ABpQ9QOUnOZq/izBr08vhQ6UHxsFiAX/zHtyKnLv896VNylNMtFpw0vTkTETtwU6p/vSo634rJ5UmcYB3AZ9kTO6lScqjAdRh9mbWcI6MyZs+hiWUJjdwCW6d5JWMyb00q98L3+xpbVuQTKVu5mY9D928h6Qvs5dp3Clqc38QeDfwi0TZ0iTlKbqcbE3e77c09feqA7bk2F+OwLEaEvrI0ZTCKPF7jweZc8ktIm8EzqnqbsFzps1oaE75bKHqlqpepqqXFbxmEBi72JqZy46nyOPBYPzQVclapVFK4VRZXweaRTT3q4A3icj3gNuBK0XkE8CTIrIfwH+e8/ufBZ6fOP4S4HFffklKuWEK89IfT8ejVFqovF1g6n37RPnzd4kFI/QOM7a5LSlPE7ZjpWPzbf5Mibpvd+V+q+ItSSG3JeVhNh66an5sRacyxKovb7IdrrfnDYhDlsrk7kK67rRCHTuRoKccQdzrRiaP1axrxSSMKmjrtHO7S44SeGbVs+v+NnJ32tGzZWVymky8VKrgMU0QsA9+byN3Zx2fUZ6hdWNyxv8IMblqW5o2p16z9SxW7277eJY7FjjVEFSrrW5R3KoZZRwLMvKSVD15mfrUvmqoXqhFBZZAB0TKh3rCj9/M0EbOwvhhCxVp5LaowDnQrzXg442qHCwTJFs/eCrfL17pWomrRpIyuRM2THMPBMl47NiUqR0RwS5lM83dEYo8ncl9JvWN/xJlH7vGxkQ8dt3EHr3M9X3j7V6ga2/IMnhLpqVspF32+26yPRtZZbCoyzA8D8lEG80VGLhEKTctniDSRcjVABGn6pZ0I4bgGjRytySdz+KpZmrrUhLNv17XfW7kblPmEYICkywT37skUPq1Q5q5NHL3TOogTx3rLvuwGieNO+YKDABKVxMkStK30l09qkPNFRgmuiNU4rqaQewo/wwhuwZNc7eASf2Ys59OzvQHEc8xXammjyuJ0ppbRL4nIveLyH3xol0RuVhEjovIw/5zX2L/m0XkjIicFpGrEuUH/XnOiMitIkEFbTSHgo/nRG9EARAbQJy5ouh4rWZU5LgAbm3BQd33gOdMlS1N3hKTWenc3Tklpb0lpJN7afKW1NL5IdSh1IRQitdGNYj2TLQthTtFB5QKfFlEdkVk3Zc1lrdkiJh8cXVHdYjNnAVNhmRI7dbuChzth0+lKLlfpaqrwBuAPxKR1+bsWzlvSWhJeeqEIE6PL4BSieinUCW3yeh4VUf0yNerpnM3hb1FdlLVx/3nORH5PHAIn7dEVZ+oO2+Jqm4BWzAMb0kSur3DouNo+UnFa26CVFjQoJvA6inWp+odvPYuYPM+C7ggsf1vwNVY3pLy9uHE9mL2b56kRglmzVBGKfXKm4lM2T8kKTWgBF6EI+s3gAeB9/hyy1tSWgoSOo9Qi5BNtZYp9Dj9RPIhCsVrUtpb0qV03Wl9ldy4FH/nK10j4XEJIfakirfE0BHKDtbW78yxscWtqq8EkdG4ONRMsDb9biiFEbEDGVWqBU6FiXlZXduFzvdURoB4YkfM378jmOYeKiLqXakeAe9TVLqMYsxGmuY2coeCiHDTJvTgNX5G7h5AtZ40bHWiap3agNncPcAkieY917O/L0ZspUisS1/1i5E7aMwjalV1KoUehq3d1YrXSaK9B8XI3UMs6vuO9y8bkbh+8N5Sx6WjPfvGbO7AoSjHdldbyeQ6e21GitZsbkPtEGQizLRR+GussTExSRM6sbNgmtsAOELfzXv9t/6x2TR3z1EojcIoJsp/LWhn380t6PYJ1nj//J17AtPcPUPyTQrN+J+VVM2debGM/VuGae4BYP3gKafBo3xi66ZX4gt6SOK93b/E7LFpnvU8dLoEret4bYvnLiv58dhxjLXGHJ86JnUFUJTcl5mFCSFL6XhuEblIRD4jIt8WkYdE5JWWlKcrqP/MX2gcx4IISd2a3Jrs+jgyUWB0XkHy48LzalnDoubKKKg9bwP+wG+fB1yEJeXpVlNR7wqYCU2uqrrpl5UFsMqmUP3TuFOAXBcCjzK13hFLytO6lDERiq5xTH8FSYLwUfftz5Oy5H45TuP+NXAK+AhuFfyPpvb7of/8EPDWRPlHgTfjFgd/JVH+GuCujGuuAye9dN5xQ5BUG3uB30OXsjb3XmAV+CtVXQF+gjNDslA5KY+qbqnqZap6WYH6GXJs7xgS5Q9vhuhwLULus8BZVb3Hf/8MjuxP+mQ81J2Ux7AgiozLo/yf814mOw686hfmkltVfwA8JiIv8UW/BXwLuAO4wZfdAHzBb98BXC8i54vIpcCLgR11+QSfEpErvJfkbYljDLlogVYnDmU+APGkUdEc48GgoLfk5Tj795vA3wP7sKQ8vREtul8N78/prI1lBpRdS9ed1ukNq5Fs81x6setvolx1oXOkntcPVLXpvkrhjsWWLBnyX0Wi1BknorT3Aim12JIeIUPpLKKL0qIIs4itKHVrka5TQBi5Q4XIOA+2ThYXPsVUOoZc8mqx9ZR9gpE7ZETuQ2QOMQsi+cbJJFTh2L3O3RdW9quK6HrAaAPKMqK1nUuVian1ugZ+bUcTmrdkKDLlxci96RkPwigoKuqWlHWJkbuH0mRyd2WSzLHrMX4gsh6MWEJ6EIzcA5K6SN/ExE0Xk0FG7iWSogRTdRp4ev+0hyfk2O407tgkzhJDUdg+Aa871NvcJDHUJnGWDFH+z4K4gCkIY1lYzTByB4QieUnW2ChOxKjAPkdAjsKxK9NXqWf9bU6vatftHTdCDekh6dqmNpt7MdnaXVGt+7zR2PORZqvnDl7dTepczOYeLJQqAU9buyvcuHpqxu6udtZ2YTb3YFGNgut3nuIamTV3+kLsLBi5DeMYlorvmA8OBWzelwD3JeTHwDuBi4HjuJU4x4F9iWNuBs7gVuJclSg/CNzvf7sVW4lTk2jliRNFg3nVdan6V53EAfYAPwBegCXlCV8WGOzpJi7OJCp+/smHoduBZR3k/h3gX/22JeUJSKrOHqona9F4kdDWW9ZB7o8Bb/fblpSnqRuVQxxtSEOusaG66c8fiHtvoT5L4VHhAaWInAe8Cfi7ebumlGlO+WzhkiflyV7j2NzSrbu5BTnikvfMOtX6iUW8JW8A7lXVJ/13S8ozREQAwiLTH6Gu3lmE3G8BPpX4bkl5BopjuysImv63moK7uaXR+pRGQVv7mcB/Ab+cKLOkPD0RnbtPio0dj8g0Zf8AbfLKA8oupOtOG7RkkHRrd8WROmLEkAkiBeYpASO3ybSkkDvW1Kqobu84X3ZivxCJ7epr5F4aiUm6sETuc2JtpYa1XjK9vbPcsdiSgWKRlTWxt0PRzBjwG3Nez60L1KtNWMhrD6Da7iuqVQGZnJjQ7Z2ZDFYhQS3ktZ9okthpPmoRYHpV0Il0Ynf6nsl56NqmNpu7gD3J7MCv6fNP5zMZDTQD6I/UNpjN3U9UnXKfp12T5481+Y2rYxv72JUrXCMboNqrBQxmc/cUSvXc11m2/Khc3VWIzZTDh4Ilt9ncfUc03hSEeU/9PL2VpTdiwqsIInANG3D48lk7PHCY5u458t+UMOfYDM29xgZf3L4a2fYvgYoSP0az+4eANM3d+YDRBpQVRJsdaI7SPYDq9k6wg0lIH1Ca5h4q6nCO+3OMbkDAvm6zuYeGKOe3FGIvlg1KQXBZsNR9P3bBHy5Uvc7RtdlhZkmFv+IS/u9k4NO8IChlcl2lBtDmzLqmcadr8hq5uyBCwf02E+nbVJtJ5VZbm0qSG/gT4EHcQoNPAU/H8pb0VkYaW3XOqvkxZ8Yk6r7+qW0qQ27gAPAo8Az//dPA72F5S3olyRwjs+Gr2eaNTu0f35mu2zNTzxTuFB1Q7gWeISJ7cUvOHgeuBW7zv98GXOe3rwVuV9WfquqjOC19yC8ivlBVv66OtR9PHGNoEesz4auSPdjc3uHGpz6c2LP7l6cWxVxyq+q/A38OfB94AvgfVf0y8Dx1i37xn8/1hxwAHkuc4qwvO+C3p8tnICLrInJSRE4u1hxDFuYt4pUjoOikB0bHHhId75nYDhtzyS0i+3Da+FLg14Bnichb8w5JKdOc8tnCJc9bUhaaw7qkZs4KpBLcW4t10xM9Lt8+RPJWHQs5zDWBImbJbwOPqup/qOrPgc8Bv4nlLQkOeXM2csRH/EVpZsnsvoKACOur97oyZfTmh2TEYNAoMKB7Bc5T8kyc9r0N+GNgk8kB5Qf89suYHFA+wnhAeQK4gvGAcs0GlP0Q3d5RVT+4jMJbU1nFFXgU+DbOFfg3OOJa3pJARePPJt4xSfqr/bqWNO5YbImhFFSBr4UTa6IWWzJ8tKUJREAOH0K3d4LNFVhqSrxNIYC/vD5JF29HUHShpPWN1MHMEkNTiF2HXU3wmFmyBOjqJacxqZP+8a5hmntoUG03g08gVTDNvQw42n3ch0gYCelNcxsagUKry9LSNPfeVq5sWDoIwHa3PnDT3IZBwGxuw1LByG0YLIzchsHCyD0gaM9y+TUNI/eQkJEgfllh3hLDINBXP/f/4hY99BnPAf6z60pURMhteEFaYR/IfbrvC4VF5KS1oX2YzW0YLIzchsGiD+Te6roCNcDa0AGC95YYDGXRB81tMJSCkdswWARLbhG5WkROi8gZEbmp6/okISLPF5GvishDIvKgiLzDl18sIsdF5GH/uS9xzM2+LadF5KpE+UERud//dqtIuwu0RGSPiJwSkbv62oZMdJ26ISOdwx5cVqoXAefh0rO9tOt6Jeq3H1j12xcA38HlJW8lZ3nNbXkX8EngLv+9d23IklA19yHgjKo+oqo/A27HZZoNAqr6hKre67efAh7CpWPuVc5yEbkEuAb4SKK4V23IQ6jkzsrxHRxE5IXACnAPDeYsbwgfBN4N/CJR1rc2ZCJUchfO5d0lROTZwGeBd6rqj/N2TSnTnPLGISJvBM6p6m7RQ1LKOm3DPIQaW5KV4zsYiMjTcMT+W1X9nC9+UkT2q+oTPchZ/irgTSKyhnuB14Ui8gn61YZ8dG30Zwxy9uLyel/KeED5sq7rlaif4GzLD06Vt5KzvIH2HGY8oOxlG1Lb1XUFcjp8DeeF+C7wnq7rM1W3V+P+er8J3OdljZZyljfQniS5e9mGNLHpd8NgEeqA0mCoDCO3YbAwchsGCyO3YbAwchsGCyO3YbAwchsGi/8HwAe7C/a3EasAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "segments_rgb = label2rgb(segments, bg_label=0)\n",
    "plt.imshow(segments_rgb, interpolation='none')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We save both the `segments` and the `segments_rgb` to disk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = profile.copy()\n",
    "p['dtype'] = 'int32'\n",
    "with rasterio.open(site_dir/f'segments_{site_name}.tif', 'w', **p) as ds:\n",
    "    ds.write(segments.astype(np.int32), 1)\n",
    "\n",
    "p = profile.copy()\n",
    "p['dtype'] = 'float32'\n",
    "p['count'] = 3\n",
    "with rasterio.open(site_dir/f'segments_rgb_{site_name}.tif', 'w', **p) as ds:\n",
    "    ds.write(segments_rgb.transpose([2, 0, 1]).astype(np.float32))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Obtaining the Undirected Network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "rag neighbors: 100%|██████████| 29857/29857 [00:01<00:00, 22482.19it/s]\n",
      "adding edges: 100%|██████████| 29856/29856 [00:00<00:00, 180876.24it/s]\n"
     ]
    }
   ],
   "source": [
    "chanG_undirected = get_undirected_channel_network(segments, \n",
    "                                                  dist, \n",
    "                                                  profile, \n",
    "                                                  interface_adj_segments,\n",
    "                                                  connectivity=8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Obtaining the Directed Network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "chanG = direct_channel_network(chanG_undirected, \n",
    "                               # The keywords below are how the pruning occurs\n",
    "                               # Specifies pruning will be done\n",
    "                               remove_dangling=True, \n",
    "                               # Do not prune nodes within 1 km of interface\n",
    "                               interface_buffer_m=1_000,     \n",
    "                               # Remove edge groups with an degree 1 endpoint and size <=3\n",
    "                               group_min_size=5,\n",
    "                               # How many times to do this pruning\n",
    "                               dangling_iterations=1\n",
    "                               )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do a quick sanity check."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fc8726807f0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALwAAAD4CAYAAABWgIzRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd5xU1fXAv296252dnZ3tvcEuIFU6iCgWNNg1lthii0ZTTNToL4lJrFFjYkti7LGBFYgCAlJl6XV777Nldnpv7/fHrMIKoggIynw/n/184L53z73vzZk755ZzjiCKInHinChIjnUH4sT5LokrfJwTirjCxzmhiCt8nBOKuMLHOaGQHesOfJekpKSI+fn5x7obcY4S27Zts4iiaDrYPSeUwufn57N169Zj3Y04RwlBENq+7p64SRPnhCKu8HFOKOIKH+eEIq7wcU4o4gof54QirvBxTijiCh/nhCKu8Ccg0eiJeyT8hNp4OpFZuKOTcCjCZ039iEhISVRx3znlx7pb3znxEf4HhC8Q5k+LKmnocQIQCkeweYIAbG6xsaHZQqvVy6d1fdhcXoLhyLHs7jEhrvA/EMx2Hyuqe2jsddJmcfNZQx+XPb+BtfV9ADh8AfzBMPkpOohEWFFr4ePd5mPc6++euEnzA+H5dU1sbugjTa8hQS1jSVUvzX1u0hMV9Dm8NPe6cfmCPHH5WHa02VBIBWaUHvSc1Q+Srx3hBUEYJgjCzn3+nIIg/FIQhPn7lLUKgrBznzq/EwShURCEOkEQztynfLwgCHsGrz0lCIIwWK4clNcoCMImQRDy96lzjSAIDYN/1+xTXjB4b8NgXcWReinHO6Io0mXz4A2EvyibWWSkzxtiV5ed97d10tbrYHh6InqVArlcikwmoFIK1JmdaBQyQtEonIBz169VeFEU60RRHCOK4hhgPOAFPhBF8bJ9yt8D3gcQBKEc+DEwAjgLeE4QBOmguH8CNwElg39nDZb/FLCJolgMPAk8OigrGfgjMAmYCPxREATDYJ1HgSdFUSwBbIMyTghWVJq56oVN3PL6Vl5c20BTv5tZZekUpyaQplOyo8NOo8XL7OEmSjP0GDRKxubo8QTh2VWNVPW4aR4IcN1LG7C4/cf6cb5TDtWGPw1oEkXxi2OYg6P0pcBbg0XnAW+LohgQRbEFaAQmCoKQASSKolghxkIlvAacv0+dVwf//S5w2qDcM4HloihaRVG0AcuBswavzR68l8G6n8v6wTMyN4kLxmShV8pZvKOTW17dxIJNrfiCYUozE7lwXDYlaYlcN70IiUQAIFmnQUKUHncYkwpSlVBt9vLYkhpquuzH+Im+Ow5V4X/MXsX+nBlAryiKDYP/zwI69rneOViWNfjvL5cPqSOKYhhwAMaDyDIC9sF7vyxrCIIg3CQIwlZBELb29/d/w8c8TnDv7W+f00ed2cn7W9rI0Gu4Y84wnrpyHOeOycYbivLC2kZ67H4ae91sabIgiiKdVvcX9a+elo9MIiFVI+OjX53Kr88eTqpeSX23lb99UofNGzgWT/id840nrYM28jzgd1+6dDlDvwTCAaqLByn/NnUOJmtooSg+DzwPMGHChO+P1Wpth8U/J1R+Cf1F53HnO3sIREQ83iDLqnu5cnIBM4eZuH5GMZOKUui2+yhO0fKvtU1023wIAjy8pIa/XDCatEQ1cgFsviAqmYRknYoJ+ancOBOcbh/LqvuJRr4/r+ZwOJQR/mxguyiKvZ8XCIIgAy4E5u9zXyeQs8//s4HuwfLsA5QPqTMoUw9YDyLLAiQN3vtlWT8MVAmgzSCy+mFee/dDNHIpckGk2+rD7PRjSlACIJEIjMo2cObITIrS9Vw9NR+ZTEq2XonbG6ayzUavw8ejS6rQqRRIJVFqzHaK0xM4f0wOH1dZ0Chk6FQnxoLdoSj8l0dygNOBWlEU9zVVFgE/Hlx5KSA2Od0siqIZcAmCMHnQBr8aWLhPnc9XYC4GPh2085cBZwiCYBicrJ4BLBu8tmrwXgbrfi7rh4HGANN/QShtLLsdOuaNyWJERhKCXCBZIyMtQX7Aahl6DSadggFPAE8oxL9WN3DWk2upaLYxpTAJEQV//jC2OaWUSbjjtBJ+dcZwFDLpAeX90PhGCi8IggaYw+BKzD7sZ9OLolgFLACqgaXAbaIofr6l9zPgBWIT2SZgyWD5i4BREIRG4NfAPYOyrMBfgC2Df38eLAO4G/j1YB3joIwfFunlvJL1BxQpmUwsMFKYrmN0dhL+YIjnVjfhDYT4cqjElAQVp5el4QnD2FwjHc4gLn8Yk07GqSVpJKplnDchm5QEFRqljHNGZzGtxMTgCvERxd7noWp9Jy6b74jL/rYIJ1JsyQkTJojHvRN3awX07ITsaZBRxn0Lq5FLJNx//igAOgY8PLSkmh2tVlRyKb+bO4IzR2UMEbGiqpMnljdx3eRsNrc5yTaoOGtkBh/u6MLs8PPkj8d9sXpzNGna0cOnb9SiUMr4yZ+nIpEe3Y19QRC2iaI44WD3xI8WHG8MNEDdClj9F3revI0PNjcgle41N3KMWkpSdIQjYRQyKSr5UMWt63byl8XVtPS7KUlLRCoVMOpUbG61srTSTHlW4nei7ACNG5oJuqNEwlG8juNjFejEmKl8Hwj5QaqA8VdDxmh4+ypSnR0USkZjduZhcflJSVABoFJKyUrWMS5bj1YmYX29mY4BPxZPgCqzi2A4ilIKaoUcfziK3RMkEo0SFUXmjTng6u0Rx9/djbeyAYkkB0OGFl2y+jtp9+uIK/x3ScgPXivoM/eW7XgDPFbY/G8omA2n3g2WRiI+K7UU4FTmEgpFefijah66cBRKhZxAKEq71UsgGKG210MoLBIB0hLlJKgU/HLOMD7e1cOuLjtyQaCh343D46M0PRGDRvmdPGr/q6+R1rSb0MTryR+VEnv83j58DfVoTjoJWWLid9KPLxNX+O+K5ffD7vmgTIBTfgOjLo2VN64EUYBQABo/gda1cPHLcNGLLNolY377H+kuuY9napN5fEU9P5mUz+q6PkZnJ/Hkj8eyp8tOkkZOSYqOYEREo5Ahl8tYXtNLbbeDym47XVYvJ+cn89Tl41DJj9xqjCiKB5zs+rq7cb/zDoZIhNIR/aScdjoAltdexb9jJ4pHHj5mCh+34b8rdOmxEd7ngOUPga09Vp5SAiEfZIwFbTrMuBPkcqTDz+IXmQ1keOoZX/dPTi0xUtPl5IaXN7Kry4VKLsWgVTGzNJ2Tso2oVUr0WhVyeWwMUytkrGvop9PqxRUCjVKORnngpcxvi3XBAjr/eD8Rj2dIuX/rVqS5uUhLSnCtXoMYDuPcsgVvVTX+5mZ8tXVHtB+HQnyE/66YcgtMvhk+vANq/wfv3Qiz/w+0WSCpgowJ0LEBlBpIHwmAZsxFoNFCzSIud7+MbPiVLNgZAvzs7rAdsJloVGTJ7k46LF4s7gAGjZSIOwLikXf2kKamEly0mPY770SelgaAdsYMhIREFBkZ+HbuJFJVRduvfk3UYSdo7oFAAKnm2NnzcYX/LhEEmH0PqBPB2gZrHwN1Mnh6oPQuIAo734KUUkgrh4TU2CQ2KiBvXs4V6v8x5ye38dCyOnodflbXmJlVNnRJcsmebv62op6xWXoEIcJ5o7OZv72LgrSEI/44SaeeCl4v/qpqghYLfq+T/scfRKs2EDGbkWq1SAsKkCjkBN0e5MZkEi65hITJk494X74p8XX4Y8VAC/zv12CpA7UBxlwVG9l3vgF+F2SNg5m/jn1JAHa/B/osyJmI3RvkkSW1dNu9DM9M5KczikhLjI2aTX0umi1uPL4wb21t54qJOXywo4tZpWlcO73wqDyKKIp0ubv4c8Wf8W6s4NerElErE9FOm4oQCuHetBmJ14vx9ttJufSSo9IHiK/DHx/0VMGBBhVjAWRPgpRyMA2P2fQF0+GshyAhC7p3QNPavfefdBHkTQaJhCSdirvnDuf0EWlUddt5cV3zF7cVpSYwpzwDqzdIW5+Ll9Y3opBIKU3THbVHdAfdvLT7JcSoSHT8KHZOMSE3JIFlADESQaZWozv9NIyXXPz1wo4ycZPmaBKNwqcPQu5kmH7H/tfL5kJ/FYgScLZDbx2kD4fZd8OKP8aWKg05YNx/ZDZoVUwqMvH+DjOBcJRoNIpEsnf80iokeMNRWvp8nDMqmynFR8+dr85ah9ljZm7hXGbnzcY71YWh100kHEE1rJSoy4U8OfmotX8oxEf4o4lEAilF0F8LkQNMGjNGgSoZXD1gN0P7Z7FyrRGm/AKc3dCx5SvFr6jupbXPRbFJN0TZY01LyNCr+fmsIm6YWXxUzsp8TqO9EVvARqY2E71ST4YhG9Xw4WhHjkAqlx83yg5xhT/65EwDRzfsemf/a4IA466ApMzYBLV/n+U6XUps8tpT85WiL5uQy7CMRNY39fNmReuQa1Z3mOLUBG46bfh+X4YjSb21ns3mzVi8FtZ0rCEajR61to4EcYU/2gw7I2ajV38AzZ/tfz13Epx8CzjaYkr/OeokSEgDSy1sfBFc+3trpSQoefPGyRQk61jb0M9HO/c6h72+uZUtbdb96hxpVrWtwhv0cmHxhXijXuz+49tdMK7wRwgxHPM29Ia8Qy9IJLGd1YRU2PoitG7av7KjOebO17Zz7wRXECBvesy3q34pfPwr+PAX0PbZkEmwVCrljtNLSNcr+WiPmY5Bt7455amcNiztaDzqEHp9vRQaCun19dJoa0TK8X2uPq7wRwDPtu30PPIo7+54nXfqD2C66FJg6i9BpoStL4Gleej1cBBMpRDxwNb/xkygj38HFc9ANAJBFwy0Q90nBFY+Ct6BIdVlUiklJg31fR4qOx1AbL6skB+9jzcQDvBe7Xt83PIxMkHGnt49NNubkcmO73WQuMJ/CwZefoXex5/AU1kJQNDnJWg2o6xsYUPnBpx+5/6VTEUw6jIIeWHz8zGN/BxPP+jSILUMahbif/FcOhp2sTNhBu6pv4UJ10L5ufgzx7DIdxJL6lxDRM/f3MJza1uZVmxCK5fwn9WNrG/oxxM4OqH07F47d668k4c3PYwn7GFR0yKy9FlcMewKtArtUWnzSHF8fx2PQ/qefoaBZ58FiQTP7t1kPfk3JCYTEa+XYn0BMkkPm3s2c3r+6UMrdm6Fzf8CRyf0VsVWZYy5kDYS+htApYfT/w92vol8+3+x9w/wZ0cm6T4JxWmjyUyawpSzfsbb83dS1uLgzDF7lyFbBnyEgkEsTi//WGWnfcCDRiHB6vHTYXOTYzhya/AVnRUsblpMg70BhVSBGjXFhmJ+Mf4XFBuKj1g7R4u4wh8izoULUY0ahahQ4OhrI3DP3WizcyESISM5H3VgJ9t7t++v8FnjIX8GNH8KLWsg4IDgWGjbAI4uKLsAFBqYeAPSUZdS+tr5/Eayjc9SxmJxBdnVbuPtTa3YvEEaLRI6Btx8vKeHfk+Qm2cUYkpQ4vCGKEwVyE5SYXWH8ATCbKq3kDPp8BT+6a1P44v4iIgR9vTvIVWTymOzHuOTlk/o8nRx57g7yTZkf72g44C4wh8Cod5eQp2dKCdPQjd9Opubl+Gr+IyR29twpOvInzAB9ZZl6OQHUDBBgGm3gc8GfXWgNcG4y0GdCgtvhtY1scmoIIAqAWVyIVNNqUydNQKAtgEXT3xSz0CLBSEa5WevbSE7RcfoLD1pSRp+Nqvki6bC4TA/fXULXl+YitYBLpqY963X4UVRxO63IwoiGdoMut3dRMQIBfoCegO9ZCdkf2+UHeIKf0jITCZUJ0/A/dHHhMxmCoMOapJNPH+6HKc8zFsrb8MT9pClO4hXkVIP2lSwNcLuD+D8f0ByIfidMecQrTGm9Bc8B7K94TLzjAncPC2fR1x+kjRyJpekMrs0lVG5hv2akMlkJKnlpGiV9DkCvFHRwlVTv905GkEQ+P303wPQaG1kIDDAmfln8uzOZ4lEIpyWd9q3knusiE9aDwFBIsFw440Iw4YRaO8gWWeifEBJViSJ6lArlbZKwoSZVzLvq4Uk5YBnAMJh6NgIi38FyUWgS42t0LRWQDgwRNk/Z0SukfEFKTgCEcZkJx1Q2T8nKkJhmpYMg4odHXa8vtBhP78z6KRuoI7Xdr9Gt7ubaZnTGJM25rDlfpfEFf4QSZo5k4SRIxDNZhaO9KHV6rlBmMHyi5dzUclFGJVGHAHHVwtQ60EmAbkytrHUtjZ2wMzWAs2r4dOHYeHt0L37gIfObpxZiEYu419rG/nfzk5cB1DkQCCAVikjUaHg8pNzcQfCLNzdud99h0qZsYxQyE2jo56Z2TO5aPhFhy3zuyau8N+CzHvuwfH8H6gtlLNoqgzVZeeRrk3nhpNuQKVQ8eLuF2l1tO5fURRh478guQTKfgQKNWRNgUgA7J0xpScMvfXw+mXw38tg5YOw8QXoqQZAq5RzzeRcehwBHl1ay22vb+HxZbVDmqloHKDT7iU9Sc2o7CSyk9R8WtPHlpZDiK0ZjUDLBtj0H9j2OlQvQS1X8yd1MS+Ek7h4+LE/+fhtiNvw3wJBJmPK1EtwtiQyPzqfVkcrKZoUktXJTE2fyt93/p2mVU2cU3AON4y+YZ+KAoy4AHImgMcC6/4WOwNfOB12vBo7NSkoQOqDSBB6doFECr1vx7yhyi+Ezi2MD4e4Jm8C73cb8IWiqAdDdYTCERbt6GRFTR+5Bi3jcpKQyaRcO72AR5fU8M9VTYzKTv56v9ZoFDY8B1teg7ALZCoIeqFrG0VduyBr4lF8u0eX+Ah/GIwxjUGj0rC4efEXZRMyJ+AP+xmbOpbcxNz9K439MaQUx862Z06MrcnnTQdDYWzSqtbHXP/m3A8FMwARJt4cU8KN/wJ3H4rOTVxt+ycP524iXRnE4gxS12Xh7vd3s6y6h4wkFbfMKsKYGAvrkZOs4+7pSVyU0okq5IDFv4aqxfv37XN6KmHdk+Bqh0gI0oZDxkgYqAFNKmSPP7Iv8jskPsIfBmm6NJLkSSxqWsRY01jOLz2fAn0BWoUWX8THGQVnHFzAyVfD0lrY9RYkZoGwHTq3Q+lZMOEnUHgK/O83YGuGmXeB1wLpI6BxBZJVf6W8/j+clyXlxd6pXP+aGalE4PfnjmTOiKFuf5iryPnsAXJa10OTCYQIZJ984D5ZW2HTc5A6DLJOjrkYJmbGfG09AzGH86TvzzLkl4kr/GESjobJ1mbzafunKOVKxqaMRS/Xs828DVfARYLyIL6kCalQNAvqloHHFjszH3aBf3DSm5wLoy+J+blGZVA0M1ZeMBPMlWxs6mGDdBx3njGcpZXd7Gx3IBHE/cNnVC4ApxXkhtgZe2PRAVeBANj2asx8OfW+mKm1L1rjt35Pxwtxk+YwuWrkVfz79H+TrErmjao3uOPTO7AH7GRqM7n/s/vpdfceXMDYK2MeUSptLPJY0A9tW+l0dtDmaEM0lsTK1z8MjavA0gQ1i2D4OSxOuYG1rW5eWNtEKByloc/Fg/+rwe7+Ulg7rw2yyuGyF6D4FEgphOID/PpEo+DuA5Vhf2X/gRBX+MNkZMpIMvWZnFdyHo6Qg2A0iDPkZG7RXGqsNdyy8ha63QcJXf/5Duwpd8HIC4jKFbxm3cSfVv+Op7Y+xSJvG5z7OIS9sP5pWHgHLH8YPn2YX3fdydOyZ8ikF4snTKFRgz8Y4tY3t/HPT+vpc/rBbY2dztSYIHcijLoEImH44GdQ98nQvkgkMZv9e2yyfB1xhT9C6JV6+jx9XD3yamQSGa3OVk7JPYUeTw9PbX2KDV0bDi4gZzzMfQjJ+J9Qq5TjiwYJRoJs7t1MSJsGs34Hnt6Ys7chAxRqUrKHo4j6yLdtIBAK4faHSAjbqem2UdXQSM2bv4O3fgyuLmhaDrs/BDEM5fNAqYWdr0HVor19CHoh6OaHnN4vbsMfJmIwyMA779BQImWMcQwGmYFbTrqFiq4KPEEP7pAbmUTGgvoFOAIOzi48++AC82dzU+921GXXsFupYFXHqlgM+KwxkJQLxaXg7o2F7Jh9H4qXr6Wwr5ocoZ9Rum46EBDSykhJzyWtcnlMwU1jQCmF2kWxSWfJmTDvH/DuTbD1ZUgtB1MxbH0t9otTMOs7eXfHgrjCHyahvj7633uftM522m9R8IH8A+6bfB9EoMHRwPml53P58Mu5fsn1PLHlCbZ2b6U8pZxzi89FKdsnsGnrBmhcAR1byVck4dakka3Wcv/U+1FIFbFNK0N+7DhC0WzoqwSJhJwZPyHlf3ch8degLJxLitcJgVYY6AKNEUrOgLb1IEuCU+6EDU9Dw+rYIbaADQaaY2Wn3h3zpjKVx0yfHyhxhT9MpHo91NUhTdSQa5MjL5CTpkvj5nE34w15UUgUSCQStEotlgELa81rqbZVk6xJ5tTcU2NCWjdCxbOxyaVMRdTWzB+2PIzMUMQDMx+I3SMIUH4+9FaCuwdmxw50UTgLdfmPQGeECT8FlQ7WPwe1i2OeVJ1bQKKKnccPuKDsHGhYBuZtkDkBSs6FppWw+nEIOmNt/ICJK/xhIiiVZDz2V+SZWZSHVtPuaicSiSCVStHINUDsiG3VQBUyiYxpWdMY8A7wdu3bqKVqjN1VlNR8BPpcuOg/EHDj2f0Wnv71KEK+oeZ0/iQYNhfaN8Z2PwFkcjjrz0M7NfUWsNTEJqUD9bE1/oAL7K0w5sexpdCgDzTJsPVVsDdD22owlsQCRP2A+dap5wev3T6YXr5KEIS/DpblC4Lg2+f+f+0j6weXel6iUJA0dy7aMaMpNhTjDXtZ27YWq9fKwvqFPLn1SV7Y/QJ2v51AJECiIpHT8k5DKSh5tepVftv0Oh9rNTDzNzEHkIRUEqb9grnlVxESQ6zpWDO0wczxsR3Zpk8P0ikJzPg1ZE2I2f26nFhAp22vxs7sKBNiewBd26FxOWRNAqkK1MbYBtMPmK8d4UVRrAPGAAymkO8CPhAE4VRiGbRPEkUxIAjCPjEmaBpMSf9lPk89vxH4mFjq+SXsk3peEIQfE0srf9k+qecnEBvrtgmCsGgwK/fnqeffHvxS/XRQ/jFjVu4s3qx+gwc3P0ipoRRBEFDL1LS72jEqjRjkBppsTdw46kZGGUdh9VvZ1beL6WWXgnJovPRsbTZhMYxKrhraSNYo0GfETlP6HLGjCAfCWABXvh5z+NamxML2rfgjbHgGznwApLLYvCDghfxZ4O0FnzV2xkebcnRe0HHA4aSe/xnwiCiKAQBRFPsOVvFESD2vVWi5aVsS134SYUrmZO6YcAePz36c3538O0anjabL3YXNZ0MQBIqSizg582RuGHMDiV9S9m0923ir/i1yEnMYYRwxtBG5OjYZbVgK71wH9cu+ukOCANoUGvvrcAtRGHM1dG6Gba+AuQayxkJSHrSujm04+ayw+eUj/l6OJw7Vht83TWUpMEMQhAcBP/AbURQ/jwtXIAjCDsAJ/J8oius4hNTzgiAc0dTzxH5VyM09wGGuI0SwuxvLggWkbm4l1eojaa2L9J/F3O7GZo4lQACn34lBafjK6FzRaJQVrSt4r/E9DHIDN5x0A8nqA4SpG3M5pA6H6sWw7TUoPj12qvIADHgHeHTjnxjbsZNbNSWxCAnrnwTjElCqY7a92hg7mx8xxibFju6haXl+QHzjEX6f1POfB16RAQZgMvBbYMHgyGsGckVRHEss5+qbgiAkcgxTz4uiOEEUxQkm09ELKBru68Px2n+hrw/CYewfLqTt9jsw/+Mpuh59lKJPavlF9FR8ES/za+cfUMbqjtW8Xfs29dZ6pmRO+WpXwYyTYPw1MbtfFGHNYweOUAwkq5MZnTmV6tRhWEUfONrB2Qm2NujeGZsTTLop5mYo18R2Wju3HanXctxxOKnnO4H3xRibgSiQIopiQBTFAQBRFLcRS0Bcyg889by6vBztjOmg14NSCcEgod5ePHV1iB4P3h070H20gal1EmoGarB79w9Jt713O0qJkglpE6izf3VaGJvXxlbzVmo0esTMiTE73D3o3OEe6uQhCAJGpZF2IcJ2Q9HeC/aGmJNHxmhwdceOJ/vs4BuArrjCw/6p5z8kZkcjCEIpoAAsgiCYBie3CIJQSCz1fPMPPfW8oFCQetddaKcPHrrSasGUglwqRRBBM2UKETHKmH41gijw3M7ncAWHBlSSCBIkUgmZ2ky6Pd28W/fukOsev4f5NfN5aNNDvLDrBf6959+8qRLBvAOa1sU8pd7/GWx/Y0i9fEM+KcpktI4WUKbE1uUVybEzNR/dBRX/jgV7Sh2Gx9rCW66G/TJ8/1D4Rjb8Pqnnb96n+CXgJUEQKoEgcI0oiqIgCDOBPwuCEAYiwC37pIv/GfAKoCa2OrNv6vn/DqaRtxKbKyCKolUQhM9Tz8P+qeffFgThAWAHx0HqeVVWFln3/g7rsGH4tm0j7HAQQSDc3Y2vrxeJKCJv6WXOmXNY6P6M/+75L7eOv/WL+ucUnsO9a++lorsCtUxNo62R7b3b0Sv0OP1O2pxtmBJMZGoymVc0j1Vtq9ju6WSGKonciidjHlMBe2wpch/kgpyAuwOZ1BBLmpY5Ajp3wK53QSLCnD/GfGzPfJDlhkxe6PiIni1/55cn//Kohtk+FnwjhRdF0UtsorhvWRC46gD3vge89xVytgIjD1DuBw6YC0UUxZeIfbm+XN4MHHd74FKVCtOVV+CfOQPXypX4KquIAn7rAKLLDV4vxW9s5IzxJja4d/K+9n0uHH4hAKXJpeTqc6l31iMTZQSiAaosVQSjQaKRKGqZmouzL2Ze0TwkEglKmZLXq1+na/J15DasBUsjjLkbCmd80R9RFFlvXo9RnUFOymQo/1FsTT59VCzrd8ATm7QCKNSMKDsPdc9qVrauJFWTypUjrzwGb/HocULvtHrCET7ssTJBr2ZYwpFNCaPKyUF17bWEQyEcCxchLlmCv74BfD5sjj50ZhOlmVL+E/0PlQOVjEsdR4+nhyZbE+WGckanjsYZdKKT69BINaRp0yhNLmVM6pgvQuyFo2FERJTKJDaUz2FzdyJSXzNZ9e8zr2geMqmMDlcHLe0VjLG2kT7giq3CNK+D2z6DKbfC2r/Cm5fAVe+Bqc34LkwAACAASURBVIQSYwnziuaxsHEhS5uXUpZSxrj0cUf03RxLTmiFtwZCPN3WS5ZcwhSDnksyk8jVaKh2+aj3eBmv15GrUX29oIMgk8sxXnwR0qxMBp4LEOzrofm0QhYGKjh1Z5hL3lWz6dJGKgcqsfgsANx98t1My512ULmugIsdvTvodHXy/I7n6Qn2oJQoydPn0WxrpsHawI0n3ciq1lXInN1Mc3lgysXQvycWjfi9m0CXAV57bLL6/g1wzSJQ6bm8/HKaHE1UmCtY377+B6XwJ3QWv/5AkFsqW5AhIgH8UZBKBAwyCbZQBINSzjxTEj9KPzKubf7ubhzvvU/KddcSVSmpnzwZiUJJwhWXsWKmnlUdq7iq/CpOzz/9oLZzKBLiLxv+QoO1gbAQZphhGGa3GY1Uw7ySeYTCIVZ0rKDf148n5CEQ9vDg+LsZUzA75kZobYPNz8WUHSEWqrt+acwLavbdAOzu28196+9DKpHy6tmvold+xY7uccQ3yeJ3Qo/wJqWCv5Xnk6dW8kZnLy92DuDyh4kq5dyRbWKpzcWb3QO4IxEuy0g57NQxqsxMVLf/HAApUL5pM67NW9CUl1HurmN3/26kgvRrJ4oyiYx0XTqnK8aQv7oRyVmnoh1XzDM7n+F/Tf/j0tJLOTPvTFa1r6JQX0itvZa32pehNxRQkFQApkKY+xhEwyAdzM7ttceCQnWcDjnjKU8pJ12dzh7rHhbULuDG0Tce1rMfL5zwHk956tiZ9Cuz01g+cTi356UhSCSsc/u4f1gOw7UqlvU7eabFTORAickOA0EqJXHKZGR6PWPSxuAOuum2dVDx3nxqN6w7YB1Hfx/2XjNjW03waRX+Fcux/vwX3PXeT6m317OxZyN/3/53lrYupcnRRK2jlstKLsPhd7C9Z/s+jQsxZRfFWCoeSw30VceiGQNrO9bSYGvApDBRb62n03X4kcuOB07oEf7LSCUSrs5No87jZZ3Nxa/J5Pcl2TzdYmaDw4uuq5/rctKOylKdQqbgzKIzyVRl4NyyDUt7Cxklw9Dok2jY/Bktu3bgdzrp72hDjEZQJSRQNvUUMq68GkfVTrLlm1hv2YJRbuTk1JOZlTcLiShhfsN8Xql6hXA0zO7+3fuHx6tfBttfB4U25gnVVxUrttYDMCljEmafma3dW8ke9v33dY0r/AHIV6mp8froCwQo0Gq4ozATc00rb3ZaWNXZwyNjysjS7PVWCkejzDdbOTtFT7JSTq3Lx/IBB9dnm9DKvnnOo4tKL0IURRoTI9TuWM26t17F73RiNXeSMXwkUUFB3uhTSMrMZvj4MgyZsaMHCbl5/DY4h/plN9Jka6LIWMTEwehgSoWS/1b9l3p7PRaPhbXta5mZO3Nvo6Vnxs7G500huPE53ujfhnXTE+y07OTaUddyzchr+L+1/8eLlS9SnFzMSNN+q8rfK054k+ZAXJlj4sFhBaQqY0otCAJXZZuoDoRZH4AnGtrwhPYGMd1gdbCkz8ZWWyyhWIXDwXzzALscrgPKPxiCIJAxuQzdSWm07tqGY6AfVUIiAZcHl2UMPc2FdNfqGeiTsWVJC3vWteN1B9AoNJxbdC4lySUMNw7/Qt6o1FE8MOMBpmVOIygGabQ2frlBGHUBJKbzut7Ay1InFr+FM/LO4PyS8xEEAYlEQoe7g9Xtqw/9ZR5nxEf4A6CRSSnTqYeUJcqknKJX44tGWO/08ceGLu4szCBDpWSCXscb3VbqvD7OwMBZhiTWWFystbqZakw65PZ1+Smc+rOb0C42YjN3UTZtFhs/eBcxuoHknHlYO3189m4DepMKQRQwNzgom5rJZcMv46z8szBohobRVsgU3DbuNjrdnbgj7q9sV6pJRilVcd3I6yg1ln5RXqwvpjipmHMKzznkZzneiCv8NyRXo+aVk4pxR6P8vamDzU43lkCYDJUSjVyOSiKwuN/O9TkmMrQq5qUm81J3H2929nFeejLab5HdbuKP9trbzoF+GjdVkDPcidepJRqNTaADfj+2XhcC0NviIK0wCVmqnwRDbP8g4vEgCAIajYaixCJ29OzA6rOSrE4mFAkhIFDRVcHKtpVU9FTQ5+3jxT0v8uisR79o26g2kqHNOO6TDn8T4ibNIaCUSTEq5GhkSlLkMkp1ezelJup1+CMRflndSp8vwAyDjrEJWp5r72WT7atH1W9K+fRTiRLF2dtG0TgTYkSCRBDQpySSmp+E1xOkbksP/3tmJ2versJu8RDwhPBWVtJ57334KisREWmwN/DM9mdY27aWe9bcw73r7uXlqpf5rOszrH4rUaKIoogrsNccKzeWIxNkbOnbcpAefj+Ij/DfAnc4TKpSgVK6d0J6ZU4qvnCIh1v7uXZ3I+ekGRilUtCtVbPT4WKNzck2h5d7ijKYnnzomzgSiYTCMRPoqKpi3NxMJs4tIhgIU722mxEzM1Go5Lz3RAWuAWjb5cDavQNTTgInT0gk6vOy4tk/8t6UXlK0KZjdZv5T+R+i0Si+oI/ZebO5tPRS2pxtLG9bTpuzjRd2vsAdJ9+BVCKlwFCAVq6lxlJDNBo9qqnsjzZxhf8WjE3Usthio8buoixpb7DUS9IM/M9iZ7vbRW97CIlUQo5SyWcWF1GplHSljO12N1P0OqTSb756Y+lw0bSzF6clFbdLz44Vrcz6cTlee5C2KgvtNf3MuLSM5PREAt4IDmsQXZIMQYBgUhqKYcMwfrCVWxu85Jbmk5iagPVHZzKibAYCAgmqBOQSOeFImJzEHFx+Fxt6NvDAhge4YfQNZCVkIZPI2NG3g7eq3/peHyj7/n5VjyHTDDpUCCzsH5raRq/R8Jv0AOezhnxplBy5nEg4hCUU4i6jknmt1Zg2reWll15iyZIl9Pb2Iopfbxd31tvorLUjkaqJRiV01zmo3dRFoknN9ItLsPf62fPkWyiadhDst5HoMzP3glTOunE0mYXJaE+ZhcmYw2hpLunoUVS1MNqixqgxkqxJRi6J7bbKpDLmFs3lkrJLuKDoArrcXbxW+Ro2l41eby8KiYKV7St5ZdcruAOHb6YdC+Ij/LcgQ6MiUSbDFo7s9xNfohjgCs0uEnLn4FFlcX99J5ZwmF3r1xGwWZHL5aSlpWG321m1ahUq1QoMhgROOeXJr2wvb4SRkadkIZNJ8dh9bF/ZRv3mXqIREWOWFq1WRON2kiGzoA7sRtFcjXOlj4jLhaq4BMd77yFLMpD9zDOosrIQAwEExVdHNZFIJMzKn4Uz5OTjlo9Z2rqUDlcHD0x9gH/v+Tf/2vUvdvbvZK49mylpk1GNGok86dBXo44FcYX/lhRoFGxyeIYoezDooKPjTRJ0OYzIiJ0wPDs5SpfnIzyOAcaPnseIESMwmUwMDAzw4Ycf4nAa8XjUrFmzhmnTpiE7wGqOVq9k60ethEJhxp6ey4yLhrPy1UpadvfR16pAZ0qk8NobMGYlkhmN0v2Pf+Bas5qBBQtQ5ucTdTohFCLY0YHgD6DIzUH4BnZ4aVIpq6SrcEfciBGRECFuHXsrD254kJa23QhLd9PtW4Nm5EgMV/8EXXn5kXvBR4m4SfMt0UtliKLIu61ddHv9BMMR+vo+wR8wk56+dznx2nQNZ0nXUFa8h8rKLbS2tgKxDSaPx0OCbjRZWbNoaWn54tqXkaukJBiVOPq8rHytmm3LWkEm0NfmpnttNdLaCsKfLsTy4osEu7tJmjWLqNcLnZ0EOjtIvv46lMOGYfnb3+h59FHMjzxKwGL52mcsMhQRIYJSqmR8xngi0QijTKN4cMaDeLUC4UQdvoATv7mb/mefJdDefgTe7NElPsJ/S+ZlJFPl9fGO1c0nrgCCGCbJ2c3pohpBsnfTKhgcQK3Wk5Y6iVCwhIGBAQAUCgVarZZgMMioUaNobW3F5TrwzqwgCIyYnk3WMAPmegedDQP4PRESkhSkurtIrVyKda0VWV4evppaIhKBiAjk5yNJ0mP5578QjMmox08gLJUS7u6m94knSL35ZlT5+V/5jHKZHJ1ChzPk5Fcn/wqDKrahpVfpGZs+nk9O72PmO2Yy9VKS27uxr11H2lXH94Q2rvDfkgS5jN8UZlDp9NLsCWALh/nIcRLmsI7EllfIC5yHyTQJpdKIVKoAVOj1enp6egiHw+h0Os4991yWLl3K4sWLCQQCeDye/dPV7IMmQUkwEGbmpaV4nCEECRhSJ+OuPAPzE08geDwEBAF5OILhjDkkX3kltnfexfrOO4ht7bh6ekmcMR3NhRdgf/MtBhYsIOM3vznoMmOiPJFmRzMJir2rUSmaFB6a8RA2n43H0x5HawtwzSdypImJRDwepFrtkX7dR4y4SXMYJCsUzExJ4tq8NH5VlEVn1MAqxvOAYwqP1W/mk+pnAQX6xEn4/DvQ69uIRCJs3rwZgLS0NMaPH4/T6UQqlVJXV0d9ff1XtmdustFZZ2XRUzvZurSVDx7fzqrXa6jrUJPw87uRmkyEKioQfV5Srr8eeWIiqT+9npyn/oG6rAypWo1n126Cbe0knXcegcoqLK+8SjgcHtJOKBLC7DYTjUbJVGeyvn09N39yM1af9Yt7ZFIZJp2J0WmjGdCG8apl2BcvovPe+xh4++2j8r6PBPER/gjy++Isnu/owx4pYVUwjRaLGb/vGUYYcnG7O0lL7cNgGEF7ezsjR44kMTGRsrIyGhoaaGlpQS6XY7fvH68GIBoVsff4CPpD2Cwe1ElKfL4wvW0O2mvC9OUlMP2ue/Bv2oB7+Qq8O3aQMHkyANrSUsQbfkr/3/+BKJPh+mQZKddcg3LcWFwVFXRGB9gyIQFXwMUE0wTW96ynylKFRq4hGA0ik8losDVw76p7OTX/VKblTCM7IXZU+LS809jau5XaFCcjGnrw+v2EzWY0kyahLjj+IhHHFf4Ickm6kTyFgk8GbDR5VTQ4ElgUVpLpXI4Y9ROJBJg9ezaBQACNJhalVyqV4nA4CAQCSKVSRo0a9YW8+k2f0dfawqjZc+htg92fdqLQyxBECdYONzKJQNAfRp2goKN2gM1KGalFE+nX9tP5QSXJ4SyyhieTlKJFPXw4aNTITSYQRQZefRX95ZfjXrOWxk0tVOWXUWmpZEHtAqZmTWV48nDCYphJ6ZNI06Tx/J7naXe382b1m/S4e7hp9E2s7lxNjbWGmoEamlPcPCw3gESD6PEQbG6OK/wPHYVUwnSTnukmPfM7enjM7aNPKGK+Yjin614nKdSPVBpGr997tKC+vh5BEBg3bhxZWVmo1bEJb9OOTXz27lu4+vtQaNSk5E5HqZMz87IivPYAFR+0EImKnHH9CFRaGXtWdeGwBmjdZUPMGAnqBNzVNjprbUhkEjRaGZ6BJNKqN5J3/73Y3ngD69tvo0w1MXLm6UyacA7P7HyGVmsrP8+7mhyVCWVODiG7HUGh4N5J9zK/dj411hp29O3grvV3ERWjaKVa0rRppNZ4wedHd9HFRKtr8FZVoz/ttGP1UXwlcYU/Cph9QQaCEZIUCvJUCvpDIkslF3JF+AU6Ot6ksPAmAAKBABs3bkSr1TJnzpwvJo+hQIDa9euQyeSoDAZyykaSmKpHpZNRt7GXmZcOZ/cqM2IENn6wjUhwF8ZsAyVjRqPPKCHJpEShktPX4WTPmg666m1IlCHU6lwyrJ9gfuwx1OPGoc/PQ3/GGXT86U9s2ryM2ok+QgoJa168nxltGjJ//0fsr71GJBREIsINEybgOOsSagdq8Uf9mDQmTk47mTprHa/IX6avy4x0zRoikQjhgJ9AczPKwsJj+VHsR1zhDwNRFPlDdSsBQSRfo+LqLBPbrA6e6RwgFIkwSa/l6qwU+gNBXjEL/M01l3sk7zFgnUiyYTSbNm1CIpEwevToISsle1Ytw9LeRmpJKdaONmor1jPpvHRKJ6ZT9VknFYv3EA460CbL0CS6qduwjdbaEIbtWymdNJWxZ56LRKskPT8JpVqOx16LKIYouHg6STndeHfthL5eXFVVWN95F9paGRYVuTdzHr3nTUYVrEbuMtN1551IUk1oxk9AKggIchlFyUUUJRcNeQ9j0sdgakllReourtjkALmcYDBI1+//QP5/XzuuDpvFFf5b4AlFqBiw8XafnQ1WNzIR3DhZ3GunxxfELMJIAf5v3DDUMinD0GKLiFQrS/H4czCb30IqSaO7uxuj0UhJSckXsusq1rFz+TKUajWm7DzqN6yjr7UFgJwR5+Bzhtm1bCkBTz1yRSJphSaS0tLpb2smEg7SuG0T3U116BKSKJl8DnmjizjlsuGse6eW5no/w++9j0hTI7L0dDpuuhGcTmQTJ0FlJYp11Yy9+GZUV16Ao2QLA88+S8TrRYaIPCubpB+d+5XvZGrGVJ4ev42NiYWcVzwP69tvEWhpoe/V10i/7tqj/ZF8Y+IK/y14qauP51p7mZik5bnyXMr1WpaZbWx0e7gjP43dLh/TkxNQ7+PP+qP0ZH6UnkxPz0zMPR8QCnWhUChIS0sbIrt553ai4SCn33AnNrOZgrHj6aqpomZjBfb+XAzpKsSgB60+g4DHTl9zI/r0DJLSM4gEQ/g9bpBIsXV2YutPwNwIJ52WxcS5hWxf0c7q/1Yz7swCUo16dNOn46mrJ/vBB+h97DGCbe2EtbHJtH7iyahS/oBj5Up81TW4N25CqtOhP2POkP4GLRZsH3xAqUZJdkoGyrnTSB5xPvLhw+i89TZsCxZgmHs2yi8957Hi+Pmt+R7R4QtQpFHxt7I8ZqUmk6pU8pP8dJ4dWcTZaUbuLs5m2gHOvIuiiMtdiyCEEASBSCSC2+0ecr2/pYmcEWMIBwNkl42gYNwEdAYjfo8DV99aOqtfx9rTglwtRZ+WQVJOLh6HHUEi4Zxf3sW0i69ALpFgKiohozgDa5+HFa9U43IGOOmUbIK+KLs+7SAajaKcPJmopZ9wt5n0229HYTLheOllosEgAMrCQlJvvJGsP/+JqFqNd89uIl7vkGeyfvgh9g8+xLN8Fae+shv9ko1EvV4SysrQnX0WghjFNn/B0f1ADoG4wh8i3lCIDTYXJyWoMCoPLY9aZ+fbuN2VJCacRHLyWEKhEHK5/IvrIV+Q6bN/jMQFS/75FJ88/zS1a1cDIlJAwIVCoaBw3CQSk1MQBIG8spF47TZSsvNIyy9i2JQZFIwbj7Ovl6bti8koCtHf4WHjwkbUCQpKTs4gHIji7Pchl8uRyBUIGjWKrCzUc+YQaGqi/Z576H3xJYI9PUAsU6EqPR3Hp6vovPsewoPHI8J+P6GGRhJmzyb5vHmUG8dQuLIR14YKIqEQWbffjuHSS/E31OMc3Gw71sQV/hD5d1svOpmMM1MOzWspELDR07MERDV5eTfT1NSEp9eBZIcL145eel6txLKgDkWLnFxPEQViGbYeM2UzTiUYCBAOhUjOziUc8JOWl0dW8XCioSDtNVWodXp8LicBjxuFSkUkBEmpmUgkMlKyMzHmaAh6o2z5uAV7nxOiItEoCHIF0hQj7g0VAJh+dC4JEycSHRjA+vrr9Pzr34iiiFShQJaXhxAK4d2+ne6/PoZj/XrMjzxCxO1GWlhIsKYWwWyGzk76//53ep98kkgggPGii5CoNTgWLybs9x+Nj+SQiCv8IbDW4qDa6+dMYyIzUg7t/LcgSIhGvbjcOwmHA2zfvh11oo6MRBOuii7CHTairS5wB9AWp1F49mSSMzKwdraj0moRJBK0BgMBr4f80eMomjiFga52LE31KHQ6Bro62PXpcga6OrB0NtNRtQOJVEptxXYCngjjz8rB0eejcWcf6iQZSalqtKNPQj28DO+a1bgGR+DkK68g65FHSL3tVoLbttH5lwfwNTRguuxS0u76LYopUxBEEe+GDYh2B9qJJyP29uKtqUYikyE1JKGePg3PmrV0/f4PCCoViefNI9Lfj2/7jqPxsRwScYU/BPzRKMlyKbfmpSM5xOhjCoUemexi7LbJOJ1WZsyYwUVXXIzxzGIIREEmRz7SACiQ2KOkZheQO3Is5pYG+ttbSTClMdDZhi41lUgwyCfPP03Q6yMsClg72gj4/DRv2UBSWjrn/vJuzrjp5xgyMumuX4lS7cDvCaM3KXH2h7B0uhjojs0djNddizQ9Hcsrr+JYtw5RpUaZkUHyRReRdNGFhHt76X/+P/j7LYS7upC6XSCKqKZMIe23v0E7cyberVuRGJJRmEwYb7+d9N/+FlQqfDU1uNavJ2naNIREPc5VB8kt+x0RV/hD4IxUA4+WFaA8hGhi+2I263A4ylCpksjMzEStVqM1JaCblIGQpABXiMSzchHkUgKNdsbMOYuMgmLUCQmE3G7aK/eQnJaBKa+QSDgCMjmjTj0DlU5Lan4+Uy6+AjEq4FnbjcmaTl75aMSoD5XWhaXdTUSAhGQ5A10+Pn29hv52BzKNhvRf/gpJQgLmZ57nnfvXYOlyIwgCKddei/GWW0Aioe+RRxAyM0GXgHvbNgb+9iTmhx6m769/RW4yIZVJEbRaki64AKlUSupNN6IpL8e+cBFBmw0iEfxNzYScziP8qRwacYX/Dunp6aG7u5vdu3cPCcyqm5RB0oUlCCo5kS4vjqwQne2dOCq6/5+9946T8yzvvb9PmWd639kys72q92LLkruNjcE2xWCKMQFCSAgckpM3vCQ5gRBS4A34JDmHnHCAUGICscENbMvGtmzJlqwuS6tdbe+zs7M7vc9T3j9mtdKyAmywbMvZ7+ejj1Z3m9lb19xzP9d93dePbVe9m3f+v19g/Y034/BXUUinOPDAfRQyaUyyidnRAaqb2vBa66hWGph5pI/MyQi5yTQzpwdJz0QJdbmx+RUiAylCnR7WXhNCEAVeemaCXK6AqcpP6K+/iKWzlUI8x/RwfOG9udauoequDyL7vMzecw+loSHI5zE8Hqzr1+G8/nrc73kPgm5gW78eaf7Glvstb8H/sY9i5PNMfelvKIWnKB0+TN8nPkH6+PHXfO7P8l86P/xrTTqd5ic/+QmxWIzW1lZ27tyJ31/JPW8YBtPfOoE55ERotDH+6Gl8mg3nrgbcV1QiE2OTEwweO0T3nqexxa1oVgGLzUFADFLjb8bu8aPP5LB0+fG+rZW5yCQP3/MtmjffTiknk4hmkU0Svjo7ikUiOpbC6jTTuNrPqh31xIYjHHh0ilyuRH1nFY2rfQTbPGiagagWyQ8OIvt8GGUVSyiIMO9hiv/0p6Qe3433dz+Ga/36Rb/zwHcfYvLxF2lZYWFwRGdcbcZhz9ESlFlx1w2vaoDZq5IfXhCELuB8YdFW4C8Nw/ifgiB8CvhDQAV+ZhjGn873+RwVKXgN+LRhGLvnyzdzTtTsUeC/zQuhmakoc28G5oD3GoYxMt/nbuAv5l/7S4ZhfHe+vAX4IeADjgJ3zetOvWFxOp1s2bKFffv2EQ6Hef7557nllluQJAktpyLZFNRYnqq3tKCP5cgPzlFOF1CzRWS7GV+oHldNLaMnjtPl3ITX8FUSzZtk5KATx7Ya1Lki9lVVSBYT1U3NBFpuY+x0mmCrlTW7QvTun+LEEz9GQKJ5y41Ex7NoqsHgkVnaN1ex8sogiekCE30xohMpDutgcyms3FFLoHUFxVyJglpCy2k43BWDLwwNI9ls2FesWPI7y+u2kEiGsLyjleBYGqF3irm+GQanYuT/+QEu+5vfR7Ral/S7WPxagzcM4wywAWBejnISeEAQhGuoSMavMwyjKAhC9XybVVRU+FYDQeDngiB0GoahAf9CRRX7ABWDv4mKkt9HgbhhGO2CINwJfBl4ryAIPuDzwBYqwsNHBEF4eF6G/svAPYZh/FAQhP8zP8a/vCqzchHp6OggEokgyzLhcJiJiQmampoQrRJywEI5nKU4msJzfTNSlYXCsVmSe8bxvbUNQRAoZtIoNit6Qq8sMy4TUsCKliogWk14LlssvlzT7GZ2bISJU89RztZTyiSx2DJomo/4VJ6urTW0balh+PgMJ56eQJJFbv3UBtbsDPHEd08Sj2SY7FGJRzJYHAqJcA6bR8FqlfEE7azeWk2prw/J50Mym5f8vt56DyIzvPDjYUIdbjp2rWDnh7byn194juRsgvgDD+B///tfo9l/5Xv464BBwzBGqUhQ/r1hGEUAwzBm5tvcBvxwXqB4GBgAtgmCUAe4DMPYP6+z+j3g9vP6fHf+5/uB6+a1XN8CPGkYRmzeyJ8Ebpqvu3a+LfN9z471hkZRFK699lo2bdqEYRiEw2GgkhrDeUUDhqpTHE8hmSQ8W0PY1vopz+aY/Wk/pVgem9uDyWwD1QCXjOx3YN8RAiSyfbElr7f26ka2334ZTRuuwBMMYnWtx1F1Ff6G67E7FRJzBfJ5lctu66RjWy1mm4mff/c0vQen6dxci1YGREhFckz2JWhZ5yPU4Wa8P8GJJyZ5+J+PMdoTQzzPaVUslCkVKtmV3QE7njoryZk8/ScidO8d5+SzE1irXNgtAoWBwddg1s/xSmNp7uScOHEnsEsQhL8BCsCfGIZxCAhRWcHPMjFfVp7/+RfLmf97HMAwDFUQhCQVmcyF8l/o4wcShmGoFxhrEYIgfJzKtwqNjY0XavK64HA4sFgsxOPnHhBFs4hokhD1c89Vju1BDN0gfTxMSgPfLW2IokDZokJJRQ2nOXXfCH6rC2e/QNofwdzqRjabEBUJu9vMml31rNn1XgDUssbcZJp8pozTZ+bJfzvNwQcGiaxNsvnGZppW+xjrjtG9b4ItNzZz8yfWY3PIHH1qlNmxNMMn5thwfSOdW2tQSwbxcJr+DR/Htr2OyFMjxKaylPIasiLRuMJH12V11LZ6Sc8V6bg8wNDRWfoOzVDf5aZZENHG4pSTSUzu10ZD6mUbvCAICnAr8Lnz+nqBy4CtwH/OK29fyEFt/IpyfoM+v2qsxYWG8Q3gG1B5aL1Qm9eLQCDA+Pg4yWQSt9tNtjeGaF4hpgAAIABJREFUoRuIdbaFNqIoIigS5HVKE2mmR4fJpmK4W2uRSk60sTQ9xSGyhRI32y7D2DsBPx9DtEhUf2QtktW06DVlk0RNc+XQLJfNo1hFZidGmOjN07LWT22LF5tTYao3ztRwgo5ttQiCwNV3rqJcVHnuR328+MgQoihw/e+sxmwN0vPCNLPjBRALyIqEN2glPp1n4FgEp9+C1SmTjhUwCSZS0SIWm8DG65vJPeAme6bwml76fiUr/M3AUcMwIvP/ngB+Mr89OSgIgg5UzZc3nNevHpiaL6+/QDnn9ZkQBEEG3FQUuSeAq3+hzx5gFvAIgiDPr/Lnj3XJ0N7ezvj4OMeOHaOlpYXRo71YkrCx5dzDX2LfGLlTUYiX0Y0y06fG0Moq62+9iezeMNpImis9G3iqeIzZbSLrQ83khhMUnp0ifF833stC2NqqwDAQ5HM72PDwAM/82zdIzeYoZuMU01Uce2KQy26/BZPZRlkbYLo/QCLSiLe2YpAms8zl72inKmhnvG+OvfedwWxVWHF5LQ63mapmFy5P5QE0OZfhoa+dYN99fZTVErmkiqbpbLqxkd4Xwuy77wxdhRKCxVzRmXqNeCUG/z7ObWcAHqSyj94jCEInoFAxxIeBHwiC8DUqD60dwEHDMDRBENKCIFwGvAh8CPjn+bEeBu4G9gPvBp6e997sBv5WEISzGf5vBD43X/fMfNsfzvd96BX+7q87NTU11NXV0d3dTXd3N8VMAbPNgnHMw/bt2ynnS+SORBHrbRhtTsTBHOXcOIrdgdlup2i3gFuh4aNraT1QYHR6jO1XXoZS6yBeUCmOJMjsmyb+00EknxXfdS2YG1wAVIUacPh9ZBNxbC4LxVKM0ZP7GTv5AorFgiBKGLqNWNiN9zydVptTYf31TQSaXcTCWfoPT3PksVHMdpmNb2nEtS1EPl3k2O5xXH4z8ViG3KyO3WPCXWWlqsGFbJY4+MggvSMF1gRc8BrmnX9ZBi8Igg24Afi984q/DXxbEIRTQAm4e3617xYE4T+B01T8CJ+c99BA5UH3O1Tcko/N/wH4FvB9QRAGqKzsdwIYhhETBOGvgbOJyb9oGMbZJ7PPAj8UBOFLwLH5MS45tm/fTjAY5Pjx4+QcOQqFAqlUqnJNLpJDEMG1PUh5JEk6NoOvqYHIwCDRsVGsXjsUSsz+tJ98Mk7GlCGdTuN2uwnc2kW5XKbcn6QwmaI4nSP+6BDOnSHsqwOYFDM3f+Iz6JrGsz/4Dqeefw50nWI2y9Uf+ii++kZOPvkYU/09tG1cLEys6zrekI1gu5eW9dWMn55ltCfCSz8fJzaRJTVXQFcNVu4KMnkmztDxKCabiGKvbK+aVlbR//wo0SE3hcwU4gW8OxeLl+WlMQwjZxiG3zCM5HllJcMwPmgYxhrDMDYZhvH0eXV/YxhGm2EYXYZhPHZe+eH59m2GYfzh/AcEwzAKhmHcYRhGu2EY2wzDGDqvz7fny9sNw/i388qH5tu2z/ct/raT8XpgsVjo6OjA6/WSyWRIJBJ0d3dz+vRpZI8ZQzNIvziJqc4BOjgyXtYVd1A+lcHe5gGLDIMZatzVxOIxjh2rBGip6RKF7jnUgoqlqwrfuzswfDKZ4xG0eQ+KyWzBbLPTuX0ntY1NKFY7ZquN2pY2Wtasx1FVzdjJk/QdHiYxkyE8lGB6NMnub53i/r87xOnnp7C7FFZcFmTD1c0ko0UGj8xgcypsfEsDXVuDeAJ2zDYZdDjwkz7OHnSuurYF1eRmKOtn5rvfWyi/2CzfeHqD0NDQwMDAAJIkoSgKDz/8MN6776745k/OEetPQdnAGM2gVNkRRsvkpzOAiNzpZusdXQz8xySz8zkjM0enybwwXglME8B+bSOO1TVk90+S7ZnFtbEOgNPPP8vAoQN4qutoWLGa/sMHGTh8kDMHnic6Nk45DwNHEhx9PEIhryKJoJZ1ijmV/kNhqlscuKtsDBydoZzXaF7r58r3di1kT2tYaefEnhfIp1dSLmlMDsSo7/BT2+yl+cpORveUmH7wfrw33oCpru6iz/Oywb9B6OzspLe3l0KhwJo1a3jwwQfZu3cvV/k3IqXsKA1O8mfmUOpdGFUK5UNREvf1IbkVfNc2IkkSoVCInp4eDhw4wJr1q6ndVouWLDJ7fz/ZZyYQAxb0okr+TBznhlp0XWPqzGlSc3Ncd/dHySSTDJ08zkvPPImmlrG73Nz8yf8HtehBR2D05DTDJ+OYzBJdl9WQiBR44f5+EEQMQ8fqNCHJIqWCinneO5RLjJONHENxqGSTlzFyYpb6Dj+iKHDZbR0kogUGfL/DqqrAr5mhV4dlg3+DIAgCLS0tHDt2jPHxcRoaGhBiJSSfiFrScGwJ4rqiAcEkkHh2lHJWhUYL1XevQ1Iq/41bt25F13V6e3uZnJzk6quvxl/rJ/jJjRRnciT3T6Afi1K25dFUDdkkk5ybo5BM8Mz3voVsUvD6qimm0xiGgmKzY+hlmtdVA5CK5pjoTdG1rYbVOxvIJfNM9MdBFymXVXxBO0OHIzz1vTLe6iSBRoWa9hVYXW7q1zQz2i0wevIkLu8Ia6+5AZMis/aqerqfFRg8OkPX9uBFn+dlg38D0dLSwujoKHNzc9glKx0pD+WTceQ6B7LdhOSsXCnUYyXQQFIsCPq5IwmHw8F1111Hd3c3Tz31FPl8Hqh8mCw1diy3dzEny+QH55j5XjcJYxq9WMRbF6R9++X4G5qprm9k3333MtHdjcPnw+nzLYzvDlhwVFno2FyH3W3G7jbjDTr5ydcOoZV1aho9uINOMrNjDB++n5qWFrSfPYKBjlG2YXbJzE0e5uDPkqy55noEoHl1gJ4DYcKDqWWD/6+G0+nkhhtuqBxEudwk//0MotuE5+bWBWMH8L+jk1j9BMWDEYqTKaxtXvSyhmiqxOnncjlkWV5I53c+vltayfY6STw6ghEv4W+rZ9Mtt+KpPrd/vuz2O9gdDaOWVDw158olxUQppxKfTaKpGoMnIvS9OEMuWeSGj6wi0OhGliWszk0MvxQiPNjHySd34w2G6NhxOeV9k+SiaQxNpZDJYHNVTleDrR6Gj82QjOVw+5a+51eTZYN/g2GxWLBYLJTyJUSbCcmlYHJaFrWRHAr+Xc1M98RIPjmMDuSPRnC/rRWTVaGzs5OeU6d58f49bA+tw311I5K98oERBAFbVxWJh4dxu3207LocW3UVAJqmIUkSqWiE6cFhMHTC/WeoaW1HlCQcbgWESfb9YBiLqwXdgFCXh9Z1NdSv8CKbzl2MaVm3CbWsYXHa6bp8F2qhn+kzP8Qf9NK2aTuK9ZxhN67yM3AkwuTpGO6dywb/XxLFqiB6FLRECbVYRjYvDhEQZRH/bV3EHukn/vQoxHIop92YNtfhdrvxuL1M9U2SS8+SH0oS+NBqFHflg5N8dBByZSga5AcTWNv8FKczzP18gPF4DwORw5X7tJEwT37r67zrc1/A7vbiq3NQHczSM7KXYOcuuq64nvq2WoTzIsc0TWPizGkmTp0kOjLI3Pg4A08dwBKTyRUivP2P/pjalsXp97w1dgJNLkZ65li1s56LybLBv4GxdHrJ7J8iczCMZ9fSwDfRLFeiiqJ5MISKr36emlANk9OTaLsCyC8kSDwzihJyQFEjeypaibYE8vunyXb4SNzXRz6fYih8mKyc4M7Pf5YDD9yHWsxjcbgWxr3iPe9GN/JM9RylbcNqBLGy5TEMg+NPPs7UvhOk4rPY2wIEV6yiHCvRll+Jy+OjytmyxNjPUtvmomdfloFjEdo3XrykTcsG/wbG3uGnOJSkOJgk5Z/Gtap2UX1hKoWWKON+aysmrwVL8JxKR3t7O8899xyT6QgrW+vIHpyidGIG8lTOuc1AyAaxMsmfDUK6jNXrwO2qJzucJjU7g9luZ3qwj1wijtNf2fYoFivX3f1xSoU8iuXcxY3jTzzKyEtHaM2tpLylnVXvuAGLw0mbaQP5vmkkh43O5sU5Kc+neXUVE6fj9B2M4Kux4Tvvd3k1Wb7T+gbHuSOE7pD58aMP0/fwUdLdUbRSJVKjnC6AScAUcmBpWhxe63A42LVrFw6HA/euBtxva8f1znYIKSBJSE0uBF3A+74VOC4LITQ6QBdZv/lqbF4fBx68j+rmVtKxWcJDA0ve1/nGnozOMNb9Ep66BlZ+7hY23nkbNpcbNaOixYuYmwOYfDayByPEfzZIObk4P42qamSTJbbf0ooAHHtqDK2scTFYNvg3OCanmerbOgm1hMhGk+SPRJj7zx5iTwyh2M1gEsgdnEbXFgdgjY6OcuTIEVasWIFokXGur8HRVgWiBEEbvlvakWUT2WencG2pw9zphZKBMmfjyvd/iGQkwsDhg6iGzrHHf0oueWFlEoBEZBKxAC0bN6L4bUjzzxsmpwnJLCEYBlKdHYoq2SNTlCbOibflM0We/Y8zdD87gd1roWt7NaWsykvPjl6U+Vw2+EsAySRx47tuYd2HrsS2tQal3kVpJEVhOAW1FtS5PHq2chcmGo1y4sQJPB4P0WiU4eHhhXEEi4xokZARUfw27FtqKYXTRO/vRY8VQNCRvGb0TAlNU5noOYmhqkiyCVX75RGN0pzMmsx2XIZ/UbkgCAh+K6Vojuz+cGUr5bAgVlcennXd4PiT4+TTZYIrKjH6LetrMNsFho7NEh6ae5VnctngLykkk4RjZQDP1Y2Y2t0Y2RJmj51SJEfmaJhyusDzzz/PoUOHEEURs9lMuVxe6F9KFhByOoZNQhAEHGsDOK9uoHR0llLPHCgy+T0TjN67n0w8zjs++5eEOlaRTszyyNf+jkf/91c5/NMHlgR6+TtbMDltSEtvGGLr9CDJEuSLgIHsUDC7K65HVVV56dkxJvrm8FY7MAyD3hcnyac14lN5Tu0Lv+pzuGzwlyiSJFFOFlHcZqQaK5kDUxT6EtTV1ZFOp/nmN7+JYRhUVVUt9FE8Fiwr/ZBXKYxVAl+dG2uR25wgSWAWoATWai871t1KoL6JW//4s1z1/t+hZeNGpoeHOPLYIxQy57YkekEl9dNBMIuItsU+EEM3ENxmVFUHuwUyoJZVRKXirxcFkdYN1VisEs/9qJeDjw8wdDSK3WehcY2X8ECc+HSGV5NlL80lijnooHAmTrEnjhbJgUVGMEtsWrcJj8fDQw89hMViobq6elE/QxHRJtOUZ/NYGt2Ikoi5y4uamca5s57MwWnqwvWgQ3E4jb3LR+v6TZhMCgNPH2DrineQ+tEw+XYvoqhT6E+iDqbAIWNtP5dvM304TPJgGNEqYkylQYNIfozsZJbpR8ZZe8NNmC0Wrn7/Cp7+fjfjPQkEQaCmxc3Wt7YgIBCfziArr66JLhv8JYql3YvkNlNK5tCfngSLhGNNJeKwra2NlpYWEolERd7eWXHxCYJAoXsOfDasK87FyEhWCUpgb6/C2uQl/uIo5f1R4o8NYG3ZjKhIiLNwlfOdSGGJsruIIcSRzDKSx4y6xoOcUdEL896jeIHkkyOQUdFdEggixbLOydkXkepEtL0vkUrMctX77mZ2Kk0ho+Kps9KxpYaOrbULEjkXwzW5vKW5RBEEAaXajqMjgGVjAHRI91Ye8hKJxMJ2RpKkhT23qqoUfCDYRWSbglZSyfTPkXsphrnVjSAKzD3ST/nZKJSotJl/WLUkrUiSBDYJx9UhfHd04bqtncC7VuCod2IoIlKg4qqM7R4Ci4yyuRrRbQUMzApUueq48gMfJdi5ht69z5CenSUTL5OeKxBq81Dd5EL+DfN2vlyWV/g3AY41AdTRNPljEUw+M7IiE4lEqKmp4fHHH8dkMmEymcjn85Qmk5SsBp6Hxwkm3VRlLFgcNqzrAwgWqZIOQgCMin/8bNSAoapozRZeCvchPzHCprqrqWmsRDeWzXAk20vDYY2uzpUIgKXZge+2DgoDMRLPTaAPpale346vPkgiPEpN+woESaSu2U3bxgDDJ6Pk0yq73tuJ6VXexpzP8gr/JkC2KTh2hkA3SD09Tv9LZ8jlciSTSTRNQ9d1MpkMiqIQ2tZO67pOBocGeWbyRQ5yBs/bW3B0+BElEceVwUrCExmYzpEerHxrKCE7qq6TbTB4SRzjaPcJUvOZgFNelZ7CGPueeo7E08NoyRJKyIUoitg6qxDdJhBh9VXX4XD7cPqrKaRTTJzpxe4xE+zy4/RYEWXhoiv+La/wbxIsdU70XfVk900wfrAfs9nMpk2b2Lx5M7quUyqVUBQFURQZHh7G4XCQNtKM56M8+swT1NbWUi6XWde6GhwC5AxQoTSZgdUguhRsJgsbxHYG9GFOnTxFNBolFAqxYcMGTCYTiWyGyPgMVS4vpfEUbAtRzpVQRxJggCEbDJ04zBXv+SDP/cd3Of74w0z2TROdqEct6Ky7th5JvrgGv7zCv4mwNbpxXtuIvdqJxWJZyLQmiiIWi2Vh9dQ0jXg8jsVioVwuMzk5ycmTJzlw4ACHeo4it57ztlhrKwFpkt2CNpXCGRW4TbqctwV30tLcwvT0NPv378dpc+IXHXirPBi5Ipgra6k2mQVNABl6vrObh77692STKa666yMkpqc5uftB4uE8oklgZix70edo2eDfJOi6zpEjR/jPJx6g7JNwOp0cOXIEVVWXtG1ubmbjxo0UCgUCgQBvf/vb+chHPkJbWxvZbBbv1U3gNoFXQZ/Pxyy6TAg2C7LfRt3mFlxJmY2udpqbm+nu7iaRjLFOaEUsGmDAgBgmGo2iNLlw39oOdVZKhQJ2jwenz4vN4aR989WIzltpXO2ncZWf1MyywS/zMkkkEjz33HOMjo4yMzPD9PQ0J0+e5MiRI+R+QWpSlmVUVSWZTJLJZEilUsiyjMlkIpPJINfYCPz+RpRmB+lnx8meiFCaSGMKWHHsqMNzbRP71TN878n7OXbsGPl8HkGUCK1sojSVZlCd4fGDT/P9738fUZFwdAUgX0QHunZejcPrQ5JNbLzlPSjWGgxBw9ANeIUyQr8Jywb/JmFgYACv18sNN9zALbfcgsPhQNM0XnjhBSKRyJL227dv56677sJms5HNVlbWQqHAzMwM6XQas8uM963tmGqtpF+cIvnYCOigNLmRTTJpe5mZ7Bxzc3PYbDbqgnX0e6I86z9Dn2mCYDDITTfdBIAgCRAyY8OCyTjndqxudLPjna1MnErRd3Aah/fiJ2RaNvg3AaVSiYmJCfx+P9u3b6empobbbrsNu92OIAg8/vjjC2m5z+L1enE6K3v9mZkZjh8/zsjICMXiuXxWJqeZwJ2rkDwWSrNZsAqYLJVISMWqoCgKNpsNURSZnZ1l8Hgfdf5a3vO+O/ngBz/IivMEEuwddVglN5FDfTzyT1+pKIYD669qpG6Fi1CXjw3XN130uVo2+DcBZ/PZdHZ2LiRAOnny5EIWg7m5OQ4dOrSkn91ux+l0Mj09zZNPPkl1dTWqqtLb27vQRpREpFUuHlWO0J2vhOxmMhmmpqZwOBz4/X6uuuoqspksbtHG2vEazBlx4YNwFn0mj6chSJdpC/EjI4ycPAFAuaiRnimQmsujWH610zA2NcGT3/w68fDkbzxXywZ/iZPNZjlz5gx2u5329nYAxsfHGRgYYOXKlQwPD2O1Wlm7du2Svna7nWuuuYYtW7bg8XjYsGEDfr+fPXv2LIqyjMoZ0tYyvZlR8vk8DoeDz3zmM2zatAmr1UpTUxPBUJDVl6/HkA3yU4tj58vJAoUzSQSfCTduQGR2tBK2LEoCvgYrsklidvKXB4pl02n23X8vkdFhssnkL23361g2+EsYXdfZt28fmqaxadMmJElC13UOHjyI3W5n48aNbNiwAUmSOH78+MID5vm4XC4KhQLyvPqe3+/H5/Mt/BugvqEem81GMp1kel6OPpVKMTU1xdzcHE6nkzvuuIO2tV0oNQ7UqcUPyZJFBkVAqXVg/6OVmJocpGOVlIAmRaZxhQ8Bg1915nTyyUcpZjLc8LufpH7Fqt94zpYN/hKmr6+PeDxOe3s7DQ2VlPy5XI5MJkNbWxtms5n6+npaWlro6enh0UcfrWQp/gWvzdTUFLOzs+zevXuRB+csdrsdr9eLKIocPnyYsbExnnrqKSYnJ6mpqUFVVaxWKzaHDanKhhorUvqFa3xGXkeSJSRRZ3qkn/Tc7EJddq6M2WHCd94l9PMZOPIiQ0cP0bpxKzWNzb/VnC0b/CWCrutksxOUSjFyuRwnTpzg0KFDWCwWNm7cCEChECYef5orrmhny5aKeqMoiuzYsYNgMEi5XGb37t089NBDC1I7fX19RCIRSqUS5XKZgYEBEonEIoOHSi77puanyeUe4uGHH6anpwebzcb4+DhHjx5deI9atojslJGk8/bvAiDqqDN5XIFqWjduI5dMMDMxPl9fJDV9lD33/huj3SeY7Du9EPA20dfL0ScexVlby7rr3vJbz+NyaMElQLEYpb//H1C1NKVSjplINanUGlwuN9u3b69EMQIDg98gEvkhba2fw2y+YqG/z+dbyCrs9/uZmJjg/vvvx+l0ksvlMJlMC4dRp0+/gMFxNO3nRCLt1NRcha7r1NY+CcIkzoYAp7st2Gw2nE4nNTU1bNu2DQCjoENeRzcqyaLOkuueg9kiWn0lfPjauz7CD77wWQ4/8hOueNedhM88RHToGNODBkcf/yl2r4+GrpVoqkopl8VitbLmquswKb+923LZ4C8BZueeJZXuprr6djQ1SiH/PM3N2/B6ZXL5H5NOX4/T2UWhMAlIeLwbl4yxbt06MpkM1157LTMzM4yNjXHq1CkaGxu5/vrrkWUZSSpSLp8glZ4gFk+TSL5IPH4Aw1DRjSwe92bs9lo++MEPks/ncTqdSJK04BkSLRJatowpULnCZxgGubEE6ecnwGvBs6OSZMnhr8JdXUs6Os2+H/07kqBx4+9/GqNU5sWHf0z7tu30HdiPWsyz9tqbWHPNdbgDr06ummWDvwQwDBG3ew21NTdjtdaQzY6STH2fcrmWROI0hqHhdHZhGHlsttXYbYv92ZIkceONN6LrOqIo4vF4CAaDRKNRXC4XLpcLXdc5c+Yf0XWdNau/hNPZzuTkgwyP3IOmCdTVvY1S6TSF4jiCkMXr9S95n6XZPNpUBsyV7Uzu9Bzxh/rAqWB9ayMm//xd1nKRxHSY7FwUu7+ane/7EKuuuAqAlVddA8CuO+8mGQnjqX11E6wu7+EvAUrFWcqlFIIAgiDi8a7BbPZgd3TicLYTi50gFnsRk+zBYq7DZHJdcJzz/eKCIJBMJkmnE5RKeTKZYSIzD1EqzWC315PPTxEMvg0MGUFQ0LQSAiU0LUP/wFfQ9aWi54YIFHS0sQSxp0bI9s1CWoNUnvwP+sgOVp4bYtMRNLVMTWcXgaYmVu3Ydd57lBDFyrfGq23s8FtIzwMe4HeB6Hz5nxmG8aggCM1AD3BmvvyAYRifmB/rv7T0/G+K1VrHdORB8vlh7PYQjQ0fobbmVkAkHj9MJnOE8fFaVDWJYq5C0wpIkuVXjtnT04OiKHi8pzl85DuIogdJMmFQ5sRLn0KSLEiSF1m2UyzNUCwOU1v7DrK5KcLhH+H1XE4otFgL2uyz4npHG8XJDKXBBGqmBC4JclpFpdfQySYTPPK1v8Ppq6KmpYNCKoUoXdxbTufza1d4wzDOGIaxwTCMDVQMMgc8MF99z9k6wzAePa/b4Hnlnziv/Kz0fMf8n5vmyxek54F7qMjKc570/HZgG/D58xT9zkrPdwDx+THedMTiR5ib208+30d+/qRTlm3YbI3YbPUEg+/G5dpCdfVVJFMniEYfYmj4Xy84lq7rlMtlMpkM/f2VmHldP4OAC497JatW/iOh4Huorb0FRWkmkeihWJoGcmQyo7g9O6iteQ+lop+52OLkrpqmIYgCrq1BArd3UvORdVTduQLzCl/FynwmCrk0u//1f2G22XjX5/6KmuYW0rFZho4dvsizeI7fRnr+FbEsPf+bkUoeJTLzCIpSSyBw9aK6QmGaaPRpLOYQXu9Wtmy+F7O5i3jsMPl8dMlYhw4d4p577mF4eJjx8XGy2W7S6Tlcrp2sWPHn+P2baGh4Pw31d1Ll34qqTgFlFFM9ilLFwMBXOX36GL291zM4UCKdTlMonMuFcz6CLGIJurBvC2Fa6QcBBp7dj66rXPX+D2N3u7F5fKRjszz3H99FVV+bL+dXavDnS88D/KEgCC8JgvDt81ZegBZBEI4JgvCsIAhnN2ghXqb0PPCqSc9f6uQLUUCmve1PsFgW/4pT4YdIJY+hqiVi8WcBgRUr/hIEkcnJHy0ZS5ZlCoUC+/fvx2LJ09Y+SlvrWlpabl3SVlE8KKYQYGA2B1i18q/J58cYHn4MXZfo6+vjO9/5Dt/97nfp6+sjk1kaFlAKZ8junUAvlLF2VNHytsvZcccHaVyzjrmpCXqee4qa1jY23fR2ZFlZ0v9i8LIN/jzp+fvmi/4FaAM2AGHgq/PlYaDRMIyNwB9TESl28TpJzwuC8HFBEA4LgnA4Gl266r3RKZdnUUwtBAI3LalzOjpxey5HFEsMDX2dcPhBfN7NOByNRGYeJJPpX9S+urqajo4OANasncPn89Dc/HvYbEtzsvt826iu3gmUKZeT2Gz1FIuT1AWHsdvtKCaF8lwOb8bKlVt2cM01Fe9KcTJN4pkx4s+MEHtsgOJUBvumWjw3NlHb1UawrRLvk4pEyMTjrL3+ray79sZXedZ+Oa9khV8kPW8YRsQwDM0wDB34v1T22BiGUTQMY27+5yPAINDJy5Oe5wLS8xeSsV+Qnr/AWIswDOMbhmFsMQxjSyDw2ijFvZqEgh9ANqn0938ZVV18IygQuIaNG76K07kGwzBQFA+CINPR/idgWBgZ+eai9g0NDbz3ve/ljjsux+Mp4XZvxuls/6WvXSxVLnALIuzddyOQIFB1HU6nE9kkI1n1FFveAAAbvklEQVRNtOdqqM445zOh5UnuHSdzMkJxJAmSiGNdNY5VASTb4j2/tz6EoWuUfyHM4WLzSgx+kfT8/J78LO8ATs2XBwRBkOZ/bqXycDpkGEYYSAuCcNn8HvxDnJOLPys9D+dJzwO7gRsFQfDOb5luBHbP152VnodLVHr+5eDzbaKt9VMkU0cYn/jhBduUSlHs9lU0NX0EQRAwDBlJklHMSz/goihSLE6SSZ8mlx2hrKYuOObE5AOk06dwODYiy14UxYzFvBqPZwNvf/vb6ejowOdwozhNiN7KdqQwlKI0koByCUM3cFwRxHVlPeIFLmanxmdITE+TCF9wnbpo/DbS818RBGEDla3EyHl1VwJfFARBBTTgE+fJxS9Lz79CBEGguvoGJiZ/iKrmL9jGbl9FLvc46fQIHs/KSvIlQDHVXrC9270Bl/ty4ol9pJI34vfvXFSvqkUS8RcRBFBMfgzDSrp0EJdzE3v2jKHrYURRxOfw4BqWSb8wia5pyB4zoteCns6jlVUUv30hj+RZtJKGaBKxTCu8fesncV/Z+arM08vlZRm8YRg5Kg+K55fd9Uva/hj48S+pOwysuUB5Abjjl/T5NvDtC5QPMb+NerOjaUWKxVlczgt/IfurriSeOMTA4Jex29qYjjyAogSYie7G6ezEam3EZjtn/LJsB0p4PJtwOJaG2oqiiXJ5Do9nJ5o2w8xMkFTyRhz264lEeqiqqmLlypWs7lqJHi2SeHyY1FPjuK5uxCjqkNChSiTx9Ciey+tRgpUoyOJ4mvTBKUq5IqRKOFcFcfmXntheTJZPWi8Bstk+VDVNVdWuC9bbrLU0N92NyeQjMvMIhqHhcKxEFBV6ej9HOLx4/Umle8jl+jGbuzCbfUvGK5fzCKKC1eKjq/Mv8PtuYGSkk+ee68HpdFIqlaitrcXmclDWNIjkMcJ5SlOpyve9zwyzBUqHosQeG2T2wTMURxLkx5IUhhJYq+04NgdxXnbxdVl/keVYmkuAYmkOXS+TyY7jdq9GEJauU5pWIpPuxWbtoLX19/D7dzE0/L8ZGXmJWHwfdblbsdkqz/9u13qczrWk00fo7U3Q0PBe7PbG88aqhAAYhoyi+AmFdHzzAsXpdBpd18nn8+h5ldxjo+Ay4dwRojCSwOS1YN4QgLxK+tA0UpWF0kSGjCagmzXMdU7s2+pQvBdXnvKXsbzCXwJ43Jvw+bYyOPj3lEoXVsVwONqxWKvx+3dRVXXl/MOrDpRIpQaIRp9eaCvLZmzWJlKpo0yFH6B/4CuUSufi300mHxiQSB5FVYsEAgE8Hg9Op5Pbb7/9XDZis4SRK0O6TG4ggZHWsG+uxb2uBvf2EKE/2ETgti7sawPke2cpvpTAkAUk22vjc78QywZ/CWAyufC4t6CqZYwLnjaAIJgACV0/57psbPgA69f9L6zWauRfCCgTBAFFqWHN6q9iGDp9ff9ANjs5/3o2DGSSicP09H6ebHacjo4OisUiuVyOu+66q3JhXBQwd/jBLqP1JjDV2bCvPifAcFa/tSKKLEBKRfIqCKbXz+yWDf4SITKzB5PJgSguvQRRLmc43fNFdF2lquqGhXJF8SObvICBIAjoeuVgWlVzzER/TrEwgyw7Cda9C1VLMDj0dQqFJOHwTykURnA4V1AuRRkb/zatrQ1omsYTT+xm7949hMNhSvE86lwOpd6NvDmAY8vSPXlpLkviaBhsErTYyHfPkjk5c9Hm6dexbPCXAKqaIZXqo1xOksmeWVRnGDovnfwUsbnnUNU0srx4Jde1HKVSkjNnvkwy+RIA5XIaw1AxyDM88nUs1gaqA28nn+vh6NEPMzT8dXK5GZqa7iYYvINsdownnvgsU1PHaGr6FsnUP/H8C/cTGQ2jZcqUeudQJAmlzr7otdOnIgz/6ARMZpCsMtY2HyRL6NHF911fS5YN/hJA1TK4nM0IgsDo6L8vOnHV9RI2azOablAoTFMqLT7I8Xg209jwIXQ9zlxsP5pWoFCIYJKrqGyBCoyPfwe/fxutrZ/B59tCXe078XjWMTZ2L4q5ilDonXi802zeMovT5cLpnEZTw4xlw7iuCEJARvCbF24+AWjZEuG9Qzw+9wK7LSepes9KtGQemu24rry48vK/imWDvwSwmGtZu/Ye3O7VxGLPMTp670KdJFloaPgAimJB1wXGJ+5F18/llBFFiWSqF0mqpVCY5OjRj3Gq+7Ok072EQu/FZHKTy41TLIYJBK5kxYo/p7X1YzQ3341isjM29m3MSjWBwCYCVT50fSPgQJI2kEwmEc0yosmCbcO5U91yLM/MI/2IkTySRabogoJJxburiarr2hEvouDBr2PZ4C8RzGYvLc2/h9lcj8m02KVnszWybu0/YpJdaFoJUTT9Ql83suQgmXwJl2stTY2/w4YN91BXdzuFwih2exsu18pFffy+Haxe/VUymRFmZp6iWJwjk0kj8DYmJz5ELJYgP5EksWcUk9uyIEYMkOuNofUlcLZV89abb8HlcrF3715EvwVL44VvY71WLPvhLyGKxRk0LYPN1rGoXBQVLJYgHk8XqXQ/3ae/gNUapK72bQiChWx2gGJpGo9nC52df4phGExPP8rs3H7yuVG6Ov/8gq/X1//3FIvDzETLlMtlpsNbKRbHmZqaQhRF3FYLSsCOdV1gkdqeGs6CV6HqnV3U2E1k9DynT5+mu7ubdevWXdQ5+nUsr/CXENXVN+JyraJ/4Ivk89OL6nRDo1CMUCzOks32MTb27co9V5Mdv+9yLJY6DCNPoRBmbPzfGRn9DtHoz7DZV5LLjc/77Bdjt7VgMrXjsHdit3morcvj8cisWbOGUCiEe2U1/nd34TjPFQmAQ4acTmGk4ttfs2YNZrOZ0dGLIyf/Slg2+EsISTLT2vJHlEoFCoXFBm+1VLOi63+wbev3sVob0LQy+fwYhlGmpeUTrF93D6qaIBz+GdnsKPn8ILJcg0m20Xvmf3L4yMeXvF5Dw/vZuuVfqaq6Cd0oIgiHWb9B45ZbbiGfzzM5OclcKr6kn2tnCBwSxYlKJKYoitTW1pJOp0mn00vav5YsG/wlRrk8iyjqGEZxSZ3LtQq7vbIiC4LM+MS9ZHNjCIKE1dpAlf8GZmd/js3awqqVf0djwzvJ5YeRZQvp9EH6B766aLy5uQP0nvlrpiP3YzJV4fXsJJk8yHN7f4/aunuprvkGA4NfXHRKC5A7HIFYCcl37szA6/UiCAJzcxc+KX6tWDb4SwxF8aJpJWKxpemvobKa1gXfhiy7sdvX4XRUwm8lyUxLy8cRJQcz0UeQJDtmcz3FYglFcVNbcweNDR9eNNZcbB/FYhybtYnqwI14vKvI5yeBZ3E6p7DZZsnnH2Pf8zsYGdm90E8tlCGrkT56Lie92+2ej8Vf+kF9LVk2+EsMl6sLj2cLM9En0bTyBdtYzDW43esxmx3M38UBQJIUujo/hyjA6Z7/QakUY+XKv6BcjmOxhjCbF4fqKoqXcnmKdOY4s7NP43FvIBR8D1BJASIIVkDEMAo89dRP6e7uJtcXQ0sUwSMj6BJ6SZtvKyxkTXg9WTb4SxC//wokycbMzO4L1kej+yiXIlitjej64ksjTmc7XV1fQtfLZLM9JJNHKJejjI99j3J58dYkWHc7nR2fp7Pjr1DVLFPhp0imjlLlv5bOzn/C5b4Sl2sLLtdVzM252LNnD6XJNMWBBLGgivfG5oULIJlMZj5+5/ULHINlg78kCQVvw2yuJjLzKMXi7JL6VOoEyVQ309M/o7//n5bUO53t2O3NFIsJNC0LpFC1PIax+G682RygpuY6MtmXKJZmiER+RDy+H4vrGv4qEuKP0u/GUX03ogibt+Txer2cMSb4Qfkpno4dQQ2c+3Y5e4G+tvbCt7BeK5YN/hJEFM2Egu8klx3m4MH3EIsfXVTv8W6hseGD+HzXkC+MkEicWFQvCAJu92p0I41i8gACLmcrJpNzyWtlMkN8f3KWE+pKhljBT3g3fxtr4LEUTOs+qqqvRVFqkKUT5PNT/Hz/HjRFoqqqCrO58tCqqioTExM4nU5cruWDp2V+A6qqdtLXZ6dU7iOXG8Hn3bRQ53S0YegFCoVZJqdOk0gcwuNZv6i/WQmQpZ90Zgpw4fVuWxQLc5Znp3p5thDCTpAJmljtsJMzLFjI8yH3NAHLBl6aESkUiszN5amqqsZisRAKhRBFkX379nHy5Ek8Hs8ikbPXi+UV/hJFECTa2j8N2Jme/tmCgABAdHYP3af/jKHh/4OqzmGxNCzpX1f3DhSzn4HsDPvsf0Bz8yeWtDmSSPF3c02YbWvIWtdya009O6uqUfU820zDXFf6v8zOHqNUrKe3p518Ps/KlSsXlL4ffPBB9u7dSyaTYcOGDbS2tl7MKXlZLBv8JUygajtmcyPJZDfF4rmDKE1V0XUDXS/R0vxpqquvW9LXYqki1PInPCP/N54tb+EH4RS6vvi09XvjUeKlMiHFxExBRRA1DsTmELMv8lnvacqlMP0Df0t7Rw2FQgtg5vjx40SjUcLhMIbRTVPTYd75zttZuXLlBb9BXmuWDf4SRpLM2GyNQIme3r9H10tks+OMT9wLaGhaienI/ZzLSLiY/SmZJwoyaQ3+fXKWmw+eZsfzJ/jeyDQTuTxr7RZaLQrvqPHRaTezfy7FVEFlB8ewSAJO5xpyueOEw1/H5yshCAKxWAxBELjsssu47bZPsnXrDtraOi74+q8HywZ/ibN61eepqrqaRGI/qdRpTCY3dnsrFks7kKZcLqLrSw97osUS/zYxhxf42opGVAN68yqoBn86PM1VL57hPyJxMobBk5E5hvIljpYMThbglHQVmdwEYEUQzBhGjlC9jK5XLnt/+MMfZvPmzZjNATo6/vtrPSW/kmWDv8Qxm2uorbkZXS8xO/scJpOTNau/TEfHp7Db2gkG343J5F7Up6TrfOBoH8cyeT7dVsuuGh+fbKymySzzgYYAaywmZGCkoDKVL/NQooBDhKtdFj5Q4yZvXcPfJG/ge3ErZdvliKIdu/1Rtm49xU03bae6uvr1mYyXwbLBvwnQ9QKKyUMmM8RU+McYRolU6hTFUoZM+tiS9qcSWU4WVDSgP1tANwwKhk6L1cydoQAPbl3B41s72eKy4hCgAIyVDebyBT5SX0WsqDEttfC0eCs/LN1AJclcDLtjlGDw9fWz/zqW3ZJvAgKB6ygUIiRTJ+nv//8ITz1DOnMQwxBxOJckeuMH4914SCFQjd8sIwoCPkUhXCrRk86xs8rDQFnDaZJpdZiZTBdpFsFpkrnlyCACsNFhJqKaGCvGSKglrIAsuzGZLqy1+kZheYV/EyDLNlpaPkZ722eQRBfJ1JPoegyrJUBtzS1L2v/36mn+in/gm977+Fig4pnZ5XPSoJj4+liUgqZRxMAw4MMNAdZYRNb5nLwrVI2DinrNkUyR/kIZQWnE2/wFCtYd9IqbORafXvJ6bySWDf5NhN3ezPbtP2JF11ewWNaCYCaZPEF0dg/l8rk49Lq627lu3ZcRiyP0D3wBVc3hkGUabQqj+QInE1l88xl/j8SzTBR0NjvsfKC+mj9rr0OiklFPB4bLBu8aCfAHhU/zF4X38ceDGvdOLA13eKOwbPBvIirBWT5CodtpanofGCX6+v6JM71fYWTkXxf52f3+jfz/7Z17cFTVHcc/v30n2SzZzYvNAwR8TcTqAMUHHVux4quFYn1A29HR2o4Vp1XHQVqnA63+UbGdOoxWW21n0NriY3SstKgU6HSs+ACUiEpIQgiEvElINptkN9k9/ePeyALJJlHIXbjnM7Ozd3/3nrO/vfnOzTnn3nO+Hk82kUg1iYQxUWN2XgC3UjxY3UBN7wCTfW72xwbIBd461MF7nd1cHS7gzrIC5vq9zM3xEU0aN7xiuMkV8KCo7JnYNd/Hgxb8aUoo9A1y/GeRnz+LyeGb6Il+Sk9P7ef7e3ubiESqCIdvwOs1OpqX5uVQ5nXRGh/gvUOHuTKUS20kSgTYGonz94Y2Ai4nPygvpCORpL2/n9TFsCe5nOR5PbTF4jT3W7f2TDp0p/U0JTurmJnnrQYSNDX9l7a2w8Tj7Rj+FMYDaC6Xn/b2/1ESvp4WFeTx+hZ8bg+lCFu6olT2DtA4oJic5aa+b4DYYJKVn+0j4HZzSzifR/ceaa9fXxBAlCKeTJBA2NLezdKy9NaZVqCv8KcxDocLh8NLbm4YSFBX9ww7K39KNFqPzxcikFtBX18VOz66gz/s2UpLLMHS0hBFbheC4qZwkDXnTeXyoB8P8ObhKP9sPcxLTW14BJRw5AqvFBcEsokkhWh8gD1Ra2c2jYQWvA3w+yuYUn4Hg4lDtLdvZk/1apLJASoqHmL2rLUkOIOG7jr643Gqe2LU9PSyuDDEzSUFXFUUZPn0UuYHc8gT+GogiywRVtU2E1VQ4nEyw2GM3MzLyyWWGASHsK8/Rp/Fs5uGQwveBvT3N9HV/QE+XzGh0OVEIrV0du5AKXA6/WS5W5kfmkSR38fr7Z00DCgODyY+f9gr2+3i3ulhKgLZbO7qIy6CwhBPEgh63dT1xWkeiFPg85DrEkSErZ3HW1lajRa8DeiJVtHc/AZuVxHhyd8mmUxSVbWSXZ+spGrPw/T3N3GJ4zPmxDcRkkF8QGcicVQd5wT8rJk5jWuLJqFwkO2ACLA/nmBH3wC7e2Pcvms/9X1xFhfnI0qxq2d4TyorGVXwInKOiHyU8uoWkXtEZJWIHEyJX5tS5uciUiMiVSJyVUp8toh8bO5bY7r5ISJeEXnBjL8nImeklLlVRKrN160p8WnmsdVmWWsnS2YwgdxzKSy8jEMdWwgGZzN71hMUFl5JX99uRJwEg5dSn/DycrSMDyNJctwOvhY6fmZSyOPm9xXTWHfhDO45s5QL/Fl4BdxiCKnM6+KbIT8LCoMUul00xoZ/StNKRI20wv5wBxtT4A8CFwG3AT1Kqd8ec0wFhr3lXKAE+DdwtlIqISLvAz8D3gX+BaxRSm0QkbuAryil7hSRJcBipdTNIhICtgFzMO51bAdmK6U6ReRF4BWl1DoReQrYqZR6Ml3+c+bMUdu2bRvz7z2diEb3U1m5jMLChZx55o9MC50aPO4CauPZrKlvxiXC/PwAfqeTeaEAfpdz9IqBeCLBS82dbO/q4dBAgr3RXupjCZZPnczd0yfu2RoR2a6UmpPumPEOS14B1Cql6tM8zL8IWKeMlYLqTCvKuSKyDwgopbaayT0LfAfDunIRsMos/zLwuHn1vwrYOGRVKSIbgatFZB0wH/ieWWatWT6t4O2MzzcZf24FNZ2f8mp1HQdjihJvgBsn+9nQfIjWvjiPnVvO1MD4n4XxOJ18v7SAG8MhNrR08E6XG19XDwdj/SSTSRyOzGk5jzeTJaSYEwN3i0iliPzFNA4GKAUOpBzTYMZKze1j40eVUcZshS4Mm8yR6soHDqsjMxtS6zqKU916/kThdHrIyirh9Z4Ab7c14kTYE+1nRfUB3u7u4dKCPKbk5oxeURo8DgeLwgU8cu4Ubi8tojE+yC927+ez7sy58zpmwZtt5IXAS2boSWAGcCHQBAyt0zbcpV+liX+RMunqOjp4ilvPn0i8+ddx0H0RC4rz+d3Maaw8q5SbSvJZWlLArSWhEzoFb0lpPtcXBelIJHhsXxPPHWg5bgqhFYynSXMNsEMp1QIw9A4gIk8D682PDUDqrOEyoNGMlw0TTy3TICIuYBKGI3cD8I1jyvwHaAfyRMRlXuVT69KMQMJVSq87jt+VB0B5jo/ynJNzN1REWBTOZ14whycOtLKpo5v63hg/LC8inG3dHdjxNGmWktKcEZFwyr7FwC5z+x/AEnPkZRrGvez3lVJNQERELjbb57cAr6WUGRqBuQHYrIze9JvAAhEJmk2mBcCb5r4t5rGYZYfq0oxAnsdF+2CSxAS2qQt8Pn45o4xrCvKo7Y+xem8jDdEM93gSkWzgSuCVlPBqc4ixErgcuBdAKfUJ8CLwKfAGsEwpNTSo+xPgGaAGqMXosAL8Gcg3O7j3ASvMujqAh4APzNevhzqwwAPAfWaZfLMOTRoCLicX5HhxjH1g7oTgcDi4ubSQu8oL6U0qfrP3ILu7o6MXPAmMa1jyVMfOw5IAu7p6eLSumYXFk/hu2Jr+zMeHe/jjgVYSwI3FQeYXBUctM1bGMiyZOeNFmpPOsw1t1PT2c4bPunt05+f5eWBGCTlO4fmmQ6xvmtj14rXgbUJXPM6mti4i8UEK3dY+FV6e7eNXZ5VT6HHzWlsn73Z0jV7oBKEFbxOcgDgdBL1OynOyRz3+ZJPjdrF8epgsgecb24nEJ+bJSi14m+D3eJgZyObroUBGLHkHEPJ6WFwUomswyeaO7gn5Ti14GxF0u/A5xvZ8zERxeXGI8iwv8eTEDJ7oKX42IuBy4cyQq3sqD59dPmH/dfQV3kYE3Q76hvFjtZqJbGJpwduI5v5Bkja67zIcWvA2oSsepykeZ4rPO/rBpzFa8DbhrbYuXMD5OVlWp2IpWvA2oa4vTsjj4qJ8a03FrEaP0tiE+2eU0DGQeXNMJxp9hbcJDhEKPG6r07AcLXiNrdCC19gKLXiNrdCC19gKLXiNrdCC19gKLXiNrdCC19gKW61aICJtQP2XrKYAYyGoU4FTKVf48vlOVUqlXY7BVoI/EYjIttGWgsgUTqVcYWLy1U0aja3QgtfYCi348fMnqxMYB6dSrjAB+eo2vMZW6Cu8xlZowWtsha0ELyJOEflQRNabny11Ihwl133m93wkItvMWEhENpr1bEyxGbI03xFyzcxzq5SyzQtj7fm/AevNz6uA+4c5rgLYCXiBaRhr2TvNfe8Dl2DY7mwArjHjdwFPmdtLgBfM7RCw13wPmtvBMeS6Dyg4JrYaWGFurwAeyYR8R8g1I8+tba7wIlIGXIdhyDAanzsRKqXqMAwc5pquJwGl1FZlnPEhJ8KhMmvN7ZeBK451IlRKdQIbgau/4M9I/Y61x3x3JuY70m+wLFfbCB54DFiO4ZaeilVOhKOhgLdEZLuI/NiMFSvDOgjzvShD8h0uV8jAc2sLwYvIt4BWpdT2Y3ZZ6UQ4GvOUUrMwzOSWichlaY61Ot/hcs3Ic2sLwQPzgIVimCOvA+aLyF+VUi1KqYRSKgk8jeEeDl/OiRA53olwuLrSopRqNN9bgVfN3FrMf/1DpnKtmZDvcLlm7Lm1uiM50S8MG8yhTms4JX4vRtsS4DyO7ljt5UjH6gPgYo50rK4148s4umP1YkrHqg6jUxU0t0Oj5JgD5KZsv4PRNn2Uozutq63ON02umXlurRagxYJ/DvgYqMSwzkz9Iz2IMYJQhTlaYMbnYFh01gKPc+RutQ/DtLkGY7RhekqZ2814DXDbGHKcbopiJ/AJ8KAZzwc2AdXme8jqfNPkmpHnVj9aoLEVdmnDazSAFrzGZmjBa2yFFrzGVmjBa2yFFrzGVmjBa2zF/wFXU5u7cqtSXQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_edges = export_edges_to_geodataframe(chanG, profile['crs'])\n",
    "df_edges.plot(column='segment_id', categorical=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Obtaining the Widths\n",
    "\n",
    "### Using the distance transform\n",
    "\n",
    "We are going to use the distance transform within a segment. We have a function that determines the width according to $2\\cdot d - 1$ where d is the distance to the land mask according to the distance transform in a small 1-pixel buffer around the segment (if this is 0, we use the full distance transform on the channel mask)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fc87223a630>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAAD8CAYAAABAfImTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2de5RU1ZX/P7uqu4EGREUwiOJbNCpBRSRRiYrdXZ2YZJI1+IOMjk7MUifqqKDiC0UNiqImjpmZ6FKDUUbH55oZf+nqbh+Mj1+UhxLUQeMrCkoEQd6P7q7avz/uvd3V1VV1b3Xdqrp163zWuqurTt269xTUrnPOPnt/t6gqBoMhHETK3QGDweAfxqANhhBhDNpgCBHGoA2GEGEM2mAIEcagDYYQUXKDFpGYiLwvIh+KyNWlvr/BEGaklPvQIhIF/gw0AKuBJcB0Vf3fknXCYAgxpR6hJwIfqurHqtoBPA78qMR9MBhCS02J7zcaWJXyfDVwQvpJInI+cD7A4MGDjzv88MNL0ztDRbJs2bKvVHVEIddoOnWwrt+QcL/Xil2tqhor5F7FpNQGLRna+sz5VfV+4H6ACRMm6NKlS4vdL0MFIyKfFnqN9RsSLG4d43pedNQHexV6r2JSaoNeDeyX8nxf4IsS98Fg6IMCSZLl7kbBlNqglwCHisiBwOfANOCnJe6DwdAHRelU9yl30CmpQatql4hcDLQCUeAhVX23lH0wGLJhRuh+oKp/AP5Q6vsaDLlQlEQIUolLbtAGQ1BJ9vXPVhzGoA0GLKdYwhi0wRAezAhtMIQEBTrNGtpgCAeKmim3wRAaFBKVb8/GoA0GcCLFKh9j0AYDAEIiY6pBZWEM2mDAcYoZgzYYQoG1D20M2mAIDUkzQhsM4cCM0AZDiFCERAhEcI1BGww2ZsptMIQERejQaLm7UTDGoA0GnMCSyp9yu34CEXlIRNaKyDspbXuKSLuIfGD/3SPltWtsEf33RaQppf04EXnbfu2fRaTy5zeGUJGwg0tyHW6IyFgRWZ5ybBaRy0plM15+khYA6bKlVwMvqOqhwAv2c0Tkm1g6YUfa7/lXW1wf4N+wpHkPtY/ASqEaqg9VIaER18P9Ovq+qo5X1fHAccB24FlKZDOuPVTVl4ENac0/Ah62Hz8M/E1K++OquktVPwE+BCaKyChgN1X9o1qlOn6f8h6DIRAkEdcjT6YAH6nqp5TIZvq7ht5bVdcAqOoaERlpt48GXk85b7Xd1mk/Tm/PSKrQ/pgx7lrJBkOhWE4xT+awl4ikCsXfb+vIZ2Ia8Jj9uKg24+C3UyybkL4ngf3uF9KE9v3pmsGQnTycYl+p6gS3k0SkDvghcI3bqVm6k5fNOPTXrfelPSXA/rvWbs8mpL/afpzebjAEhoSK65EHzcCbqvql/bwkNtNfg/4v4Bz78TnAf6a0TxORAbaY/qHAYnuqsUVEJtmeur9PeY/BUHacSDG3Iw+m0zPdhhLZjOuUW0QeA07BWjusBm4E5gFPiMh5wGfAVABVfVdEngD+F+gCLlLtLkfwj1ge80FAi30YDIEh6cGL7QURqccqmXxBSnNJbMbVoFV1epaXpmQ5fy4wN0P7UuAot/sZDOXASs7wx6BVdTswPK1tPSWwGRMpFjJiIy9Et26jdfsj5e5KRaEInSb00xA04mt/W+4uVCSqeAocCTqV/wmqiKZxsz2d1zzmsiL3JIy4B5X0I7Ck5BiDrjAaj5vjek7LZ78mNvLC4ncmRCj4EvpZboLfQ0M3rStuIbJ1h6dzzdQ7f3zetioLwe+hoReyYxcAjXV9Nx8aIlNL3Z3QoAhJdT+CjjHoCqAheqb1NzKVls9+DUBbx2N9zmtPPlnSfoUJS8a3xvUIOsagy0zqqNo8+pKM57QnnrD+phlsplG6cdBZPvaumnDPha4EEUFj0GUm1UhbPr83L2dWplG6bcejvvSr2lCsSDG3I+gEfw5RZRhnVvmohBHYjeD/5ISU2NhZuV/f+xcl6okBLMUSM0Ib+k3nPsP45rW/IroD3v7V5X1ej3/5r2XolXdiR12HRiLwySpatywod3cKxnKKmdBPQz9oHn0J6390IDXbQQWOnvErarYrb/12RkHXbYhMLZ2nu6OT1j/fUZp7lQSpiMARNyr/E1QQTUOsdNiWz+9l+b/OYMWvL0ejULNDiXTC8efc7flajQP+rk9bSbetQibaajnFKn8f2ozQJaB59CXokHpatz7c57V37rSm28efezfRjh6FmeZ9Lqbli99kvWbbroX+d9QjTfVnExmxV9nuXywqIRLMDS+63PuJyEsislJE3hWRS+12o82dg9jw87sft3x+Lx2jd895/pIFM4h2KKc2zLPek8OYy03r9kfQ7dt7tTVMurlMvfGHaooU6wJmquoRwCTgIltL2Ghz5yCxcWOv57v2qOWkH8/PGAzi8NrTV9A1KMrJP5rPST+eX+wuuuJEqGVi5zEHdj+OjbiA9tdvKEWXikqSiOsRdLzocq9R1Tftx1uAlVhyokabOwdOdJcTuVW3uYsBGzszBoOk8sp/XYlGQSPCKY3zsp7XNOxn/nU2C85nyEbj8XNoPuRK4uvuI3b09UXvTzFRhc5kxPUIOnmtoUXkAOAY4A1KpDNc6TiRW9EdXZ4dSa8+cyUAU069Nes5rZseKrxzeTD5jDt4+bmriB12FcndB7No8Y0ATJk8l1Ma57Ho7V+WtD9+Y025g2+wbnj+BCIyBHgauExVN+c6NUNbXjrDInK+iCwVkaXr1q3z2sVA8/yr11Oz0Vvqo8MLL10LWCNhOTmlcR7b9q7h1IZ5JHcfTPSrzTSNm03TuNm88PJ1LGq7uqz98wu/YrlFZHcReUpE3rN9T98ulc/J0wgtIrVYxrxQVZ+xm78UkVH26OyrznAlCu3Hxs5i/bf3ZsmC7HvJWtO/EaBtyZx+9sofarZ1sezBcBhtNpxtK5+4B4ir6t/agvv1wLVYPqd5InI1ls9pVprPaR/geRE5zFb+dHxOrwN/wPI55VT+9OLlFuBBYKWqpm6UGm3uFNadtDeRrty/Pa1v5e8JzpWskepJLybRbbtKcp/y4k/op4jsBkzGshlUtUNVN1Iin5OXIeNE4GzgtJQSmd/D0hluEJEPsDSI59kf4F3A0RmO01dn+AG70x8RIm3ump2KRuC4n3sPDvFCzmQNTfp6r0w0H3IlrctvKfp9goBHTbG9nOWgfaT/qh4ErAN+JyJvicgDIjKYNJ8TkOpzWpXyfse3NJpi1LZS1VfJvP4Fo83dzRsLZzLx7+9iwObiG5mDdnQW9fqTf3AHkYOGu58YAiwvt6dYbrfaVjXAscAlqvqGiNyDvaWbhUDUtjJkYPHvZ5KoFU44666s5xSSRdW8z8W9nsuAAf2+lhckAYlBlZ+w4AUfA0tWA6tV9Q37+VNYBh7o2lZVi5uU7hsLZyIJmDQ9s1GnZ1E5+8le1sMtX/yG2OE9P/bx9dmqmPpDpDOJJCvCJ+kLfsj4qupfgVUiMtZumoK1/CyJz8kYtAtNQ89l8vd7sopkV4fre15/bCa125NMPqN3NpITZBIbcQHNB80EevaTU40ztufP+1zTaYu/lz3YxG9qtnTw8n9fVbL7lROfkzMuARaKyApgPHArJfI5meQMF77+m6Op3a6c9JP5VsDHxi2u72kcdBZy+jjACp90Iq6cIJP4uvtyvj++4YG+jaO/QdP42chna9BEsiSBJdJZOn9AEPArsERVlwOZ1tmmtlW5WfzITE768XzqNnfRMOlmorsPdX1P245HaYhMJXHKsa7hk+k01k3vFR7qiAiWQ9FTOhPuJ4UEVaGrmiLFqpFvT7PWwa8+eyU1r77DpsMGs+GEkS7vsmhPPsmLL17juXyNg2PMDTXTuq9TDmNuPnAGrSuqY7vKoVqyraoX7XEIte1ayOLfWw6vE866i9hhHteWKf/CXhMqYnv+nPaux/Ppqa80RM+k5RN/99ODTlgEDoxB5+CP/3EFTcfe2KvtjYUz2bx/hE+njvJ0jeSK94gdea014nlc92ZcQ+fAGc39oGn8bKKHH+Lb9SqJMBi0WUO78d7HNE68CYA2O8Oo/q/WyB07/GpXr3PqdDl25LWwboPvUr1+juays7OknvSg4OxDVzpmhHahdfsjtC2+keiXPYIFyx6cweC1XWw5Mj8Znvi7txJf+1tiu5/X3dZYN737SCfXFN2v6pJ96mGFV0TGlTCUkzUjtEdaPv1Vr+cv//dVHHxH/9aZ8Y0Pdj/OJHjgeLqzTdGbx1zmyyjfNG42Xacdx2lTbiMZFeo27oL64kafBRVV6KoAAQM3jEHnwZFX/YpkLaycezlHPDuH5Oi6otwnm6qJM8V3Ctalk77llY3TptwGwItV5sV2w0y5q4jm0ZdQt0WJdMH+C27nwOEbOHLMGg6+824OeeIWPlrlzUlWCLItt0BC95ZXDi2wKd+9FY0IL75wja99q3SqSSTQgKXcGd2Jtb/RYf2zbemwpqedGweyXaN8ttqbUWeqENlUf7Z7H1bdYznWXMgVzBLd0clL7eEWK+gvquJ6BB1j0Hmw9Hcz0ChEd7PiuYfW9ST+b9MaNiWjfP3Fvrz+6QE5r5OpQmTr9kc89SH+7q29EjQy0TT03KyvJQbVerpPNRIGp5gx6Dz539sup7bWCon8cN1efHTFjF5ZqklVRkR28dCfT+zX9Z3RO3Xa7Hiinb/OtpJTTzo9mSNXrSnpqq74bK+ohmMf2osE0UARWSwif7KF9m+y26tWaP8bDwzgy4cPQERpnHgTn15wJRe+fRZDIwk6UWoFBkonv1iWf/F1Z/ROnTY7e9npIaC6YwdN9Wf3CkTpsw2VQuzo63n+tcqW2y0eQiIZcT2Cjpce7gJOU9VvYaWCxURkElUqtN80bjYD1u5g2YMzWPnjObQtvjFjvHZHiSoZOlN1p25WOeK+w0JVrKHVYqv9tNY+lCoV2pedu2hbNqdXW+uKW3hk3AK2JKPUImxKRu1KC9m/AG7rYC+kjsxSl31t7BS2i1e4dnYxqapYbhGJishyLNmUdltepWiiZ0HW5Y5nKaF69JjV1EuCPaL1rOranS869uDVx47tfr3xuDk0j7ms5zr2OtgPw4bc8d+R3YYw5bvZRfsNgFrraLcj6HgyaFVNqOp4LF2jiSKSK+m6YNEzVb1fVSeo6oQRI0Z46WLJySQZdPB+a1jdtZWNiXq2J+uo/zLJd5tv57Qpt9G2bE7GgBDZUXyJ3Pi6+0jWBX/9V27C4OXOK1JMVTeKyCKstW/RhPYrgfj6+2ne71KS6zf02nIas+8afmo/vtXDTlR6FcdCiR1xDfGVt/U83/084hsfZMCnG3y9T9hQ2ylW6Xjxco8Qkd3tx4OA04H3MEL7tKy6p8/+sbOV5IXmg69wlSPK5bXOhA7svZbujhv/agPNB1+R17WqDb+m3CLyF3s3Z7mILLXbSrIr5OUnaRTwki14tgRrDf0cRmg/Iy2f3+vpvKbxs2n56E7X8/LxWjcNOYcvTt2jV1tsxAXWg9o6T/erZnz2cp+qquNTNLxLsivkxcu9QlWPUdVxqnqUqt5st69X1Smqeqj9d0PKe+aq6sGqOlZVW1Lal9rXOFhVL7a93RVH05BzeqU7xoafb7WlhHTGRl2U9f2xERf4Xo2iadjPaN36MKNe3cykn95Fw3es6zszAEc+uPmQK329b1iwRuCiblsFphSOIY3WrQ8TGbZb9/P4+vtp3fpwr5DO+Jp/yfp+t2l2v0hYk6C2xTdS/2UH68YP7nNKKWpKVzIet63cSuGA5extE5FlKa+XpBSOMeh+4hhlathlJpGCTLiVh+3PKNq69WGa97sUgNr1lqNt8g96b7G1bnqI5JCBeV+7WvC4hv7K2YGxj0zVDk5U1WOBZuAiEZmc47a+lsIx+dAF4uz/xnY/L2MucmzkhX3ECNqWzCE2dhad3xiGJJXopp0A6MAaIts7iH84P+v9Ml0PrCm+MytwpvPfbb69z3mty2/Jeo1qRhGSPnm5VfUL++9aEXkWmEiJdoXMCO0TqSokvdrX/jbjnrUkktR9/jXR5R8g6zbQuuIW2hbfSPydPnrrfa6XsT3DFP9/WmYxZXLf62XrU7WjHg43RGSwiAx1HgONwDuUaFeoqkbopvGzS1oadcI/3E0kofD9sZz0k/kM+utONBohuq2D1hyjcDacUbh5/8t7SSI1j7ksq4rJCy9f11P1I4Vi18WqOGynmA/sDTxr7zDVAP+uqnERWQI8ISLnAZ8BU8HaFRIRZ1eoi767QguAQVg7QqYUTiqOMadOT5vqz0YTSdp2LfT9fkt/N6NPW0P0TFrzqKbREJlKdPiexNfdh27dBvTVN9P63Oviges6aDx+Dm1L5vRqbz5oJi0fZ6+UWXX4sOeiqh8D38rQvp4SlMKRoO8cTZgwQZcuXVqUazcN+1lJakQFgYZJNwPQ/voNZe6J/4jIMpeaza4MPHi07jfvH13P+/DM2QXfq5hU9Rq6ddNDeUdi+UUu3a9MZKpImc912l+/geimvppkTePzK9UTVhRIJsX1CDpVbdBQvvzhXLpfmX5ksmVT5VUML9rz3+1ssbUuv4XG4+Z4v0ZYUUDF/Qg4VW/Q+eIICfhB4wk3Z2zv74+M22gt23d2P+61xVZjvgZQRemThh5atz7cp80REMiXtjduIHZEj5xubOSFeU/F8yE5rG/0GEB0vXvN66rAj32rMlNVXu5iUYiHPDXVMXWPuSEyNe+R2nX6nS1ZJ2mEA6EyJIbcMCN0QJEa/39rI19uYPL3+yquaA75oqoiBCO0MeiA4qWkTb60fH6vFeiSRjVWm+yDgibF9Qg6xqCrjEXxWd2plamYtEqw8iHcjmBjDDpgFOpF9+JUa/9/s/uIE7Z8ON8YdTVNuW3lz7dE5Dn7edUK7XuhIXqm53TKVBwveno5m8ZBZ3mqf+V1Xzr+3jyajukdNdby4XzfVEgrkmoyaOBSYGXK86oR2u/PVlJ74omC1sHp5WzadjzaS7/MD7GC5IC+jrf4e/Nyqq2ElmoKLBGRfYHvY+mBOVSP0L7HiIJs4ZkOmapO9hepqSk4bLX99Rsyhn6mpmL6GUgTdKopsOTXwFVA6oZl1Qjte90P1o5OIPuInqnqZD6kXje+/n6ie+yR42xvSGeiT1vqfTIF0oSWpLgfAceLjO8ZwFpVXebxmlUhtJ8J58ufV3x1HjjXdQzOyWkuZOTPJPSf3v9c5WnDhKj7EXS8jNAnAj8Ukb8AjwOnicij2JIqANUotF9MYkfnrhCZbnCFjPxa6yGApRLmmoXixSFWAf8MXmR8r1HVfVX1ACxn14uqehZGaL8oNA09t99F5bx4wdOJv387jRNvynlOdUy7PTjEwuIUy4IR2qfvejmXh9iLg0n22ZumY2/s3lJy3uPFWNOreHhFOhPEjr6+X9tsoSIEI3S+ta0WAYvsxyWRVAk67YknaIie2TMNth1jmcg10jUNPRfZZ2/i7/codZ78o/m8Yr+nuw70sJ8h0QjxDQ9YXvWamoIVPFvfstI4m8bN7h6tI5u29+qLQ2zPn+esdFnRhCBHpaoliIpFbMQFfcT0m4ae22dvORPNoy/pLqdzasM8Il1JpEt5/pXrst9v7190V8aA7FK//SU2dlZG4w4KfkgQDRizn46adZnreZ9efIWne9mxF0uBz1X1DBHZE/gP4ADgL8CZqvq1fe41wHlAAvgnVW2124+jRyTwD8ClbtVmTOhnEUg35sa66RmNOdM+csvn93ZL7Na8+CYvvHQttZ9bVYaa97k48/1SjBl60jAd4f1seN7HTnOcNe9/ubf3VRg+e7nLEohlDLpIxPb+RffjbBFj2fa34+vvp3HQWbQnnqCxbjotH99FbMQFtHzxG8DdUB1aVt2T83Wv++vxd+b2kilKVx0NDT6tocsZiGUM2oX+qpGkj5r54mxFOT8GzqjfEJlK1xdrus8zgvklx0ttq5IGYqViDNqFyOBBvl+zkJDN9uSTvfah4+vvzzoVB2iomdbve6WTLiYYtuJ3HqfcOWtblSMQKxVj0C4Uw6Prt9Joyxe/ITbywu7nqdtP7V2P+3KPtmVzINr7OxYqTXPFr9DPsgZiGYMuAV5GstSRNHUEd6b8uTK+GmqmEV/72+49a7/VTpx7JwfUhDsTy4c1dLkDsYxBl4D0kSyTcaaOpKkjuCNA6Ox3p+KMxM57W7c+nHHNny1gxOvUvz3xBI0Tb+L5V67LWfe60ilyLHdJArGMQZeBbMkbmYzRMeKGyNS+MdwZRuK2XQt7Tb+znQf5Tf0jO3oHzKQ642IjLvB8nUDjc6SYqi5S1TPsx+tVdYqqHmr/3ZBy3lxVPVhVx6pqS0r7UlU9yn7tYrc9aDAGXRT66/TKJAfsGHE+xpdYv4Hm0Zf0qw+ZiB15bZ/48lBWrwxB6Kcx6CLQnnzSU9y2W+x06g9DPt7q9sQT3dFmxcSZCaQH0lQiXqbbYUmfNPQDLxlKbs6r6F7Dux/311tdzK0lP8NLA0E1CBwYik+2kdoPg8m1tZR+32ye9Pi7txbcj0rAjNAGXyiGqH6+9+2VMZZC7LCrcl4jVNtYZg1t8EJqXHcu0kfMUuYnZ5VNcimTk20bq+Jyq6tpDS0if7H1tJeLyFK7zehyeyXZV4gvE+kjZrlG7lRky3aaD7mS2FHXETvKSuH0Uk86CH3PmxCM0PkIHJyqql+lPHfSweaJyNX281lp6WD7AM+LyGH2ZrmTDvY6Vn5njBColrjRHy9wsYQG8yVTZlXbsjm98rbDgoRA4KCQKXf16HL7gJtmdylwttLynQ6n76s3HjeH5MjCJYQN/uPVoBVoE5FlKeliJUkHCwtOkkc5M5ScrbR8p8OpQS2xo65DEolu2aJQEYIpt1eDPlFVjwWagYtEZHKOcwtOBwua0L6ftG56KP8R0ocUyP6U80ml8bg5NB1zA/F35tK63Kpe6QgX+lkRpGxUk1NMVb+w/64FngUmUsR0sEoV2vdCQ2Rq/g4jLWxxlykOPB9OP3kuHSMG9RmV+6syGliqYYQWkcEiMtR5DDQC72B0uftFe/JJYruf1/3cy2jt1RgzxZA3RM/MGAeeKREkm1Tw869cR2JA9q9KoSV+AkM1GDSwN/CqiPwJWAz8X1WNY3S5+01q0kyu0TrXVDtjCmYGw3XSLtONPVMiSK4R95X/vJKmcbNpOuaGXsEmFbffnAXB8nK7HUHHyPiGlMa66Vl/LBoiU5HauoxG7ZXYkdd2h4Q21Z9d1um3HzK+g0btpweeN8P1vJVzZxR8r2JiIsUCTn9HQE1m/6FuTz5Z8Lo8/u6t3TW4QrOWDsGUO6/KGYbS09+IK7fsrEIjuU4/8Zew28CCrhE4KsBg3TAjdADIpdpZahoiU10FGprGzSbSmeD5166ncdBZBReeDwph2LYyI3QAcAT0/aShZlq/cqizKaPEjr6e5MAatDZK+4pbuttD4+GGUIzQxqADQvMhV9Ly4XzfrueXfC9YKZTxP9/h2/UCifrjxRaRgcDLwAAs+3pKVW80ta2qDMeYY2NnFXSdQiPCMhKpkq+JP06xXcBpqvotYDwQE5FJmNpW1Ulqlcd09c50Mu1T5xMR1hCZ6mtljUrHjzW0Wmy1n9bah2JqWxncJIgKnVa3J5/0dI3OkUMLuk8qgY779jZCu9a2EpGoiCzHCoduV9U3MLWtDNnId2/aOT8I3ujAOtG8GLNl0DlrWwGoakJVx2PlK0wUkaNy3NnX2lbGKRZwYsPPR7u6eon95buH7Jyfb02tUxrnEd2Z4Osj/C/YFzQE/7elVHWjiCzCWvt+KSKjVHWNqW1VxcTX30/rpodoPmhm0e/lSAs1TrzJChwBXnj5Ot68zz0kMgz4sYYWkREisrv9eBBwOvAeJUpmMiN0hdDy8V1Fvf4psdvZMn43Tv7hfF5ZfGNR7xVY/BmhRwEP257qCPCEqj4nIn8EnhCR84DPgKlgJTOJiJPM1EXfZKYFWNtWLZjaVuGiadxs93OG/ayXKorXbaxF8VkMXdVB/aot/e5fxePDtpWqrlDVY1R1nF2X6ma73dS2MvSmdcUtNB98RffzTPnLrZse6rXezmcba1Hb1WhtNONr2ba3irLvXQ6qSbHEEBxaPrqT2KiLiO3585xZTk1DzqFx0Fl5a5g5VSabD76iW7YX6M7OSvewu/1gZBJSCCwhyLYyBl2BxNf8C0Qzj6SptO14FBHp3q7qVfwuw8gaO/p69C/W1mfLR3eidTXdKZKO4ebtYS8g57rUhEHgwKvQ/u4i8pSIvCciK0Xk20Zovzw4csDxdfflLEPjKHzGNz7YvV2Vum2VPrJOOfVWOkYMpnXLAhqPn2Nd482bSNb1z29aiVPxappy3wPEVfVw4FvASkoUm2rojSMHHBt+PnR1+Xbdjt1qefHFawDYNmYIpzbMY8rkuSR2q+vX9YJSKMAz3gNLAo0XkcDdgMnAgwCq2qGqGzFC+yWnaei53Y/j6+/3VJGjPzrgrz19BS+1X03ktT/xwkvXAgEP2fSLajBo4CBgHfA7EXlLRB6w1T+LFpsaZl3uQmjdsiD/99ge79Qfg0wM2NDRpy11lA1syKZPOJFi1TDlrgGOBf5NVY8BtmFPr7NQcGxqmHW5i4GnGO1E7oJ50U07fepN5SJJdT2CjheDXg2stjNGAJ7CMvCiCe0b8sNLjLarkF80u3/ScXC5jfIVTbWsoVX1r8AqERlrN03BClMzQvslIja8T4ae7yQH1NAwKXO9Kmfq7WXKH4SMrv5SLVNugEuAhSKyAkuF4VaM0H7JiK/vk6GXF14qX7a/fkPWKLF8yDejK1BUwwgNoKrL7TXtOFX9G1X9ulSxqdVMbMQFvlwnvuGB7lKy2Wism07nkN57zumjbbZSOblwpumVoIxSTSO0oQz0p1B8NpxAk2zbT20dj/X5wqaPtvkK6jfVn909TfdTtLBoVMsIbSg92UbDfNaomYw32/bTKbHbM7YXQkVV1NAqCv00lJ7W7Y90V6lMrVaZzxo13XhzBYfUbtzFonhhiqOVTDXtQxvKRHzjg91/+7N+Tccx8PQMqCnfvbXbIVYVEXkpcuIAAAvmSURBVGHZUHU/Ao4x6Aohdfpa6NZQagZUw3duQaPC869YqZLRkf4E8vQn5LTcmBHaEHiyrsVrplkjdVK7kzIASPqzUEwVWagIqiWwxFBeMk2B81lHZ3NMOV7n5KCerarmMZfRsuqe7h8Bt5lA2KbnfjjFRGQ/EXnJTjN+V0QutdtLkm5sDDrgZPNK+xGR1bZrIZ31PQbd8tmvgZ4fAbcfjmx9q8TpNvjm5e4CZqrqEcAk4CI7pdiUwjFkx83YvDrRuuojnPST+X1K2mZUNEnxtuei4qbbYE+pC3eKqeoaVX3TfrwFSztgNKYUjqEQvOwBNw74O7Z9I8qWfWv6lLTNJFDgeN3DikenmGspnO7riRwAHAOYUjiGwnFb47btWogkFOlSxl90d4l6FWC8OcVcS+EAiMgQ4GngMlXdnOOuvpbCMQYdILzocDUN+5lnZ5QXUQKNiHVk8bdki8FO74OjROp1Wh40/AwsEZFaLGNeqKrP2M0lSTc2Bh0gvOhwJbdtdxUryIc//eZyJKlEEta3Nd3Z1t71eMYfmu4gFVvWN77mX4hveKByp+XqLm7gReDA9kQ/CKxU1dRpT0nSjY1BVxjtXY/3kdItNIos2gERW4Eo1dnmjM65fmjylfUNNP7sQ58InA2cJiLL7eN7lCjd2NS2CgGFJkFIElRgwj/czdLf9RSmq4gMKR/xIxJMVV8l8/oXLHGQTO+ZC8zN0L4UyFWKtg9mhDaw9KEZiIJW87dBgaS6HwHHi4zv2JSpw3IR2Swilxmh/eDQOOisgqO2EnWQ7J8Ed3iohtBPVX1fVcfbFemPA7YDz2KE9gND245HM3q081EJiXRCdJc17fZK6tq9krXEHKoxOWMK8JGqfooR2g8UGWO+81gDj3j6XSKd6jmJv6FmWq+1e0VridlUi4xvKtMAx61phPaLRK796GyjbsFC+LU1aFSQpHLC37kXlw+dw6zasq1EpA74IeD2U2yE9gsk1zZRsQwpvu4+liyYQd3WJHVbK0Brx2eswBJ1PYJOPiN0M/Cmqn5pPzdC+yHk1WevJLojweTv3+H5PaXQDS8JSQ9HwMnHoKfTM90GI7QfWnYOr2XrPjV865JfeTq/UN3woFA1I7SI1GNFtzyT0myE9isENw90ev7y64/NZNdwoasejj+3r9e7EjS286aa1tCqul1Vh6vqppQ2I7RfIWTzQDvbTun5y7EjrqH+r8qQ1Ulqtyc58W/vLHofy48/sdzlpppjg6qeTL+nzftfTnzlbezx9mb2iL+PJGHwZ9u6X2+Inhk+D7eDUf00BBW30jeQZasran0lIrs6ia+/n9eeugLpTPDtadZWlpeMsIpEjdC+IcA4pW+80LzfpYC1lm752DJc2bGr+/XNh+2Wc7rppRheRWBGaEMpKLYTqmXVPUDPWvqU2O1sO2Jk9+t//I8riO7K/mWOb3igqP0rGdXiFDOUl2KuWZsPmtmnbVF8FgPX7ujVVrexM+P7YyMvLEq/yoEkk65H0DEGXQEUc4Ru+fiujNtasrOr+/GUyXOJ2s/Tw1Lja39btL6VFKXqAksMZaLQEdpN0SR1W8uRFEodjeq+2EjnbnVM/sEdoXWKCe5BJaEJLDEED8fwvJBN0SR2xDV92to6HqNpyDnE35lL07E3WkXnd3WQGBCpiDVkQRinmKGUpE532zoec1UJdc1RzrImdDzkKkJ83X2s+cH+bB5Tw6YDavPrcKVhDNpQStKnu+2JJ3IatVuOsuzsyNjeUDON00+e21P3SkCjEu5vi49raBF5SETWisg7KW2mtpUhN4110wta07Z8mjn5or3rcZ5/5Tqim6296G37wPZvKNv2Df4IVQg+erkX0FeNx9S2MuQml4SuH5JAiWEDAUgOsI/aMBu0h+m2xym3qr4MbEhrNrWtDLnJJQyYabqdj5Dg5DPuoGNYLZN+ehcj3lT2eksZ05J5LzoUKF4N2nNtqzRMbSuDv7TteLSXUccOuyrruS8/dxU7htfQNVBYd6zwxsKZ1G7OvOYODd7W0J5qW+WBqW1lsOiPjljbjkd7ytf8ObcqyV7tHxPdpYxckuwObmmYdHP+Ha0QirwPHZzaViJyuV2N/h0ReUxEBhpd7solUl8PWMkYsbGzsp63rukgEgOEzvoI7V2P0/76DewYNahU3Sw9xd22CkZtKxEZDfwTMEFVjwKiWF45o8tdQaSGj2qHNXVu3fQQ8fdvJ3b41Rnfs/R3M1CBrhQbjnQpp025rah9LQuqkEi6Hx4QkceAPwJjRWS1iJxHwGpb1QCDRKQTqMca+q8BTrFffxhYBMwixWsHfCIijtfuL9heO/tDO147I0NUCrTny5geORZ/bx6xw68m/t68Pm8TtQT4HQau3UnbGzcUrZtlxafAEVXNFsZX/tpWqvo5cCfwGbAG2KSqbRhd7orCbb/aMerY0dd3tzXWTSdZY5XIcQrCt71xA03HhNigwx4pZq+NfwQcCOwDDBaRXPsfRpe7TOSK707dl86mZhJ/bx7xt39J7IhriI2dxQd3TmDncGHncGHzwSn/VWF0fSihKFbnZcp9OvCJqq4DEJFngO9ge+1UdY3R5Q4GuQJN2pNP0jT0XMBdzSS+smeNfOitd1OzXaj/IsJxP7+bZQ/MyPHOSkZ7LUsqFS9e7s+ASSJSb3vbpgArMbrcFUfrlgW0blmQ13s+uHYGK395OfVfJhn4dZKTfjKfD6cP49DbvBe1qwgU35xi5cTLGvoN4CngTeBt+z33Y3S5A0vz/pcD/mp9LX5kJq89fQXrj6hh4FfCwK/COO2u/DW0BF0ae8KECbp06dJyd8OQxkk/mU/N1gSL2jJveZUSEVmmqhMKucawur31OyP/j+t58c/vLfhexcTrtpWhQmg+cAYtnxR/OvzqM1cCcPrJc9GI8ML/XFv0exaXyhiB3TChnyGj67PV7if5yPOvXIcklcYTKjwkVLEEH9yOgGMMOmSUQ/Pr+VeuA1UaJ95U8nv7SgjW0MagQ0bzPheX5b5ti2+0jPr4OWW5f+H4F/pZTswaOmRoGb90bUvmANB4/JzuxxWDglbJPrShgtCdO8vdBdqWzKnMNXUIIsWMQYeM9NKw5SJZE+H0k/vkGwQbs4Y2GDLz/GvXE926y/3EoKAaCi+3WUMbikbrWxU27a6AEdgNY9AGAwCKJhLupwUcY9AGA/SkT1Y4xqANBgezbWUIErG9f1HuLlQsCmhSXY+gY0boMJGs/DVg2dBwCBwYgw4R8XX3lbsLFU0YnGKBz4cWkS3A++XuR4HsBXxV7k4USJA/w/6qWpD4nIjEsT6jG1+pamDlpyvBoJcGOaHcC+YzGEqFcYoZDCHCGLTBECIqwaALre4XBMxnMJSEwK+hDQaDdyphhDYYDB4xBm0whIjAGrSIxOz60h+KSPnFn1MQkf1E5CURWWnXzb7Ubq+4mtkiEhWRt0TkuUr9DIYUVDVwB1YN6o+Ag4A64E/AN8vdr5T+jQKOtR8PBf4MfBO4A7jabr8auN1+/E37MwzAKvr3ERC1X1sMfBurmF8L0FzizzID+HfgOft5xX0Gc/QcQR2hJwIfqurHqtoBPI5VATMQqOoaVX3TfrwFq9bXaKw+OpXgHsaqfw0pNbNV9ROsUkAT7SJ/u6nqH9WyjN+nvKfoiMi+wPexyhM5VNRnMPQmqAadrcZ04BCRA4BjgDcoYs3sIvFr4CogNSuh0j6DIYWgGrTnWtLlRESGAE8Dl6nq5lynZmjLq2a234jIGcBaVV3m9S0Z2sr6GQx9CWq2VbYa04FBRGqxjHmhqj5jN1dSzewTgR+KyPeAgcBuIvIolfUZDGkEdYReAhwqIgeKSB0wDavudCCwvbgPAitVNbUyXMXUzFbVa1R1X1U9AOvf90VVPauSPoMhA+X2ymU7gO9heY8/Aq4rd3/S+nYS1rRyBbDcPr4HDAdeAD6w/+6Z8p7r7M/yPileYGAC8I792m+wo/dK/HlOocfLXZGfwRzWYUI/DYYQEdQpt8Fg6AfGoA2GEGEM2mAIEcagDYYQYQzaYAgRxqANhhBhDNpgCBH/H1h7shCvBBlBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "width_features = get_width_features_from_segments(segments, profile)\n",
    "widths = get_array_from_features(segments, width_features)\n",
    "plt.imshow(widths, interpolation='none')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now add these attributes to the graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((540462.5, 7575260.0),\n",
       " {'label': 19373,\n",
       "  'meters_to_interface': 90847.54961447642,\n",
       "  'x': 540462.5,\n",
       "  'y': 7575260.0,\n",
       "  'interface_adj': False,\n",
       "  'graph_degree': 2,\n",
       "  'cc_id': 0,\n",
       "  'width_from_segment': 54.852813742385706})"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chanG = add_width_features_to_graph(chanG, width_features.ravel())\n",
    "node_data =dict(chanG.nodes(data=True))\n",
    "random.choice(list(node_data.items()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using $\\nabla \\varphi$\n",
    "\n",
    "We compute width using the gradient of the distance function (and possibly the network structure). First we obtain the orientation of along wich to obtain these measured widths (perpendicular to the artificial flow)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "flow vector computation using gradient: 100%|██████████| 29023/29023 [00:00<00:00, 63995.67it/s]\n",
      "flow vector computation using network: 100%|██████████| 29023/29023 [00:02<00:00, 11471.03it/s]\n"
     ]
    }
   ],
   "source": [
    "chanG = add_flow_attributes(chanG, dist, profile['transform'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To perform intersection, we must extract the polygons of the segments. This takes a bit of time to polygonize the segments, so be patient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 25.7 s, sys: 369 ms, total: 26.1 s\n",
      "Wall time: 26 s\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>label</th>\n",
       "      <th>geometry</th>\n",
       "      <th>perimeter</th>\n",
       "      <th>node</th>\n",
       "      <th>flow_vector_perp</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>POLYGON ((506490 7707810, 506490 7707780, 5064...</td>\n",
       "      <td>2040.0</td>\n",
       "      <td>(506200.54054054053, 7707611.3513513515)</td>\n",
       "      <td>(0.9724516, 0.23310483)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>POLYGON ((505080 7707690, 505080 7707660, 5050...</td>\n",
       "      <td>4260.0</td>\n",
       "      <td>(504491.3095238095, 7707259.523809524)</td>\n",
       "      <td>(0.998991568985452, 0.04489816361483868)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>POLYGON ((506520 7707660, 506520 7707630, 5064...</td>\n",
       "      <td>2160.0</td>\n",
       "      <td>(506263.36363636365, 7707473.181818182)</td>\n",
       "      <td>(0.9492476, 0.31452978)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>(POLYGON ((507450 7707630, 507420 7707630, 507...</td>\n",
       "      <td>7140.0</td>\n",
       "      <td>(508104.7368421053, 7706691.804511278)</td>\n",
       "      <td>(0.44568253, -0.8951911)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>5</td>\n",
       "      <td>(POLYGON ((505140 7707510, 505140 7707450, 505...</td>\n",
       "      <td>4020.0</td>\n",
       "      <td>(504499.1304347826, 7707085.507246377)</td>\n",
       "      <td>(0.9204698, 0.39081362)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   label                                           geometry  perimeter  \\\n",
       "0      1  POLYGON ((506490 7707810, 506490 7707780, 5064...     2040.0   \n",
       "1      2  POLYGON ((505080 7707690, 505080 7707660, 5050...     4260.0   \n",
       "2      3  POLYGON ((506520 7707660, 506520 7707630, 5064...     2160.0   \n",
       "3      4  (POLYGON ((507450 7707630, 507420 7707630, 507...     7140.0   \n",
       "4      5  (POLYGON ((505140 7707510, 505140 7707450, 505...     4020.0   \n",
       "\n",
       "                                       node  \\\n",
       "0  (506200.54054054053, 7707611.3513513515)   \n",
       "1    (504491.3095238095, 7707259.523809524)   \n",
       "2   (506263.36363636365, 7707473.181818182)   \n",
       "3    (508104.7368421053, 7706691.804511278)   \n",
       "4    (504499.1304347826, 7707085.507246377)   \n",
       "\n",
       "                           flow_vector_perp  \n",
       "0                   (0.9724516, 0.23310483)  \n",
       "1  (0.998991568985452, 0.04489816361483868)  \n",
       "2                   (0.9492476, 0.31452978)  \n",
       "3                  (0.44568253, -0.8951911)  \n",
       "4                   (0.9204698, 0.39081362)  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "df_segments = get_segment_df(segments, chanG, profile)\n",
    "df_segments.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will save the geometries to a file for later inspection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "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></th>\n",
       "      <th>label</th>\n",
       "      <th>perimeter</th>\n",
       "      <th>node</th>\n",
       "      <th>flow_vector_perp</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2040.0</td>\n",
       "      <td>(506200.54054054053, 7707611.3513513515)</td>\n",
       "      <td>(0.9724516, 0.23310483)</td>\n",
       "      <td>POLYGON ((506490 7707810, 506490 7707780, 5064...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>4260.0</td>\n",
       "      <td>(504491.3095238095, 7707259.523809524)</td>\n",
       "      <td>(0.998991568985452, 0.04489816361483868)</td>\n",
       "      <td>POLYGON ((505080 7707690, 505080 7707660, 5050...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>2160.0</td>\n",
       "      <td>(506263.36363636365, 7707473.181818182)</td>\n",
       "      <td>(0.9492476, 0.31452978)</td>\n",
       "      <td>POLYGON ((506520 7707660, 506520 7707630, 5064...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td rowspan=\"2\" valign=\"top\">3</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>7140.0</td>\n",
       "      <td>(508104.7368421053, 7706691.804511278)</td>\n",
       "      <td>(0.44568253, -0.8951911)</td>\n",
       "      <td>POLYGON ((507450 7707630, 507420 7707630, 5074...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>7140.0</td>\n",
       "      <td>(508104.7368421053, 7706691.804511278)</td>\n",
       "      <td>(0.44568253, -0.8951911)</td>\n",
       "      <td>POLYGON ((507480 7707570, 507510 7707570, 5075...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     label  perimeter                                      node  \\\n",
       "0 0      1     2040.0  (506200.54054054053, 7707611.3513513515)   \n",
       "1 0      2     4260.0    (504491.3095238095, 7707259.523809524)   \n",
       "2 0      3     2160.0   (506263.36363636365, 7707473.181818182)   \n",
       "3 0      4     7140.0    (508104.7368421053, 7706691.804511278)   \n",
       "  1      4     7140.0    (508104.7368421053, 7706691.804511278)   \n",
       "\n",
       "                             flow_vector_perp  \\\n",
       "0 0                   (0.9724516, 0.23310483)   \n",
       "1 0  (0.998991568985452, 0.04489816361483868)   \n",
       "2 0                   (0.9492476, 0.31452978)   \n",
       "3 0                  (0.44568253, -0.8951911)   \n",
       "  1                  (0.44568253, -0.8951911)   \n",
       "\n",
       "                                              geometry  \n",
       "0 0  POLYGON ((506490 7707810, 506490 7707780, 5064...  \n",
       "1 0  POLYGON ((505080 7707690, 505080 7707660, 5050...  \n",
       "2 0  POLYGON ((506520 7707660, 506520 7707630, 5064...  \n",
       "3 0  POLYGON ((507450 7707630, 507420 7707630, 5074...  \n",
       "  1  POLYGON ((507480 7707570, 507510 7707570, 5075...  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# From https://gist.github.com/mhweber/cf36bb4e09df9deee5eb54dc6be74d26\n",
    "\n",
    "df_segments_out = df_segments.explode()\n",
    "df_segments_out.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# From: https://gist.github.com/mhweber/cf36bb4e09df9deee5eb54dc6be74d26\n",
    "# To make multipolygons --> polygons after merging them.\n",
    "\n",
    "def explode(gdf):\n",
    "    \"\"\" \n",
    "    Explodes a geodataframe \n",
    "    \n",
    "    Will explode muti-part geometries into single geometries. Original index is\n",
    "    stored in column level_0 and zero-based count of geometries per multi-\n",
    "    geometry is stored in level_1\n",
    "    \n",
    "    Args:\n",
    "        gdf (gpd.GeoDataFrame) : input geodataframe with multi-geometries\n",
    "        \n",
    "    Returns:\n",
    "        gdf (gpd.GeoDataFrame) : exploded geodataframe with a new index \n",
    "                                 and two new 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": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cmarshak/anaconda3/lib/python3.7/site-packages/geopandas/io/file.py:108: FionaDeprecationWarning: Use fiona.Env() instead.\n",
      "  with fiona.drivers():\n"
     ]
    }
   ],
   "source": [
    "df_segments_out = split_tuple_pairs(df_segments)\n",
    "df_segments_out = explode(df_segments_out)\n",
    "df_segments_out.to_file(site_dir/'segments.geojson', \n",
    "                        driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we extract the widths are associated geometry. This is the most expensive computation in the notebook, so be very patient. Could take 5 minutes or so."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 24s, sys: 289 ms, total: 1min 24s\n",
      "Wall time: 1min 25s\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>label</th>\n",
       "      <th>node</th>\n",
       "      <th>geometry</th>\n",
       "      <th>width_m</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>(506200.54054054053, 7707611.3513513515)</td>\n",
       "      <td>LINESTRING (506520.00001 7707687.928470092, 50...</td>\n",
       "      <td>634.600687</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>(504491.3095238095, 7707259.523809524)</td>\n",
       "      <td>LINESTRING (504960.00001 7707280.588393838, 50...</td>\n",
       "      <td>450.454272</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>(506263.36363636365, 7707473.181818182)</td>\n",
       "      <td>LINESTRING (506580.00001 7707578.09814138, 505...</td>\n",
       "      <td>632.079571</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>(508104.7368421053, 7706691.804511278)</td>\n",
       "      <td>LINESTRING (508050.00001 7706801.748071255, 50...</td>\n",
       "      <td>154.438467</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>5</td>\n",
       "      <td>(504499.1304347826, 7707085.507246377)</td>\n",
       "      <td>LINESTRING (504960.00001 7707281.183522252, 50...</td>\n",
       "      <td>1108.129767</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   label                                      node  \\\n",
       "0      1  (506200.54054054053, 7707611.3513513515)   \n",
       "1      2    (504491.3095238095, 7707259.523809524)   \n",
       "2      3   (506263.36363636365, 7707473.181818182)   \n",
       "3      4    (508104.7368421053, 7706691.804511278)   \n",
       "4      5    (504499.1304347826, 7707085.507246377)   \n",
       "\n",
       "                                            geometry      width_m  \n",
       "0  LINESTRING (506520.00001 7707687.928470092, 50...   634.600687  \n",
       "1  LINESTRING (504960.00001 7707280.588393838, 50...   450.454272  \n",
       "2  LINESTRING (506580.00001 7707578.09814138, 505...   632.079571  \n",
       "3  LINESTRING (508050.00001 7706801.748071255, 50...   154.438467  \n",
       "4  LINESTRING (504960.00001 7707281.183522252, 50...  1108.129767  "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "df_geo_widths = get_geo_width_df(df_segments, chanG, radius=1)\n",
    "df_geo_widths.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's save the widths again for later inspection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_geo_widths_out = split_tuple_pairs(df_geo_widths)\n",
    "df_geo_widths_out = explode(df_geo_widths_out)\n",
    "df_geo_widths_out.to_file(site_dir/'width_geometries.geojson', \n",
    "                          driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We update the node attributes. And we see all of the attributes we have added."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((532127.3684210526, 7598432.98245614),\n",
       " {'label': 15138,\n",
       "  'meters_to_interface': 67502.63917727376,\n",
       "  'x': 532127.3684210526,\n",
       "  'y': 7598432.98245614,\n",
       "  'interface_adj': False,\n",
       "  'graph_degree': 2,\n",
       "  'cc_id': 0,\n",
       "  'width_from_segment': 937.4709297958259,\n",
       "  'flow_vector_perp_grad': (0.24091133, 0.97054714),\n",
       "  'flow_vector_network': (0.9900715437642892, 0.14056435617964214),\n",
       "  'flow_vector_perp_network': (-0.14056435617964214, 0.9900715437642892),\n",
       "  'flow_vector_perp': (0.24091133, 0.97054714),\n",
       "  'width_m': 1020.043221184196})"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chanG = update_graph_with_widths(chanG, df_geo_widths)\n",
    "\n",
    "node_data =dict(chanG.nodes(data=True))\n",
    "random.choice(list(node_data.items()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save Network\n",
    "\n",
    "We now save the nodes and edges to `geojson` files and serialize the NetworkX `DiGraph` with `pickle`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_nodes = export_nodes_to_geodataframe(chanG, profile['crs'])\n",
    "# If df_nodes is not writing, then drop columns with tuples\n",
    "# We found that flow_vector_perp may not be defined in some \n",
    "# of the segments near channels and `split tuple pairs` only \n",
    "# inspects the first row of each column\n",
    "df_nodes = split_tuple_pairs(df_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_edges = export_edges_to_geodataframe(chanG, profile['crs'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_edges.to_file(site_dir/f'{site_name}_edges.geojson', driver='GeoJSON')\n",
    "df_nodes['interface_adj'] = df_nodes.interface_adj.astype(int)\n",
    "df_nodes.to_file(site_dir/f'{site_name}_nodes.geojson', driver='GeoJSON')\n",
    "nx.write_gpickle(chanG, site_dir/f'{site_name}_network.pkl')"
   ]
  }
 ],
 "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
}