access-contours.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Access contours\n",
    "\n",
    "Extract contours of access from a Conveyal regional analysis\n",
    "\n",
    "Input: regional analysis, [downloaded as a geotiff](https://docs.conveyal.com/analysis/regional#downloading-regional-results) and saved to the `data` directory\n",
    "\n",
    "Output: filled contours (.geojson multipolygons) corresponding to the intervals in `accessBins`. If discrete points (a \"freeform\" spatial dataset) are used as the destination in the Conveyal regional analysis, the contour containing the locations that can access at least one destination can be considered a combined catchment area or a composite \"reverse\" isochrone, of inbound access to destinations (as opposed to outward access from an origin)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import json\n",
    "import matplotlib.pyplot as plt\n",
    "from shapely.geometry import shape, MultiPolygon\n",
    "import shapely_geojson\n",
    "import geojsoncontour\n",
    "import geojson\n",
    "import rasterio\n",
    "\n",
    "\n",
    "accessBins = [1, 1000]\n",
    "inputLayers = [\n",
    "    {\n",
    "        'id': '60372f300554eb553d2be04b_603723df0554eb553d2be042',\n",
    "        'percentile': 25,\n",
    "        'cutoff': 45\n",
    "    },\n",
    "    {\n",
    "        'id': '60372f300554eb553d2be04b_603723df0554eb553d2be042',\n",
    "        'percentile': 25,\n",
    "        'cutoff': 60\n",
    "    },\n",
    "    {\n",
    "        'id': '60372f300554eb553d2be04b_603723df0554eb553d2be042',\n",
    "        'percentile': 25,\n",
    "        'cutoff': 90\n",
    "    }\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEQCAYAAACgBo8fAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVTUlEQVR4nO3df7DldX3f8ecLqMmoQTEs/sgGN5sB1FKh9Q62IVUhHSJCMT9qZq3EjqHZWScjpCES6Q/HpmObJu6UbhNDdkimcRxlaiPRMMZIpQqTSpu7ZREQRLqA7qyTvYgGMZ2kyLt/nHPleDk/vufe8/N7n4+ZO/ec7+dzznmfL+zrfO7nfL7fb6oKSdLyO2HeBUiSJsNAl6SWMNAlqSUMdElqCQNdklrCQJeklphroCf5vSTHk9zTsP/PJPlCknuTfGja9UnSMsk816EneQ3wBPCBqjp7RN8zgP8CXFhVX09yWlUdn0WdkrQM5jpCr6rbgMd6tyX54SSfTHIoye1JXtZt+nngt6rq693HGuaS1GMR59APAu+oqlcBvwy8v7v9TODMJH+a5I4kr59bhZK0gE6adwG9kjwX+BHgI0nWN39P9/dJwBnA64CdwO1Jzq6qb8y4TElaSAsV6HT+YvhGVZ3bp+0ocEdV/T/goSRfpBPwfzbD+iRpYS3UlEtVPU4nrN8EkI5zus1/CFzQ3X4qnSmYI/OoU5IW0byXLX4Y+BxwVpKjSa4A3gJckeQu4F7gjd3ufwJ8LckXgP8OvLOqvjaPuiVpEc112aIkaXIWaspFkrR5c/tS9NRTT61du3bN6+UlaSkdOnTo0ara0a9tboG+a9cuVldX5/XykrSUkjwyqM0pF0lqiUVbhy5JC2X3gf3fuX3kyqvnWMlojtAlLbTeQJ33a8+zliYMdElqCadcJC20SU5zDBph93uNRR+N92OgS2qFYQF85Mqrh7Y3De9Fn0M30CUtnPWAHSdA1/tu5rFt0ejQ/yQPA98Evg08WVUrG9rfAvxK9+4TwNur6q5hz7myslKuQ5fab94B2/v6m6ll3vVvlOTQxgxeN84I/YKqenRA20PAa7uXhruYzkUqXj1mnZI0cU3mx4eF9aIEeRMTmXKpqv/Rc/cOOhegkDSmRRsNTsIyvJfdB/YvRZ2jNF22WMCnutf53Dui7xXAH/drSLI3yWqS1bW1tXHqlNRiuw/s39Sqks0+Dpp/0Ky/xjKsemk6h/6SqjqW5DTgFjrX/LytT78L6FwD9EdHnavcOXRJ6/pNgUzrCM1+c+pNXmdR/noaNoc+9vnQk7wHeKKq3rdh+yuBm4CLq+qBUc9joEsapd+oeBqhO+hxvVMx44T/NG3pS9EkzwFOqKpvdm9fBPzqhj6nAx8FfrZJmEvaXoYF7bC2WYfmxrn0YevXF3HefeQIPcluOiNv6HwAfKiq3ptkH0BVXZ/kBuCngfXTOj5jaeNGjtAlrdsYjuNMtzTtO86BR4NqGfb4WZnolMukGOiSBpnk9MY4gbzZD5Kmj5mESa1Dl6SpmvRKknHO3bK+ff0xmz0IaZ7TMI7QJS2UaZ5/fNCoul94b2YFTD+Tfg9OuUhaaKNG5oNCcRnWhs8y0J1ykTRXTUK56Xz1qMP8By1B7G2fxrTPrKZhHKFL+o55XG5tnNH5LKY2xvkAGMek6nOELmmoRZu6GPal5brNnipgnLnwSe6XWYzUDXRpG1m0NdVb+Ytgq0eMrgdsb99JjsjnwSkXacFM85wh4y7jm4VhK0zmcaTotG31PbnKRdJMNR0h95rHXwbT+BK0ia2812GB3vT0uZLUWO/0xbA+vSZxitre5xj0fIswnTKtGpxDl7oW5fSo2pxRI/5FCvdp/T/mlIvUZaDPx6yXIvZ7/mWadnEOXdJC2fjhOWhOu0noTTKM5zGnPm6wG+iS5mLYhSOa2My6882EcpP6ph32TYPdQJe0UGYVjJNadz+rUXuTWgx0SQtrWt9dTOp85Ys2BeOh/5IW1iSCvMm68u3wZbfr0CXN3KRHvcPWvTdZEz/INE+VsJW6BmkU6EkeTnJ3ksNJnjFPkuRlST6X5K+S/PJEK5S0tAYd3DON6ZVpLzvdGMBbXfI4jf0yzpTLBVX16IC2x4ArgZ/YUjWSWmVa8+JNn7dpv81+GEw61Lf6odToS9EkDwMrQwJ9vd97gCeq6n2jntMvRSVNy6QOVlrEA462vMolyUPA14ECfqeqDg7o9x6GBHqSvcBegNNPP/1VjzzyyMjXlqR1w778HLaWfL3/pI88nccFMCaxyuX8qjqW5DTgliT3V9VtDR/7Hd0PgoPQGaGP+3hJWrcx1Dd7XdJxbXzNrZzHfdJ/ATQK9Ko61v19PMlNwHnA2IGuzdvKEXfbYbmWtofNjIiHXWd0M/82mlzLtN9fB7MwMtCTPAc4oaq+2b19EfCrU69MwGQuieVJp7RdTXqN+zj9m96fpJFz6El2Azd1754EfKiq3ptkH0BVXZ/kRcAqcDLwFPAE8IqqenzQ8/qlaHPT+B/AcNcym8fFrMc17kU0PJdLS03rCxZH6tL8TOrfn4f+t9QkL6orabpm8e/PQF9ABq+kzfBcLlLLLco1NBehjrYz0KU5mkXQ+Rff9mGgS3O0XcJ2GmcW1DMZ6NKcGXSaFANdklrCQJekljDQJaklDHRJagkDXZL6WMa18wa6JLWEgS5JPTaOzJdplO65XCSpx8YrCi3TcQKO0CVpgGUKczDQJak1DHRJagkDXZJaolGgJ3k4yd1JDid5xnXj0nEgyYNJPp/k70y+VG3WMq6nlTS+cVa5XFBVjw5ouxg4o/vzauC3u7+1gLy2qNROk1q2+EbgA9W54vQdSZ6f5MVV9dUJPf+2MY2w3fhcBrnUTk3n0Av4VJJDSfb2af8B4Cs99492t32XJHuTrCZZXVtbG79aSdJATUfo51fVsSSnAbckub+qbutpT5/H1DM2VB0EDgKsrKw8o12OniVtXqMRelUd6/4+DtwEnLehy1HgB3vu7wSOTaJASVIzIwM9yXOSfN/6beAi4J4N3T4OvLW72uXvAn/h/LkWhat8tF00mXJ5IXBTkvX+H6qqTybZB1BV1wOfAN4APAj8JfC26ZQrNWeIa7tJZ2HK7K2srNTq6jOWtEtbNizI/Y5Cyy7Joapa6dfmkaLaVpx+UZt5+ly1zpErr24U2hv7OHrXsnOErm3JUbrayBG6Wml9tN0kuB2Zqy0MdLWaYa3txCkXSWoJA12SWsIpF03FvE/R2zt37rSLtgsDXa3UphDv98Vum96fJscjRaUF1nR5pQG/fXikqLSExlkr7xGwAqdcpIXVdC29o3Otc4QuLbhhgW2Yq5cjdGkJGNxqwhG6JLWEgS5JLWGgS1JLNA70JCcmuTPJzX3aTklyU5LPJ/lfSc6ebJmSpFHGGaFfBdw3oO2fA4er6pXAW4H/uNXCJEnjaRToSXYClwA3DOjyCuDTAFV1P7AryQsnUqEkqZGmI/TrgGuApwa03wX8FECS84CXAju3WpwkqbmRgZ7kUuB4VR0a0u3XgFOSHAbeAdwJPNnnufYmWU2yura2tsmSJUn9jDw5V5J/B/wsnYD+XuBk4KNVdfmA/gEeAl5ZVY8Pel5PziVJ49vSybmq6tqq2llVu4A9wK0bwzzJ85M8q3v3nwK3DQtzSdLkbXodepJ9SfZ1774cuDfJ/cDFdFbESJJ6TPuMmJ4PXZKWiOdDl6RtwECXpJYw0CWpJQx0SWoJA12SWsJAl6SWMNAlqSUMdElqCQNdWkK7D+yf+lGHWj4nzbsAaVo2Bt6RK6+eUyWT16b3oskx0NU622Hkuv4eDXb1MtC11LZDePdjkKsf59AlqSUMdC01R6rS05xy0dIbFurONWs7MdDVaga5thOnXCSpJQx0SWqJxoGe5MQkdya5uU/b85L8UZK7ktyb5G2TLVOSNMo4I/SrgPsGtP0C8IWqOgd4HbA/ybO2WJskaQyNAj3JTuAS4IYBXQr4viQBngs8Bjw5kQolSY00HaFfB1wDPDWg/TeBlwPHgLuBq6rqGX2T7E2ymmR1bW1tE+VKkgYZGehJLgWOV9WhId1+HDgMvAQ4F/jNJCdv7FRVB6tqpapWduzYsbmKJUl9NRmhnw9cluRh4EbgwiQf3NDnbcBHq+NB4CHgZROtVJI01MhAr6prq2pnVe0C9gC3VtXlG7p9GfgxgCQvBM4Cjky4VknSEJs+UjTJPoCquh74N8B/TnI3EOBXqurRyZQoSWpirECvqs8An+nevr5n+zHgokkWJkkaj0eKSlJLGOiS1BIGuiS1hIEuSS1hoEtSSxjoktQSBroktYSBLkktYaBLUksY6JLUEga6pE3ZfWA/uw/sn3cZ6rHpk3NJ2t6OXHn1vEvQBo7QJaklDHRJagkDXZJawkCXpJYw0CWpJQx0SWqJxoGe5MQkdya5uU/bO5Mc7v7ck+TbSV4w2VIlScOMM0K/CrivX0NV/UZVnVtV5wLXAp+tqscmUJ8WjAeTSIur0YFFSXYClwDvBX5pRPc3Ax/eYl1aUB5MIi2upiP064BrgKeGdUrybOD1wB8MaN+bZDXJ6tra2jh1SpJGGBnoSS4FjlfVoQbP9w+BPx003VJVB6tqpapWduzYMWapkqRhmozQzwcuS/IwcCNwYZIPDui7B6dbJGkuRgZ6VV1bVTurahedwL61qi7f2C/J84DXAh+beJWSpJE2vQ49yb4k+3o2/STwqar61tbLkiSNK1U1lxdeWVmp1dXVuby2JC2rJIeqaqVfm0eKaiZcvy5Nn4EuSS3hFYvUSJPRtQcdSfPlHLoam/aUiR8I0mjD5tANdG1K03A/cuXVY38QGOzSYMMC3SkXbcrGoO4N4d7tg8J8Y2j7ham0dY7QNTXjjOJ7+ztClwZzhK652BjUTftL2hwDXVNnUEuz4Tp0SWoJA12SWsJAl6SWMNAlqSUMdElqCQNdklrCQJekljDQJaklDHRJaonGgZ7kxCR3Jrl5QPvrkhxOcm+Sz06uRElSE+Mc+n8VcB9w8saGJM8H3g+8vqq+nOS0yZQnSWqq0Qg9yU7gEuCGAV3+MfDRqvoyQFUdn0x5kqSmmk65XAdcAzw1oP1M4JQkn0lyKMlb+3VKsjfJapLVtbW18auVJA00MtCTXAocr6pDQ7qdBLyKzij+x4F/leTMjZ2q6mBVrVTVyo4dOzZbsySpjyZz6OcDlyV5A/C9wMlJPlhVl/f0OQo8WlXfAr6V5DbgHOCBiVcsSepr5Ai9qq6tqp1VtQvYA9y6IcwBPgb8/SQnJXk28Go6X6BKkmZk0xe4SLIPoKqur6r7knwS+DydefYbquqeCdUoSWrAa4pK0hIZdk1RjxSVpJYw0CWpJQx0SWoJA12SWsJAl6SWMNAlqSUMdElqCQNdWmC7D+xn94H98y5DS2LTR4pKmr4jV1497xK0RByhS1M2rVH2oozeF6UOOUKXltaijN4XpQ4Z6NLUGXiaFadcJKklDHRJagkDXZJawkCXpJYw0CWpJQx0SWqJxoGe5MQkdya5uU/b65L8RZLD3Z93T7ZMSdIo46xDvwq4Dzh5QPvtVXXp1kuS2qP3CErXo2vaGgV6kp3AJcB7gV+aakVSixjimqWmUy7XAdcATw3p8/eS3JXkj5P8zX4dkuxNsppkdW1tbcxSJUnDjAz0JJcCx6vq0JBu/xt4aVWdA/wn4A/7daqqg1W1UlUrO3bs2Ey9kqQBmozQzwcuS/IwcCNwYZIP9naoqser6onu7U8AfyPJqZMuVpI02MhAr6prq2pnVe0C9gC3VtXlvX2SvChJurfP6z7v16ZQryRpgE2fbTHJPoCquh74R8DbkzwJ/F9gT1XVZEqUJDWReeXuyspKra6uzuW1JWlZJTlUVSv92jxSVJJawkCXpJYw0CWpJQx0aZvzIs/tYaBLUkt4kWhpm/N8M+3hCF2SWsJAl6SWMNAlqSUMdElqCQNdklrCQJeklpjbybmSrAGPzOXFF8OpwKPzLmIJuJ9Gcx8105b99NKq6nuFoLkF+naXZHXQGdP0NPfTaO6jZrbDfnLKRZJawkCXpJYw0Ofn4LwLWBLup9HcR820fj85hy5JLeEIXZJawkCXpJYw0KcoyZuS3JvkqSQDl0sleX2SLyZ5MMm7xn38spvAfnpBkluSfKn7+5TZVD47Td9jkquS3NPdn7/Ys/2cJJ9LcneSP0py8syKn6EJ7Kdzk9yR5HCS1STnzaz4CTDQp+se4KeA2wZ1SHIi8FvAxcArgDcneUXTx7fEVvfTu4BPV9UZwKe799tm5HtMcjbw88B5wDnApUnO6DbfALyrqv4WcBPwzplUPXtb3U+/DvzrqjoXeHf3/tIw0Keoqu6rqi+O6HYe8GBVHamqvwZuBN44xuOX3lb3U/f373dv/z7wE1MpdL6avMeXA3dU1V9W1ZPAZ4Gf7LadxdMfmLcAPz29Uudqq/upgPW/Xp4HHJteqZNnoM/fDwBf6bl/tLtN323YfnphVX0VoPv7tBnXNgtN3uM9wGuSfH+SZwNvAH6wp+2y7u039Wxvm63up18EfiPJV4D3AddOv+TJ8RJ0W5TkvwEv6tP0L6rqY02eos+21q0ldT+NNmwfNXl8Vd2X5N/TGYE/AdwFPNlt/jngQJJ3Ax8H/nrrFc/HlPfT24F/VlV/kORngN8F/sHWq54NA32Lqmqr/7GP8t2jpZ0s2Z95TUx5P/15khdX1VeTvBg4vsXXmoth+yhJo/dYVb9LJ4RI8m/p7Deq6n7gou72M4FLJlz+zExzPwH/BLiqe/sjdL57WBpOuczfnwFnJPmhJM8C9tAZQem7DdtPH6fzD5Hu7yYj/mXT6D0mOa37+3Q6XzR/eMP2E4B/CVw/5XrnZUv7ic4g4bXd2xcCX5papdNQVf5M6YfOFy1Hgb8C/hz4k+72lwCf6On3BuAB4P/QmYIY+vi2/UxgP30/nRUNX+r+fsG839MU9lHf99hnH90OfIHONMKP9Wy/qrvvHgB+je5R4m37mcB++lHgUHf7/wReNe/3NM6Ph/5LUks45SJJLWGgS1JLGOiS1BIGuiS1hIEuSTOQ5PeSHE9yT4O+/6F7grDDSR5I8o1Gr+EqF0maviSvoXNk6geq6uwxHvcO4G9X1c+N6usIXZJmoKpuAx7r3Zbkh5N8MsmhJLcneVmfh76Zpw98GspD/yVpfg4C+6rqS0leDbyfzhGqACR5KfBDwK1NnsxAl6Q5SPJc4EeAjyTfOffc92zotgf4r1X17SbPaaBL0nycAHyjOhfTGGQP8AvjPKEkacaq6nHgoSRvAkjHOevtSc4CTgE+1/Q5DXRJmoEkH6YTzmclOZrkCuAtwBVJ7gLu5emrcEHny9Aba4yliC5blKSWcIQuSS1hoEtSSxjoktQSBroktYSBLkktYaBLUksY6JLUEv8fsI08az0Hvs4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEQCAYAAACgBo8fAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVkElEQVR4nO3dfbBtdX3f8fcHqMmoQTFcfMgN3twMoJYIrWewDakK6RARgnmoGazEjqG5c52MlzQ3EumDY9OxTRvvlN4mltwhmcRxlKmNRMMYI5UoTCJpzi0XAUGkF1AGJ/cgGsR0kiLf/rH3wcNhP6x99vM679fMmbP3Wr+99ncvuJ/9O7/1W2ulqpAkLb/j5l2AJGkyDHRJagkDXZJawkCXpJYw0CWpJQx0SWqJuQZ6kt9JcizJnQ3b/0ySLyS5K8mHpl2fJC2TzHMeepLXAI8DH6iqM4e0PQ3478D5VfX1JKdU1bFZ1ClJy2CuPfSquhl4dOOyJD+Y5JNJDie5JcnLuqt+HvjNqvp697WGuSRtsIhj6IeAd1TVq4BfBt7fXX46cHqSP01ya5LXz61CSVpAJ8y7gI2SPBf4YeAjSdYXf1f39wnAacDrgJ3ALUnOrKpvzLhMSVpICxXodP5i+EZVnd1j3UPArVX1/4D7k3yRTsD/xQzrk6SFtVBDLlX1GJ2wfhNAOs7qrv4D4Lzu8pPpDMEcnUedkrSI5j1t8cPA54AzkjyU5HLgLcDlSW4H7gLe2G3+x8DXknwB+BPgnVX1tXnULUmLaK7TFiVJk7NQQy6SpK2b20HRk08+uXbt2jWvt5ekpXT48OFHqmpHr3VzC/Rdu3axuro6r7eXpKWU5MF+6xxykaSWWLR56JK0UHYfPPDU46P79s+xkuHsoUtaaBsDdd7vPc9amjDQJQ21++CBhQ8zOeQiqYF5DjVsfu9xhkD6fSn12s4yfoHN7cSilZWVcpaLpHE0Dd2j+/ZPJKAXYQw9yeGqWum1ziEXSQun6RDP0X37nxay/QJ3GXvbW9Goh57kAeCbwLeBJzZ/OyR5C/Ar3aePA2+vqtsHbdMeutR+85whMijEm9TS6/WL3kMfZQz9vKp6pM+6+4HXdm8NdyGdm1S8esQ6JWliBo2Lbw7rXm0XIbxHNZGDolX1Zxue3krnBhSSRrQeNMsYJr0sy+fYffDA0tQ6SNNAL+BTSQr4rao6NKDt5cAf9VqRZA+wB+DUU08dpU5JLbbVoZlxvgCbHihdphOLmo6hv6SqHk5yCnAjnXt+3tyj3Xl07gH6I8OuVe4YuqR1vYZAphWkG78Emgy99HrdPA0aQx952mKS9wCPV9X7Ni1/JXA9cGFV3TtsOwa6pGFGPTC51dDt97qNQzGjhP80jXVQNMlzgOOq6pvdxxcAv7qpzanAR4GfbRLmkraXQUE7aN08ZsZsngbZb1hmEcfdh/bQk+ym0/OGzhfAh6rqvUn2AlTVNUmuBX4aWL+s4zOmNm5mD13Sus3hOMpwS9O2w6Yx9ttO0/nwszLRIZdJMdAl9TPJ4Y1RAnmrXyRNXzMJk5qHLklTNekzOke5dsv68lHmqvd6v3kOw9hDlzR3gw4+wuQOhPbbdq9tbGUGTC+TDniHXCQttGE982W+RouBLmnbmGQoD7sM7qC/AtbXT+NLYpKhbqBLamQeZ0U27Z2PcqncSdazaJfe9aCopL4W6RZv64bNSd/KKftNvqx6bXdS+2cWB0wNdGkbWbQ51Vt936ZTCvt9MfQ6/X/jtdWXYWy+F4dcpAUzzWuGDOsRz+N6JYNmmCzSNdQnZdzP5Bi6pKcMmo43qSBtsp1R54hP0iL0wLf6OQ10STPX73T+JmPhw8a4+42J9xtK6bWNeTPQpSlalMujbifjnj7f74SkYdMX522c/8cMdKkBA30+ZnGWZdNjB7O2lc/ntEWpAYN8dvp9eQ6b8z3KSUKj1DEvk57KaA9d0tSMO4496tzz9deMGtSLMM7eNNgdcpG0UKYZmJO+lvksw71JqDvkImmhzOK2bvO6Xvk4xh2CMdAlzc2kD3oOGm5Z9DCfBANd0sxN+mDgoG1N+25HWzWNGTbHNWmU5IEkdyQ5kuQZA99JXpbkc0n+JskvT6w6SUtt98EDMxn62Pg+0+qJb7zWyyRMY7+M0kM/r6oe6bPuUWAf8BNjVSOpVaYR3KNst2m7rXwZbGUGTr/3HaeOjRrNcknyALAyINDX270HeLyq3jdsm85ykTQtkzpZaRGnL449bTHJ/cDXgQJ+q6oO9Wn3HgYEepI9wB6AU0899VUPPvjg0PeWpHWDDn4Omku+3n4aN7/otXwrJjEPvemQy7lV9XCSU4Abk9xTVTc3fO1Tul8Eh6DTQx/19ZK0bnOoNznDdCt3ZBo0DLLVL4lpXXJg5BOLBvXCHXKZnnHOuNsO07W0vWwlCEe5TPC4Y9nD7pY0zgHRsXroSZ4DHFdV3+w+vgD41UbvrLH1O2iylW0Y7GqDWd+EYpxru/daN836h/bQk+wGru8+PQH4UFW9N8legKq6JsmLgFXgROBJ4HHgFVX1WL/t2kNvbtHvQi7N0jTPLp2UJj3yzbyWS0tN6wDLPO7oLm1Xm3v2k/r357VcWmrU/ykMcWl+ZvHvzx76knAcXBIM7qE3OvVf8zfp0461fcz7Jg7rNSxCHW1noEtzNIugsyOwfRjo0hxtl7D1L8zZMNClOTPoNCkGuiS1hIEuSS1hoEtSSxjoktQSBrok9bCMc+cNdElqCQNdkjbY3DNfpl66F+eSpA02Xwt9mc4TsIcuSX0sU5iDgS5JrWGgS1JLGOiS1BKNAj3JA0nuSHIkyTPuSpGOg0nuS/L5JH9/8qVqq5ZxPq2k0Y0yy+W8qnqkz7oLgdO6P68G/lv3txaQdz+S2mlS0xbfCHygOvezuzXJ85O8uKq+OqHtbxvTCNvN2zLIpXZqOoZewKeSHE6yp8f67wO+suH5Q91lT5NkT5LVJKtra2ujVytJ6qtpD/3cqno4ySnAjUnuqaqbN6xPj9c84+7TVXUIOASdm0SPXO02YO9Z0lY16qFX1cPd38eA64FzNjV5CPj+Dc93Ag9PokBJUjNDAz3Jc5J8z/pj4ALgzk3NPg68tTvb5R8Af+X4uRaFs3y0XTQZcnkhcH2S9fYfqqpPJtkLUFXXAJ8A3gDcB/w18LbplCv1tzG0j+7bb4hr20lnYsrsrays1OrqM6a0S1vWJMDXg95jFVpWSQ5X1UrPdQa62qZpz9xQ1zIy0LXtbGW4xYDXMhgU6F4PXa0yzri5Z9Bq2Rnoao1xD4Ia5Fp2Brpaw0DWduflcyWpJQx0SWoJh1w0FfM+wLj5JCNpOzDQ1UptCvFeB3vb9Pk0Oc5DlxbYKDN3DPntYdA8dMfQpQU16jRMr10jh1ykBbXe4x4W1PbMtc4eurTgBgW2Ya6N7KFLS8DgVhP20CWpJQx0SWoJA12SWqJxoCc5PsltSW7ose6kJNcn+XyS/5XkzMmWKUkaZpQe+hXA3X3W/UvgSFW9Engr8F/GLUySNJpGgZ5kJ3ARcG2fJq8APg1QVfcAu5K8cCIVSpIaadpDvxq4Eniyz/rbgZ8CSHIO8FJg57jFSZKaGxroSS4GjlXV4QHNfg04KckR4B3AbcATPba1J8lqktW1tbUtlixJ6mXoxbmS/AfgZ+kE9HcDJwIfrarL+rQPcD/wyqp6rN92vTiXJI1urItzVdVVVbWzqnYBlwI3bQ7zJM9P8qzu038O3DwozCVJk7fleehJ9ibZ2336cuCuJPcAF9KZESNJ2mDaV8T0euiStES8HrokbQMGuiS1hIEuSS1hoEtSSxjoktQSBroktYSBLkktYaBLUksY6NIS2n3wwNTPOtTyOWHeBUjT0Cvsju7bP4dKpqNNn0WTY6Br6W3Hnur6ZzbYtZHXctFSGzXMDUAtO6/lImGYq/0MdElqCQNdS21jr/vovv32wrWtOYYuSUvEMXRJ2gYMdElqicaBnuT4JLcluaHHuucl+cMktye5K8nbJlumJGmYUXroVwB391n3C8AXquos4HXAgSTPGrM2SdIIGgV6kp3ARcC1fZoU8D1JAjwXeBR4YiIVSpIaadpDvxq4Eniyz/rfAF4OPAzcAVxRVc9om2RPktUkq2tra1soV5LUz9BAT3IxcKyqDg9o9mPAEeAlwNnAbyQ5cXOjqjpUVStVtbJjx46tVSxJ6qlJD/1c4JIkDwDXAecn+eCmNm8DPlod9wH3Ay+baKWSpIGGBnpVXVVVO6tqF3ApcFNVXbap2ZeBHwVI8kLgDODohGuVJA2w5cvnJtkLUFXXAP8O+N0kdwABfqWqHplMiZKkJkYK9Kr6DPCZ7uNrNix/GLhgkoVJkkbjmaKS1BIGuiS1hIEuSS1hoEtSSxjoktQSBroktYSBLkktYaBLUksY6JLUEga6JLWEgS5pS3YfPMDugwfmXYY22PLFuSRtb0f37Z93CdrEHroktYSBLkktYaBLUksY6JLUEga6JLWEgS5JLdE40JMcn+S2JDf0WPfOJEe6P3cm+XaSF0y2VEnSIKP00K8A7u61oqp+varOrqqzgauAz1bVoxOoTwvGk0mkxdXoxKIkO4GLgPcCvzSk+ZuBD49ZlxaUJ5NIi6tpD/1q4ErgyUGNkjwbeD3w+33W70mymmR1bW1tlDolSUMMDfQkFwPHqupwg+39OPCn/YZbqupQVa1U1cqOHTtGLFWSNEiTHvq5wCVJHgCuA85P8sE+bS/F4RZJmouhgV5VV1XVzqraRSewb6qqyza3S/I84LXAxyZepSRpqC3PQ0+yN8neDYt+EvhUVX1r/LIkSaNKVc3ljVdWVmp1dXUu7y1JyyrJ4apa6bXOM0U1E85fl6bPG1xoagxwabYMdDXSNJyP7ttvkEtzYqBrooaF+aD1noUqjccxdA21OYSnFbyOs0vjMdDV2NF9+58K882hvr5u1LC3Vy5NjtMWNTWj9LY3jr0b8lJ/g6YtGuiaqiahboBLzQ0KdA+KaqoMa2l2HEOXpJYw0CWpJQx0SWoJA12SWsJAl6SWMNAlqSUMdElqCQNdklrCQJeklmgc6EmOT3Jbkhv6rH9dkiNJ7kry2cmVKElqYpRT/68A7gZO3LwiyfOB9wOvr6ovJzllMuVJkppq1ENPshO4CLi2T5N/Cny0qr4MUFXHJlOeJKmppkMuVwNXAk/2WX86cFKSzyQ5nOStvRol2ZNkNcnq2tra6NVKkvoaGuhJLgaOVdXhAc1OAF5Fpxf/Y8C/SXL65kZVdaiqVqpqZceOHVutWZLUQ5Mx9HOBS5K8Afhu4MQkH6yqyza0eQh4pKq+BXwryc3AWcC9E69YktTT0B56VV1VVTurahdwKXDTpjAH+Bjwj5KckOTZwKvpHECVJM3Ilm9wkWQvQFVdU1V3J/kk8Hk64+zXVtWdE6pRktSAt6CTpCUy6BZ0nikqSS1hoEtSSxjoktQSBroktYSBLkktYaBLUksY6JLUEga6tMB2HzzA7oMH5l2GlsSWzxSVNH1H9+2fdwlaIvbQpSmbVi97UXrvi1KH7KFLS2tReu+LUocMdGnqDDzNikMuktQSBroktYSBLkktYaBLUksY6JLUEga6JLVE40BPcnyS25Lc0GPd65L8VZIj3Z93T7ZMSdIwo8xDvwK4Gzixz/pbquri8UuS2mPjGZTOR9e0NQr0JDuBi4D3Ar801YqkFjHENUtNh1yuBq4EnhzQ5h8muT3JHyX5u70aJNmTZDXJ6tra2oilSpIGGRroSS4GjlXV4QHN/jfw0qo6C/ivwB/0alRVh6pqpapWduzYsZV6JUl9NOmhnwtckuQB4Drg/CQf3Nigqh6rqse7jz8B/J0kJ0+6WElSf0MDvaquqqqdVbULuBS4qaou29gmyYuSpPv4nO52vzaFeiVJfWz5aotJ9gJU1TXAPwHenuQJ4P8Cl1ZVTaZESVITmVfurqys1Orq6lzeW5KWVZLDVbXSa51nikpSSxjoktQSBroktYSBLm1z3uS5PQx0SWoJbxItbXNeb6Y97KFLUksY6JLUEga6JLWEgS5JLWGgS1JLGOiS1BJzuzhXkjXgwbm8+WI4GXhk3kUsAffTcO6jZtqyn15aVT3vEDS3QN/ukqz2u2KavsP9NJz7qJntsJ8ccpGkljDQJaklDPT5OTTvApaE+2k491Ezrd9PjqFLUkvYQ5ekljDQJaklDPQpSvKmJHcleTJJ3+lSSV6f5ItJ7kvyrlFfv+wmsJ9ekOTGJF/q/j5pNpXPTtPPmOSKJHd29+cvblh+VpLPJbkjyR8mOXFmxc/QBPbT2UluTXIkyWqSc2ZW/AQY6NN1J/BTwM39GiQ5HvhN4ELgFcCbk7yi6etbYtz99C7g01V1GvDp7vO2GfoZk5wJ/DxwDnAWcHGS07qrrwXeVVU/BFwPvHMmVc/euPvpPwH/tqrOBt7dfb40DPQpqqq7q+qLQ5qdA9xXVUer6m+B64A3jvD6pTfufur+/r3u498DfmIqhc5Xk8/4cuDWqvrrqnoC+Czwk911Z/CdL8wbgZ+eXqlzNe5+KmD9r5fnAQ9Pr9TJM9Dn7/uAr2x4/lB3mZ5u0H56YVV9FaD7+5QZ1zYLTT7jncBrknxvkmcDbwC+f8O6S7qP37RheduMu59+Efj1JF8B3gdcNf2SJ8db0I0pyf8EXtRj1b+qqo812USPZa2bS+p+Gm7QPmry+qq6O8l/pNMDfxy4HXiiu/rngINJ3g18HPjb8Suejynvp7cD/6Kqfj/JzwC/Dfzj8aueDQN9TFU17n/sh3h6b2knS/ZnXhNT3k9/meTFVfXVJC8Gjo35XnMxaB8lafQZq+q36YQQSf49nf1GVd0DXNBdfjpw0YTLn5lp7ifgnwFXdB9/hM6xh6XhkMv8/QVwWpIfSPIs4FI6PSg93aD99HE6/xDp/m7S4182jT5jklO6v0+lc6D5w5uWHwf8a+CaKdc7L2PtJzqdhNd2H58PfGlqlU5DVfkzpR86B1oeAv4G+Evgj7vLXwJ8YkO7NwD3Av+HzhDEwNe37WcC++l76cxo+FL39wvm/ZmmsI96fsYe++gW4At0hhF+dMPyK7r77l7g1+ieJd62nwnspx8BDneX/znwqnl/plF+PPVfklrCIRdJagkDXZJawkCXpJYw0CWpJQx0SZqBJL+T5FiSOxu0/c/dC4QdSXJvkm80eg9nuUjS9CV5DZ0zUz9QVWeO8Lp3AH+vqn5uWFt76JI0A1V1M/DoxmVJfjDJJ5McTnJLkpf1eOmb+c6JTwN56r8kzc8hYG9VfSnJq4H30zlDFYAkLwV+ALipycYMdEmagyTPBX4Y+Ejy1LXnvmtTs0uB/1FV326yTQNdkubjOOAb1bmZRj+XAr8wygYlSTNWVY8B9yd5E0A6zlpfn+QM4CTgc023aaBL0gwk+TCdcD4jyUNJLgfeAlye5HbgLr5zFy7oHAy9rkaYiui0RUlqCXvoktQSBroktYSBLkktYaBLUksY6JLUEga6JLWEgS5JLfH/AYyEuaFSpMs3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEQCAYAAACgBo8fAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAV1UlEQVR4nO3dfbBtdX3f8fdHqMn4gGK4oOYGb8gAaqnQegbbkPpAOkaEYh5qBiuxY2juXMcRbIhE+uDYdmzTRKb0NjE3d0imcRxlaiPRMIZIpQqTSptzy0VAEOkFlMHmHkRDMJ2kyLd/7H30cNgPa5+z9tM679fMmbPPXr+993eve89n/85v/dZvpaqQJC2/Z8y7AElSOwx0SeoIA12SOsJAl6SOMNAlqSMMdEnqiLkGepLfSXI0yZ0N2/9ski8luSvJR6ddnyQtk8xzHnqSVwOPAx+uqjPGtD0V+M/AuVX1zSQnVtXRWdQpSctgrj30qroZeHTjfUl+JMkNSQ4luSXJS/ubfgH4jar6Zv+xhrkkbbCIY+gHgXdV1SuBXwI+1L//NOC0JH+c5NYkb5hbhZK0gI6ddwEbJXkO8KPAx5Os3/19/e/HAqcCrwV2A7ckOaOqvjXjMiVpIS1UoNP7i+FbVXXWgG0PAbdW1f8D7k/yZXoB/yczrE+SFtZCDblU1WP0wvrNAOk5s7/594HX9e8/gd4QzJF51ClJi2je0xY/BnwBOD3JQ0kuAd4KXJLkduAu4E395n8EfCPJl4D/Brynqr4xj7olaRHNddqiJKk9CzXkIknaurkdFD3hhBNqz54983p5SVpKhw4deqSqdg3aNrdA37NnD6urq/N6eUlaSkkeHLbNIRdJ6ohFm4cuSQvllP1Xfff2kUsvn2Ml49lDl7TQNgbqvF97nrU0YaBLGuuU/VctfJi1bRnfr0Mukhba5mGO7QyBDAvpRR9KaWpuJxatrKyUs1wkbUfTXvSRSy9vpce9CMGf5FBVrQza5pCLpIXTdIjnyKWXPyVkhwXuMg6fbEWjHnqSB4A/B74DPLH50yHJW4Ff7v/4OPCOqrp91HPaQ5e6b54zRMaF+FbCf9F76JOMob+uqh4Zsu1+4DX9S8OdR+8iFa+asE5Jas2g8N0Y1sM+bIaN2S9CmI/TykHRqvrvG368ld4FKCRNaJnCo4lleR+n7L9qaWodpekYegGf6V/nc++YtpcAfzhoQ5K9SVaTrK6trU1Sp6QOWx8zn3SsezvTKduaIbNImo6hv7iqHk5yInAjvWt+3jyg3evoXQP0x8atVe4YuqR1m8Ny86yUNnvPG/8KGvS6TR43T6PG0Ceetpjk/cDjVfXBTfe/ArgOOK+q7h33PAa6pHEG9YqnEbrDHrdxKGaS8J+mbU1bTPLsJM9dvw28HrhzU5uTgU8AP9ckzCXtLKOGRkZtW5+WuB6e0w7RYX8pDKpvEYdgxvbQk5xCr+cNvYOoH62qDyTZB1BVB5JcA/wMsL6s49OmNm5mD13Sus0HJScZbmnadtx0xGHP03Q+/Ky0OuTSFgNd0jBtDm9MEsijPhw2D8tMOhzUFgNd0lJo+6Serazd0vQx8zoByUCXtLC22/Od5EDosOce9BzD/kpo0nbQ49pioEtaWNsZo17EA5ObGeiSdoS2A3nc6f6j/gpY3z6ND4k2Q91Al9TIPBbTarqI1iRL5U6znq1qa3+2tTiXpA5apEu8rRsUfpMG++YPpyYfVtPcF7NYL8ZAl3aQRZtTvZXXneTg57APhkGn/288gWkZxuYHcchFWjDTXDNkXI94HuuVjJo1smhrqLdhu+/JMXRJ3zXqpJ22grTJ88z6+p6jFuWah62+TwNd0kyNmlnSZP72uDHuQWPiTf/6WBTTCHTH0KW+RVketQtGHdTcymMHbRvX017kf8dp1WYPXeoz0OdjFqG8qL33rbw/e+hSAwb57Az78BzV6x53Kv526piXtqcy2kOXNDWTrnuy2VYet5WDnpMsvDUtTYPdg6KSFso0A7PttcxnGe5NQt0hF0kLZRaXdZvXeuXbsd0hGANd0ty0fdCzyfTFLjPQJc1c2wcDm051nNQshobafI2xF4kGSPJAkjuSHE7ytIHvJC9N8oUkf5nkl1qrTtJSG3aB5bZ7yxtfZ1o98Y1rvbRhGh8Wk/TQX1dVjwzZ9ihwKfCT265IUmdMI7gned6m7bbyYbDVmTuDXnfzfVs+i7TJLJckDwArIwJ9vd37gcer6oPjntNZLpKmpa2TlRZx+uK2py0muR/4JlDAb1XVwSHt3s+IQE+yF9gLcPLJJ7/ywQcfHPvakrRu1MHPUXPJ19tP4+IXbS72Natpi+dU1cNJTgRuTHJPVd08QZ0A9D8IDkKvhz7p4zWZeVx9RpqVzUE6KlTXt21lmuQ0xuYHHRBt4/knPrFoVC/cIZfpGfafalaX5ZIWSZtDIaOuQ7rV35txV0vazoHibQ25JHk28Iyq+vP+7RuBf1VVNwxo+34M9FZN+z+utGy2+jsx7CBm0zNFt7K2+yTaCPQmQy4nAdclWW//0aq6Ick+gKo6kOSFwCpwHPBkkncDL6+qxxpVqJlwCEZdtJX/y1udzdK0TdO12ts2NtCr6ghw5oD7D2y4/X+A3e2WtnNN69N+3ivLSW1Yts7ILJcg8EzRBbPVleim0V5SO2b1u+dqi0vCiy9IgtFj6I1O/df8tX3asXaORRhqG7YEgNploEtzNIugsyOwcxjo0hztlLD1L8zZMNClOTPo1BYDXZI6wkCXpI4w0CWpIwx0SeoIA12SBljGufMGuiR1hIEuSRts7pkvUy/dxbkkaYPN66Yv03kC9tAlaYhlCnMw0CWpMwx0SeoIA12SOqJRoCd5IMkdSQ4nedpVKdKzP8l9Sb6Y5G+1X6q2ahnn00qa3CSzXF5XVY8M2XYecGr/61XAb/a/awF59SOpm9qatvgm4MPVu57drUmen+RFVfX1lp5/x5hG2G5+LoNc6qamY+gFfCbJoSR7B2z/QeBrG35+qH/fUyTZm2Q1yera2trk1UqShmraQz+nqh5OciJwY5J7qurmDdsz4DFPu/p0VR0EDkLvItETV7sD2HuWtFWNeuhV9XD/+1HgOuDsTU0eAn5ow8+7gYfbKFCS1MzYQE/y7CTPXb8NvB64c1OzTwFv6892+dvAnzl+rkXhLB/tFE2GXE4Crkuy3v6jVXVDkn0AVXUA+DTwRuA+4C+At0+nXGk4Q1s73dhAr6ojwJkD7j+w4XYB72y3NKldmwPf4xXqGs8UVWccufTyiULaHr26Jr3O9eytrKzU6urTTjqVWrGVsLbHrmWQ5FBVrQza5nro6oQ2etueQatl55CLtIFhrmVmD12dYBBL9tAlqTMMdEnqCIdcNBXzPsC48SCpwzHaKQx0dVKXQnzQDJ4uvT+1x3no0oKadCqmIb8zjJqH7hi61BGe+SqHXKQFtd7jHhfU9sy1zh66tOBGBbZhro3soUtLwOBWE/bQJakjDHRJ6ggDXZI6onGgJzkmyW1Jrh+w7fgk1yX5YpL/meSMdsuUJI0zSQ/9MuDuIdv+KXC4ql4BvA34D9stTJI0mUaBnmQ3cD5wzZAmLwc+C1BV9wB7kpzUSoWSpEaa9tCvBq4Anhyy/XbgpwGSnA28BNi93eIkSc2NDfQkFwBHq+rQiGa/Ahyf5DDwLuA24IkBz7U3yWqS1bW1tS2WLEkaZOziXEn+LfBz9AL6+4HjgE9U1cVD2ge4H3hFVT027HldnEuSJretxbmq6sqq2l1Ve4CLgJs2h3mS5yd5Zv/HfwzcPCrMJUnt2/I89CT7kuzr//gy4K4k9wDn0ZsRI0naYNorYroeuiQtEddDl6QdwECXpI4w0CWpIwx0SeoIA12SOsJAl6SOMNAlqSMMdEnqCANdWkKn7L9q6mcdavkcO+8CpGkYFnZHLr18xpVMR1feh9ploGvp7cSe6vp7Nti1kYGupbYTwxwMcg1moGspbSXIDUF1nYGupbJTe+RSEwa6lsqgXvbm8WRDXzuV66GrszxwqC4atR66PXR1lkGuncYTiySpIxoHepJjktyW5PoB256X5A+S3J7kriRvb7dMSdI4k/TQLwPuHrLtncCXqupM4LXAVUmeuc3aJEkTaBToSXYD5wPXDGlSwHOTBHgO8CjwRCsVSpIaadpDvxq4AnhyyPZfB14GPAzcAVxWVU9rm2RvktUkq2tra1soV5I0zNhAT3IBcLSqDo1o9hPAYeDFwFnAryc5bnOjqjpYVStVtbJr166tVSxJGqhJD/0c4MIkDwDXAucm+cimNm8HPlE99wH3Ay9ttVJJ0khjA72qrqyq3VW1B7gIuKmqLt7U7KvAjwMkOQk4HTjScq2SpBG2fGJRkn0AVXUA+NfAf0pyBxDgl6vqkXZKlCQ1MVGgV9XngM/1bx/YcP/DwOvbLEySNBnPFJWkjjDQJakjDHRJ6ggDXZI6wkCXpI4w0CWpIwx0SeoIA12SOsJAl6SOMNAlqSMMdElbcsr+qzhl/1XzLkMbbHlxLkk725FLL593CdrEHrokdYSBLkkdYaBLUkcY6JLUEQa6JHWEgS5JHdE40JMck+S2JNcP2PaeJIf7X3cm+U6SF7RbqiRplEl66JcBdw/aUFW/VlVnVdVZwJXA56vq0Rbq04LxZBJpcTU6sSjJbuB84APAL45p/hbgY9usSwvKk0mkxdW0h341cAXw5KhGSZ4FvAH4vSHb9yZZTbK6trY2SZ2SpDHGBnqSC4CjVXWowfP9feCPhw23VNXBqlqpqpVdu3ZNWKokaZQmPfRzgAuTPABcC5yb5CND2l6Ewy2SNBdjA72qrqyq3VW1h15g31RVF29ul+R5wGuAT7ZepSRprC3PQ0+yL8m+DXf9FPCZqvr29suSJE0qVTWXF15ZWanV1dW5vLYkLaskh6pqZdA2zxTVTDh/XZo+L3ChsdaDeNI56IMC/JT9VzmXXZoSh1zUSNPe9ZFLLx/Z1jCXtmfUkIs9dLVqXPAb9tL0OIausWY19u04u7Q99tA10iRDLZM+RlK7DHSNtJVhkPXHGOzSbBnomppxB0g3tpO0fQa6psqwlmbHg6KS1BEGuiR1hIEuSR1hoEtSRxjoktQRBrokdYSBLkkdYaBLUkcY6JLUEY0DPckxSW5Lcv2Q7a9NcjjJXUk+316JkqQmJjn1/zLgbuC4zRuSPB/4EPCGqvpqkhPbKU+S1FSjHnqS3cD5wDVDmvxD4BNV9VWAqjraTnmSpKaaDrlcDVwBPDlk+2nA8Uk+l+RQkrcNapRkb5LVJKtra2uTVytJGmpsoCe5ADhaVYdGNDsWeCW9XvxPAP8iyWmbG1XVwapaqaqVXbt2bbVmSdIATcbQzwEuTPJG4PuB45J8pKou3tDmIeCRqvo28O0kNwNnAve2XrEkaaCxPfSqurKqdlfVHuAi4KZNYQ7wSeDvJjk2ybOAV9E7gCpJmpEtX+AiyT6AqjpQVXcnuQH4Ir1x9muq6s6WapQkNZCqmssLr6ys1Orq6lxeW5KWVZJDVbUyaJtnikpSRxjoktQRBrokdYSBLkkdYaBLUkcY6JLUEQa6JHWEgS4tsFP2X8Up+6+adxlaEls+U1TS9B259PJ5l6AlYg9dmrJp9bIXpfe+KHXIHrq0tBal974odchAl6bOwNOsOOQiSR1hoEtSRxjoktQRBrokdYSBLkkdYaBLUkc0DvQkxyS5Lcn1A7a9NsmfJTnc/3pfu2VKksaZZB76ZcDdwHFDtt9SVRdsvySpOzaeQel8dE1bo0BPshs4H/gA8ItTrUjqEENcs9R0yOVq4ArgyRFt/k6S25P8YZK/PqhBkr1JVpOsrq2tTViqJGmUsYGe5ALgaFUdGtHsfwEvqaozgf8I/P6gRlV1sKpWqmpl165dW6lXkjREkx76OcCFSR4ArgXOTfKRjQ2q6rGqerx/+9PAX0tyQtvFSpKGGxvoVXVlVe2uqj3ARcBNVXXxxjZJXpgk/dtn95/3G1OoV5I0xJZXW0yyD6CqDgD/AHhHkieA/wtcVFXVTomSpCYyr9xdWVmp1dXVuby2JC2rJIeqamXQNs8UlaSOMNAlqSMMdEnqCANd2uG8yHN3GOiS1BFeJFra4VxvpjvsoUtSRxjoktQRBrokdYSBLkkdYaBLUkcY6JLUEXNbnCvJGvDgXF58MZwAPDLvIpaA+2k891EzXdlPL6mqgVcImlug73RJVoetmKbvcT+N5z5qZifsJ4dcJKkjDHRJ6ggDfX4OzruAJeF+Gs991Ezn95Nj6JLUEfbQJakjDHRJ6ggDfYqSvDnJXUmeTDJ0ulSSNyT5cpL7krx30scvuxb20wuS3JjkK/3vx8+m8tlp+h6TXJbkzv7+fPeG+89M8oUkdyT5gyTHzaz4GWphP52V5NYkh5OsJjl7ZsW3wECfrjuBnwZuHtYgyTHAbwDnAS8H3pLk5U0f3xHb3U/vBT5bVacCn+3/3DVj32OSM4BfAM4GzgQuSHJqf/M1wHur6m8A1wHvmUnVs7fd/fSrwL+sqrOA9/V/XhoG+hRV1d1V9eUxzc4G7quqI1X1V8C1wJsmePzS2+5+6n//3f7t3wV+ciqFzleT9/gy4Naq+ouqegL4PPBT/W2n870PzBuBn5leqXO13f1UwPpfL88DHp5eqe0z0OfvB4Gvbfj5of59eqpR++mkqvo6QP/7iTOubRaavMc7gVcn+YEkzwLeCPzQhm0X9m+/ecP9XbPd/fRu4NeSfA34IHDl9Etuj5eg26Yk/xV44YBN/6yqPtnkKQbc17m5pO6n8UbtoyaPr6q7k/w7ej3wx4HbgSf6m38e2J/kfcCngL/afsXzMeX99A7gn1TV7yX5WeC3gb+3/apnw0Dfpqra7j/2Qzy1t7SbJfszr4kp76c/TfKiqvp6khcBR7f5WnMxah8lafQeq+q36YUQSf4Nvf1GVd0DvL5//2nA+S2XPzPT3E/APwIu69/+OL1jD0vDIZf5+xPg1CQ/nOSZwEX0elB6qlH76VP0fhHpf2/S4182jd5jkhP730+md6D5Y5vufwbwz4EDU653Xra1n+h1El7Tv30u8JWpVToNVeXXlL7oHWh5CPhL4E+BP+rf/2Lg0xvavRG4F/jf9IYgRj6+a18t7KcfoDej4Sv97y+Y93uawj4a+B4H7KNbgC/RG0b48Q33X9bfd/cCv0L/LPGufbWwn34MONS//38Ar5z3e5rky1P/JakjHHKRpI4w0CWpIwx0SeoIA12SOsJAl6QZSPI7SY4mubNB23/fXyDscJJ7k3yr0Ws4y0WSpi/Jq+mdmfrhqjpjgse9C/ibVfXz49raQ5ekGaiqm4FHN96X5EeS3JDkUJJbkrx0wEPfwvdOfBrJU/8laX4OAvuq6itJXgV8iN4ZqgAkeQnww8BNTZ7MQJekOUjyHOBHgY8n31177vs2NbsI+C9V9Z0mz2mgS9J8PAP4VvUupjHMRcA7J3lCSdKMVdVjwP1J3gyQnjPXtyc5HTge+ELT5zTQJWkGknyMXjifnuShJJcAbwUuSXI7cBffuwoX9A6GXlsTTEV02qIkdYQ9dEnqCANdkjrCQJekjjDQJakjDHRJ6ggDXZI6wkCXpI74/wpWntDpN+gnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for layer in inputLayers:\n",
    "    tiff = rasterio.open('data/' + layer['id'] + '_P' + str(layer['percentile']) + '_C' + str(layer['cutoff']) + '.tiff')\n",
    "    crsEpsg = tiff.crs.to_string()\n",
    "\n",
    "    # Set up tiff processing        \n",
    "    latRange = numpy.arange(tiff.bounds.top, tiff.bounds.bottom, -tiff.res[1])\n",
    "    lonRange = numpy.arange(tiff.bounds.left, tiff.bounds.right, tiff.res[0])\n",
    "\n",
    "    combinedGeojson = {\n",
    "        'crs': {'type':'name','properties':{'name': crsEpsg}},\n",
    "        'type': 'FeatureCollection',\n",
    "        'features': []\n",
    "    }    \n",
    "\n",
    "    figure = plt.figure().add_subplot(111)\n",
    "    values =  tiff.read(1)\n",
    "    contourf = figure.contourf(lonRange, latRange, values, levels = accessBins)\n",
    "    contourJson = geojson.loads(geojsoncontour.contourf_to_geojson(contourf = contourf, ndigits = 3))\n",
    "\n",
    "    for i in range(len(contourJson['features'])):\n",
    "        feature = contourJson.features[i]\n",
    "        minReachable = list(accessBins)[i]\n",
    "        \n",
    "        # Replace default properties from matplotlib contourf\n",
    "        feature.properties = {\n",
    "            'minReachable': minReachable,\n",
    "        }\n",
    "        \n",
    "\n",
    "        # The next two steps remove shapes with no coordinates, zero-area artifacts, \n",
    "        # and self-intersections\n",
    "        nextContour = shape({\n",
    "            'type': 'MultiPolygon', \n",
    "            'coordinates': list(filter(\n",
    "                lambda coords: (len(coords) > 0), \n",
    "                feature.geometry.coordinates))\n",
    "        })\n",
    "        \n",
    "        nextContour = MultiPolygon(list(filter(lambda poly: (poly.area > 0), nextContour))).buffer(0)\n",
    "\n",
    "        # countourf returns the area between the requested intervals (e.g. a bandreachable in \n",
    "        # between 15 and 30 minutes). Union this band with the previous ones to show the cumulative\n",
    "        # area reachable.\n",
    "        if i > 0:\n",
    "            nextContour = nextContour.union(shape(isochroneJson['features'][i - 1]['geometry']))\n",
    "\n",
    "        feature.geometry = json.loads(shapely_geojson.dumps(nextContour))\n",
    "        combinedGeojson['features'].append(feature)\n",
    "        \n",
    "    with open('data/' + layer['id'] + '_P' + str(layer['percentile']) + '_C' + str(layer['cutoff']) + '.geojson', 'w') as f:\n",
    "        json.dump(combinedGeojson, f)\n",
    "    f.close()"
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}