Finding_contours.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finding closed contours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Construct some test data\n",
    "x, y = np.ogrid[-np.pi:np.pi:100j, -np.pi:np.pi:100j]\n",
    "r = np.sin(np.exp((np.sin(x)**3 + np.cos(y)**2)))\n",
    "r[:30] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import measure\n",
    "\n",
    "# Find contours at a constant value of 0.8\n",
    "contours = measure.find_contours(r, -0.495)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUUAAADrCAYAAAD+FYR3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABrAElEQVR4nO1dd3wcxfX/zu41nXqz5d474IIx1WCaA4ReHEzvkASSQH4QQicEAoSQQKimBELv2IAJvRnjSgzGGPcmS7LV2+na7vz+mDe7e6uTdCfLtuTM9/O5z97u7c7Oztub+c57b95jnHMoKCgoKAhou7sCCgoKCt0JqlNUUFBQcEB1igoKCgoOqE5RQUFBwQHVKSooKCg4oDpFBQUFBQc86ZxcVFTEBw8evJOqoqCgsHTp0irOefHurkdH+Nnhmby6xkjp3KXfRz7gnB+zk6vUZUirUxw8eDCWLFmys+qioPA/D8bYpt1dh1RQVWNg4Qf9UzrX22dd0U6uTpcirU5RQUFBQYDD4ObursROgeoUFRQU0gYHYGLPXA2nOkUFBYVOwYRiigoKCgoAAA6OmJo+KygoKAhwAIaaPisoKCjYUDpFBQUFBQIHYOyhYQdVp6igoNAp7JkaRdUpKigodAIcXOkUFRQUFCQ4B2J7Zp+oOkUFBYXOgMEA292V2ClQnaKCgkLa4ABMxRQVFBQUbCimqKCgoEAQztuqU1RQUFAAIDrFGN8zY1SrTlFBQSFtcDAYe2jgftUpKigodAomV9NnBQUFBQBKp6igoKDgAoOhdIoKCgoKAiLytuoUFRQUFAAAnDNEub67q7FTsGd29QoKCjsdJlhKn1TAGDuGMbaKMbaWMXZ9kt+vZYwto88PjDGDMVZAv21kjC2n33Y43ahiigoKCmlDGFq6hlMxxnQADwM4GkApgMWMsTmc8x+t+3H+VwB/pfNPAHA157zGUczhnPOqrqiP6hQVFBQ6gS41tEwBsJZzvh4AGGMvAzgJwI9tnD8TwEtddXM31PRZQUEhbUhDSyqfFNAPwBbHfikdawXGWBDAMQDecFXnQ8bYUsbYZZ17IhuKKSooKHQKRurO20UuXd8szvksx36ygtqKwXMCgK9dU+eDOedljLFeAD5ijP3EOf8y1cq5oTpFBQWFtMHBEOMpdx9VnPPJ7fxeCmCAY78/gLI2zj0Trqkz57yMttsZY29BTMc73Smq6bOCgkLakIaWVD4pYDGAEYyxIYwxH0THN8d9EmMsF8BhAGY7jmUyxrLldwDTAfywI8+mmKKCgkLa4GDpTJ/bL4vzOGPsSgAfANABPM05X8EYu4J+f4xOPQXAh5zzZsflvQG8xRgDRH/2Iuf8PztSH9UpKigodApduaKFcz4XwFzXscdc+88AeMZ1bD2A8V1WEahOUUFBoRPgHGrts4KCgoKEMLTsmcv8VKeooKDQKaggswoKCgoEDqaCzCooKCg4oZiigoKCAkHkfVadooKCggKBqXQECgoKChIixamyPisoKCgAEJG31fRZQUFBwQHlvK2goKBAEPEUlU5RQUFBgaBSnCooKChYEC45iikqKCgoAFBrnxUUFBRaoStDh3UnqE5RQUEhbYjQYWr6rKCgoGBB6RQVFBQUCCJKjpo+KygoKACQy/z2zE5xz3wqBQWFnQzBFFP5pFQaY8cwxlYxxtYyxq5P8vs0xlg9Y2wZfW5J9dp0oZiigoJCp9BVK1oYYzqAhwEcDZEDejFjbA7n/EfXqV9xzo/v5LUpQzFFBQWFtCGtz6l8UsAUAGs55+s551EALwM4KcWq7Mi1SaE6RQUFhU4hjelzEWNsieNzmauofgC2OPZL6ZgbBzLGvmOMvc8YG5fmtSlDTZ8VFBTSRpo5Wqo455Pb+T1ZQdy1/y2AQZzzJsbYcQDeBjAixWvTgmKKCgoKaYMDiHMtpU8KKAUwwLHfH0BZwv04b+CcN9H3uQC8jLGiVK5NF6pTVFBQ6BS60Pq8GMAIxtgQxpgPwJkA5jhPYIyVMMYYfZ8C0XdVp3JtulDTZwUFhfTBuy7FKec8zhi7EsAHAHQAT3POVzDGrqDfHwNwOoBfMsbiAFoAnMk55wCSXrsj9VGdooKCQtro6iCzNCWe6zr2mOP7QwAeSvXaHYHqFBUUFDoFtfZZQUFBgaCCzCooKCg4wMEQN/dMO63qFBUUFDoFlbgKQCz6HbaW9rH2DcdvpjyH3CajZIqPUMjyMPfQ1mtd08x9AIBGI0Psm36xbwYAACHDb50bMsW5LYa4voX2I4YoN2KKrdMvKmqIe0u3APlbKrRfYzxh62FGwj4A+DRxzK/HaV9s/bTN0GMAgKAWta4J6hEAQLYWpt/Efiadk621WOcGNHF9JhO/BagOXiZa20s+ql7H43QUIN4pMymrGL3cMupJmGQm5RM2bZk1mlJW4rdQGjKTMorSVsouatq1lrKJcz1hPx2Zeah9NNr6dPup5W9SRlJ2GdT+qcgsk2SWrQtZSfkAQIDFaEvlk8x8Umb0GE6O5ZTZPnt79u3wQbsDuJo+KygoKFhQOkUFBQUFF1SnCEBjDEFNh8lbLy00aCoXs6Z0Yrqg07matTXti+TXHUgKpiFxmhtxTMXkb3LarLkUw+1522vWFIymZJrYl9MvAPDS9LmjabOcfiUca2PaLKdmAJBJUzE5BQu4pmA+4eDvmorRMZb4wkqZGY5loT46Jcql7EhFQOfocqs5ZNYGjB1YHOVUSdiy4lRvOX3uWFYSbll5NHv67JcqD0tmoo3bk1lH0+agQ2YBS1ZyS/8H+t1LctGT6OM0xlrJrbuCg8FQhhYFBQUFG8rQAoCBIcA8CWxDwoRkMcQ6uMsYADEy685rNeti0I9pwzaISGZnP1JUsjypvGekvLeEKUZz5zTAKs/FQDVLQW+zjraYopttSKYhfhPHcnQytDDJGGlL7FCcK8qTbCNALMLrYoPeJA2nuxhHMplJZihlFKP76DyR5cNJxFzkQHexNB0ds0rZprrLMALYTN/DEo1iUmZy38ku7apJVm8klC/lAzgZYuoy64ghZjK7/KDFEMW+ZIZeeh6LySdh1jpjYD2ko+HK0KKgoKCQCK46RRG4zAMdniRtYdIoLRmjTixEMgfdYiE205J6I11Lom9Ea12RKE+yDFmuj85NZHaArdeUriZx2pesI1lUYL0NpihZoXTNEeWL71IfJV1o3HpDp85JMkRLh2ixDul+42A1VL2Ai21IZqiz1qxDa4NpSJmZDsYoWYtkkZqlUzQT2kDXbPYqdcOWHHjHzNAtMy8jtypI3Z/NdL2a+B4zXS45achMMkVvUp0ise9WskqUmWSFgIMpSt2ilBnJKtOhc3WzeUtm9IxSVm3Jqed0M10XEKK7QTFFBQWFTkExRQidos6SW5zkWG9YOh85Qkq2JlmIgylqUnclRm2dmKMkPnoSPZjucszVLPYh2YHtaBwh/aJHOgITS5COwe3BZoqJ5XucrECT1mHSS+mSGQomYTOM9lgHXevSRYnvGpXvoTpItkEWzDZk0R70JHuG1AfTfsxi+VIPbD+zJo+1xe4d5+ou9m7JzMUcnc7hcVMy8kTH+1SU+m4Heykzv0OnaMnMxRQzXaw+J5kemMlzpMwSdb3iu073Ts4MO5JZT9IpGmbPqGu6UExRQUGhU1DW5xShW1Y2ATdzdOq/3ExE6q5svZWToUhdJXftu5miU+cnGEiMmGHMWjoWp20yX7FEdmOVr7UuX7INuaTL9j3smHVkk5VZ6qMk2/A7mERAMl1Lh7hz/MLcMrNZDcmO20zLtkwTK0uF3btkJc/xcsmA7TaNuXSK0v+xveV+bpm534WAQydqy0wyxY79RbMlm7RYPaiMRCYv7rnjbL4ngENNnxUUFBQcUIYWAAAHt3SDQHILmntkdLOQZNAkUzBl2XF5A7scFwOxrc/EDsxkTFF8dzNFIwU9ldR7WayD9EjO8t1swwrg0Ia1EmhtsQwSG5O6KL9DJ5oq2zAcej0zCVNLfC77mdvUD8vjXNbZye6JzUvdsNSxmqQX1uz7a5aM3PpfaksKJpHAFKWsiH1Z7D4FWemt9MtSp2gzRStgg0tW0l80xwrUYV/jZohBYrNSB+t1Ws9Z+/pqI4m13ikzvmOJ6HYpkixs2yOwZ3J7BQWFnQ7OWUqfVMAYO4YxtooxtpYxdn2S389mjH1Pn/mMsfGO3zYyxpYzxpYxxpbs6HOp6bOCgkLaENbnruFUjDEdwMMAjoZIWbqYMTaHc/6j47QNAA7jnNcyxo4FMAvA/o7fD+ecV3VFfdLqFE1whLg9FXQuapcGFDmlaGvq55xeWOfIcqSynqYYmulQ8Guu6bM8B+4pmX2NnD5FyOXDsOIGduySo7tccuS0K2H67FLatzKwuFw4gNZuHEGaLntZ66lYW5AqDOkobzjmMdYxt7ohyfIy9zH3vdtTfdjl0HTa5V4FtO1iZTtzJ6o+xLNR/E1qF+mSY6Rg6ZSykeXbRhVHvMM2ZCZVHNKoEnTIOVNLNIKlYwCT02U5RXbLTjybU349Z07ahdPnKQDWcs7XAwBj7GUAJwGwOkXO+XzH+Qsg8jvvFKjps4KCQqeQxvS5iDG2xPG5zFVUPwBbHPuldKwtXAzgfWdVAHzIGFuapOy0kRZTNDhHvek0tNiQYaz8VpCBtpXREm6Fvu24LA0tzjBjUrEvRnLpNOxebuY0yMiABj5ikVFyAZFWs2ThrqTbSCuXHJa4PEzUty3WIQ0wiU6+4rfkDtntMUTJNuKQhiO5Jedn7mQdAm51vnxS510CloO9NDzRsrh2GJC7nkGrXMmGHXeWDJkOyRBkFrvn0nhlM2npyO2lZ3S75KQiMylvyyXHEWQjYLnXJA/qIBmiZIdA5xiiLSPakmQiJKtoAru3YfQQ6wVH6vpCAFWc88nt/J6soKQNwRg7HKJTPMRx+GDOeRljrBeAjxhjP3HOv0y1cm4opqigoNAp8BQ/KaAUwADHfn8AZe6TGGP7AHgSwEmc82qrHpyX0XY7gLcgpuOdRno5WqCh0vBZ+15HwAZ3cM2g1MW5OEunGCPgcP1wZhmB1a1rVqgyh0sOMRGZFyZA+q6oS6doJgRUcAWwcDmHB7Rk+ThkMIfkDNG5DKwzDDFCzxGxnkccD1F7hR3PE2tDB2c5TCeRmSUrGR6N6u2n1yMVxmjVGTHHntQpJjJGt+7YGVRCd+neJDOUiwCkrLRWXBjwWcv8EvXAiUzRLbPkDLGzTvSy3hF618K0HyIGGKbnkDmLxDV2ebGewlM4wLtumd9iACMYY0MAbAVwJoCznCcwxgYCeBPAuZzz1Y7jmQA0znkjfZ8O4E87UhllfVZQUOgUumpFC+c8zhi7EsAHEBqepznnKxhjV9DvjwG4BUAhgEeYIBlxmpL3BvAWHfMAeJFz/p8dqU9anWKca6gwchysw7HQ3sWWYjQix+gcg6zHAcctnY7KQHuMEZCs0bCsgomM0dIlOh2+XUwkJpeVWSwkifO5yxFYMhJfMv2UK3NbW0FhA0mce9tkWg7WJJlhiFiHZBuNpsy2J55HZtQDbBYcc9mMJYP2OayqUp8Wpvpna7JdYgnX+h0ya8+bAHDPDKSlNTljdFujAdvLIGo5bUum2LYVunVoMimH1h4D7syI7kAcbv0h0DFDjDgs7mEpM1MyRHG8kXSlzTRrcWa1lO+lAYY43570Ht0RXan+5JzPBTDXdewxx/dLAFyS5Lr1AMa7j+8IFFNUUFBIG2rtMyHKPdgaK7D0bE5LrLS8hmlBfRRiuZQhdXA0cuoOq15H/nHOBFnSii1ZmAx3FXAxRsDW+YWlToyYRFRab1PwBXSzYck2fI77uXP7tkobwFo/X1vBRd36Q+d3yRDrKF9yI+VTrjOF7bfZdOZa9lN57uWWtKzQmRiL8jrHtJC4RradJgN1OHS6hLb0jPK5nM8q/e9s9kiyY26ZOcqRbM9amkksKoWETm42bDNG54wmUfctk4BZlnjLz9bhz9kGQ3RbmAGbITZSh1FHspJ5se282bbMwpYfrWZ5SHR7cACqU1RQUFCw0UO8h9KG6hQVFBQ6AdaV1uduhTQNLTq2x3KsqCTO6XMW5R4JazQV0OUSMZnnQ0zNEnK0dBCVODGii8x9QYpwSzmfOJ1OVMQnTv+sqC2uSCXOqaYVfce1zE9OmwMJOVRc+ZjpuJWfg6ZiHofRo62pmOWY7ZieN1M9pWGljqZgdYaYNtcYWeJ3I8O6JkTTtbanz7Z6QeYhkbIyEKJnF8c1Tbq42HXSraWZyaMhOR3TZJZA05KRVC+0k9dFXs/k/eS1MkeLlvA8TtjL/KRLjlR9mI5zpNsR6BxXtj0r703bf3jDMtxJY4pj+uwyrLhl1UDT5yYjYF0jp88mWEpR4bsNFFNUUFBQIHBlaAEgXHKqYlnWSOzMfWGNdh4ZaTsxlp7XUqC3Zlodsg/AGpUMy43DFcXZxRgBOBYPiXvKYAOyToZl6Gmt+Hc7bfuQ+BwA4LOeDbRNzjrag2QdcnlXxMFiQ/TSNbqU9ZJ11MQlU7RZR5PhT3hWCfkcWbptaHGfo7tYsZcJ9u9z1Mndzu05MrvZveFijHJZoZmE3UtXqJjL8NGerGzjWGI2Py9a11+3rk0MipGKg7ab1YcdyjXpJlVHjFAawyzGSDKTcgKAiCldcjTEedv37XZQTFFBQUHBCcUUYXAN9bEMS6co8x0DNutyRz+2nboTnWYBIGAtzet69gE4cwbTkCbvLRkjWjtxu3WJ9n7rZXLuIAtu1qFZbdL281ghpaguzkABcilYM7nOSN2h3NbLbdzWKbYYgrFHzUTR+rREtizuLZltoqysiOLWskib3cckC5OMyvU8zmeVLlVaG+xenulsUykLGcxDykzuW78noSnu0GFy60vQKcpt17H6sGMaGZJO2uRyY8sqSFtyyYnbTFHKzARrpQvu1mhHNdyToZiigoJC+lB+igIGZ2iO+yzWlCxvhrRISyfisGtZU8wRpNaQATdZYiDOZPY3t0N3R+wDsBmIxQQty2YiQ0w2NrsDKHhd9xG/ybpJPaorC14KFkzJbGOWRdM+RwZ6kG0n2Ye0MIcoOEcy1hE1E1tRZskzHUxEMquQ5ksoP9uSmXSgdjBFWW/ZtrztPDKtsgK6ggpb1ukEJ32iH7LpZH2tZZet4WaGblavJZwrt51n9e7gvk727X7fQy6ZSVk1OwKrWEyRM6s9ewKUn6KCgoKCE6pTFIElw4YXniQWQA+NytIKLfVhcj8q8/kmhEmSjG3HW1eyESeTk6zLOkaDsMVuXD5wzmPJdFZAcqZi6xbTH+Ul25DqmZijjJhrqZsVqp/aVFotIw5WGDbEsbiLKUr9qubQr8nrZXlt3c9ZJ3m1rLcnjUfW22Dm7WnR3HJoT6doBwZO/E1nzu87zsQMSw8s4Hynoy4Zuf8Hss0lOwSAqCGuiXMdvCcZL3oQq00HiikqKCh0CkwxReGsGTc1O7CrabeK9K+Km1L3J0PcJ26ThX6SOiXTYjFte/W79XRSj5dMwWEzweSWy/ZYh7sMPckbsDPWHjh1StKiH22DwcnVD05WKL+38ncz0fpcVz5sWb68nzvAq0DX/RNkTZyBymQ7W+9EG8Fy2zumdcD2Acd7Y13TMeux/Cp5cvYqznG9967/hft/Io7p1rEe4xDNGaCW+SkoKCg4oJiigoKCggOqUxRtYHJmO0U7p3qubGtWpOR2IiYbrsX/XYEddX11u3Wkdo17KpZ6LVrnbm7HjQeJ7W7nLXFO31rLBgBM1vpcwyXHVHIrG61UHR3Dbo/EZZ2xFP5VlhysZZ7pvyw7yx3aSFJ9d+AK+/135ZpxyMf5f+pR/UyPqmzq6EHu8woKCt0G0nk7lU8KYIwdwxhbxRhbyxi7PsnvjDH2IP3+PWNsUqrXpotOdYpOtri7oDPWJe4VeyqkjLqDrLozusN71FNlxHhqnw7LYUwH8DCAYwGMBTCTMTbWddqxAEbQ5zIAj6ZxbVpQTFFBQaFz4Cl+OsYUAGs55+s551EALwM4yXXOSQD+zQUWAMhjjPVJ8dq0oDpFBQWFTqGrmCKAfgC2OPZL6Vgq56RybVpQ1mcFBYXOIfUpfxFjbIljfxbnfJZjP1lB7u60rXNSuTYtdHmnaCaxiCrsHJg7GGZqR69X6BhteQP0eKQ+NQaAKkpc3xZKAQxw7PcHUJbiOb4Urk0L6l+hoKDQOXSdTnExgBGMsSGMMR+AMwHMcZ0zB8B5ZIU+AEA957w8xWvTQtpMkXNmBQVwMo22WIfh6nc742fWVZA+iGYby/1SKiOdAAgphKFqfdzpv6Yl3brbNJnPm3u5WDLfUgkDye9jL9Vs+6FTC/fWdp5ncZ79PZYi+0hlud+uQlIf3Db4hrttxffkMuvuSMNVtV1wzuOMsSsBfADxKj3NOV/BGLuCfn8MwFwAxwFYCyAE4ML2rt2R+iidooKCQufQhWMQ53wuRMfnPPaY4zsH8OtUr90RqE5RQUEhbaRhWe5xSDueYpxrVhQSZyQWK/qHK/KKnBrICCwJ0wbauuMpyqjU7U0/3Wg3zqEVeUXGFEzcbzeaSjsRV9xOv+ks73PDsLbOKDmupZO0327EFddvEowiGiVGZ0m8vq37Jdap87DbR8pXvht227YV2cj9e0K5rn+nlHd770Q6MNx5wi01jIDznbbfc/lsrohGvLV8DGtKzVQ8xW4AxRQVFBQ6B8UURcjCSNyTNFiCjLzdolPEYTMxV4WM0Rd1qOStAAdyfweiOcesfRsy3JtU3mtuKbZzn7biKCZlHzuwTKw163BG3nZHbU7cyujNYUcU5whF3o4ZrsjbeiJLBgCfIXhfRE8s147Ench6nPXckWjp7hiYSfOuuOIqpiIrt3ylASepEagLoqTL99f5TkddeXXCblnFqa0N+68Xidvt3JPynqjps4KCgoIE7zrrc3dDWp2iyRlaYh5orlwngD1K+3TKWqZTFjPXVo6cABCmvC0yr3BXZPVzBgOWI7vMGmdniOv4Wd26RG8S1uHOzdIVWfxkTg/xPTEjXKMRAGBnggvFaRuzM8OF45IpJvIvqTd0un14NHHvZsrfLcvP1d2ZGO06xSgX9I5k9bN0u/S7s029lowSdYvtQXOxei8SGahT99sVWfxkbm53/hXAkb3PoKx9MnsfyaqFmGIo5sjRYjHFHubkrZiigoKCggOqUwRMU0Mo7IdGDCMat8d4OcIxGq19mtBXZegiz3O2HgYAhPSQdU2Yi+8yr7BkH940crVIFuC1gnja9lF5teliH7BYX9tSdeeR7irWIRGjesYs1iFq28xt1tdAzK3eyAAANMZpPyq2DRGxbYzY10SIgZiu/BkRTdwn5m1dN5mdMYMYY70u7pdL98/WWqxzsylvd9Cqv0FltC2rzrB7WUufO3d3O/BajJHuS8cTc4FLGXXM6iXaZvXiDiFu591uJFnVxoMAgLqY2G+QMosSk3TILE6s3jQ1mGZn7OO7B3uqTrHnSEBBQUFhFyA967PJEGn2AcQ6NI+taZWjneEa6fyaYIFBYoxBLWL9FmTieyYT5wSIdXgl+3AwRTf7aot9eBOGL6qfPMfNGNuBZBs2Y5GM1GYWnWEdcUiGRXmyqSqNpmAOdUamdW5VPEdsY1liGxG/1YTFtjYkWEg4Yuun4jHRZtxIrAvTJVN0ZPNzWailjtEnZaYJmWXSFgAyGcmRCfYYsNi9TDXQ9jjbFruH0w/QqnaipVdLQWZuva9X3sfx7qTK6p2+iZINh+lYiBhinSHYYHU8yzq3ir7XRIWMqkhWNS3i3OawkHPEITMzTnUwGXhPypC3hzJFpVNUUFBIH8r6TDAY0OCx6JPhsYeKkJ9WshBTMcxEvzgP6Rj9xAoBIJNYo9wGmNAxSqYo2QcAa1Ryj+ySfXilTss5elmDLq2goD171UTbcLMNHa1Zh7ynZLRtsY6YQ88pvzcT66gny2WdKZhEZTzbOndbTDDFijBtm8W2plmc2xISrMNscYgxRrpdF1M0iSlGvfa58QwpK3tFBWDLzGL5DnZvyypObSD0kFJmzlU9Xpee0WofKaOkkfBMqgM9jpvlJ4ElV4uxy/uLb17HjEPWqS1WLxli3PF2hOnZQlTvOrIw1xiCFW6L5VrnVkSEjMpbxLayWTDFxmahU4y1EEOM2HVicaqLCfEf6ylQTFFBQUFBgGHPNbSoTlFBQaFzUJ0iwAzAX6ODk6HF4bMKwy+mKjGaFtSTEt/tYuBcNhfQaOrFEqfLGiK0b0+1g9Ko4aqTPSWTCn/HGdaUWy7OT5yaOarfCpa7Tavps11+qsp65/Q5RN8bqV2qadq8NZYPACiNFljnbgnRscY8cW6jmIqFG8X0jTXT9D3kMCTQekf39JnrSWQWFNc1R2kabbim0Zbxyq6/LSshGz9rBAD4WOvps21Yad9I5k0yY3Q7z6fixO1Wcci6JMqsbSdzwHbUDnP73WukYBrVNG2uiOcBAMqkzCL51rlbmsX38gYxfW5sEsYwo1E0vEay0iOOhQ9SZiYD25GIG7sSKkqOgoKCggvK0CJGtEAlALlcy8E64gEKa5VJrCwiRsT6mHSqbh35WbIBjUmWJ0Zn3VNH+7YrCMgtJAgf/daWEt+WlDzHdDlZmylI02aBbgdth/N2BwwxQmwjxGPWb3V0622kpJcMcXO0EACwvrnIOndTo/itslYYX2L1gql4GoghNom6eGzfauhkE3EzDk4+RobtZ2wZWuIt4jnCUbHdTi4iyYiA7paV5TAtGKOmxVpd46fXzN1eUj4Jy0VZoowM1xK79rAjMotZxhTxnjWadgPWkLtURVwYVLbEBJtf31Isto2F1rllxBCb6sQMAPXiT+JrJFk10zbsqLdsMu743gOgmKKCgoKCE6pTBAoyMnDSkOFoDEWwrboRm7bVwpDBS4kpxrIonBKNiFFyVm0kxrjFEaiAd5D5T0etY4+G0FQZI2z9o+FyNZG/tMc+krEMd/lutMUQaxy3qSDn7M3ENtZFegMAVjWJ7fo6m3VU1wo2yWvEs/prRb19DeJ3byPp5EKScTP0zc5Gn145KMwOgoGBc47KmiZsrapHZW0TDL8jNFmQZJZNsiLXnmiUUV1JT+xw8m5LVprltNbsOJhIe9pijAn7Uu9LMrICArfjGO8O7tuR3tCJthhilSMc21bSIa6P9gIArA6VAADWNQhWX1pnu+SEagRD9NSIZ/XVibr4BJGGt4lc1MIcXo+Owf0K0KswE1mZfmRnBHDb8x1WuXuAY5d0ioyxAgCvABgMYCOAGZzzWtc5AwD8G0AJxKR+Fuf8AfrtNgCXAqik02+g9AVtIq1OsaQ4B3+65gRrPxyNYdWG7fj6v+vx2Q/rsKGsJp3iFLoABdlBHHvwSOy/1yCMH9kPOZmBNs+tqm3Ct6tK8c33G/D5krWoQw+aq+0hGDe0BIdNGIaDxw3GsAFF8HgSB/Xb/rCbKtYJ7KLp8/UAPuGc380Yu5723a0UB/B7zvm3jLFsAEsZYx9xzn+k3//OOb8v1Rum1SnW14XwxccrkJWTgb7989Gnbz7Gj+qH8aP64VeYihU/bsVbH3+Hj+etRJPUobSIbSQibhWK2cvYNsYTQ7Qb7YRNMlAHADBBekZijAF6hPYcqPU2LNft3S+dVAgxV3CEJmKIUn9Y4Vi6tzEq2MXqcB8AwE+NgiFuqBPMsbbadt7WqgRbCVSLugRqyKm6jmPyuIE4e8Yk7DdpCDy6XdfqykaUba1FTVUjTMOEx6OjqFcO+g0sQFF+FqYfMBrTDxiN68+P4ctv1uDlNxdhRWk1ACCSRwwrLFoqQsELquJ2+dIp30phkO1qJ59zR7BG08UYpTO1m+UDSVikvG8XycrtnC2tzPXEECspLJtkh4DNEFc2C5mtbRC6xNJqcU602h6I/FWixv4aUV9/LUdmwIdTJ4/BST+bgMH97ZmAaXJs3lSFbWV1aGoMo6kxDC3aU8zP2FXT55MATKPvzwL4HK5OkVKdltP3RsbYSgD9APyITiCtTnFbeR3uvOF1cPJpyczPxPiJg3DgISNwyGGjMG5sP4wb2w+XzDwYT723CO/MW4HInqp42E04ZK8h+PWxB2H0UNGZxuMG5n+9Gl99tQrL/rsJVVvEzILFE1UD3KNhwKBCTDhgGA6bNgbjxw/E0dPG4uhpY7Hw+4147JV5WFazfZc/z56MoN+L84+fgLN/Nhm5WaLjrKpuxBdfrcLCL1bjxx9KEQ7HwGKiI5R5dHoKdtEyv97U6YFzXs4Y69VunRgbDGAigIWOw1cyxs4DsASCUdYmu9Yqg6cR/zzX14sfVDTDsj7DY/epvvxMTPv5BJxy0aEYPEzUe8uWajzw/BeY/+16tBSK0byl2B7xw71o+V2xMMX1KawHAIzKE3/OkZkV1rkj/NsAAH094nmKablZnibK9TNRFyf7SMZEdhTOQAHSpy1CzLDRFKyjxhT3rSAL88ZosXXN6rDQR61qIIZYKxhiQ6U411tpt2lGJaOtiZGDeuH/ZkzDhL0HAABqa5ow54UFeO/VRWh0qi2ixMoMF+PQqS18gn326pePky6dhmNPm4yMoGCEn837CY88+Tk2mYLhtRRLmdnvSKxYPGNOcRMAYEi+uPeoHCGfkQFbZoN9Qo1TootzC2ipZ7YmntHPyHcvBevwjiDZMkup960zhTwryQexLC4s/mtI1wsAq5tJZnXivS6vFjpEo1J0dIHtdp0zKjk0xnDaPmNwxcypKMgTs4Tl/92EN19YgEXvL4MRN8HCDs+KOPlEco75Va+iPrq926/1y+g9gA8/+5qUzv3h79dsAlDlODSLcz5L7jDGPobQB7pxI4BnOed5jnNrOef5Sc4FYywLwBcA7uScv0nHetO9OYA7APThnF/UXn27zPocaYnhg9cX44MPV2DqkWNx3q+PwIABhbjvj6fi66XrcMfrn2BbbVNX3e5/BkG/F78760DMmD4RuqahvqEFLz7zFea8tRS8SnRenfkHbd9ai1n3vo8XH/sMp11xOE6dMQWHHzIaB0weiife/gYvffBt1z7I/whG9ivCLedOx7jB4j++YsVWPP30F1j++WoAgBbfM5z7GNJ676o455Pb+pFzflSb92FsG2OsD7HEPgCSTmcYY14AbwB4QXaIVPY2xzlPAHi3o8p2uUuOaXJ88dEKfPHtepx88r449/xDcPC+w/D62P544M2v8NxPy7v6lnss9h8xAH+aMR19C3JgmCZefXsJ/vXCPMTKxODSFcJragjjqcc+w9uvL8Hlv/8Zjjh0NH5z5mE4ev/RuOH1D7C2oroL7rLnw6tr+NW0A3DJ1Mnw6jq2VTXgkee/xFdzxPu+R/q+7ZrZ/hwA5wO4m7az3ScwxhiApwCs5Jzf7/qtj5x+AzgFwA8d3TCt6XMOK+D769PB5Do53Z6eMppKswwxpWDZYjqYO7IPrrruOBx8+BgAwFfL1uNPT36AusYWhEpEOS29RB14iZhG9yoSPicj8iqt8kdnig5/qF8MFAO84s9arAvP5Tyawfgd068ATand7jWdUczLqbJzKianYI10Tg25cVQYwoFXGlXWtNhTsZUNgkFsqhEzgFClmF4FKkRdM7Zx+Dw6rj7lUMycNhEAsHpVOf5273vYtGA9AIA1CIbIm2kbtqPYcEPGkHQxEnpm5jDKsAAtF8wUdeA5mdj3kBG46rZTUNInD7G4gVkvzcOTC5daIShbeos2DJdQBJ1iUYdBBUKtMSbHofLIEDIb7BMzpxJdyLWAInxns0TVB9A6is2OyEq6XDmX7EXoHGkEq6RI2Vti5DwfEVPkn5ptma2pE+qP7VVCrqxCvOMZNMsdzXPx51/9HGMG94Zpcrz97rd49m8foCUUhV4jBjDeSNsW8Y7zuF0nqergJsdC40M08JpuP30O9h7AR8xMbfr8/QPXLG2PKbYHxlghgFcBDASwGcAZnPMaxlhfAE9yzo9jjB0C4CsAy2Gvs7mBcz6XMfYcgAkQXfhGAJc7Osmk2OkDWE1VI26/7hUcdOpEXHPldEydMBQv/flc3PTo+/iqrnRn377HYVCvfNx78c8xqn8vxOIG/vXaN3h91lcwDLPdtdpdhaXz1uCSC57A5b86EiecNAm/PvcwTNh3EG569n3UNrV0XMD/GE6cNAa3nHwkMvxebK2sx933vIvvV5TCH4p2fHFPxy5gipzzagBHJjleBuA4+j4PbczmOefnpnvP9DtFboLHqTUcox2PEFsJiZiIWlMzbcX+gn8vwK8+XIn/u38mJkwchIevOw1Pvj4fz769EBotL2uhnCAV5ALijAwdNUVVIzKDmsw0R4wxrIv75Wl2nWLWkrTEuHru3CpOuB26bXcbyTrs3xvJTaSG2EaFIRTwkiGuDZFjdoNtMNsi3Ti2CSffjApRp8xyjiMmjcDtF/wMwaAPpWW1uOvaV7H2p3LoVfXQAfB6wbQMauMEtpEiuNM7JixYC5OyahR1Mppa8MhvnsPCT1fi2ttPwcFjB+PV687BzXe9jcWbBXtnMSGHFooKvl46epu2zCKu3N9hn9wXBrUCytuT7VjOGZCyShILEWhfZjK/S8wV3TzkmA3V0XtUSW5SkiFuiAg2KJ3o19Xayy2rKoWblHebcNfJKGfweXTceuzhOPGIvQEAH3/0Ax78+weIbqhCBgDUiGc0G4XXtin/H+3NzHYgf/guxx4cZHaX5mipqqjHtde8iOeenQfGGC7/xSH4y9UnIujfFRyo+4Ix4FenHIy//uoEBIM+fPLlSlzy22ex9qd2Wf5Ox6J5a/DLmY/huxVbUFyYjQfvnomTp+61W+vUHdArNxPPXDkDJx6xNyLRGO675z385c9z0Nwc6fjiPQk8xU8PQ3qhwzwe6AVFgFwsH3MwxagY7eWIaBILkS4iGu17whG8eNubWP3Zj/jDP8/DYVNG4KWSfPzh7reweovQu8h8xlVxe/mUzBwYlnlzs8U54YDYb/YKnVYvvdG6JteK7E3hriDDXiWGA3PCDkBA9+XS7UbsNzqiYMiwXzKU1KaIYBfrQoJ1yGVgW6vyrGvM7YINByvEeFRU5cEdVxyLwyYNh2GYePK+9/HWc/NFm1WKZzIaBEO02DjB0uP67SgPzEfe07prvCNdo5STszzJOOV9GB3XwhHUlG7HDWfPwmXXHouTzj4IN58/HWOy8/HXOfNgcg6N2HwLMa+NjgyPMes3khm5/jRTfZt5HQCgULMzPGaTo7efaIgdkoweqx2ZxSymKNBMDub1pu1cvd0QrE+G/VoXFix+TaPYul2kAMC3TbRzsJxh78El+MdvT0BRfhbKy2px+x9exaZ5qyyfdbNGPJPZTM8k/ysySEWAdO4+h5e7jIau6WA1Pccks6cGhNht2fyWfPETrrzoSWzeWImhA4vw+F1nYeyA3h1fuAehKCuIx2+YgcMmDUdDQwuu++MrVofYnWDETTz6l/dw/81vIhYzcPqp++He35wAv6/n/IG7AofvMwxPXnUGivKzsHT5Zlx54ZNYt2ZbxxfuqVBMEeABL+Ij+wHkea+32Aoq1kxMsJGWdjWQLqVFKOelbsXJVLbNj+HqY+7Fjc9cgUkHDMNTV52OW+6dg49LN4syHYvyG0wxcq+l0T9qShZCASfI6t3ozbCuKSZrZ44m6haUQW0hWUhrpYjMLxyjpYfufMwyLwcAbKdse5sjiWG/ZMiviirBdNk2m8lJhjjGyMWDl5yGvr1ysXVrDW656Cls3VgFVNouMEYD+XUS22BeUQctS7AylkMBI7KC1jVmkPK2uPI7azEKtOAwAEh9L6f7mKRblAzSqBaO2XpMXPPhw3NRtn47bn34PBw2aThmXX0afvfgbNQ3tYAZiXphANhiinaImfI3YowUiksyxkZPg3VNATl6y6yBASuorcyz01pmUZJVjPSPIWLzDcQQK8kbAADKYnkAgE0tQlbrm4TsNteKujZXirb0V9jvXrCc44yp++D6GUdA1zXMfWsp/nn3e0BpJTwAjBp7gYTF5iUzDIrytBzBUHk2Wfoda9SNDLqXxsCX9RxVkmKKAHJzgzjk0FGYuO9g9OufD9YFiuFQYxg3/fZFfPTuMgQzfLj7xlNw9MQRO1xud8bIkiLMuvlM9O2Vi59WluE3V/5bdIg7CMYYinplY9z4AZh62GhMnTYahxw6CqPH9kVeQWbHBaSAH5ZswNUzH8G2inqMHdsPj900A4V5XVN2d8WF0/fDjTOPgq5rePKFefjHn9+B0QVO2LpHw6BhvTBpvyE49PAxOO6ECTte2V0FDqFjSuXTw5CWn+LkyZP5kiVLrP1QOIof11fgk8Wr8fXnq1FT1wx/tRjhvZXEcqrFKGrWCWuc02JqMZ+CPADAxfechTPOOxiGYeKuB97Hm2tW2ffqk+gfl1EkWE6/fFHuwExxn/4Z9qjdx1sHAMjTxbnZ5NMoWYiWhH0YNE6EiW2EuGAzMsevM9veVgpDXxoS9d/cIParKKiDtp2slRX22DPZX4wHrzsNudkZWLxoPe68+ElEWmIwiSGaTY5VPzKTHfkRagWifKNIMNBocQYCfi8O/tkYHDRxKCaNG4D8HJs1ulFR3YAlyzfjqyXr8PV/18NDsvJVinbRq8hiSszHJD9Ip1VUyxLstGjcENz5zKUYPLIEpaU1+OW9r2N7TSOa+trP2lJCfoO9xH2KCsVsYWAOySpYBwDo57dlVuwR50iZydzgMnWFlJXhGM9NYoqSzTeSN4CUWTmxQwAobaGgvpQ2YGutaMuWKnGu9BcNlou2/+XxB+Gynx8A0+R48M538P5bS8HKq6idRP15zGbfUs+r5ZE+vJDYcrFot8JxvXDEQSNxyEEjMXJgMXzexMlasNcAhLZv6fZm6MziAXz0Kan5KX77ROf9FHcH0po+1zW34NNla5GT4ceAXnnonZ+NyWMHYvLYgTDOOQJfLlyD2S8uxPLvt3SqMk/+40OEW6I49/LDcePvjkPoX8B/lqzq+MIegtF9ivHQZacjOzOAr+etxh23vQXW0rnwXX365+OMK6bhiMNGIxi0p+e1jSGUbq9HZWMzTNOER9fQKz8bA4vyUFKYg+On7YXjp+2Fmrpm/OfD5XjljUVoqUzf/7BmewP+cPaj+PMLv8SIkSV49MYZuOyOl9GEPceX8fKfH4jLfn4A4oaJu/4xF1+/tbRT5Wgaw0GHjMSpp++HiZMGJ/y2ZVstymsb0dAcRmOkh/k27qHT57Q6xS3NDbjo67etpEjFLIipQwbjuNEjceiQwTj8oFE4/KBR+Or79fjnQx+jfFs9MsrEH1b3CuZlVNnBC+QIK0dcnTG8ePsbMAyOC351BO44/xgYW8P4Zul6cFrpwGn1QwvEyL6FQlm5LZ0A0EQrNgq9gvHkegT7kLmLZRImJ1OMkg4xTKxD6r+qYoL9VUZtneJWYogyBH1dLTG6SmKI2wSDySw3MbAkH7MuPhXZmQF8+eVPuOeKZ4C4CU7tkcAQCXoO6cJ6C/1XpF8ecnIycN5Vh+PEw/ayYvF9u6UMc75fifkbNmNtlMrzJ76xeoRhdHExpvUZgpP3GYsRvQpx1owDcNKJk/Dy3KV4ds4ieDaL9vIFRP31bYIRGY22RV/WkzGG+qoa/PGsR3HXvy7ByL0H4OHfn4bL7n0V9U3keUC+hi1km63iog3jpH8Mk6ycMqv3CbkWecU9gzJEHLF7aY12MsWY1PtSUIf6uCijmsLUlYdtL4ayZvG9ok7UJVIjWKWfLMwZFaLdLpg0HlccfyAMw8TdN76BeR+tgE5LHg0XQ5QzHgDQi4T12ugrdJWDDh2Ga644CmNGiLBjzZEoPl25Du+tX40Fm7egKRqFEaAVXTqHQyXb7cHSmGX2JOyQ+bCmpQWzf1yJ2T+uRF8zCzP33wfnHDARU/cZiin/vBCPPPMF3n9yXtrlvvjUl8gI+vCLCw7BHX84Eb+56RUsbO65Ya0KcoL45x9OR0FOEIsXr8edd85pFdorFUyZMhTXXfdzFBRkwTBNvLXwBzz9yRKsjjsiIWUlv5YDWFlZibUbqvHk/CXYp18Jrjxgfxw2biguPvVAHLH/SPzl9tlYuzY9a2pTQwtuuuQp3PvKrzFkSDHuv+YU/PIvryIa60FxAV04dvJo/PYisYji3oc/wJcfrUi7DF3XcNGvj8RpZx8ETWPYVtuI5z9aipd/WoGmSBSx3B7eofRQy3Iq6DKXnMrGZjz48Tc49ZZ/Ye7ClfD7vbj68qNw3e2nwN8J5+yn/vkx5nzwHfx+L+7648noldfGv72bw6vruOe3J6JPUQ6Wry3DrbcKt5Z0wBjD+b+bjrvv/gUKCrLw7cotmHnds7jl5Y+wsbLd0HBt4vutFbjyydk478FXsGFrNYb0K8TDD5+PY47ZO+2yGmpDuO66l7FtWz32HtEXN1x0dKfq1B0wblBv3HbudADAP5/6FHM/6TB+QCvkFWXj7ofPwxnnHgyTc7z01iKceuszeOGTb9HU06bI7YDx1D49DWkxRW8ghv4jt1suFo0tNtdvrhff11dyXPXV+/h09Xr8acZ0HHXcePQd1Rs3XPosQs2RhEVbciptTaPJGKOTsvqRP72DAYU5mDh5CB646ERccfNLoFV+4BSMIgJaYkdlxh15psOGOFlOo/PIXceakiVJnWaQ0l5Om5so/V1VRHTKlS0Ol5wm8b2+XkzXWDXlUqFI2RmVHDddcDjGj+yH7ZUNuPXGN2BuqhGTyeo6cT85NaWpiDVlBoCSYmgaw+/+eTamH7kX4oaJh9/+Go/9sBRmFreMTt5826k7N0vo9DK8iUsAW2hZnsxDDAANtX58gc04+ZkXcd3RU3H2lAm47rrjUViSi+denI+A5eTuaB9y8Jb11kkODSs8uO2Xz+L+F6/AcYeMw4YVFfjXwu8AAKaHpsuaaB9KW2Ll6HHmgGnOEO3dRO46WZSeUMpMT+JGZRnFSGZ1MfGMlWEhn+0h2zhW1SCm1HLa7KNI2RmVQEFWBv7x6xPh93rw7ptLMOeJeQgA0CvrxHNIY6Fr2iynzABQvN8I3PfY+ehVkouq2iZc/fx7WLaxDKGBlMuomBYU5Nrp/LIzxHevZqIq0HNSRKhlfmnig2WrcfaDL2HrtjqMHVaCW++ZAa83vaCvhmHijhteR3lZLcYML8FFZxy4k2q7c3D4pOE4+dC9EYnEcMMtb6Cmtrnji1y46rZTMP3IvRAKRfDbf76FZ/+z2Eoi31UIx+L409zPcOdzH8EwTVx83lSccep+aZezYVUF7rpHhKu77KJpGN6/qIMruhduPvMoFBdkY9nKUjx8/3/Svj63IAt3PXg2epXkYvnqMlxw/XNYtrFsJ9S0m0A5bwM53jCOLFmFiFxUH7VH4M25wvVgU6bYNvoz8V9U4qKn3sDzV56JifsNwe/uPhX3/eZl6xqNHLml86tcGsjI8VvP8CNU24z7/u9l/PX5K3DOSVPwxffr8cO6chh+yg9MxoaILphqtcOZQcpDBpNo9gn2kekRo7Vfax1QQbp3SOV/fUyUWxsRbLCm2XZ5aWgk1lVHuX3r7FwqBdlB3HiuiJ355AMfYfOiTfABYLXEtMhR2nK7kcu/yD0JAE679hgcO2MKwpEYrnrwLXwVLwMGAZH+ot2KewvmMjTPdvjuF6gDAGR5EpcENsXFs28N2+WvzxHGgMqAMD48u2k5aufEcd/Jx+LXlx2B7Vtr8dW81fBHbfYiZSZlJR2+9VrBmhbO+R7vjBuAE07fD3dedCwu/f1zMMjoIx3KIx5yrtZs1ioh8/VIlp/rFffJoHBjWhJ6It/HZnrGuqgotzYstjWNth9lhGY0njrKpVIrXphTxo7GEfsMR1NTGHffPhu8stFiyFwuRJBLV6VjNjnPm70L4PN7cNusC9F/YCHWbNiOS/79FppzomgaSAsF+oh2Gloo3vWBWQ43JEr159fiWO91JITuzuihU+NUsNOX+W2tbcAVT72F5lAERx8yBkcct0/aZaxYuhGvvroQuq7hmnMP3wm17HpccfyByMsNYvGyjZjz6uK0rx86pi8uvOBQmCbHjU/MxbK1u4ZxzPlhJR58/SsAwDW/Pw65ua07ro7wxN8/xNbN1Rg2uBgn/ix9ee9qZPi8uGrmoQCAhx/+GNu3N3RwRWucfelhGL13f5Rtq8fv73wDzXuQ7rBNKKYoXFn2z1xnhYKqjtv6tWKfWFDv9wj2tYr0RS2RTCwztuOe1z/Dn847BldceyzmrypFfX0LMsiVRaPF8xb7oECceqM93Xzxn59g+lF7Ya9hfTB93HC8t3EdACCeQbmj/cQYHe4R9Xoiq5D6xjA5zPq01gYPm6mIZ2ykjHYNYQpiELJ9As0mcY63SVwjc/sODebilIP2gmGYePgfH4I1hqwwBjwkdH6Ww68m6s1o6Z5RmA1NY7jq3jOh6xpeef9bzC1fD/RmiJQIttS7pA4AsHehiKIzJtOOptOXAmNIR3UJ6dDcz2+nt8gmNrmchvxtlJ/k0f8uwb7jB+LgEYNw2dXT8fdrXrWukSHhZKAP+RzyubSmFkSbWvDE3e/itkfOx3mnH4i3F/6IlkgMccozbQRoWR6FEmt2yEnTyOGb3h/J8gPEFD1JmKJc8tkUo2WDJLP6kHjmSLP9TmhNlEmwSZTvbeS48LiJKMzNxIpVZfjk1SXwAGCOd89oSWRvGuk7Wa54f/tNGYQzzj8Epslx/b/nYrOvBTKusLdEvNsje4uQa3vliMFtWMD2pij0CDenAIvhRa1nRNphUExxhzFnwY9Y8NMm5OYGO6WvCodjeOHJLwAAZ3Xi+l2JM382CR6Pjg8/XI7Nm9MP53/QYaMxelw/bK9uxKxXv94JNewYt8/5BJFoDMdOHYuBQ4o7vsCFBZ/8iBUrtqKgIBPHTR27E2rYNfDoGs46YhIA4LF/f9mpMi449xDouobZ73yL7zfs3nBvuxLM5Cl9ehrSYoo+ZmCAp84KlpCttV69IC2A1blC91aWK0bVaK4XD3+xEAeMHoTp0/fGY29/A5PSPrIMmqJZoayIhTimIFo4io9eXYiLf3MUxo7si8HZeSjdXgeKUwo9QpbMsCPIKTGRkE4BFei4tHp66bjmGPIkm4yQTqs5Kp6nJSK28YjdZFqY8kmHZb5kDo+m4aj9RwEAZr+0EHooBtZij/5mNHFaxbwyjYNog3i2H8eeITr9Zz9eglpfDFHyPfbniYftn10HABgeFGxDpmgAgBJd6BmDLst6SBPMx2lxl07PtaSDq8sT22iuB+uMOry7YCVOO3QfHDNjPzzyyCeiDlRP5qXw+hEKpU/PpTmedc6rCzHu9lNxzAFjMPudZdDDPKG9DGq/eMBu0xZp0dVI10qykXLxuJgkYFuvm2PkcE+yioSF/LnjnfBEpKzE/tThg5GXlYF167fjh0Ub4acse853T76PUpcoZWVmBZCbn4kD9x+GuGHiqXcXIlpEebFzxTV9c0U7Dc0SjvCSIcpMh4AdOs3LTMs5vdujh06NU8EuDR22dONWlFbUoldhNsaP7p/29dFIHPPmrwEATNt3eFdXr0swYVg/5OcGsWFLNdatruj4Ahfy8oLYd98hiEbjeG/Byp1Qw9Tx5lffAwCOPHJcp66f/+UqtLREsc/oft02aIR8jz79rFN50zH16HHweHQs+O8G1NSFOr5gDwIzU/vs0D0YK2CMfcQYW0PbttKbbmSMLWeMLWOMLUn3eifSYooaOILMsML8R2HrWuQC/lyPYI9ZXvLlCoiRzyTWtmD5Jpxeko8RI3pjjVf86Zk7IKp0OXHmLo6L78u/24zpR+2F4X0LocU5tDiN3tKQ7BCCSUsAJcuIu8KOSTAHU5R+ipJ9yCVpBm254QiOQNWTRmwtLiLgAMD3K0rBZBIpx3NYiaWse1N5xBgHjegFTWNYtaYCtfEI4GMw/eKaoE/cSFpkc2WgCwdjz7aWxSUyDsm8mrl9rrxeluen8pvIWvzDlko0hMLIz89EXu9s4VIkmS3VW7YcT/KskaYI1q4sx96TBmFo30KUbi9NaC9ZRdPRprKdZbtLOUjGKBkidzBFKU+3nKX8ne8Ec8gKAIb3I3kt2wxmcOs9S3j35PsoZUXvK/fqGD62HwDgmx82wvAx0ETJeu/l/0D+L6zgJJr93wlSZbzg0HoS/do1Vb0ewCec87sZY9fT/h/aOPdwzrk73FQ61wPYDUFm120RdR4+oHM+bFu2CIfvQX077PB3CwaXiHpt2VrTwZnJMXAAxffb0j1Si66vEM8h65UutmwW8h7Yr6CDM3cPBvYheW3pnLwGDxMGRvle/y9hF61oOQnAs/T9WQAn7+zr0wsyC4YYGAwapaVPH2Drp2TiImnF5TTiy9E6EqEQUJpmj+BuZ2Q5Imus1bFoVIyqmqYBDJBVsKrS+hKLCcqt1sa+qErycyw26fSDZK23mkbW61jc8RyOscf5TE7QfT10TSQct5TUcgpiEAOSOt0YLe8xEuQgEz+5E3BpSc71JP5G5ctLmckRoZQTOiCYVFuO41qSZ2UM0SgF3fCypO0lznNeklw2kkFZjDHh1sllxZKU735fdKpvNBYnk2rb755dCLWBCXhotU4kHBPtJjPMSsYrZUX/C/k/cf535P8JTPzHegQ42k/C1XXoLVOScs7LGWO92jiPA/iQCeE/zjmfleb1FnZ5PPksuYwr1DnXgwyKLB1q6Z5+YCFS0GcGfR2c2cb19FxZmf4Oztw1yKKIOS2dbO9gd5dXOIqMgBeZmX6Ew+kvsWtuEu9xtiOS9v8K0tAXFjn1fABmOTotMMY+BlCS5Lob06jOwZzzMur0PmKM/cQ575Q7wS7vFMcOF89etr2+U9cPGyo6+orq9B1sdwXKqF7DBqXvxgIA5eWiXYYN63BA2+nwenQM7SOmzRXbOtfeQ4YLh73yTjhE7wqUVzWgMC8TQ4YUo7q6dfi2jlBWWoMJ+w3B2GElWLR8006oYfdEmn6KVe0FmeWcH9XmfRjbxhjrQyyvD4Ck4bIoDzQ459sZY28BmALgSwApXe9EWp1iHAw1RgBRWgDlzFdSQZn3yiNiK5fDsWZylg0BGmM4YO/BAIBF89dBk+4PchmZpOPk0AyPXT1OhprxFKRz8ZpSxP22YttKsue1hy/dI6ZuftoGyLE8wyPu5yMFd6JLTqIRJkZL06Jxmmp6HUYTH02F5DI2H8PC9SLA7sS9B8Kk33VHdGVGARQsA4WMRB4WjGPVt5vQ0hLFsKG90Nebhcq6JniayRGenJC3hYTTcDktz8vVW6+pzuGJDscyX0mFI0NiOfn6yPJk+V663359+yLD58W6tdvQUCZyT8t6unNOy+eC41mzirIwbGQJItE4vl9fDrNAOtpT/EAfTYkdbeql7wEKaOHXE2XmIeuWmWBo8SQck/nCIyR3w7Hm3qTv8r1ZsqYUew3vgwmTh2Dh8s3We8Yc756MCylz5XArQ2UUSz//Ccedsi8O2nswXnjhG3hDpD6h917+D8qzRFv38onBwekaFYXojH0wEO8x0+d2VCldizkAzgdwN21nu09gjGUC0DjnjfR9OoA/pXq9G7vU0HL0+BHIyw2itLwWmzthiCgoysKUfYfANDkWrdy8E2q441hfUYOKqgYUF2Rh3/2Gpn19LGZgyaL1AIBfHD2xq6uXFmb+TDg1LyA3qHRx1M/2hqYxLFtZaukmuxsW/LARAHDUtDHQ3V4QKeC/C9YiEolhn7H9MWxw52YHPRW7yNByN4CjGWNrABxN+2CM9WWMzaVzegOYxxj7DsAiAO9xzv/T3vXtIS2mGDa9+CnaxwrVVOXIV7IhJKzJP9WKaV9dtWCRAQrNFKwGfnn4/gCA156bD39NBKyBMv+FE/WLmhytM+x1t0aWD8efdxC8Xh2fLViNsmgzkMkQDxJLC1K+laD958vMEOXmBMiFxSfcIoLEOvzJmCIpwAMUtEDXEhUn0mUEAFqiZFTKIGU6EefXPv8OV50+FaedcyAWrtwCT9DWNzHKOuhmXJyWOnrqWvDqI59i6mGj8YujJuCt1xYjXC3OjWULPePmgLCYyiWVUnkPADV+UQmnywcANBJTlHllAGBVo5jabq6hY1WifH81w/jBfTB1wjCEWqKY8+SX8NS1JNTTYopymSI9F6dn1XQNp/xiCgDglS+WIZapWe0TJ7HyDGKFGTZryqEwWrl+sc3xJQaEkMv8nExRBoRwM3+DXHKcMjMjJLMWcWxBaSnWVVRjWEkhDj1qDL5+5VtxH8e7p/nke0pMkWSnNTSjpaEZ787+FqfN2B8Xn7o/VvzrfQBANEfcpy5TPPRPtAxWGoyaHO9EEeWlCWgxhM0eFEx5FxBFznk1gCOTHC8DcBx9Xw9gfDrXt4ddxhQvOml/DOtfhG0V9fhg7ndpX19UnI1TZ4hO9eV3l3Rw9u7FW18sR0skhin7DsG+Ewenff2aFVsx76MfEPB7ccNvj4PWBVkT00HQ78UdM0Wg1dffWYr6mvRDnp1w5v7oU5KH0q01+Oq79V1dxS4D58Bzn4mO8NJfHAxfJ3JZv/LCN4hEYpg6bQyOmjKyq6vYbaGCzAJoNAL4tHa0pcOpDdthtMobBWtsrBQjo79cnJO5leOgcYNx6UkHwjQ5/v6HV4CyOnjQOminOziCmWvrLC+/9lhkBH34ZOkaLN5WgUgJLevLIfeLHFFGdqbNkAozBaspDIg/dYFP7GeQg3MypihZlwwdJhlKsuAR0jk4Qg7kUQrkWtkSweMfLcTvjj8Ev/vNdFy2aANi5JriDYln0ijIgBkKJWw1Cj776HUvYq9Pb8SUiYNx5Zr98PQb34BTmLRmiPZZGRP71QW2HHoHBWPP8rpCh1GwhG2OgKvba4Qu0dwmWEtmqRgj7zj5SAzpVYD1a7bhtX9+agXEddZTQqMAvixbPFcsP4jCwiyc9xsRffuBF79EOIfah2Rl5JC+MFvUMT/LLrM4Q8gq309Ozp5Epijz6iRjinKJoM+1fNPp6N1AjuIxag8twvD6Tz/i7NKJGNG/CGdcfhiefn4etCb73WOUu9wdBEO+vw0ry/Dk3e/i17eegpsvPBpbvy3HD1spz7ku3qMKQ/hpNtMSxNLsPKv8/IB4Vp8WR6PRQxK1cQBGD+zxUsBOZ4qTRvTDvZcfD01jeOa5eVj2zdq0yzjqhAmYdsgohFqiuO+Vz3ZCLbsez36+FOu3VGFAn3xccvkRaV9fXVGPe255E6bJcemMg3HGMTtfv8gY8MeTpmH6wWPQ3BLFn//4GiKR9HSBjAHX3nA8MjP9+HLJWsxbum4n1bbrEDdN3PX8xwCAc87YH3uN6Zt2Ge+88A0+f3cZgkE//nb3LzCgOK+La9n9oJgigOaoHws2DoFJTr5GyL5cr5e6Q5nBjuPQCUNx56XHISPgw8dvf4vX7p4LbLdXalisQ4a9J4YICrQa7ZWJwYOLcNXNJwIA7n/pc2zSm4EihmiB0C3xfDFq5+WKskqy7cxzfTKEpa/IL6x7uRROK0gh7gMsSToCuNMREIvSKZucbncSOklcaoEiXOihtLiOCAzc/PKHeObqX+C0GVOwdEMZvly4BlpEsDOd/Bk5LSfjFAxDZjYEgGVvLsLjeRn45R2n45oLj0Teqxl46tX50EkvFqLUAuW1tk9jeZZoO93r0oXGaPxrsvPleGvEsaxtDB5dw19OnY7ph4xBNBrH3Te+gfJvVkGHnb0OsHWJjMJnaRQ+yygS2zMuPxSTpwxDbX0Id7z1KUJFGiKFop1iBeQNkC/YXy8KllCSabvr9CJZFfrENosiN1jpCJJkV5eh7EKUOkKySykzr4PlS8fuOmL5UYOC424vxzNfLcUFU/fFrX88CZetfxoNDeJ9kUF2ddrKVAzy/dXpnf77ZY+jcPZ12Hu/IXj6mhm44ebXsHS9COXW3EJ1rBdW6NWOWZDUg2uaieZo9/BPTQl7aDa/ncIUPZqGK045CPddeRIyMnyY+/53uP+G16zVIqkiNzcDd9x2KgIBL96dtwKzv0g/idDuxA+btuH+t0XA1huvOhbDO2GdnPP0F7jvjtkwDBMXzTgIf7nuJBRkpR/4tT0MKs7DM7+dgelTx6A5FMGNv3sRC79O3+J8wCEjcfGZBwMA7njgPWyvT18XuTtx/3/m4fs1ZehdmI3bbjk57fQZkVAUN1/+LyxbsA6FhVl48O/n4ORD008E1lOwpzJFlk5HFeg7gA+67BprMb3uMHD6GkQ5U/v1x+8uOhJDBhTCMEw8948P8Mqjn8KsFKNpsvzGWpYYNRnlN471y4Pf78E9/zwHe4/oi1Ubt+PsJ15BSyyOlhKyPvYSLCA/X5Q3IEfod/oH66xyS/zimLTuZVOF7bzPraeGchmcxT4ol3ANBdStitkjfCmF9t/SJKy3FXWCLYWrRacVqBBM+oFjj8Ux+43G9rom/PKa57G9qhHBUsqfXC7WzJquXMIArMRIWkEeDpi+N6596AJkZgdQU9uMR5/8DO8sXw2Tc0TyHUv3qHqmTQhFGUSKvY7mz25iOOPoibj8lIMQCHixvbwOt1/zEtZ9/kNKdQIA3kfIbNDhw/HgXWci4Pfi0XfnY9Z7CxDqI9iYTLAVKBTMqyRPMMMBFJK/P6VQAIAiqmABBV4NaomsPlniKrlcUeZ9biR2L70jKiK2b2Yp5ere0iCO1daSjne7eK4hLTl4+aqzUJyTif98vRK3PTYXmaX0vmyl9qB82G29yx6vjisfvgTH/kIYBr/77ybc/9znWLupEi2F5Pea4wiCQYZo0wNsmnU/wmVbur2zYnZOfz75gKtSOvfzj65f2p7zdndDl61o2X/MQFx4zBRMGTMQAFC6pQb/+OtcrOiEpdnr1XHLradg7xF9UV7VgKvvexMtwfR0W90Jt/77AxTnZmHfkf3xt9tOx1U3voJoaXqrJxZ8uByXn/EQ/u9Pp2LClKG48drjMWNzJZ55ZyE+XLce0XhrQ1BbyAn68dDlJ2P8EFt39sHHP+Dx22ejuTGctvvwoKHFuPeW0xDwe/H2/B8w670FaZbQfVBR34QrnnoL/75iBo45eAzqGkN44p4P0iojHjPw4I2vY/mi9bjsphMxfuIgPL3Pufh8wWo89dkS/LS5B7ndtAEGWgu/B6JTnaLGGPrn52B0QRH2GzYA08YNRf8CMfI2hyJ47s2FmP3U14jFDKSb8dnj1fHHW0/BAQeOQF1jC66+9w1U1zUDwV0e0KfLEI0buPrx2Xjqml9gxMAiPPDnGbj+wqdRn2b8vcqKelx/+TOYds4BuOCcgzFiYDHu/PXxuK45jM+Wr8XStVuxrKoCZbUNaIY9iGT5fRhYlIvxA/vggAH98bNJo6zfGprDuPPud7Bw8Xr4G9NPmjRwSDHuefwC5OUG8c2S9bjzhY/TLqO74afySvzxgTm47/cn48xj9oVWH8Pjj32adjmfzf4W83/civMuOhQnnLwvjjx4NI48eDTWbK3CZyvX4b8by7CmvApl0aYeqZ5jPbHSKSCt6fP4CRP5N98sRDCjdbCD6qpGvP36Ysx9Zp5gG5S1TrotuDOhAYCeTe4hvSgs1aBC3HzrKTj4kJGob2zBRU++gZ/KRYRiOW3mvcVUprhQTIkH51IosaDY9vfVWuUXe0QdCnTByjKtqZiM/NL2s4dpSian0XKZXGXczsu8LSYGgk1h4W6xsVE8R2mdOB6qFIYjOY0uzAzi+bNOx9A+hVi3sRLX3PIawj9sAwBoFUK9YNTWWeVL44t0VdIpexzLy4XX58HRFxyKY06ehBGj+7Sqf3MoAsPk8OhaUnkBwOKv1+CO615FrJzyb5OsjAZisbSsTRpVAEDPzxM/lRRi0NBi/OWJC1GYn4lvftyIqx+bg9oim7HKaXOwWOgW++eJ8gdni2cdFBD37e2118FLmeWQ87mcNgeSqDokDOK2Ydc0Wi5DdcqsNEpZJ0Mks3qxrawW7yLbJq7NqNAwbcxQ/OPc4+H16Hjx/SV46u4PAQDe8jpRGBlYrNzdQKvsjFqeeBcKRvXHqecfjGNO3w+ZWYnBI0yTo7klgsamCEaO3gtNtd1/+pyT3Z/vN/nXKZ376ec37LnTZ69Ht/5g26sasWFDJVatrcA3i9Zh7eJN4BzwdIJtAEAw04+b7/kFJk4ajPqmFlx1z+v4Kbpnxairbg7hsn+8hlm/OwPDBhfj4Xtm4oZLn8G2srq0y4pF43jntcV457XFGDixP/bddwj23nsAhgwrRnFRNjKDdkfWEo5he1UDVq/bhhX/3YxFi9ajoqIennrR6aZnThAYs09/3PHA2cjODWLJsk24+uk5tJSv2/+fU8bnK9fj6uffxd/POR5nHTsZBfDib3+d2/GFSVC9vQFP/PV9PPGvr7DPxEGYfOhIjBjWC0MGFyMvN4jszEAPi7Szy9Y+73KkxRSLcvrzQ0dejJbGMDjnYBHbpUUufeLNrqV7km3QAnvpwiEKFKN05rg++MtdZ2DkyD6orG3CJc+8ibXbqhEqsZXqrER0tiUFgkkMzaWcF0GxHeQX2xKPzTryNBnlODEatVcGu07yjJLnyDuHyfE3RA7CdaZt+a0mJrIlRqwjLIwO65vEdmOtYCWN223jjL/Cg/xgBv51+qkYO7A3qmqa8H93vI6yr0WEFb3CZrqmzIftcpi2jB2ZlC/FsSQNAT8YY8jIC4IxBtMwEWqOgEmdo2NJJSenZLPZlWGQoAWFU7hWYC8NNErysd+Bw3HTPWcgI+DDZ9+txR+efg91xaL8SInN6LJ7CcY5OF88h8xTMjggtgO8gikW6rZ+NY+iiMto1AErnqJAsg5cvoUxepXD5IDfSG5Vdabt3C4DYmyKCBmto+Wp6+vFtoIc2nmF3UEdnT0E/zj3BAR9Xnzx7Vrcc/PbCIdj8JWRK1GVvY7frBfH3MsgLSf3TDF7YAF70GJBP4JZfmTnZeLNr/6O+nB5tx9ZcrL78SmTUmOKn3x5Y49iimkp6mJRA6GmSNquNe1h0IjeeOTh8zFyZB9s2VaLS//8MtZu6x5Rp3cWakMtuPSB17B49RYUFWThkbvOwpSDR3RZ+ZxzNDeG0dTQglBz16bMPPGM/fCn+2ciI+DDe58sx/89+U5aRp6eiK/XbMIlT76B+qYWHDZpOO7/29koKOi6nDOmYaKpvgXlnYz+vdsgI+V09OlhSM/QEo2Cb95qd4oJuUfouyufhcU2iCGaDk//iSeNxy1/OAHBoB/fry3DFa/MQW2gxdIfenrbDGlAUR0AYHiO0DHKTHaDfIkMsTAhXwkFHKC6eCGjUoutlmRMMIkjGqRvzKbnCZNuK9sR8kkGXZC6SrkvHY1lgIL1jqAStbrQXW3zAOfPfhN/PXw6fj55DG7/25l4+LkvMOcJO6WpVwbGqCW9LOmuJKMz6ojZNdhMywrh5Y7wTVG8uTP3iJnYmVlsnnS9LF+wqlhJLnRdwy9/czROPk1kGnzsPwvwyPvfINSPHLN7i3bJL7b1a0PzxeA2MkvISmay60e5qYt1cW6eZjNUSg2NgJQRbXWalrcnsxinLe3L3DXZpqN8kpW1JUd+uYxQBtnYouVZ1zTrQcw3tmLmE6/g8fNPwejRffDQ4xfiplvewJp12xDwOxYx+AQ7lYzRlGxcsn3at+QEAPSdMQa4sj12W/A91/q8W0y6jDGcdek03HXraQgG/fhg4U/45V9fR22odcrUPRlRw8Afn/sPHp47H7qm4TfnH44bbzsFgYx0bfY7FwWFWfjbP8/Fyafth2gsjtsffA+PvP/N7q7WLseGqlr84tGXsGzNVpQUZOOh+87G9CM6l+lwjwBP8dPDkF6OFs5tK7IbUnfiJ3ZDS/ZYjmAd8SKx9Q7Lx81X/xwHTh4K0+R4aPY8PPTtYqAECJfQMrBeYlQdVGhPJ0bmCJYxKijShsq8uSU6WZiJbWQ6GFKQicfzMqobJEPsWGVjkjRNchb2E2MMcFtn5mdRuo9gPpajsZYYmswZTGKtLsqr8oj2aNb9uH/FAqytqsGfZ0zH4UeNw6DRvXHDg+9g27cUEovCTOm1QnfIGwUzNCmohJWXGACPdTCVdVj/Ld0khf2SQR3MfMHqI70yMWH8QNx080kozMvE9ppGXP3Ce/huczla+lNwhxLyBigSrG94nm0cG5UlLOvDA2I7wCuYY7Gl6xVtkcnssdlPMtNZIjNMRWZ+K38LZT/kcqZgt0mQkScCk++LZI6JQUKcSwM36aLcFk8QLWjCuS+9gVuPPwKnT9oLN157PMbt1R8PPvkpYnEDgYB4/z2yTRuI3TeRrl0GlXAG6ZXh44AuVU3tbOypLjm7lClO2Hcw/vXA+Thw8lDUN7TgyofexFPvL9qVVei2+Oj7NZj54EvYsLUaQ/sX4Zk7zsHPj0saIm6XQPdouPCCqfjbX2eiMC8TS3/cjPNveh7fbSrfbXXqLogZBm6a/RH+9PxHiMbiOPnYCXj8r+dgyMDOZajssVA6RYBpGrRgppWxzRmynfnJFy4o2IyZI3SJkYIM+P0eXPSraZj5c2GA+m5DGa57di7WexuAAUC0lxg9s3qJ0XRIgWCIo7K3WeWPyBDfB3sFQ+xH/my5NKJnE1MNMLtOHrJV6iz9vl93fTOIfUjWCQBeLn3oRB18zA4rL47T7w49pAyvv4bCW1V4yOnd68dyVOGMJ17CTScegZMnjcX/XXMsphw5En/72weoqWuGn4LMeuoEC9cbKeCrk71LJuJ6Ga380j57as7Il45nC1nF8sS2ZHwf3Hj1cRgzog9Mk+ORTxfg0U8XwMjjaOkl2oGRv2jfQqHvHJEn5DI6s8Iqf7hfyKyfR+oQhXokl9i8lJWf2XWSjLArZCbl73UwRT/5O3qZ9EgQdfPCSDhXMkbADkW2wSO8DJo8ov2f27Qc3/57Gx488XiMGNoLT9x/Lh57eR5eef9b+LMojFkN6dQbKDQcqYh4xNYdWqzRNMFCPWSRAgeSxObYI7DTJTBh/EA8NetizPz5ZMQNEw/PnY8LHnwV5bWNHV/8P4hQNIYbXv8A173yPhpDYRw6fhiefegCHHfkXq0ybXY1vF4dM887GE/943yMGdEH5dUNuOJvr+Ghj7+BYfa8EX9XYEXFdsz46/N485vl8Hs9+O250zDr9jMxePCezRoZOBhP7dPTkJ712esB69vbYorckaTIJONAPEts84bm41fnHYajpo4BAKwur8LNr3+IJS0VAC3AiBcLVpNXJBjWELJWjiaGODLDZh2DvUJX1Zt82gpJzxN0sQ0nk+tK6JYV1IZkNZJd6JCBUFsSjzuCGMhj0jLtIb1amUarYLyCvb1atQILHyvFnSdNxyHDBuGPvz0Wxx2zD/7+78+w4b9loqwGwUL0kM1Epe8oM1zDOOUf4X6blcWD4nssx4f9Jw/Br644EoP7CjY0e8EK/Omzz9EUiSLSz2ZaniLBSiVDHJUnZDUqKLZD/fa63gEewfiLKYRXNrVhUCO92w4w+VRgy8wuX+oodWLvPpKHhvqEfZ9jBY2UmYe2G0hmdR5aMeMFrv3qI7xftg53HHUk9hrRF088diFe+XQZnvv312gKReCrF3L1NJF8Whwyi9lMEZt2eYLNzsPc+VSRMVYA4BUAgwFsBDCDc17rOmcUnSMxFMAtnPN/MMZuA3ApgEr67QbOebse+F0ugfz8IGacsT9OOnlfBPxehKNx/Ou9hXh88RLEDBMo7Oo77rmoaGjCxc+9iVOHjMbvTzkU40f1wzN3noNPPl+Jf73wNSoaKjoupAOMGdcPF/7qCEwaPwgAsLGsBne+8SkWrd6CsJJVWvh03Xp8//1WXH3MIZix/z44e/q+OG7/MXh5zmK88/LiTuWV7rbYddPn6wF8wjm/mzF2Pe3/IaEqnK8CMAEAGGM6gK0A3nKc8nfO+X2p3rDLOsXho0pw/Kn74ohj9rEscJ/NX4X73v4SFdWNiBV1eyf9bou5S37CVys24Iqp++GMn03EkdPG4MhpYzD/q1V4+40lWD5vDcw0prdZ2QHse+AwnHjWgdhrnwEAgIbGFjwzZxFe++C/aMjteVOe7oKmcBR3vP0p5nz4Pf5v5uGYNLI/fnnuYTjrxP0w9/3v8J9Xl6B8a23HBfUA7KKp8UkAptH3ZwF8Dlen6MKRANZxzjudhDutTtH0eRAZmA/OgNycDAwdVYJ9JwzCwVOGY8ggW4fyxbdr8cgnC7GydDuieQD6M0RkpOxCW8FcVEDT5rxEJ9+RAWHhHOyz3TtKKLdxAc2EgjRdtl04Op6CGTxxaDOTOFG5XT/aK9c9PdNIraDxKJUlo0XbK3Q0Gl41mlLL6N1eUuaXevKsc8NeYbQyAjqaEcZf5s3DMz98h18dMgU/P2AsDpo6CgdNHYWqmiZ8sWA1vl2+BSu/L0V1TRO0uP1smsZQ1CcX+4zrjxuvPT7hGRqaw3j9y+/w7IdLUOWLAn2ASH5ipGwZBxEA+hfUAQBGkBP9SHKRGuoTshvosf/wRbp0eKeo7IzySqeh4uhqmdnuWdLlh6Je6zJeo9B1aw6Vhy0zcW+p8thA7lXVXnsZZ4gMjt/WVOGsF17DtJJBuOKYAzBhSF/MPPNAzDzzQKxYU46v56/Bdyu2YN2GSoTrhUqCccAs70HT513TKfbmnJeL2/FyxlivDs4/E8BLrmNXMsbOA7AEwO/d02830pLAoEGFeOHfVyA7K4Ds7MTF6w1NYbz/xQq8Nu97bCqvRSRXMcOdgfK6Rvz5+Y/xyOz5mLH/3jh26lgMKMnHacdNwmnHiTzNkUgMVdVNME0Or0dHUVEWPJ7WHdH9z36Kd7/4AXWUeB4+JbOuxvyfNmH+T5uwb1EJTjt8Hxw9ZSTGjeiDcSPsyEb19SE0NUXQ1BzG5MkP7cbapoO03G2KGGPOFJyzOOez5A5j7GMAJUmuuzGdGjHGfABOBPBHx+FHAdwhKow7APwNwEXtlZNWpxjwe9G3Tx4AIBSOYs3WKqzYWIEvv1+PRVu2Im6YIvJzP4ZYtmiweC6NsvnChaNXvm11HpQjOmzp5DuMXDgkQ5TsEOgcQ5Qsw3LEdilBjHaEKp2HTRkKKgVXEclCgiD3JGv5mlOXlHyQkvmAnUaZUj0PANBMjNGk5WTxoIZmtODBRYvw4KJFmNirBAeOGYTJI/pjZL9i5GVloF/f/ITyq+qb8ePGCmyorMX3G8rxybI1iAcZUAjEsmRmRFrimCc6ycx8wRD759VZ5QzPFrIZLg0rLoYo2SHQOYbYlsxSkZVEKjKTx/zuv4DlPuVIpeCKvCZlZDFH3TZEVZJTftgn3KdkTvAFjRVY8H4FcmZ/igNGDcJho4ZgzMBeGNa3ELm5QeTmCqOZ2VMGJo50svlVtRcQgnN+VFu/Mca2Mcb6EEvsAzslUjIcC+Bbzrnly+f8zhh7AsC7HVU2rU5xQ1kNTr72KTSGI2gIhWE4ro4Heogw90As31iO5RvLMet9QI8AmQEfCrOCYAAMk6OqrhlhU3R0Rg/Ki7QnoiUax2fL1+GrxSKrpa4x5PkCyA76kR3oSaHDdplOcQ6A8wHcTdvZ7Zw7E66ps+xQafcUAB0mekqrUwzxGNboDeDZALIZTMfVRoCYYZBG+GxybcgRDLGIMrcNyraZ0tBMwTqG+IV+aiCFkpLLwPIcA3yqDNGpg4qTQ65kGe5gD2YSocrE8zIMlW653eitrmmLibTNGAGTWKOBOtrXqNzWzyOZyFbSYTUSY4yQIUuyEE/YHpC0CBBCDFVmvXwAoBDg5DBtOjrFOMnMyKY2yxZ1y8kRDLFfrihjaJatEx2aIWQl2Xw/j3iOApf+EOiYIbpZofgugzuQ7Fw6xKQyk3Kw9ISSdVO7OC5xy6otxmg6HO4NTrm56f2UMpP5fDTHDaTMqnShZ4wSYzQCog3imVRXh8ya42EAYbCW+la5dbo1dk2neDeAVxljFwPYDOAMAGCM9QXwJOf8ONoPAjgawOWu6+9ljE2AeAs2Jvm9FXqQVldBQaHbgMOKvLRTb8N5NYRF2X28DMBxjv0Qkjj8cc7PTfee6VmfdSBcyCFJDfc6Rng/Lf+iBFPBLMEQC7PEKNs3i0LRB23WMdAvvve1loEJPU4eMaOgY8leRwxRMgun3lAekyGlLIboutapr7LZBj0X7ctrvQ73bYsBcWlJTs4YA45mNolJGMQeY7pol6hPLie0GYTpsqpWSMboEYwxTkv2jLDDOTkqrmFG4rWcnN1Nn/2sPEAh2jKJ5WVRtj3KnT04i1I9ZNheAP19QmYlVG8ZiEM6ZjuXWabKEKWcACBG7F4yQrfMpKwS9IiWDlFA5oaWoeKczSjLddfNzRid75FBjvYm6RmjFMQ2Sjpek7eWmfQqqPaI52khxhjNoPBsEVtmLEYyM8V/rGegZ65rTgWKKSooKHQOqlME4OGI58XtIdljj6YevxgRM4KCIeZlUPqATBG4oR/l9u3js9MFlJA+qpAYogwKGyS9lDNQQEc6RLcuSnyXQUeJdVjsA3RNa8RasQ7JEIkNOAkYvRMdMUaPg13K8FaZpLPK5WSVp4CrYY/9zBFakhdz0QddE2U0eYl9hG0xGlE618UUQUxR8zmCIwQEA8rKoPBfmULv2y8oZNTXXwcA6O1I8SDrmUsht2SoNsnkPe1kfGmLIcaQRGYdyCqWoNulR5TnykfnrSUsU1EYPLk+2GaMthwMam9ZT/m+yqRmEUeQjYjZmj0CQB35NLaQzOIRRzvFqQ4mAE8P6Wg4APdS0j0EiikqKCh0AjzpoLMnQHWKCgoKnYOaPgPQODzZMTDptOq1pz0BH7lzBMgFJ0NMxXr5aesT0y6Z1xcA8nQZgVlcG7RyqSQuxUoGORWTbjduowpgT5ujXE7XxPFUxjfL0GJNxeQL4Lja+k1s5DTaPYF0TtHks/mp3pmkMsjjwsjR7LHzrYTIfyZEeT8ipnQLoqkfGV6aPbZ3cTRK0zcjcVqo0fTN57Ojv2T6hZEkLyDuXRQQ9y4mWRV5xLbQUSeZbS+T7u1nMmZhxxFv2po2J8hMyor2jU7IzCsF0s402i07d72d756UWZDJ3C8U3Yne35DjnbZlJabUUlbSVaeB2i3smHLHY6J8zhmg9ZCOZhdZn3cHFFNUUFDoHBRTBDSNIyMYscZQn8dmipk+Yh1+wSTyfWJb6BMsI99DxhTdDi6QwyiniQyK0AnWYbhcN5zOvm0xROn2YvC2mah0qbCchVsxRtsII/O4mK5yk9XfisEIGSmc3I/I7SOb21G0JRNppBiLIcpjHDcdzAKwmDsARGjJWdx0GXsks/PaTDHHJ9q/wC/uU+AT23xvc8L9ZZZCZz0DnWD10hhmbS352PV353BOR1ZAorzlvjOxoXS01qRBja5Nhd17rYjqYivf30bHOy3fc5sx0tLMNoxlABCldekc4j/WY6A6RQUFBQUC5wkpjvckpM0Ugz5bp+h3LITPItaR5RXbHI8YPbP0xNzIMosaAPhpxJXaFdtxumPWIdEe6zBd2yh5ndsO2a3vI9mfRTjpFJtZOOoiWYflttMx7HBjxLS4ZFyiljKvC2BnB5RMLZfatMUr7hSVOkbHc0jW1BZTzPA6Ajb4RLnZ3sTy7fzVksnb13itZ6XnYOnnVHGzeudfy2iDIbYnM7ejt89i7gKxdti9lYOnDXcqoPUyQiln+f4632n5nsv3Pscls3iS5Zwe3Wb+iinufiimqKCg0DmoTlEwooAn7mCKtn4qQAEBMlzbAJM6KJm/xL7GCrDqyvGbClqFAUsSfNTNOtxsw+1gCwCtiIhlnWy9RNBioi6m0paeygkrn7HUzdF9nTmK3dkAZSY6q22p/aOe1mHuY0bi3WUQ2wzHufJ6WZ4sX97Pur+jTrKeWhoyc+t/ZXu5mbzze1sMMRWZyXM1a0mf87dEVi/fI61daYHO0ai8xKV8znfabrPEtpXbsE46Ro9jmR+VwzlLCC7RvcGV9VlBQUHBAge4ct4WI5pXN6zRzOfQKfq0xHy5MgOa18qzm3gcsPU7uks/lQ7cDDGZ6rc9y6UbFhORjMjSe0l9kn0/gycy3c5At/RVct8u32e1IVk9tVjCvmxrnyNHcVxLzty8mtHqXLes3OX7rCyFTp1cYr13BMkCx7p1hu0yxDbLley7tc5SMkQ3Y0wFlv5UehdYucCdTDr5+27LitrW8d/hDlbs9CTo9lDL/BQUFBQInO+SFKe7A2l3ihq4tZIiWXBNOyFT4igq952soy1tVGfyACcTT6sQYW2wDic70V3M0woFZenF7N+8LvJi6zlT0U+lznwkI9FcYfDdW+d3Z1qDjs5tJbs0clem8xxuPbBEsqj27rBp1rkueQBts/tU3olUIN/HOE+ch7iDhohzk7/3rWTnvIb+TzDTf+93K5ShRUFBQcEGV0xRQUFBQYIrpijBEqZfZqvvWg9SFCdzBO7S8ttxCG4LiSqJxKmXRFvTaff3hHKTTp+TT5fd9+2sm4jb0b4r0Z7KozvCraJIlKlQt/QoI4sKCKGgoKBggwPgaplfeug5TqgKCjsPe+z/gHMVZFZBQUHBCb6HTp8ZT0NZyhirBLBp51VHQeF/HoM458W7uxIdgTH2HwBFKZ5exTk/ZmfWpyuRVqeooKCgsKejh3mLKigoKOxcqE5RQUFBwQHVKSooKCg4oDpFBQUFBQdUp6igoKDggOoUFRQUFBxQnaKCgoKCA6pTVFBQUHBAdYoKCgoKDvw/j+MrohGJ9mgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "im = ax.imshow(r, interpolation='nearest')\n",
    "plt.colorbar(im)\n",
    "\n",
    "for n, contour in enumerate(contours):\n",
    "    ax.plot(contour[:, 1], contour[:, 0], linewidth=2, color='white')\n",
    "\n",
    "ax.axis('image')\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determine if high"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.ndimage as ndimage    \n",
    "\n",
    "\n",
    "def analyse_contour(r, c):\n",
    "\n",
    "    # Create an empty image to store the masked array\n",
    "    r_mask = np.zeros_like(r, dtype='bool')\n",
    "\n",
    "    # Create a contour image by using the contour coordinates rounded to their nearest integer value\n",
    "    r_mask[np.round(c[:, 0]).astype('int'), np.round(c[:, 1]).astype('int')] = 1\n",
    "\n",
    "    # Fill in the hole created by the contour boundary\n",
    "    r_mask = ndimage.binary_fill_holes(r_mask)\n",
    "\n",
    "    # Find the change in mean value of inside cf edge to determine if high.\n",
    "    big = ndimage.binary_erosion(r_mask, structure=np.ones((1,1)))\n",
    "    sml = ndimage.binary_erosion(r_mask, structure=np.ones((2,2)))\n",
    "    \n",
    "    hilo = 'lo' if np.mean(r[big]) - np.mean(r[sml]) > 0 else 'hi'\n",
    "    area = np.sum(r_mask)\n",
    "\n",
    "    style = 'closed' if all(c[0] == c[-1]) else 'open'\n",
    "\n",
    "    return area, hilo, style"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/matt/anaconda3/envs/geocomp/lib/python3.7/site-packages/numpy/core/fromnumeric.py:3335: RuntimeWarning: Mean of empty slice.\n",
      "  out=out, **kwargs)\n",
      "/home/matt/anaconda3/envs/geocomp/lib/python3.7/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[(38, 'hi', 'open'),\n",
       " (614, 'hi', 'closed'),\n",
       " (38, 'hi', 'open'),\n",
       " (22, 'hi', 'open'),\n",
       " (196, 'hi', 'closed'),\n",
       " (22, 'hi', 'open')]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contours = [c for c in contours if ~np.isnan(c).any()]\n",
    "[analyse_contour(r, contour) for contour in contours]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding all contours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_largest_closure(r):\n",
    "\n",
    "    levels = np.linspace(r.min(), r.max(), 100)\n",
    "\n",
    "    all_largest = [0]\n",
    "    largest_contours = []\n",
    "    for level in levels:\n",
    "        contours = measure.find_contours(r, level)\n",
    "        \n",
    "        contours = [c for c in contours if ~np.isnan(c).any()]\n",
    "\n",
    "        largest = []\n",
    "        for contour in contours:\n",
    "            area, hilo, style = analyse_contour(r, contour)\n",
    "\n",
    "            if (hilo == 'lo') or (style == 'open'):\n",
    "                continue\n",
    "\n",
    "            largest.append(area)\n",
    "\n",
    "        try:\n",
    "            argmax = np.argmax(largest)\n",
    "            all_largest.append(largest[argmax])\n",
    "            largest_contours.append(contours[argmax])\n",
    "        except:\n",
    "            continue\n",
    "\n",
    "    print(all_largest)\n",
    "    record = np.argmax(all_largest)\n",
    "    record_level = levels[record]\n",
    "    record_contour = largest_contours[record]\n",
    "    record_area = max(all_largest)\n",
    "    \n",
    "    return record_contour, record_area, record_level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0]\n"
     ]
    },
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-7-5b8a49582583>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mfind_largest_closure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-6-bd30c6a06ffc>\u001b[0m in \u001b[0;36mfind_largest_closure\u001b[0;34m(r)\u001b[0m\n\u001b[1;32m     29\u001b[0m     \u001b[0mrecord\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mall_largest\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     30\u001b[0m     \u001b[0mrecord_level\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlevels\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mrecord\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 31\u001b[0;31m     \u001b[0mrecord_contour\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlargest_contours\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mrecord\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     32\u001b[0m     \u001b[0mrecord_area\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mall_largest\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     33\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: list index out of range"
     ]
    }
   ],
   "source": [
    "find_largest_closure(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geocomp",
   "language": "python",
   "name": "geocomp"
  },
  "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
}