2_homogeneous_permeability_zone.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "placed-letters",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import libraries\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "formed-allen",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load sota basin shapefile\n",
    "sota = gpd.read_file('./data/input/sota_basin.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fantastic-immunology",
   "metadata": {},
   "outputs": [],
   "source": [
    "# change projection system of sota to EPSG:4326\n",
    "sota = sota.to_crs('EPSG:4326')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "numerous-plain",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load soil data\n",
    "soil = gpd.read_file('./data/soilmap_benin.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "hourly-happening",
   "metadata": {},
   "outputs": [],
   "source": [
    "# change projection system of soil to EPSG:4326\n",
    "soil = soil.to_crs('EPSG:4326')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "second-official",
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract legend information from avl file\n",
    "avl_data = [['Rivière', 0, 0],\n",
    "        ['Sols minéraux bruts', 1, 2],\n",
    "        ['Sols peu évolués lithiques', 3, 4],\n",
    "        ['Sols peu évolués modaux', 5, 5],\n",
    "        ['Sols peu évolués hydromorphes', 6, 9],\n",
    "        ['Vertisols hydromorphes', 10, 10],\n",
    "        ['Vertisols modaux', 11, 13],\n",
    "        ['Sols bruns eutrophes hydromorphes', 14, 14],\n",
    "        ['Sols bruns eutrophes ferruginisés', 15, 15],\n",
    "        ['Sols ferrugineux tropicaux peu lessivés, peu lessivés en argile, lessivés en sesquioxydes', 16, 27],\n",
    "        ['Sols ferrugineux tropicaux peu lessivés hydromorphes', 28, 29],\n",
    "        ['Sols ferrugineux tropicaux peu lessivés jeunes', 30, 30],\n",
    "        ['Sols ferrugineux tropicaux lessivés sans concrétions', 31, 44],\n",
    "        ['Sols ferrugineux tropicaux lessivés à concrétions', 45, 59],\n",
    "        ['Sols ferrugineux tropicaux lessivés indurés', 60, 72],\n",
    "        ['Sols ferrugineux tropicaux lessivés hydromorphes', 73, 79],\n",
    "        ['Sols ferrugineux tropicaux appauvris sans concrétions', 80, 80],\n",
    "        ['Sols ferrugineux tropicaux appauvris à concrétions', 81, 81],\n",
    "        ['Sols ferrugineux tropicaux appauvris peu ferruginisés', 82, 82],\n",
    "        ['Sols ferrallitiques faiblement désaturés appauvris modaux', 83, 87],\n",
    "        ['Sols ferrallitiques faiblement désaturés appauvris hydromorphes', 88, 88],\n",
    "        ['Sols ferralitiques faiblement désaturés rajeunis ou pénévolués avec érosion et remaniement', 89, 93],\n",
    "        ['Sols ferrallitiques moyennement désaturés typiques', 94, 95],\n",
    "        ['Sols ferrallitiques moyennement désaturés appauvris', 96, 96],\n",
    "        ['Sols ferrallitiques moyennement désaturés rajeunis ou pénévolués', 97, 97],\n",
    "        ['Sols hydromorphes moyennement organiques humiques à gley', 98, 99],\n",
    "        ['Sols hydromorphes minéraux ou peu humifères à gley de profondeur', 100, 102],\n",
    "        ['Sols hydromorphes  minéraux ou peu humifères à gley lessivés', 103, 104],\n",
    "        ['Sols hydromorphes minéraux ou peu humifères à pseudo-gley', 105, 107],\n",
    "        ['No Data', -1, -1]\n",
    "]\n",
    "# create avl dataframe\n",
    "avl = pd.DataFrame(avl_data, columns=['Label', 'MinNum', 'MaxNum'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "sound-grammar",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add legend data to soil data for further analysis\n",
    "df = pd.DataFrame(columns=[\"SOIL_ID\", \"label\"])\n",
    "for i, row in avl.iterrows():\n",
    "    for j in range(row[\"MinNum\"], row[\"MaxNum\"]+1):\n",
    "        df = df.append({\"SOIL_ID\": j, \"label\": row[\"Label\"]}, ignore_index=True)\n",
    "soil = soil.merge(df, on=\"SOIL_ID\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "organized-master",
   "metadata": {},
   "outputs": [],
   "source": [
    "# clip soil with sota\n",
    "soil_sota = gpd.clip(soil, sota)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "monetary-rings",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Rivière', 'Sols minéraux bruts', 'Sols peu évolués lithiques',\n",
       "       'Sols ferrugineux tropicaux peu lessivés, peu lessivés en argile, lessivés en sesquioxydes',\n",
       "       'Sols ferrugineux tropicaux peu lessivés hydromorphes',\n",
       "       'Sols ferrugineux tropicaux lessivés sans concrétions',\n",
       "       'Sols ferrugineux tropicaux lessivés à concrétions',\n",
       "       'Sols ferrugineux tropicaux lessivés indurés',\n",
       "       'Sols ferrugineux tropicaux lessivés hydromorphes',\n",
       "       'Sols ferralitiques faiblement désaturés rajeunis ou pénévolués avec érosion et remaniement',\n",
       "       'Sols hydromorphes minéraux ou peu humifères à gley de profondeur'],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# show unique values of soil classification in sota\n",
    "soil_sota['label'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "involved-extraction",
   "metadata": {},
   "outputs": [],
   "source": [
    "# group soil types by main categories according to name\n",
    "sl_conditions = [(soil_sota['label'] == 'Sols ferrugineux tropicaux peu lessivés, peu lessivés en argile, lessivés en sesquioxydes') | (soil_sota['label'] == 'Sols ferrugineux tropicaux peu lessivés hydromorphes') | (soil_sota['label'] == 'Sols ferrugineux tropicaux lessivés sans concrétions') | (soil_sota['label'] == 'Sols ferrugineux tropicaux lessivés indurés') | (soil_sota['label'] == 'Sols ferrugineux tropicaux lessivés hydromorphes') | (soil_sota['label'] == 'Sols ferrugineux tropicaux lessivés à concrétions'),\n",
    "                 (soil_sota['label'] == 'Sols ferralitiques faiblement désaturés rajeunis ou pénévolués avec érosion et remaniement' ),\n",
    "                 (soil_sota['label'] == 'Sols hydromorphes minéraux ou peu humifères à gley de profondeur'),\n",
    "                 (soil_sota['label'] == 'Sols peu évolués lithiques') | (soil_sota['label'] == 'Sols peu évolués modaux') | (soil_sota['label'] == 'Sols peu évolués hydromorphes'),\n",
    "                 (soil_sota['label'] == 'Sols minéraux bruts'),\n",
    "                 (soil_sota['label'] == 'Rivière')\n",
    "                ]\n",
    "sl_choices = ['ferruginous', 'ferralitic', 'hydromorphic', 'slightly developed', 'vertisols', 'river']\n",
    "soil_sota['G_Type'] = np.select(sl_conditions, sl_choices, default=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "collaborative-boulder",
   "metadata": {},
   "outputs": [],
   "source": [
    "# another categorization by sensibility of soils to water\n",
    "sw_conditions_ = [(soil_sota['G_Type'] == 'hydromorphic') | (soil_sota['G_Type'] == 'river'),       # grey soil, recharges local groundwater flow\n",
    "                  (soil_sota['G_Type'] == 'vertisols'),         # clay-rich dry soil, expends when wet\n",
    "                 (soil_sota['G_Type'] == 'ferralitic'),         # red and patchy yellow-red humid tropic soil\n",
    "                 (soil_sota['G_Type'] == 'slightly developed'), # slightly altered soil, erosive\n",
    "                 (soil_sota['G_Type'] == 'ferruginous'),        # dark red-brown savanna soil / ultisol\n",
    "                 ]\n",
    "sw_choices_ = [5, 4, 3, 2, 1]\n",
    "soil_sota['W_Type'] = np.select(sw_conditions_, sw_choices_, default=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "separated-invite",
   "metadata": {},
   "outputs": [],
   "source": [
    "# dissolve soil data with water type column\n",
    "soil_sota_class = soil_sota.dissolve(by='W_Type').reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "intensive-chest",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALAAAAD4CAYAAABMvGxfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6h0lEQVR4nO2dd3hUVfrHP2dakklvpFcIXQgdRBFEEBFFVteyKurPta24uqur7tp1XV17WdfuqmtvuBYsKCK9hN4DJCGV9J5Mpp3fHzPpbUKmZML9PA9Pcu895573hnfOnHvOe76vkFKioOCtqDxtgIJCf1AcWMGrURxYwatRHFjBq1EcWMGr0XjagL4QEREhk5OTPW2GggfYtm1bmZQysuN5r3Lg5ORkMjIyPG2GggcQQhzr6rwyhFDwahQHVvBqFAdW8GoUB1bwahQHVvBqFAdW8GoUB1bwahQHVvBqFAc+yagx1vD63rdoaGjwtClOwatW4hT6R01FLY9e8Tz7txzCdIXgglvOIXZotKfN6heKAw8CqstqKM0vJzw2jF8/2YChzsChjKPEDYsmdVwScy47DSEE7z7wCdu/3w3A8hdW8M2rK7nnw9uYecFUDz/BiaM4sJdjMpp45/6P+fqVH7sts+2n3ZiaTPzy4fr2dZtMPHfDq0w/bxJqtdrVprqEXsfAQoi3hBAlQoi9bc79VgixTwhhFUJM7qbeCCHEzjb/aoQQt9mvPSiEKGhzbaHTnugk4/U73+vReQF+fHt1J+dtpqq0hj9MuoucfXmuMM/lOPIS9zawoMO5vcBvgDXdVZJSHpJSpksp04FJQAOwvE2RZ5uvSylX9MlqhRYytx3t9z2ydh/joQufdII17qfXIYSUco0QIrnDuQMAQghH25kLHJVSdhkSp+A4Ukp2r9nPN6+u5Ni+PHL2OqfnjLthCMfqc0nyT3TK/dyFu8bAlwIfdji3TAixFMgAbpdSVnZVUQhxPXA9QGKid/1xncVfz/k7xTmlqDVqLBYreQcLnN6GsdLE5sItJA5L6EvH5HFcPg8shNAB5wOftjn9MjAUSAeKgKe7qy+lfE1KOVlKOTkyslNA/qDHYrZwaMsR8g4VkrMvzyXOC7Dj7/v59tpf2FCyySX3dxXuWMg4B9gupSxuPiGlLJZSWqSUVuB1wHvncVyEVVppqG9k49cZ1FbWu6XNwnUlFBUdx5vEbtwxhLiMDsMHIUSMlLLIfrgE20uhQhs+PPYxv9ywiWMrC93a7tpXtqCfFMiCa89EpRr4C7W9OrAQ4kNgNhAhhMgHHgAqgBeBSOBbIcROKeXZQohY4A0p5UJ7XT0wD7ihw22fEEKkAxLI6eL6SU8Y4ZQfqHb6fVU6FVajtdvr+187wn6OMGbmCJJGJzi9fWfjyCzEZd1cWt7xhJSyEFjY5rgBCO+i3JV9sPGkoqa8lhVv/MzanRtQ+zi3BwxJC2LU60lsnL2n17KmJrNT23YVykrcAMAqrVQbagjxDeb1u97j+7dWuaQdfZQvNfoafMN8MFQ0AeAf44dGr6H6aG1LOR8/HcGRQS6xwdkoDuxhyg0VfPLpcn59ahOjh49g7eebXdaWscZEnbQwakUiPlv11EyrxIoFHT4YCOTIzCLUPmqGnp9EYFiAy+xwJgN/lD6IsVgs/PzeWr646nvK91S51HkBqrPriM9JpEHVQOHUPHRSR+DqMAILgrE+ocFqkugCtcTeEM6WQ9swWU0utccZKA7sQb55cyWvX/+e29oz1ZpYf8kOQhtDCRKBlPypmq137kH6W4laantVqS9q5KczNvHK796hurjGbbadKIoDexB1Csx6o8tYKJeiytZwYFY+hRtKAaj+soHdyzLblUn+vzg2l2yhyljldvv6guLAHmTRvAXMPm8mU/55ilvb3XzjHiwGS8tx4KU+qH3VqLStS8hb79rD1vztvHzkNbfa1lcUB/YgFcZKdm3ch7HUNWNNXbAWlU6Fb5gPo38/FF2wFqCd86bfPQrtYR+i3wtixv/Gt6tf+3Mj5jwreQ35LrHPGSgO7CGsViufvfgVK29Yz66nDrqkjcmfjWb6K+MwVDZRl99I+oo0gocGtlyf+XU6JRcU0DSykVprLdb89kvIZVvKEf5wtC7LJfY5A8WBPcSbf/2AXTv2EjW10zqP0zD6N2EebkQ/xI/CNcVsmbm/3Xyv+ahtRS5fnY9EUri8rF39hN/EYdpuZd+Wg1hl96t3nkRxYA8x+48zUN9mRmhcE7qYuiSBIlUR2modE18eyaTHxrS7nnJePMbpjfgIHUGqQEJUIRSuLW5XRhOnQn+WjkjzwI0CVBYyPMSXD31PWUYDpTsrHCqv9lW3G7v2Rum2CsblDMcSY0alU5H5ck6767HXRUKVxBhmJORwGPnvlWCqa10+Th6bwLwz5zA1fDLa0VqH23U3igN7CK3QOOy8I65KISBVz7YH9gEQNiaEpPNj2PHYgW7r1ObWs/7SHd1eX3++7dqEe0bRIJposi8tA6hUgqsfvpSZkQM/ylVxYA9gNpnZ/nPvATUAp76RTt6Hxzn0TnZr/XozWZ84Z2Zgx6OdPwSnzBrtNVvtFQf2APs2HKIoq7jXcvooP2p21GNptBA/O5r81ccBqMmpc5otUZMjMNYZiZwUjk+IluPry5i2cKLT7u9qFAd2M8YmE2/97QOHyjYUNxK0xI/gMwLITcyG6c63Jzw9GP21WnT1OgwBBi6+4HwWzpzn/IZchOLAbuaDv3/O/o2ZvRe0Y9KbKAjIQY0KXbAWY7VzFz32v3EU3gSkbeGjPL2etEeHMebUEU5tx1Uo02hu5soHf4vWx/G3+tJ/1BCqCgUhmPLhGHxCdWgDnNzv2NcvjNUmjv1awG2n3cuDFz7J8ZwS57bjAlymzGMvlyOE2GNX38locz5MCLFSCHHY/jO0/4/iHZiazASG+jtcPuebAiz/tsk+5YXnMuKHeJJWRTDt6XGuMhGA9cu38NcFf3eKcIorcZkyTxvm2NV32jr63cDPUso04Gf78UnB1y//SMXxql7LqX3VzHjNFptw8J0sImQEYYTBx1r8/hvEvmdd71j5mUXcd/4/yXXRVn5n4C5lno4sxrZRFOAdYDVw14nezJsoPtb713L8mdFEPRqMrlxL7GlD0AZqqVfV03C3lfxVx91gZSsVRZW8dc8H3Pfxn1FrBp4AoKvHwBL4UQixza6w00xU87Z6+88h3d1ACHG9ECJDCJFRWlrqYnNdy8FjmezYv7vXcqYaM77H9ZQOKUE81YT5gTo0aIh/IILTvpzAKVtSCBsT4nqD7axfvoXnb3rdbe31BVc78Ewp5URs4iY3CyFm9fUGg0WZR0rJN/Xf4hut63Rt3O0jOO2zCWgDbS93xRllrF+yg6wFxQS+Ygv2ybuonJy/HCcr+iiB+4Op2FflTvPJz3SvPoWjuNSB7dvskVKWYNuG37y8UyyEiAGbyAkw8F93+4kQgjEbx5H5QU6na5n/ycEaZMVqah/xZahoomx7JcnZqdQVNHB8Sxk+jwax5Q+OreI5k/P/cLbb23QEl80DCyH8AZWUstb++3zgYfvlr4CrgMftP//nKjsGEgv+70zyDxWSe7CA/Rsz0Qf6UlfVQF1FPRvm7+qyTtnuSsoua9U9zP7a/cHlGq2aaYsmub1dR3ClMk8UsNz+oqcBPpBSfm+/7ePAJ0KIa4Fc4LfOfayBiRCC655or+lSePQ4/1z6Isf25+M/zI/kF6IwfGZBE6hm+8P7PWRpeybOG4efv6+nzegS4U1CbpMnT5YZGRm9F/QyDm/P4t7zHmfUk8nUpdWgM+sQCLbM7JsDD70wkfAZwWy5w7lDjKdXP8S4WaOdes++IoTY1mEqFlBW4gYENeW1VBRVIk0SlVSxd0k22Tf3fbrMVGOiclYZ09aO6b1wH0gaHe/U+zkTxYE9TMGRIl6/y6YNoQoV6I8FEDMzstdY4dQlCS2bNJsp211F+JYhlD1X202tvjNyWhrBEQNXZkoJ5vEw/7n3Q47uzAFg3QXdB6B3JHxaCNGLw9lw9c6Wcw3FjWxa1vXL4Iky7vRRTr2fs1F6YA9TmlcOwNib00g6J67l/IxfTmHKY2O7rbf9wf0ceirH1eYRGB6I1TowN3SC0gN7lPzMQg5vyyJ5YRy+l6s5cn7rruB6fT1l33YfuG4xWCjfW+VwWxq9mslPjqV0VSU+4Tr8EnzYdt++Xuu9+df3yT2Qz5QFE5hz6UyH23MXigN7kPzMIkxGM3E3DuHon/IwlNn2pQ3/XTLiMx2F65y3vmNusGCZYCQmJhzjMTOaQMfjGla++yurP1rPzCVT0fUhFNQdKA7sQaYunEBQeGDLBstmAm/xwSos4OTUbZq9OkrTSyi+u5aqzL4J95mMZsryywdcbmXFgT1Izt48aspbZwy0gVqSz43FZ48eoXd+exuv798Lnko98F6ZBp5FJxE//OeXlt99QnVM+mYk9bdVYiwxsf4Kx2ck3IGPn46wmIG370BxYA+yY1XritnwpcmU+5UR+k4UGfcOvKRNRoOJ8kLHdCzcieLAHqK+poH8zKKW44KfShiSFc2+lw97zKaglEAmPTIGvyGd4x6klOxZ072QiqdQHNhD5B8qxNTUusO4Yl9Vj0o6rkYboCH8Ix/8RuloLDF0WebTp7+iumxgqbYrDuwhjAYT+iA/T5vRQvr9o1AJNRk3dh9AlLM3jz+fcT/GpoGTO0NxYA9xyumjPB7h1Zaq2WVENka2zEV3R+6BAu5d9Bgbv87AYnFcbNBVKA7sAYwGI2WFFVSVev7rWKVTkbwwDoO1EXWBY4sUO37ew/2L/9myDO5JlHlgN1OaX8695z1G1q5jnjYFgGk/juGYr82W2i0NDteLHx5DZLzrxLkdRemB3czH//yyX86r0qmY8dp4pr3QXthkyhOnMPG+vg1JEufHtjgvQPE6x3vU0OiQAbHN3mXKPEKIBCHEL0KIA/ayt7a59qAQosCu2LNTCLGwq3sMNr55dSXfvfnzCdf3DfNh+oqxVKdXULGhNRH4iKUpFMw6RvGavn2lt036HaIKoTijrIfS7Zkwx72ZlbrDlco8ZuB2KeUobLqKNwsh2nYRz9oVe9KllCv6YLPXUnikCKPBRFByQLvQSUcIGxPCxCdGkxOQQ1BjMIc/svWcqRckEDUvnOhv4yn4tXfJ1rYk3dka1xB61PHhQMKIWM6+Znaf2nIVvTqwlHINtk2cbc8dkFIe6qVekZRyu/33WuAA0Lf/tUHG3Ctmcdp/JxL5iT+mB2qJmhzhUL2QtCDSHxhB1qu5pFSmorKoOOWW4S3XrSFmtj/St/1zM14bT3aYLfuQRmg48PfsXmrY650/mUe+vpshiQNDo8MtL3F2aaoJQNtkwMuEEEuBDGw9dWU3da8HrgdITEx0saWuw2q18vxNr1EfUk/CbTFkh2SRcH4UAUl6jn6R26IQ2RVpv0+iOqEKlY+ajefvwmpqLZz1ZR5ZX+Y5bEfYqGBCRgZxbFyrwyYUJrF+X8+LKAkj45h14XSueviS/kiKOR2Xv8QJIQKAz4HbpJTN80YvA0OBdKAIeLq7+oNFmSc/s4gDmw6T+30he246SqIxkfIFJTT+pYqRS1N7rLvtvn00qhvgKQMzPkxvkVf1DfPpkw0qnYrhLycS+ecgkmtSiBARBKoC2f/X3oUCb3jySq5+5NIB5bzgYgcWQmixOe/7Usovms9LKYullBYppRV4nVbFnkFJdVkNG77c0nJck13Lpll7EY/6Eq4KR3czTH2me7lUq9FKaKltjFqTWIlQ25xo0nujGHZxUs+NC1pSyE57aRzm9ZKC+8rZMH8Xu6dnUbK0lsqDvc9HB4YF9FrGE7jMgYXto/omcEBK+UyHazFtDpdgeykclFgsFv50+n282SatwORHbXvdsv6Xh/oLX8qsZRyfmc/Mr9O7vc/O6zNJ2pOK7y+BLSrt6xfu5MgnPU/Jzfx4AjO+ssm0ChVoThMUb7PPVkgccl4hBIcbBqZOsCuVeWYCVwJ7hBA77bf7m33G4QkhRDq2kV8OcIMzH2ogkbM3r13U2fQXx2OZZMQnVEdTpZEDrx4lekkgJmnq8X+jvqCBvY8cpv54o8Ntp16QgDHOgO6YL9OeHkfT6EYs/1NhqnUsliEw1J/7Pr2dhJFx1ARW917BAziiD3xZN5eWd1G2EFho/30d0OWASUp5ZVfnByPbinfQVv0o4669pN8zivH3jWTLn3cTmOCPRdpiCqrCK4lMD+tWE6I2t97hdlVawZA/B1OuLachtYHAYYHU3G7qlI0TAf6xeuoLOq/C1VbWE50yhIjYMCIIc7htd6KsxLmYEUPTEEKg0atR6VSYGyxk3LMXU5ktK2ZjeRP+/w4lTBWGESOx87qVSu4Tw69IIVeXS721HonE76sgCtcWE3t6FBp96wrapAdGE/6FLyFpNvGSsDEhDLvINtvjF+6Lob7n4B5Poziwixk/dCwjZgxl2henoNa1/rl3/MM2b1tf0MDhD3PYP/MYfu8F0ljUdSxuX8l8P4dYa2zLsX6ajnF/GoH6KSMxMyIJSg5gyhOnYG2SBIgAjDUmhFow6t4UGu+oZtyWVKasHkX86JgeWvE8igO7gaGzkqFakHh2q0P5x7TGAltNEqtJUrSmjLocxwNqesJqtFJ4VyUhqhAARJPAelkT8Q3xFGeUM+KGZApmHaN6SRmhmeH4x/ox/V/jKP2miqSGJMqsZVwQdz5a1cDaRt8RxYFdzP6NhzAWmVh/yQ6Ofp7bcn70n4aRvKj9wmTpzgqnakEUri2m5m4jfio/suKPEtAYgHkLTP50NCq9CoEg+KtI9tx3hKqjtey47yARVwdRFlTGwpgFpAUOc5otrkIJp3Qhv366kb9f8kyn8zPfn0DJV5X4RvRtIeJEyF91nMCL/En7JA5hVWGd1YTVqEWTosVP5YfvXDXJ50UScCAIVZBg43m7WfTKmVyy1DskmxUHdgGGhiaydh/jX7e82eX1jGX7aKo0us2e2tx6VD/pKJlfTKO1kQR1AkXR+RisTdT611L3sInK+46RWJyMtFjxSRrYw4a2KEMIF/DR48u59dR7qCrpeu7Unc7bTMa9ezk6tZjg1yPJ0+RhsNpmF6yvawi71z6VZ4BpK8cSPzy2l7sNHJQeuBe+PfgdC0cscCgGIHPbUZ657hWydw+M3RZd4RutI+z4UOpiajAKI/W/r8QotGCF8tRSbht+C4n6BE+b6TCKA3dASkltZR21FXV89szXrHjzJ3KfKmLRb84mIa7naNDnbnytRet3oLL76UwsTbaFE6ESDF2SQMif9YTKUMICQ0nRJ6NSec8Xs+LAHVj98Qbee+RT/IP1HNhkExn54tbv2PnWfi7683kEhQVg9DVyZF0OMxdNZfikoTTWNfLkNS9xeFuWh63vHYvB0pIJ9OBvcjnyWS58BglzoxELNGzI2sJp5073tJkOozhwB9Z+sYncA51zA2ftOsYTV/2r3bmvX/iB659cyuHtWaz9fHOnOgOVpNtiMFuMiDZifXk/Hyfv5+Oc8UWfc1F6FMWBOxASGdzj9XG3j+DQW9k0VRqprazn6d+/7CbLnMeWy/d2GdATMjyIYyOyaE1jPfDxnsGOm2hOqTrt6a7jc7M+zvfILIIz6S4aTR/lx7CwoW62pn8oDtyGV/78NrkHbJkwN9/edVLuunzHI8K8DX2sL6n1Pe8OGWgoQwg721bu4vPnvvW0GR5FHFYRGz+wFNh7Q+mB7WgHWO4HTzDrohlodd71d1Ac2E7h0b5pKgxGPn/uW4qyvOvv4DJlHnu5BUKIQ0KII0KIu9ucDxNCrBRCHLb/9Lh2fe7+1q3p8bO962vUWVSVVHP77Ad45c9vYzF7XnnSEVymzCOEUAMvAecAo4HL2ijz3A38LKVMA362H3uU8iKbLIVKKxhy+sDcPuMOSvPL+fy5b9slnxnIuEyZB9tW+SNSyiwppRH4CFhsv7YYeMf++zvABX0x2hWMnjECsAWX91XlZrCh1WmwWntQWhlAuHIMHAe0lYzJp1VaKkpKWQQ2CSrAORvB+oE3LAO7iwv/tIjwAZiRqCtc6cBdhW/1+WMthLheCJEhhMgoLS11glldkzQ63mX39jamLpzoaRMcxpUOnA+0jcuLBwrtvxc3i5vYf3a7j8Yd0lK1lXX8+M5ql9zbG/non50UEwYsrnTgrUCaECJFCKEDLgW+sl/7CrjK/vtVwP9caEevVJVUk7PPcYG8wc6RHTmeNsFhHJlG+xDYCIwQQuQLIa4VQiyxq/TMwKbM84O9bKwQYgWAlNIMLAN+wCat+omUsjk9+uPAPCHEYWCe/dhjnCzj30c+Cmf+5b2Pbc+6/HQ3WOMcXKbMYz9eAXQSr5ZSlgNzHTfTtcSP6PsWGv/IIGJPH8XhL7wjjPLyuyLwD7Lw4/tdqti2kDw2gTmXneYmq/qPshIHDJ80lAlz+yaZn3DOBA4v39J7wQHAuNMCmDHfxJ8X9p4qdtZFMxg2IcUNVjkHxYHt3P/p7USnODabp48IQvjpQA7suVJff8Gtz0byyHvFPH5DE91I1bVj/Zfe8aFsRnFgOwEh/vz36EtcfMf5vZeNCeHwu7+6war+Mfe3oZy2sIR9GYnkH3FMsurozhxe+8u7lA3AxN5doThwBy6+czFRSZ2n65pFpQH848IwNw78oPZv367gpjN1BIf3LTXsp09/zeVJN/H5s9+4yDLnoThwB4Ijgnjg8zsYktg+Acvwi2YA4Busxyc8yBOmnRBlBSY2/6gjZUznDPQ9YbVYaaxzjtCgK1EC2rsgODKImrI2wSxCcHyzbYdyaFoMR77Y5CHLTox3H3M8/1tbfPx0TrbE+Sg9cBeEx4Ry3k3zW09ISXWObRm7KOOoVwwfnMFQL5iNUBy4C4qyijm45YinzfA4ufvzqa2so6Zi4IZWKkOIDlSWVHPTxDsxNAxsZXJ3EBQRyONXvkBxTikPfH4HCSMGXp5KpQfuQOiQYJY+dImnzRgQvHrHu2xZsYNj+/O5a94jFB493nLNYrHQUOt4whlXoThwF5Tktg/bDE6OJHKM9wjeOYuKotZl59L8cq4Z8Udeuf0dVv73V84PvJIrkm/iq3//4EELlSFEJ4qPlbLms/azDLFnjSf7s40esmjgYLXKdnPDRoOJF5e9wdGd2fzx39eh1qh7qO0aFAe2U1/TwBt3vcfe9QepKKpEHxFEQ5ktCWD+d9sxVA1eQZP+suKNn9mxai/DJqbwl7f+gF+AX++VnIQyhACMTbae5JtXV5KzN4+g+DAS5rUG90RNHfi5IjxNUVYxaz/bxNbvd7q13ZPegY1NJm4/435+fm9ty7m4+elkfrKh5bhiT25XVRW64JGLn+HxpS+w85e97RI8uoqT2oEtFgtv3v1+uznfgOgQDPnlSIv9jy8EtV4S2DJQ+Pm9tfxl7kPs+HmPy9s6qR141Qfr+OL59npoal8thRszW46HnjsRnX/f4ggUbHSls+xsTuqXuH3rDnY617xk3IzaT0d9ae8Z3RU6Exod4vI2TlRaqldpKCHECCHEzjb/aoQQt9mvPSiEKGhzbWHH+q7EarWSvecY9TWds2JqfNuL25mVFbkTxtff9XnwTlRaqldpKCnlISllupQyHZgENNB+H92zzdfte+fcgsVi4XeJN3L9+DtY/fGGdtd8Av0YunhKu3M1XiZ2N5DIO+j6IcQJSUvRd2moucBRKaXH809lfL+T8sKuNzYmzh1LzZHjSGwKLCMuP50yN4zjBiMarZopC9Jd3s6JvsT1VRrqUuDDDueWCSF224co3e71drYyT8Hh491e04QFognwJXZSKjETUjBXK4sXJ4rZZOGW6X/j8HbXSha4fBbCLmpyPvBpm9MvA0OBdKAIeLq7+s5W5rFarV0mLRQaNahUaMYlU7gti+BTR1CyPbvf7Z3MNNYZyPhhF0aDkdL8cpcE/5yoAzssDYVNXnW7lLJlMCmlLJZSWqSUVuB1bEqWbmHJrQvRB3Ve6kz7zVSackvRFle17N1trKhzl1mDEq1Ow6Gth/nnVf/id4k3cvfZj2A2mZ3axok6cF+koS6jw/Ch2fntLMGmN+wWSnLLaOpiZuHI8q00lNVw8JMNRI5JQG00kTBrdBd3GJxMf3F8u2ONvv+BOSajmfVfbmXNp7ZAqAObDvPXcx7lpT++RWO9c/bbnZC0FN1IQ7WVlrIf6+3Xv+hw2yeEEHuEELuBOcCfnPI0vWA0GHn7vo8wmzqrj1tNZnxmjkStURM2ewwanZbq7J6+WDyHf5zeEYmHPiE6rAhMf7frNGP9ZeeqvXz5r+8oPNL9u0hf6I+0VCdpqC6kpRqA8C7KXdkHG53Ca395l18+Wk9ZQdfLwkPGJaGvb0KaLahqGzHpfag4XORmK3tnyuOnUDKnkMBb9Bzf0vtmzfBTQlDr1JRsK++xnCqmQ1+mApVORUhaEBX7qvphcdf8tPYg8SPj8Olncp2TZilZpVZ167wA4aeOoGTdAQJjQ0GnJffzgbnzOG95EU1/FA45L8CoR1JIuy2xxzJqXzXGyNavdG2gFikg5bx4Rj3c+8bO6f8a32uZjny3M5t//O55fnxnNY/+7jlK8k5s5/RJ48Dzlp7R4/Xqndnoo0KIGJeE1Uc7YJePj28u4/hmx/+zTeFGsO+O9wnVdRrrAkx/Yxy6gtZ4j/R7R1L0bhlqPzX49R5RZqm3cuob6Q7bBDDE35cNyzfz5DUvsfqj9Vw76rZOcSmOcNI4cE+at/qIIELHJKAN8kPtp6NmU2/pP7wHyxqB+biF0dcPY9T3iVTvbD+zotIKCJPkvWGbJIqeGoFunJqs5XmYa81IB14Dtt61BznMwim3DnfYrsb1B9odGxqayNx21OH6zZwUwTxGg5FfPlrX5TX/yCDiLpmJMBg59tNuLE3OnebxNBn32CZ4Yk+PQv+ML0c+Pkb6XSPxn+RHfVItKqmiRlVN6LggEm4ZQk5kNk1W23Di+KYyYi6O6On2LRy6OY9hr8aSmptA1vLexcKPdzHHPnxS3/M0nxQ9sMloJn3O2E7nhVow/Nq5qAFTg3HQOW9bCtcWc+Rj20q++cJGypNL8a8MQGvSom/SEzDbF2u+JH59EkIKpr0wjuFXJ1GVVsGwS5J6vX/Fvip2LMgk5rIITn073dazd0PqkoROQ5m0Samcfc2cPj/XSeHA/kF6zrtpPpHxrRMiQi2IGp9MU6MRa1k1havdNhXtEJMeGsO0Z10zlRVYH8TRM49jPQabT9/H5tP3sX7RTjbetIvNt++m4X0z2lg1+jRftFJL+NRgh+5rqjOz/uIdNB4xMv3N7l/sov4YgiauvestvnkB/kH6Pj/LSeHAAD5+Pjz+433ED48hYcowUs+ZiG9kEHmfb0QV5E9tNwE+nqI+qxEmW5w+3wuw6ay9mOrMiKFdv6CFnOvPuot2sOH3O/E/Fsjm23f36f7Zn+RTPrKU0b/vekhQ6V9J9erWOJMzLp7B/Ktm96mNZk4aBwZIHBnHm/uf499rHsanzkDI6AQSZo/FkOu69F0nysF3sii4r5wp/+g89OkvzdulSl7p5kPbYPvUTHpoDE3xfY9fqMqsISwjgvrrqpj86NiWD6FQC079Tzo11hqqD7TKVY2cmtZlfIojCHdsvHMWkydPlhkZGU65l9lk5pUXv2Pbl1vJX39wwKute4LUCxLwi/Vl378P97luUHIA8Z+EUWWtIlgVRIApkBptDbVWm+PWXWqlJqcOIQTvHv0X0ck9BzQKIbZJKTvl5T4pZiG6QqPVULcvj+PbjrY4b9LcU2iqrOvyDflkJOvLE089VpNTR9lNPsS9FEWxtZhqdQ1YbdcSmxLZlGN754gfHtOr8/bESTWE6IhJStJ+15qRR58YgdnQNzXzwcyoa4cy9cm+Jb9pS8m2cnZMP4zm/gDi1iQRtzqJ8Pej2Xq2bQ54+OShPLj8zn7ZeNL2wAAaoxlDGyHrA2+vVoYSbQj9nT8Hbs7p931yfywk98fCdudGTh3GTc9eTeLI/ilenrQ98O6tR9j+0x60IW2mbpzkvNHTIpi+xvkvX+7G4G+g8qDzl9T9In257OEljJ4xot/3OikduMlo5qN//0BlSTVWFyxeJF0cS9VbnXc8exMqraBW5Xxha58QHfO+OZUpZ01wyv1OSgcWArI3H8YnwJcTnL3pEc0IFYc/8vj+1X4RPSWCIIvzk9kIjSBbZPPzJuekKTspx8ANtY2U2rd8Z/+wy+n3X79op9Pv6W4KN5RScZ7zhw+xM4dgEkZSfZ2Tf+Ok7IGzDhQgsM2vW4yOzzpMuHe0bTeEk4no54uMqzCUOU/UZeZ7E5j+r/FUHqomMiOGoeOTnXJflynz2Mvl2LcO7RRCZPS1vqs4tM82vymBkdee5XA9v9ka9EOcq5PmHxmEqotAo8HG1pv3suX2PSy+fCFXzuxuk0/fcZkyTxvm2NV32q6i9KW+0wnw0aLz1RIzMYW+DII3nLWb0h3OVaoMTYsh0APK5u7GWGPi9F8ncPb1s0kZm4ha7ZxndpcyT0f6W79fnHfVbKadO5GQ6cOp2Z3jzqY74T8xBfXxgRVI5ArG3JhG0pAEIiI7bZHsF65W5pHAj0KIbUKI60+gvtOVeZox1DdRuyObgjZSqo4y9uY0Upf0P+mLBMpMVrK+dk58x0BFqAU+VwlOCXb+UMnVL3EzpZQTsYmb3CyEmNXXGzhbmaeZcqP5hJxXpRUIBDG/632nwpQnel6GjZmYQqRaDPrl65RF8YRlRRBdG9N74T7iUmUe+zZ7pJQl2JQpp/alvqsoK66i/gST9vkN8aMqswbDvt7f0MVMC/FnRnd7vaGkmrxPB3/2o6z/5eG/IYjoxBMP2ukOlynzCCH8hRCBzb8D82lV4OmLso9TsVqtfPVVBmFdiJu0Q4A2QIMuWMvMDycw862JRIwLpb6ggbyVx9n24P5e21Lt1BKU6t/ltaQ5Y4iblz5gdz87i2EXJ3H6+5P4w3PXuOT+rlTmiQLWCSF2AVuAb6WU39uvdVnfFVgsrY6ae7CAK1NvZv2eXHL/t6XbOtNXn0Lq5kimfD2ayZ+OpvTrSiyVZvSxjqePEmoBBsnBd7pWZ9SMTaJie9934Q50ks6JIyDe9qFVaQVR54fBFxpUKteMVl2mzCOlzAK63BglpSzvqr4ruH/xP8k/ZIuEqiqpQRseQJSPBmNcWLe93+5LjhA3ZwjmSZLCb0vIX30cPuhbu9Ii2fbgAaymzgFCWr2OcoORyl3HXLFjyKPEL4skdE0gO/95kNkfTGOaz1TGPec6jblBvZS89fsdbP1uZ7t0TyMvmErZmv0c39V9rEJDcaMtluGj/rVvqu365WzYJbYY5Kr+3X5AodKpQEq2X3OQpkrb+8GZ+jlMO2eia9t16d09zPM3vd7OeYVaIPx0BKXFeDTutz7Ah7xvBtfU2eR/jCFl/RAmvjAKq0mSMDKOcW5Q9xzUPXBdh/SwyXPH4RvgS8muHM8YZKfsw/UtaWwHC1vu2EPq4gR2bjnIiClDeejLu/BzQ3qyQdsDSykZf8mpAATGhpI0ZywNxVUY6hop2Nz3TYrOxF3Om3ZpElM2jCJ5kXuChbL+l0d9USNL/ngu4THuCW8ZlA5cW1nHnWc9RJU9ziF2zliKt2dRlVWCtdoWaK720eAb7PzIsoFExBkhFKgKiPxrECOvSnVbu87S/nWEQenADTWNVEiJRqu2BetIMFQ3kHrJqRyzJ6MevnR2nwJ5vBFTqW23SdE9lYSeFsTMD52zC6InNFo1kYmO6ak5pT23teQmcuvzePLGl/CLiqKhoJKRNy/g4L++I2xYNLLBiMHeA5vLajFUDd4sRD6hOpoW1JPUkETd+CY2XreT8FNCXNLW5AfOwdRkZt+LvzJr8VQWnIDG2Yky6By4fq+BrF9zGfrbeCx1BmRpNQBxs0Zz7MedLeUOL9/sIQvdg7HGRPU9TRw7Uo42wBa6WL6nyuntxF4wlouuO5vEkBDksgvB4t7ZnUHnwJEJ4ZiazGh9NVSX1GDJt0nrq3w01OSfPFnnpUWSu9I5KRK0PlouvO1cIuLDSRgZh8Vs4dt131OdFMrp08YzKdb+kqjvetnclQw6Bw6JDGLC3LE0SoEuRM/xbUcJHxGLtUPsg2+YPwHRoZTtz/eQpd7D+TfN59rHLm93bsK8sZQaKojROz9Apy8Mupc4jVZD+txxSIMRjZ8OQ0U9PsF6zNb2X20qjYaoFNfPU3ozoVHBzL96NksfuqTTNY1K43HnhUHYAzfWG/jg0c+RKkFgfDjRE1PwCQsAbftHrS+p5vd3BqK5I5SAAAMPXiE5dsg5ucsGA9POnchtr1xPRJxzd1A4m0HnwL56H2YunsLP76+laX8+Iy+fBSqBubj9th0BPLN5NKG+OqbV76C0sHW86J8ajvbOaWgsUHbzCk4mtD5a/vDcNSy6YZ6nTXGIQTeEEEJQX2ObKkuZPx6roYkyvQ/Z9vnfZgLSItkfbKbyQBlvPVhEQ23rGNl/WDjbDBXoC52f23egc+riyV7jvDAIHRhs+90koB0RS8X+fAJ1mk7bdvxvmYyUIML8uOIvYbz0SzBC2IWff8wkxi8Ac2bPyQG9leiUIdz5zjIeWn4n8cNbt/lofbTMv8p9c7jOYNANIaxWK/mHCkmdPx6ZV0bE5GGIutaeVKVTEX5qCjtUtUw7CJZYf+qCEzAbs5HStjKn8tVQbmxkRKQ/frFB+EYHUrm9wFOPdELc9OzVxKRGsfrj9az6oDVDU+KoOP763q0Mm2BTxlFrVGz6ZhtTFkzA0NDE1HNcv1rnTAadA9dXN1BWUEHob6ZR+dNuaqek4ftTa46H8fdPIjMtkAnFPuQ+9x1BwyOJvFXP+09ZW8pIi2S0LpiSCRIxYSb+ZRavcuDTfjONC245B5VKxfRFk5h/1Wy0di2MgNAA4tNae91p505i2rmTPGht/3CZMo8QIkEI8YsQ4oAQYp8Q4tY21x4UQhTYFXt2CiEWdqx/ovgH61l067lYdBoizjyFVLWK6rzWoYDQqIkNF/g1mYl6cSHyjil8uzuaLStbI8SkyULtDSsIaRTsq69E1HvXruHq0pqWnBNCCCbNG8+4WaMZOTWtnfMOBhzpgd8G/gW82+Zcs7LO40KIu+3Hd3WoZwZul1Jut2/u3CaEWCmlbN4N+ayU8qn+md8ZlUrFsqeW8s6Ha9m+5SgH/rOqZduOSqfCd7gPuYd9KPvLCoRGEDIuFhbaMkyqfDXE3ToTa6gPTYEaLICfSoMl2MfZZrqUgsNFSClPOHGKN+HInrg1QojkDqcXA7Ptv78DrKaDA9sFS5rFS2qFEAeAOKD37bz9xNBkYs3WLIrf+KndnrNJX85hQ6WZfHM1Q16bhwqBvkpg3mzbMxc6OQEhBKriBmqf3UXw+Bjif5tC0yfelXp21PQ0l22iHGic6Bi4nbKOEKLHJRn7B2AC0DaCZpkQYimQga2n7lJfya7ocz1AYmLPWdebOXS0GP/iKqxI4h+dR/49K4lbOok9BTqKjSZmVPpR+PhqLA0mmke2Kl8NdVeNpBQLSQVq/FPDKPkxE37MxNvmIiae1ffs8d6Kyz+mQogA4HPgNill80DzZWAokI6tl366u/onoszz3ep9mOsNxF0zmV/9awl5ZSHrJunYV1/J5AI15l3FRD9wJgk3TSfuStsLjNVgRvf0NmosJtaHNeI30nkqQO5k7hWns+hG75nH7S8n2gMXCyFi7L1vt8o6QggtNud9X0r5RfN5KWVxmzKvA9+coB2dKK+sZ/XmTNKmRlLy/jZm+4/GmAiRvn4UNzawekgjZ/hHY1UJZHwgjd/Y5KVUPhrUATpCbvyFYMB75hxsBIYF8MDndzD+jDGeNsWtuFKZRwBvAgeklM90uNb2VXgJrYo9/WbTjixGChV1+47jc9tUTEfKOf7HFfjfsIpTV1YzVRdGZoQVYbFiWptLxUbb9npLkxnz9eMImRzvlVoNNz591UnnvOBaZZ6ZwJXAmV1Mlz1hF77eDcwB/tTfBzE2mfjmx13s25jJwffXUPTJLtTLD1OzMImgVxeS+PBZYJWU/2EFSR/mYAjSoJvSutkxfukkTEJy6JpU4m+cbnt2jXe48pDECGYumdp7wUHIoEk1u+qT9byz9iCVH62jsbyu5bxaryXqkbls1tVitlqZlael+D8ZGMvbZBFSCeKumEjVmizCzhvF/rG+jKzRUhqiwnjdd0SckUpdZhmGooG5Ff65dX9nzKn9T1k1kBn0qWYP7cmj8KXvSbjrDGg0Y04JoslXYH1zN4W3f09aTBCRV08kd4QPOfePJ903FL8fcjn+xV5iL00n+7QwhkboaRgWSK2phjXaBjSNKpJUAsulI4ksNJD3yCqn2526OIGgEf7sfOLgCdW//N4LB73z9sSgmCy0WKxk2JVu6lZlsXaolXWikq1NFWRek0rC3+ZgrGwg77HVcN0PnLq6Ft8miTYh2FZndxEpGytpGBaI/kgNw3VBBGh1mK1Wgl8+h73mGtCq8I0JIv7ReQQOd3CGQsDMr9OZ+sy4bov43aPCf5zjooERcWG8sPEfXLDsHOZfPZurH77U4bqDkUHRA5cUVnDMLpRXuTWPtOvHENagwvq/THxGRGBNCiLyifnoqozk3v8TRcv3wvLW90b9sHBUei0Vt/5AuX3rURQwNC0SeUcYBosZShsouH8ipSoLQxodXFqW0BBZj74xoNsiwTXB7LrbcaHt+uoGhk9KZdS0NIfrDGYGhQMHBvmhjwumscC2A1nzwAZKKxqwGsyw6gjaIF8sRjNWg7nTDIM21A/dyEjyn12H7LBvru5wKdq/rMLn0UmsHmqBJgMTXs2mJq/KYdtyryvFL7z7jJcbztrd7bWueHLVg6hPgqQwjjIohhA+/j4Ejh6CPRqSxsIam/PaMdUY2h23xVxnJHeoD/4vzkcTqOt0PXRaAsY2GsONfXyRK99TZZNndRI15c5P/+rNDAoHVqlU+FwyBlWbnkkCQ15cSMLfeg7QliYLYW8dJNtcT8jjZ7V8CAA0gTp2LY6m7TxN9LLpjhklbMn9uhMTUWn7NkUXFh3C3z64jdEzhvep3mBnUDiwWqViQtgQYi9NJ2HZDBKWzSDqRfuUs8VK1AvntCsfNjUBbVDrjuSKjcdIeGwXFTor8f88u+W8udbIOALx02iY4hPGTEI5kCCoeeVM4h/pOUFiSFoQhEsq9ld3ujblsbFMWNM3R0wZl8ScS2fiHzS49dz6yqAYAwPsNlZhOFZB2dpsouYNx1evwxLmg7WkntLn17crW3vdGKIPNZD71BoS/zCD3Jc3EvGbMVT8Ywtrbx3OxHExVO22bfIsXbaCGFoT5QX764i5dDx1u3oWDanKrGH9OTtbjmNPj6KpqonyPVVYDJLjoojkhXHkrOh90drX34cLbzu3L3+Ok4ZB0QMDxFh06KcmkPjXOWjmpVK14hD5f/0RVXQAoc/MbykXMi6Gg/VV4K8leGw0q4dLfMP9ERoVwdMSOSNfh35E+2mytl/2lnoj+W9upSqjb4IoAQl+DH0lFpVORVOJEYCgUd3PToAtGP2Rr+7mo4LXmLLAu7b6uItB0wM3rsom791NYLIQOjGOxqIagkdG2eZ+2yAlaFUqMFkIOTuNSb6BVJXVk//SRhpfmMOwzCYKP+3bzIAjZH6Q05JnY98rhzlFN5yjn+X1WEdKychpw5RhQw8MCgfeX1rCt+mCKaPnUR+gwiqhTmUiOrOJmgPF7cpW7yliZtZQdo4Ag1niKxuI9NFgNZgJuHMd+XWuEzfxj/GjvqgRJOx5vvPcrxAC/2A9dVX1BIb6c+1jlxMSGewyewYDXu/AxiYT7975Hsn/20mpgKBH53BENBAmfCgd4Y/aX4el3tiuTsErm9C/MBuD2cz4zXXkN9mm2ExVrtGB0AZqMTeYSX96JFtu2Ntt8pfgyCDey36J9V9uJSQyiIlndb+Cp2DD6x04t6iUpiAd0X+cATo1qtwGxlUZaNidTfmGrjMRWQ1m4t87yvC5qeS9va7LMs4k7owhxCyKoG57I1aTtdtysUOj8PHz4czLTnO5TYMFr3fgYcmxPPv8MqoMjTQ2mQjz11NaX89PRw9TkV3G7n+spGBHbkv5wOGRqIN8KF+TTfmabLfYmPNNATnf9D7bMGFuz7mVFTrj9Q7cTIivHyG+tqCY+OBgrp44GSbC/tMm8Na3v1Lw1lbKNuQQcvFYSpN8oY+zCK4iIMSmqTs0PZmpC12bU20wMmgcuDtGR0Xx1P9dzLq5U3n+5185YGkk5cPDDBQdysd/uJekMQn46r1r6/5AYdA7cDOnJSVzyu+i+TbzIDtEGLq4YAr+u43QKQnoRw6h4L/b3G7TvKVnkDo+Ca1O6/a2Bwu97sgQQrwFLAJKpJRj7efCgI+BZCAHuLirbfFCiAXA84AaeENK+Xhf6nekpx0ZfWVjXi7/Wb+JLVUlTC5Q0bAxD+tFwyn/g+vlVE9dPIVFN85n8vzxJ4X4iDPobkeGIytxbwMLOpxrVuZJA362H3dsUA28BJwDjAYuE0KMdrS+q5mRkMjzFy7hl6uu5arfL0S7bBK+m4rRJ4ag9u8cleZMIuLCmHJ2uuK8TqBXB5ZSrqE1FKCZxdgUebD/vKCLqlOBI1LKLCmlEVvq7MV9qO9y/LRaQv38OD0pmdcWXYBlQTJhN0xBuHiBPf+wc5KvKLhWmScOaLtWmg9M60N94MSUeU4EH42GlxYv5rfB12A1W3qv0AVpz5xFpNrAhlu7nlsOiw7htldvIHZYdH9MVWiDK/uarr4f+7wF+kSUeU4UP7UGVQ/f6lJAwr2d44ubhxyNe6tZf1v3CyMX/2UxM86bTNKo+H7bqmDjRB24uFmcpAdlnnwgoc1xPFDYh/puJ7uiApOp+943eGQUpUm+qPWtswZCI9A+NxcJ5L+1FdHhIxqdHMkfX/o9Z105iwtuaR+XrNB/XKbMA2wF0oQQKUIIHXCpvZ6j9d1O9Z7j7b42gkZFtRM3qTlQTEg96MJbE/pJsyRiWyXxj8+nK5Y+dAnn3XQ2d71zi7KXzQW4TJlHSmkGlgE/AAeAT6SU++y37bK+pxkxMZXJZ9umtuKumIi4dRJhpya3K1P7whZ87jmVoNFRRMwZCoDQqvBZldduMSJtYgoJI+OISfF8LrXBzKBR5nEmR3Zmk1FdwtqXV3Fs9SGaStonBZcCfCP84bEzuLg0mNqpkVw5Lp38zELumPMgJqOZl7c9Qer4JICTRqvXlQx6ZR5nMiw9hWGkUGBqoGRBHOY/fIfFvqtZ5aMhaHQU9cvGceecM5idnNpSb/SMEXxa8ia/fryBiPgwxXHdgNIDO8D+TZn88NYqsvQmxOnxzB81kuFRkQwPj3C7LScrSg/cD0ZPH87o6cp29oGI8h2n4NUoDqzg1SgOrODVKA6s4NUoDqzg1SgOrODVKA6s4NUoDqzg1SgOrODVeNVSshCiFOhabmfgEgGUedoIJ+HJZ0mSUnba0eBVDuyNCCEyulrD90YG4rMoQwgFr0ZxYAWvRnFg1/Oapw1wIgPuWZQxsIJXo/TACl6N4sAKXo3iwE5ACOErhNgihNglhNgnhHioizJCCPGCEOKIEGK3EGJAigE7+CyzhRDVQoid9n/3e8JWULYUOYsm4EwpZZ0QQgusE0J8J6Xc1KbMOUCa/d804GVapbYGEo48C8BaKeUiD9jXDqUHdgLSRp39UGv/1/HteDHwrr3sJiCkWZ1oIOHgswwYFAd2EkIItRBiJzaZrJVSys0dinQldhjnJvP6hAPPAjDDPsz4Tggxxr0WtqI4sJOQUlqklOnYNOCmCiHGdijiFLFDd+DAs2zHFpswHngR+NK9FraiOLCTkVJWAavpLArek9jhgKS7Z5FS1jQPM6SUKwCtEMIjIhmKAzsBIUSkECLE/rsfcBZwsEOxr4Cl9tmI6UB1s0byQMKRZxFCRAu7vLwQYio2Pyp3s6mAMgvhLGKAd+xpFVTYhAy/EULcCCClfAVYASwEjgANwDWeMrYXHHmWi4CbhBBmoBG4VHpoSVdZSlbwapQhhIJXoziwglejOLCCV6M4sIJXoziwglejOLCCV6M4sIJX8/89JOKHeE4LcwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot categorized soil\n",
    "soil_sota_class.plot(column='W_Type')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "early-forth",
   "metadata": {},
   "outputs": [],
   "source": [
    "# export to homogeneous permeability zones of sota\n",
    "#soil_sota_class.to_file('./data/output/homogeneous_permeability_zone.shp')"
   ]
  }
 ],
 "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}