edge_consolidation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import LineString, Polygon, MultiLineString\n",
    "from shapely.ops import triangulate, polygonize\n",
    "import shapely\n",
    "import geopandas as gpd\n",
    "from geopandas.array import from_shapely\n",
    "import numpy as np\n",
    "from scipy.spatial import Voronoi\n",
    "import pygeos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "realdf = gpd.read_file('/Users/martin/Dropbox/Python/Jupyter/edges.geojson')\n",
    "realdf.to_file('/Users/martin/Dropbox/Python/Jupyter/edges2.geojson', driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 1), (2, 1), (3, 1), (4, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<shapely.geometry.linestring.LineString at 0x128a234f0>",
      "image/svg+xml": "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"-0.2 -0.2 5.4 1.4\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,1.0)\"><polyline fill=\"none\" stroke=\"#66cc99\" stroke-width=\"0.10800000000000001\" points=\"0.0,0.0 1.0,1.0 2.0,1.0 3.0,1.0 4.0,1.0 5.0,0.0\" opacity=\"0.8\" /></g></svg>"
     },
     "metadata": {},
     "execution_count": 95
    }
   ],
   "source": [
    "top"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<shapely.geometry.linestring.LineString at 0x128a23cd0>",
      "image/svg+xml": "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"-0.2 -1.2 5.4 1.4\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,-1.0)\"><polyline fill=\"none\" stroke=\"#66cc99\" stroke-width=\"0.10800000000000001\" points=\"0.0,0.0 1.0,-1.0 2.0,-1.0 3.0,-1.0 4.0,-1.0 5.0,0.0\" opacity=\"0.8\" /></g></svg>"
     },
     "metadata": {},
     "execution_count": 96
    }
   ],
   "source": [
    "bottom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "'POINT (1.585786437626905 1)'"
     },
     "metadata": {},
     "execution_count": 12
    }
   ],
   "source": [
    "pt = top.interpolate(2)\n",
    "pt.wkt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "pt_b = bottom.interpolate(bottom.project(pt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "'POINT (0.2928932188134524 -0.2928932188134524)'"
     },
     "metadata": {},
     "execution_count": 15
    }
   ],
   "source": [
    "pt_b.wkt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = gpd.GeoSeries([top, bottom, pt, pt_b])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<matplotlib.axes._subplots.AxesSubplot at 0x1178a4b20>"
     },
     "metadata": {},
     "execution_count": 18
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAClCAYAAABfoUaZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAPaUlEQVR4nO3da4wd9X2H8edrY8QGElYUA/baxkR1rZIiMDnYQVQtaSm3RLGhTWWaW1OiFW2QUrVyC6WqlL4IUaNUFQUlbCmlpASCGkzcxsVAYkpRAmHN1caYbpCB9bqsMbVjFxPA/vXFjmG9Puu9zJy5fj/Symcue/7/Y8vPjmbnnFFEYGZm9Tej6AmYmVk+HHwzs4Zw8M3MGsLBNzNrCAffzKwhHHwzs4Y4qugJHMmJJ54YCxcuLHoaZmaVsWHDhtciYna7baUO/sKFC+nv7y96GmZmlSHppfG2ZXJKR9KtkoYlbRxnuyTdIGlA0jOSzs5iXDMzm7ysjvBvA24Ebh9n+yXAouRrGfCN5E+rmXuf3MbX1m1haNc+5nZ3seqixaxY0lP0tMyMjI7wI+Jh4PUj7LIcuD1GPAp0S5qTxdhWHvc+uY1r73mWbbv2EcC2Xfu49p5nuffJbUVPzczI7xx+D/DKqOXBZN32nMavteE9b3Ld6o3s2PPzQuexaWg3b+8/9LOZ9r29n1X/+jS3/WhrMZOyypvbfQxfuewMut93dNFTqby8gq8269p+apukXqAXYMGCBZ2cUy0M73mTK/oeZWjXm5xz2gmFzmVs7Eev/0DXrJxnY3UQETy4eZiXdj7GHV9Y5uinlFfwB4H5o5bnAUPtdoyIPqAPoNVq+aM8j2B07G/7/Dks++AvFDqf8776Q7bt2nfY+p7uLm7/g6UFzMjq4KEtw/R+awOfusXRTyuvN16tAT6bXK3zEWB3RPh0Tgpliz3AqosW0zVr5iHrumbNZNVFiwuakdXB+YtPou8zH+a/h/fyqVseY9cbbxU9pcrK6rLMO4EfA4slDUq6UtJVkq5KdlkLvAgMAP8A/FEW4zZVGWMPsGJJD9dffgY93V2IkSP76y8/w1fpWGqOfjZU5hugtFqt8BuvDlXW2Jvl4eDpnUUnHefTO+OQtCEiWu22+bN0KsSxt6bzkX46Dn5FOPZmIxz96XPwK8CxNzuUoz89Dn7JOfZm7Tn6U+fgl5hjb3Zkjv7UOPgl5dibTY6jP3kOfgk59mZT4+hPjoNfMo692fQ4+hNz8EvEsTdLx9E/Mge/JBx7s2w4+uNz8EvAsTfLlqPfnoNfMMferDMc/cM5+AVy7M06y9E/lINfEMfeLB+O/nsc/AI49mb5cvRHOPg5c+zNiuHoO/i5cuzNitX06Dv4OXHszcqhydF38HPg2JuVS1Oj7+B3mGNvVk5NjL6D30GOvVm5NS36Dn6HOPZm1dCk6GcSfEkXS9oiaUDSNW22ny9pt6Snkq+/ymLcsnLszaqlKdFPHXxJM4GbgEuA04ErJJ3eZtf/ioizkq+/TjtuWTn2ZtXUhOhncYS/FBiIiBcj4i3gLmB5Bs9bOY69WbXVPfpZBL8HeGXU8mCybqxzJT0t6T8kfWi8J5PUK6lfUv+OHTsymF4+HHuzeqhz9LMIvtqsizHLTwCnRsSZwN8D9473ZBHRFxGtiGjNnj07g+l1nmNvVi91jX4WwR8E5o9angcMjd4hIn4WEXuTx2uBWZJOzGDswjn2ZvVUx+hnEfzHgUWSTpN0NLASWDN6B0mnSFLyeGky7s4Mxi6UY29Wb3WLfurgR8Q7wNXAOmAzcHdEbJJ0laSrkt1+B9go6WngBmBlRIw97VMpjr1ZM9Qp+ipzd1utVvT39xc9jcM49mbN89CWYXq/tYFFJx3HHV9YRvf7ji56Sm1J2hARrXbb/E7bKXLszZqpDkf6Dv4UOPZmzVb16Dv4k+TYmxlUO/oO/iQ49mY2WlWj7+BPwLE3s3aqGH0H/wgcezM7kqpF38Efh2NvZpNRpeg7+G049mY2FVWJvoM/hmNvZtNRheg7+KM49maWRtmj7+AnHHszy0KZo+/g49ibWbbKGv3GB9+xN7NOKGP0Gx18x97MOqls0W9s8B17M8tDmaLfyOA79maWp7JEv3HBd+zNrAhliH6jgu/Ym1mRio5+Y4Lv2JtZGRQZ/UYE37E3szIpKvq1D75jb2ZlVET0Mwm+pIslbZE0IOmaNtsl6YZk+zOSzs5i3Ik49mZWZnlHP3XwJc0EbgIuAU4HrpB0+pjdLgEWJV+9wDfSjjsRx97MqiDP6GdxhL8UGIiIFyPiLeAuYPmYfZYDt8eIR4FuSXMyGLstx97MqiSv6GcR/B7glVHLg8m6qe6TmTsfe5mXd77BP/1+y7E3s0o4f/FJ3Pzps3n+f/ZwzxODHRkji+CrzbqYxj4jO0q9kvol9e/YsWNaEzpqpnj7QPDIwE4i2g5jZlYqEcF/vvAa+w8EM2d05nqaozJ4jkFg/qjlecDQNPYBICL6gD6AVqs1rVr/4a//Iq+8vo8b1w8A8KcX/hJSu585ZmbFiwi+/G/PcduPtnLlr57GZ889tSPjZBH8x4FFkk4DtgErgd8bs88a4GpJdwHLgN0RsT2DsduaMUN85bIzABx9Myu1sbH/y4/9csdalTr4EfGOpKuBdcBM4NaI2CTpqmT7N4G1wKXAAPAG8Pm0407E0Tezsssz9pDNET4RsZaRqI9e981RjwP4YhZjTYWjb2ZllXfsIaPgl5mjb2ZlU0TsoQHBB0ffzMqjqNhDQ4IPjr6ZFa/I2EODgg+OvpkVp+jYQ8OCD46+meWvDLGHBgYfHH0zy09ZYg8NDT44+mbWeWWKPTQ4+ODom1nnlC320PDgg6NvZtkrY+zBwQccfTPLTlljDw7+u8ZGf8ure9i0bTfbd7/J3O4uVl20mBVLOvYR/mZWA2WOPTj4hzgY/a07/48Hnnv13fXbdu3j2nueBXD0zaytssceMrqJeZ3MmCFe3vnGYev3vb2fr63bUsCMzKzsqhB7cPDb2r77zbbrh3bty3kmZlZ2VYk9OPhtze3uart+zvHH5DwTMyuzKsUeHPy2Vl20mK5ZMw9b/6Ge432PXDMDqhd7cPDbWrGkh+svP4Oe7i4EzD3+GD7ywRN44LlX+fr9Lzj6Zg1XxdiDr9IZ14olPYdckXPgQPAXq5/1dfpmDVfV2IODP2l+c5aZVTn24OBPiaNv1lxVjz04+FPm6Js1Tx1iDw7+tDj6Zs1Rl9hDyuBLOgH4DrAQ2Ar8bkT8b5v9tgJ7gP3AOxHRSjNuGTj6ZvVXp9hD+ssyrwF+EBGLgB8ky+P5aEScVYfYH3Qw+ivPmc+N6wd8yaZZjdQt9pD+lM5y4Pzk8T8DDwF/nvI5K8VH+mb1U8fYQ/rgnxwR2wEiYrukk8bZL4D7JQVwc0T0jfeEknqBXoAFCxaknF4+HH2z+qhr7GESwZf0IHBKm03XTWGc8yJiKPmB8ICk5yPi4XY7Jj8M+gBarVZlzo84+mbVV+fYwySCHxEXjLdN0quS5iRH93OA4XGeYyj5c1jSamAp0Db4Vebom1VX3WMP6X9puwb4XPL4c8D3xu4g6VhJ7z/4GLgQ2Jhy3NLyL3LNqqcJsYf05/C/Ctwt6UrgZeCTAJLmArdExKXAycDq5C/vKODbEXFfynFLzUf6ZtXRlNhDyuBHxE7gN9usHwIuTR6/CJyZZpwqcvTNyq9JsQe/07ajHH2z8mpa7MHB7zhH36x8mhh7cPBz4eiblUdTYw8Ofm4cfbPiNTn24ODnytE3K07TYw8Ofu4cfbP8OfYjHPwCOPpm+XHs3+PgF8TRN+s8x/5QDn6BHH2zznHsD+fgF8zRN8ueY9+eg18Cjr5Zdhz78Tn4JeHom6Xn2B+Zg18ijr7Z9Dn2E3PwS8bRN5s6x35yHPwScvTNJs+xnzwHv6QcfbOJOfZT4+CXmKNvNj7Hfuoc/JJz9M0O59hPj4NfAY6+2Xsc++lz8CvC0Tdz7NNy8CvE0bcmc+zTm5HmmyV9UtImSQcktY6w38WStkgakHRNmjGb7mD0V54znxvXD/D1+18gIoqelllHOfbZSHuEvxG4HLh5vB0kzQRuAn4LGAQel7QmIp5LOXZj+UjfmsSxz06q4EfEZmCiv/ylwEBEvJjsexewHHDwUxgb/SC4+qOLCp6VWbaC4G/u2+LYZySPc/g9wCujlgeBZTmMW3ujo3/T+p9y0/qfFjwjs85w7LMxYfAlPQic0mbTdRHxvUmM0e5faNyTzpJ6gV6ABQsWTOLpm+1g9M9ZeAI79v686OmYZW7O8cfwiTPnOvYZmDD4EXFByjEGgfmjlucBQ0cYrw/oA2i1Wv5t5CTMmCF++8Pzip6GmZVcqqt0JulxYJGk0yQdDawE1uQwrpmZjZL2sszLJA0C5wLfl7QuWT9X0lqAiHgHuBpYB2wG7o6ITemmbWZmU5X2Kp3VwOo264eAS0ctrwXWphnLzMzSyeOUjpmZlYDK/C5NSTuAl6b57ScCr2U4nSrwa66/pr1e8GueqlMjYna7DaUOfhqS+iNi3I97qCO/5vpr2usFv+Ys+ZSOmVlDOPhmZg1R5+D3FT2BAvg111/TXi/4NWemtufwzczsUHU+wjczs1FqF/wm3mxF0q2ShiVtLHoueZA0X9J6SZuTG/B8qeg5dZqkYyT9RNLTyWv+ctFzyoukmZKelPTvRc8lD5K2SnpW0lOS+jN97jqd0klutvICo262AlxR95utSPo1YC9we0T8StHz6TRJc4A5EfGEpPcDG4AVdf531shHRR4bEXslzQIeAb4UEY8WPLWOk/QnQAv4QER8vOj5dJqkrUArIjJ/70HdjvDfvdlKRLwFHLzZSq1FxMPA60XPIy8RsT0inkge72HkM5p6ip1VZ8WIvcnirOSrPkdr45A0D/gYcEvRc6mDugW/3c1Wah2CppO0EFgCPFbsTDovObXxFDAMPBARtX/NwN8BfwYcKHoiOQrgfkkbkvuDZKZuwZ/SzVas2iQdB3wX+OOI+FnR8+m0iNgfEWcxck+JpZJqffpO0seB4YjYUPRccnZeRJwNXAJ8MTllm4m6BX9KN1ux6krOY38XuCMi7il6PnmKiF3AQ8DFBU+l084DPpGc074L+A1J/1LslDov+bRhImKYkU8jXprVc9ct+L7ZSgMkv8D8R2BzRPxt0fPJg6TZkrqTx13ABcDzxc6qsyLi2oiYFxELGfm//MOI+HTB0+ooSccmFyIg6VjgQiCzq+9qFfym3mxF0p3Aj4HFkgYlXVn0nDrsPOAzjBzxPZV8XTrRN1XcHGC9pGcYObB5ICIacZliw5wMPCLpaeAnwPcj4r6snrxWl2Wamdn4anWEb2Zm43PwzcwawsE3M2sIB9/MrCEcfDOzhnDwzcwawsE3M2sIB9/MrCH+HzVCjXg6wJ34AAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "source": [
    "s.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = LineString([pt, pt_b]).centroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = gpd.GeoSeries([top, bottom, pt, pt_b, c])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<matplotlib.axes._subplots.AxesSubplot at 0x129613a30>"
     },
     "metadata": {},
     "execution_count": 24
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAClCAYAAABfoUaZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAP70lEQVR4nO3dbZBVhX3H8e+PFceNJu5YUWEBIVPK1JRRzBXi0GlNm4qaTEDadLB5amqGsY0z6bRDK7XTmfRFdOqk07E6iVtrranVOg0QWqmoEWudROMiKiBiNw7q7lIXsUugrhHh3xd70WW5C7t7zz2Pv8/MDveec7jnf2H4cubs2XsUEZiZWflNyXoAMzNLh4NvZlYRDr6ZWUU4+GZmFeHgm5lVhINvZlYRp2Q9wImcffbZMWfOnKzHMDMrjC1btrwZEdMarct18OfMmUN3d3fWY5iZFYakV8dal8gpHUl3SRqQtH2M9ZJ0q6QeSS9IujiJ/ZqZ2fgldYR/N3AbcM8Y668E5tW/FgPfrv9qJbN+ax+3bNpF/+AQMzraWb10PssXdmY9lpmR0BF+RDwBvHWCTZYB98Swp4AOSdOT2Lflx/qtfaxZu42+wSEC6BscYs3abazf2pf1aGZGeufwO4HXRzzvrS/bk9L+S23gwDvcuG47ew/8LNM5dvTv59DhYz+baejQYVb/6/Pc/cPd2QxlhTej4zS+efUCOj50atajFF5awVeDZQ0/tU3SKmAVwOzZs1s5UykMHHiHa7qeon/wHS6Ze1ams4yO/cjlH2mfmvI0VgYRwaM7B3h139Pc+9XFjn6T0gp+LzBrxPOZQH+jDSOiC+gCqNVq/ijPExgZ+7u/cgmLP/pzmc6z5ObH6BscOm55Z0c79/zeogwmsjJ4fNcAq767hc/f6eg3K60fvNoAfKl+tc4ngP0R4dM5Tchb7AFWL51P+9S2Y5a1T21j9dL5GU1kZXDZ/HPo+uLH+e+Bg3z+zqcZfPvdrEcqrKQuy7wP+BEwX1KvpGslXSfpuvomG4FXgB7g74A/SGK/VZXH2AMsX9jJTSsW0NnRjhg+sr9pxQJfpWNNc/SToTzfAKVWq4V/8OpYeY29WRqOnt6Zd84ZPr0zBklbIqLWaJ0/S6dAHHurOh/pN8fBLwjH3myYoz95Dn4BOPZmx3L0J8fBzznH3qwxR3/iHPwcc+zNTszRnxgHP6cce7PxcfTHz8HPIcfebGIc/fFx8HPGsTebHEf/5Bz8HHHszZrj6J+Yg58Tjr1ZMhz9sTn4OeDYmyXL0W/Mwc+YY2/WGo7+8Rz8DDn2Zq3l6B/Lwc+IY2+WDkf/Aw5+Bhx7s3Q5+sMc/JQ59mbZcPQd/FQ59mbZqnr0HfyUOPZm+VDl6Dv4KXDszfKlqtF38FvMsTfLpypG38FvIcfeLN+qFn0Hv0Uce7NiqFL0Ewm+pCsk7ZLUI+mGBusvk7Rf0nP1r79IYr955dibFUtVot908CW1AbcDVwIXANdIuqDBpv8VERfVv/6y2f3mlWNvVkxViH4SR/iLgJ6IeCUi3gXuB5Yl8LqF49ibFVvZo59E8DuB10c8760vG+1SSc9L+g9JHxvrxSStktQtqXvv3r0JjJcOx96sHMoc/SSCrwbLYtTzZ4HzI+JC4G+B9WO9WER0RUQtImrTpk1LYLzWc+zNyqWs0U8i+L3ArBHPZwL9IzeIiJ9GxMH6443AVElnJ7DvzDn2ZuVUxugnEfxngHmS5ko6FVgJbBi5gaTzJKn+eFF9v/sS2HemHHuzcitb9JsOfkS8B1wPbAJ2Ag9ExA5J10m6rr7ZbwHbJT0P3AqsjIjRp30KxbE3q4YyRV957m6tVovu7u6sxziOY29WPY/vGmDVd7cw75wzuPeri+n40KlZj9SQpC0RUWu0zj9pO0F5iP36rX0sufkx5t7wIEtufoz1W/tSn8GsaspwpO/gT0BeYr9m7Tb6BocIoG9wiDVrtzn6ZikoevQd/HHKQ+wBbtm0i6FDh49ZNnToMLds2pXJPGZVU+ToO/jjkJfYA/QPDk1ouZklr6jRd/BPIk+xB5jR0T6h5WbWGkWMvoN/AnmLPcDqpfNpn9p2zLL2qW2sXjo/o4nMqqto0Xfwx5DH2AMsX9jJTSsW0NnRjoDOjnZuWrGA5QsbfXyRmbVakaLv6/AbyGvszSy/8nKdvq/DnwDH3swmowhH+g7+CI69mTUj79F38OscezNLQp6j7+Dj2JtZsvIa/coH37E3s1bIY/QrHXzH3sxaKW/Rr2zwHXszS0Oeol/J4Dv2ZpamvES/csF37M0sC3mIfqWC79ibWZayjn5lgu/Ym1keZBn9SgTfsTezPMkq+qUPvmNvZnmURfQTCb6kKyTtktQj6YYG6yXp1vr6FyRdnMR+T8axN7M8Szv6TQdfUhtwO3AlcAFwjaQLRm12JTCv/rUK+Haz+z0Zx97MiiDN6CdxhL8I6ImIVyLiXeB+YNmobZYB98Swp4AOSdMT2HdDjr2ZFUla0U8i+J3A6yOe99aXTXSbxNz39Gu8tu9t/uF3a469mRXCZfPP4Y4vXMxL/3OAtc/2tmQfSQRfDZaNvo3WeLYZ3lBaJalbUvfevXsnNdApbeLQkeDJnn3k+Y5eZmZHRQT/+fKbHD4StE1pzfU0pyTwGr3ArBHPZwL9k9gGgIjoArpg+BaHkxno93/153n9rSFu29wDwB9f/gtIjf7PMTPLXkTwjX97kbt/uJtrf3kuX7r0/JbsJ4ngPwPMkzQX6ANWAr8zapsNwPWS7gcWA/sjYk8C+25oyhTxzasXADj6ZpZro2P/55/+xZa1qungR8R7kq4HNgFtwF0RsUPSdfX13wE2AlcBPcDbwFea3e/JOPpmlndpxh6SOcInIjYyHPWRy74z4nEAX0tiXxPh6JtZXqUde0go+Hnm6JtZ3mQRe6hA8MHRN7P8yCr2UJHgg6NvZtnLMvZQoeCDo29m2ck69lCx4IOjb2bpy0PsoYLBB0ffzNKTl9hDRYMPjr6ZtV6eYg8VDj44+mbWOnmLPVQ8+ODom1ny8hh7cPABR9/MkpPX2IOD/77R0d/1xgF29O1nz/53mNHRzuql81m+sGUf4W9mJZDn2IODf4yj0d+97/945MU33l/eNzjEmrXbABx9M2so77GHhG5iXiZTpojX9r193PKhQ4e5ZdOuDCYys7wrQuzBwW9oz/53Gi7vHxxKeRIzy7uixB4c/IZmdLQ3XD79zNNSnsTM8qxIsQcHv6HVS+fTPrXtuOUf6zzT98g1M6B4sQcHv6HlCzu5acUCOjvaETDjzNP4xEfP4pEX3+BbD7/s6JtVXBFjD75KZ0zLF3Yec0XOkSPBn63b5uv0zSquqLEHB3/c/MNZZlbk2IODPyGOvll1FT324OBPmKNvVj1liD04+JPi6JtVR1liD00GX9JZwL8Ac4DdwG9HxP822G43cAA4DLwXEbVm9psHjr5Z+ZUp9tD8ZZk3AD+IiHnAD+rPx/LJiLioDLE/6mj0V14yi9s29/iSTbMSKVvsoflTOsuAy+qP/xF4HPjTJl+zUHykb1Y+ZYw9NB/8cyNiD0BE7JF0zhjbBfCwpADuiIiusV5Q0ipgFcDs2bObHC8djr5ZeZQ19jCO4Et6FDivwaobJ7CfJRHRX/8P4RFJL0XEE402rP9n0AVQq9UKc37E0TcrvjLHHsYR/Ij41FjrJL0haXr96H46MDDGa/TXfx2QtA5YBDQMfpE5+mbFVfbYQ/PftN0AfLn++MvA90dvIOl0SR8++hi4HNje5H5zy9/INSueKsQemj+HfzPwgKRrgdeAzwFImgHcGRFXAecC6+p/eKcA/xwRDzW531zzkb5ZcVQl9tBk8CNiH/DrDZb3A1fVH78CXNjMforI0TfLvyrFHvyTti3l6JvlV9ViDw5+yzn6ZvlTxdiDg58KR98sP6oae3DwU+Pom2WvyrEHBz9Vjr5Zdqoee3DwU+fom6XPsR/m4GfA0TdLj2P/AQc/I46+Wes59sdy8DPk6Ju1jmN/PAc/Y46+WfIc+8Yc/Bxw9M2S49iPzcHPCUffrHmO/Yk5+Dni6JtNnmN/cg5+zjj6ZhPn2I+Pg59Djr7Z+Dn24+fg55Sjb3Zyjv3EOPg55uibjc2xnzgHP+ccfbPjOfaT4+AXgKNv9gHHfvIc/IJw9M0c+2Y5+AXi6FuVOfbNm9LMb5b0OUk7JB2RVDvBdldI2iWpR9INzeyz6o5Gf+Uls7htcw/fevhlIiLrscxayrFPRrNH+NuBFcAdY20gqQ24HfgNoBd4RtKGiHixyX1Xlo/0rUoc++Q0FfyI2Amc7A9/EdATEa/Ut70fWAY4+E0YHf0guP6T8zKeyixZQfBXD+1y7BOSxjn8TuD1Ec97gcUp7Lf0Rkb/9s0/4fbNP8l4IrPWcOyTcdLgS3oUOK/Bqhsj4vvj2Eejv6ExTzpLWgWsApg9e/Y4Xr7ajkb/kjlnsffgz7Iexyxx0888jc9eOMOxT8BJgx8Rn2pyH73ArBHPZwL9J9hfF9AFUKvV/N3IcZgyRfzmx2dmPYaZ5VxTV+mM0zPAPElzJZ0KrAQ2pLBfMzMbodnLMq+W1AtcCjwoaVN9+QxJGwEi4j3gemATsBN4ICJ2NDe2mZlNVLNX6awD1jVY3g9cNeL5RmBjM/syM7PmpHFKx8zMckB5/ilNSXuBVyf5288G3kxwnCLwey6/qr1f8HueqPMjYlqjFbkOfjMkdUfEmB/3UEZ+z+VXtfcLfs9J8ikdM7OKcPDNzCqizMHvynqADPg9l1/V3i/4PSemtOfwzczsWGU+wjczsxFKF/wq3mxF0l2SBiRtz3qWNEiaJWmzpJ31G/B8PeuZWk3SaZJ+LOn5+nv+RtYzpUVSm6Stkv4961nSIGm3pG2SnpPUnehrl+mUTv1mKy8z4mYrwDVlv9mKpF8BDgL3RMQvZT1Pq0maDkyPiGclfRjYAiwv89+zhj8q8vSIOChpKvAk8PWIeCrj0VpO0h8BNeAjEfGZrOdpNUm7gVpEJP6zB2U7wn//ZisR8S5w9GYrpRYRTwBvZT1HWiJiT0Q8W398gOHPaOrMdqrWimEH60+n1r/Kc7Q2BkkzgU8Dd2Y9SxmULfiNbrZS6hBUnaQ5wELg6Wwnab36qY3ngAHgkYgo/XsG/gb4E+BI1oOkKICHJW2p3x8kMWUL/oRutmLFJukM4HvAH0bET7Oep9Ui4nBEXMTwPSUWSSr16TtJnwEGImJL1rOkbElEXAxcCXytfso2EWUL/oRutmLFVT+P/T3g3ohYm/U8aYqIQeBx4IqMR2m1JcBn6+e07wd+TdI/ZTtS69U/bZiIGGD404gXJfXaZQu+b7ZSAfVvYP49sDMi/jrredIgaZqkjvrjduBTwEvZTtVaEbEmImZGxByG/y0/FhFfyHislpJ0ev1CBCSdDlwOJHb1XamCX9WbrUi6D/gRMF9Sr6Rrs56pxZYAX2T4iO+5+tdVJ/tNBTcd2CzpBYYPbB6JiEpcplgx5wJPSnoe+DHwYEQ8lNSLl+qyTDMzG1upjvDNzGxsDr6ZWUU4+GZmFeHgm5lVhINvZlYRDr6ZWUU4+GZmFeHgm5lVxP8DEnjDjCv2UJEAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "source": [
    "s.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "edges = gpd.GeoSeries(triangulate(top.union(bottom), edges=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "orig = gpd.GeoSeries([top, bottom])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<matplotlib.axes._subplots.AxesSubplot at 0x12b093970>"
     },
     "metadata": {},
     "execution_count": 40
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAClCAYAAABfoUaZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAASIUlEQVR4nO3dfZBddX3H8fdnNwlBCAYlQCAJQYwoyEBhDVCq9QmHB6eRTq0wKjqjk1Kh2vaPDh07nbEztJ3+4Tg+jCGCClNH6oxS0pqC+NBa5HGDREHApCnKNsENIHmQhJDk2z/2XLNJ7u7ee8+555x7fp/XzE727l72/C4M7/lwuburiMDMzJpvqOoDmJlZORx8M7NEOPhmZolw8M3MEuHgm5klwsE3M0vErKoPMJ3jjjsuli5dWvUxzMwGxrp1656NiAXtPlfr4C9dupTR0dGqj2FmNjAk/WKqzxXylI6kL0sal/ToFJ+XpM9K2ijpJ5LOLeK6ZmbWuaKew/8qcMk0n78UWJa9rQS+WNB1zcysQ4UEPyJ+CDw/zV1WALfGhPuB+ZIWFnFtMzPrTFnP4Z8MPD3p9lj2sS0lXb/RxnfsZsOvdvBPdz7J3NnDVR+nNC/v208EzJmVzovNXtq7n2HBrOF0HvNfXPw6Xn/iPOa/Yk7VRxl4ZQVfbT7W9qe2SVrJxNM+LFmypJ9naoTxHbu58sb72PTsiwwJRk45Fqnd3+5m2R/Bj3/5AsND4rxTjq36OKV4bPM2dr60j5Pmz2Xxsa+o+jileGzzdt7/pft5/cJj+NpHz3f0cyor+GPA4km3FwGb290xIlYDqwFGRkb8ozynMb5jN1etvp9ntr/EicfM5Zntu/mzdyzjzcvaviKrcS78h++xZdtubrjiLF57/NFVH6fv3nfjfWwY38n2XXu58ePnJRG/9914Hy/sepkN4zt5/00POPo5lfXfhWuAq7NX61wAbIsIP52TQyv2W7bt5isffhNLXnUkc4bFZ767gVR+5PVJr5zLkODz399Q9VFKs2j+kfxmz15uvud/qz5KaeYfOZvVHzzvt9F/4cU9VR9pYBX1ssyvA/cBp0sak/QRSddIuia7y1pgE7AR+BLwsSKum6pDY3/+a16NJE6afyTrfvFr7tn4bNVHLMWs4SFOOGYua9ZvZuP4zqqPU4oj5wxz2VkL+cqPnkoqfG89/XhHvwBFvUrnqohYGBGzI2JRRNwcEasiYlX2+YiIayPitIg4KyL83VQ9ahf7luPnHcFJr5yb3Mo/YtZwUiv/429fltzKB0e/COn8r/4GmC72AJL42Ntem9zKv/p3T0lq5Z9+4rwkVz44+nk5+ANipti3vHdkUXIrf+WbX+OVnxBHv3cO/gDoNPYAR8waTm7lv/roI7zyE+Po98bBr7luYt/ilZ+GlFc+OPq9cPBrrJfYg1e+V346HP3uOPg11WvsW7zy05D6ygdHvxsOfg3ljT145Xvlp8XR74yDXzNFxL7FKz8NXvkTHP2ZOfg1UmTswSvfKz89jv70HPyaKDr2LV75afDKP8DRn5qDXwP9ij145Xvlp8nRb8/Br1g/Y9/ilZ8Gr/yDOfqHc/ArVEbswSvfKz9djv7BHPyKlBX7Fq/8NHjlH87RP8DBr0DZsQevfK/8tDn6Exz8klUR+xav/DR45bfn6Dv4paoy9uCV75VvqUffwS9J1bFv8cpPg1f+1FKOvoNfgrrEHrzyvfIN0o2+g99ndYp9i1d+Grzyp5di9B38Pqpj7MEr3yvfWlKLvoPfJ3WNfYtXfhq88meWUvQLCb6kSyQ9KWmjpOvbfP6tkrZJeiR7+9sirltXdY89eOV75dtkqUQ/d/AlDQNfAC4FzgCuknRGm7v+d0Sck739Xd7r1tUgxL7FKz8NXvmdSSH6RSz85cDGiNgUEXuA24AVBXzdgTNIsQevfK98O1TTo19E8E8Gnp50eyz72KEulLRe0n9IOnOqLyZppaRRSaNbt24t4HjlGLTYt3jlp8Erv3NNjn4RwVebjx1ajoeBUyLibOBzwL9O9cUiYnVEjETEyIIFCwo4Xv8NauzBK98r39ppavSLCP4YsHjS7UXA5sl3iIjtEbEze38tMFvScQVcu3KDHPsWr/w0eOV3p4nRLyL4DwHLJJ0qaQ5wJbBm8h0knShJ2fvLs+s+V8C1K9WE2INXvle+TaVp0c8d/IjYC1wH3AU8DnwjIh6TdI2ka7K7/RHwqKT1wGeBK2PAp2RTYt/ilZ8Gr/zuNSn6hbwOPyLWRsTrIuK0iLgh+9iqiFiVvf/5iDgzIs6OiAsi4t4irluVpsUevPK98m06TYm+v9O2S02MfYtXfhq88nvThOg7+F1ocuzBK98r32Yy6NF38DvU9Ni3eOWnwSu/d4McfQe/A6nEHrzyvfKtE4MafQd/BinFvsUrPw1e+fkMYvQd/GmkGHvwyvfKt04NWvQd/CmkGvsWr/w0eOXnN0jRd/DbSD324JXvlW/dGJToO/iHcOwP8MpPg1d+MQYh+g7+JI79wbzyvfKtO3WPvoOfcezb88pPg1d+ceocfQcfx346Xvle+da9ukY/+eA79jPzyk+DV36x6hj9pIPv2HfGK98r33pTt+gnG3zHvjte+Wnwyi9enaKfZPAd++555XvlW+/qEv3kgu/Y984rPw1e+f1Rh+gnFXzHPh+vfK98y6fq6CcTfMe+GF75afDK758qo59E8B374njle+VbflVFv/HBd+yL55WfBq/8/qoi+oUEX9Ilkp6UtFHS9W0+L0mfzT7/E0nnFnHdmTj2/eGV75VvxSg7+rmDL2kY+AJwKXAGcJWkMw6526XAsuxtJfDFvNediWPfX175afDK778yo1/Ewl8ObIyITRGxB7gNWHHIfVYAt8aE+4H5khYWcO22HPv+88r3yrfilBX9IoJ/MvD0pNtj2ce6vU9hvv7AL/nl8y/y5Q+NOPZ95JWfBq/8crz19OO58QPn8sQzO/jWw2N9uUYRwVebjx36b38n95m4o7RS0qik0a1bt/Z0oFnD4uV9wY/+57lkQlQFr3yvfCtORPBfP3+WffuD4aH+vJ6miK86BiyedHsRsLmH+wAQEasjYiQiRhYsWNDTgf7091/LlW9azOe+v5FP3/1zR7+PvPLT4JXfXxHBp/7tZ3z13qf4yO+dytUXntKX6xQR/IeAZZJOlTQHuBJYc8h91gBXZ6/WuQDYFhFbCrh2W0ND4u+vOMvRL4FXvle+5XNo7P/m8jcgtXtSJL/cwY+IvcB1wF3A48A3IuIxSddIuia721pgE7AR+BLwsbzXnYmjXx6v/DR45RevzNhDQa/Dj4i1EfG6iDgtIm7IPrYqIlZl70dEXJt9/qyIGC3iujNx9Mvhle+Vb90rO/aQwHfaOvrl8MpPg1d+MaqIPSQQfHD0y+CV75Vvnakq9pBI8MHRL4NXfhq88ntXZewhoeCDo99vXvle+Ta1qmMPiQUfHP1+88pPg1d+d+oQe0gw+ODo95NXvle+HawusYdEgw+Ofj955afBK39mdYo9JBx8cPT7xSvfK9/qF3tIPPjg6PeLV34avPLbq2PswcEHHP1+8Mr3yk9VXWMPDv5vOfrF88pPg1f+AXWOPTj4B3H0i+WV75WfkrrHHhz8wzj6xfLKT0PqK38QYg8OfluOfnG88r3ym25QYg8O/pQc/eJ45achxZU/SLEHB39ajn4xvPK98pto0GIPDv6MHP1ieOWnIZWVP4ixBwe/I45+fl75XvlNMaixBwe/Y45+fl75aWjyyh/k2IOD3xVHPx+vfK/8QTbosQcHv2uOfj5e+Wlo2spvQuzBwe+Jo987r3yv/EHTlNhDzuBLepWkuyVtyP48dor7PSXpp5IekTSa55p14ej3zis/DU1Y+U2KPeRf+NcD34uIZcD3sttTeVtEnBMRIzmvWRuOfm+88r3yB0HTYg/5g78CuCV7/xbgPTm/3sBx9HvjlZ+GQV35TYw95A/+CRGxBSD78/gp7hfAdyStk7Ryui8oaaWkUUmjW7duzXm8cjj63fPK98qvq6bGHjoIvqTvSnq0zduKLq5zUUScC1wKXCvpLVPdMSJWR8RIRIwsWLCgi0tUy9Hvnld+GgZp5Tc59tBB8CPinRHxxjZvdwC/krQQIPtzfIqvsTn7cxy4HVhe3EOoD0e/O175Xvl10vTYQ/6ndNYAH8re/xBwx6F3kHSUpHmt94F3AY/mvG5tOfrd8cpPQ91Xfgqxh/zB/0fgYkkbgIuz20g6SdLa7D4nAPdIWg88CHw7Iu7Med1ac/Q755XvlV+1VGIPOYMfEc9FxDsiYln25/PZxzdHxGXZ+5si4uzs7cyIuKGIg9edo985r/w01HHlpxR78Hfa9pWj3xmvfK/8KqQWe3Dw+87R74xXfhrqsvJTjD04+KVw9Gfmle+VX5ZUYw8Ofmkc/Zl55aehypWfcuzBwS+Voz89r3yv/H5KPfbg4JfO0Z+eV34ayl75jv0EB78Cjv7UvPK98ovm2B/g4FfE0Z+aV34aylj5jv3BHPwKOfrteeV75RfBsT+cg18xR789r/w09GvlO/btOfg14OgfzivfK79Xjv3UHPyacPQP55WfhiJXvmM/PQe/Rhz9g3nle+V3w7GfmYNfM47+wbzy05B35Tv2nXHwa8jRP8Ar3yt/Jo595xz8mnL0D/DKT0MvK9+x746DX2OO/gSvfK/8dhz77jn4NefoT/DKT0OnK9+x742DPwAcfa98r/wDHPveOfgDwtH3yk/FdCvfsc/HwR8gqUc/9ZW/a8++qo9Tiskrf+++/b/9uGOfX67gS3qvpMck7Zc0Ms39LpH0pKSNkq7Pc83UpR79lFf+/72wq+qjlKa18rds3w049kXJu/AfBf4Q+OFUd5A0DHwBuBQ4A7hK0hk5r5u0dtFPxeSVv3333qqPU4rWyn/uN3uSW/nPbNvN3n37HfuC5Ap+RDweEU/OcLflwMaI2BQRe4DbgBV5rmuHR3/s1y9WfaTStFb+2K93JbXyh0RyK39/wBPP7HDsCzKrhGucDDw96fYYcH4J1228VvQB1qzfzPZd23nfjfdVfKpyzJk1xK6X9/GzLek85llDYtuul/njVfcmE71ZQ2L33n2OfUFmDL6k7wIntvnUJyPijg6u0e6f0JSzTNJKYCXAkiVLOvjyaWtF/5OXv4GP3jJa9XFKc/y8IzjhmLlVH6NU5yyen1zwzl0yn5s+/CbmHTErucfeDzMGPyLemfMaY8DiSbcXAZunud5qYDXAyMhIGv+9ntPQkJg3dzb/8icXVn0UM6uxMl6W+RCwTNKpkuYAVwJrSriumZlNkvdlmVdIGgMuBL4t6a7s4ydJWgsQEXuB64C7gMeBb0TEY/mObWZm3cr1P20j4nbg9jYf3wxcNun2WmBtnmuZmVk+/k5bM7NEqM6vY5a0FfhFj3/5cUAa339/gB9z86X2eMGPuVunRMSCdp+odfDzkDQaEVP+uIcm8mNuvtQeL/gxF8lP6ZiZJcLBNzNLRJODv7rqA1TAj7n5Unu84MdcmMY+h29mZgdr8sI3M7NJGhf8FH/ZiqQvSxqX9GjVZymDpMWSfiDp8ewX8Hyi6jP1m6S5kh6UtD57zJ+q+kxlkTQs6ceS/r3qs5RB0lOSfirpEUmF/kTERj2lk/2ylZ8DFzPxQ9seAq6KiJ9VerA+k/QWYCdwa0S8serz9JukhcDCiHhY0jxgHfCeJv9z1sSPijwqInZKmg3cA3wiIu6v+Gh9J+kvgRHgmIh4d9Xn6TdJTwEjEVH49x40beEn+ctWIuKHwPNVn6MsEbElIh7O3t/BxM9oOrnaU/VXTNiZ3ZydvTVnrU1B0iLgcuCmqs/SBE0LfrtfttLoEKRO0lLgd4AHqj1J/2VPbTwCjAN3R0TjHzPwGeCvgP0z3bFBAviOpHXZ7wcpTNOC39UvW7HBJulo4JvAn0fE9qrP028RsS8izmHid0osl9Top+8kvRsYj4h1VZ+lZBdFxLlM/B7wa7OnbAvRtOB39ctWbHBlz2N/E/haRHyr6vOUKSJeAP4TuKTio/TbRcAfZM9p3wa8XdI/V3uk/st+2jARMc7ETyNeXtTXblrw/ctWEpD9D8ybgccj4tNVn6cMkhZImp+9fyTwTuCJak/VXxHx1xGxKCKWMvHv8vcj4gMVH6uvJB2VvRABSUcB7wIKe/Vdo4Kf6i9bkfR14D7gdEljkj5S9Zn67CLgg0wsvkeyt8tm+osG3ELgB5J+wsSwuTsikniZYmJOAO6RtB54EPh2RNxZ1Bdv1Msyzcxsao1a+GZmNjUH38wsEQ6+mVkiHHwzs0Q4+GZmiXDwzcwS4eCbmSXCwTczS8T/A53FmoXSadOvAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "source": [
    "edges.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "cents = edges.centroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<matplotlib.axes._subplots.AxesSubplot at 0x129685b80>"
     },
     "metadata": {},
     "execution_count": 42
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAADGCAYAAADRwTttAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAPSUlEQVR4nO3dcYikd33H8fe3m5UuVVjaW2v2Ej2FsLRG9GQ4FaGkrXbPq/SuUksqbdV/DsVAi7DFa6G2hWBgoSRVMT3SUKXVVOhlc0hkq7bFQtFmLqe9aNySprG53ZCs2o1KBnK5fvvHzKV7m9nbeW6enWc2v/cLhp35Pb99ft/73uzn5p55dp7ITCRJL34/0XQBkqTRMPAlqRAGviQVwsCXpEIY+JJUiGuaLuBK9u3blwcOHGi6DEnaM86cOfO9zJzpt22sA//AgQO02+2my5CkPSMivrvdNg/pSFIhDHxJKoSBL0mFqOUYfkTcDbwTeCozb+yzPYA7gCPAM8D7MvPBOtbWcJbOrrK4vMLaRofZ6SkW5uc4dnB/02WNLftVjf0aL3W9wv9r4PAVtr8DuKF3Ow58qqZ1NYSls6ucOHWO1Y0OCaxudDhx6hxLZ1ebLm0s2a9q7Nf4qSXwM/OrwA+uMOUo8Jns+howHRHX1rG2rt7i8gqdCxcvG+tcuMji8kpDFY03+1WN/Ro/ozqGvx94fNPj872xF4iI4xHRjoj2+vr6SIor1dpGp9J46exXNfZr/Iwq8KPPWN/PZc7Mk5nZyszWzEzf3x1QTWanpyqNl85+VWO/xs+oAv88cP2mx9cBayNaW9tYmJ9janLisrGpyQkW5ucaqmi82a9q7Nf4GdVv2p4GbomIe4A3AU9n5hMjWlvbuHS2hGdRDMZ+VWO/xk/UccWriPgccBOwD3gS+CgwCZCZd/ZOy/wE3TN5ngHen5k7fmZCq9VKP1pBkgYXEWcys9VvWy2v8DPzt3bYnsCH6lhLknR1/E1bSSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCXFPHTiLiMHAHMAHclZm3bdl+E3Af8F+9oVOZ+Wd1rL3V0tlVFpdXWNvoMDs9xcL8HMcO7t+NpSSpVrudX0MHfkRMAJ8E3g6cBx6IiNOZ+e0tU/8lM9857HpXsnR2lROnztG5cBGA1Y0OJ06dAzD0JY21UeRXHYd0DgGPZOajmfkscA9wtIb9Vra4vPJ8sy7pXLjI4vJKE+VI0sBGkV91BP5+4PFNj8/3xrZ6S0R8MyK+GBGv3W5nEXE8ItoR0V5fX69UyNpGp9K4JI2LUeRXHYEffcZyy+MHgVdl5uuBjwNL2+0sM09mZiszWzMzM5UKmZ2eqjQuSeNiFPlVR+CfB67f9Pg6YG3zhMz8YWb+uHf/fmAyIvbVsPZlFubnmJqcuGxsanKChfm5upeSpFqNIr/qOEvnAeCGiHg1sArcDLxn84SIeAXwZGZmRByi+w/N92tY+zKX3tjwLB1Je80o8isytx59uYqdRBwBbqd7WubdmXlrRHwAIDPvjIhbgA8CzwEd4MOZ+a877bfVamW73R66PkkqRUScycxW3211BP5uMfAlqZorBb6/aStJhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVIhr6thJRBwG7gAmgLsy87Yt26O3/QjwDPC+zHywjrX3iqWzqywur7C20WF2eoqF+TmOHdzfdFnWZV3WNUZ17bahAz8iJoBPAm8HzgMPRMTpzPz2pmnvAG7o3d4EfKr3tQhLZ1c5ceocnQsXAVjd6HDi1DmARp9k1mVd1jU+dY1CHYd0DgGPZOajmfkscA9wdMuco8BnsutrwHREXFvD2nvC4vLK80+uSzoXLrK4vNJQRV3WVY11VWNd46eOwN8PPL7p8fneWNU5AETE8YhoR0R7fX29hvKat7bRqTQ+KtZVjXVVY13jp47Ajz5jeRVzuoOZJzOzlZmtmZmZoYsbB7PTU5XGR8W6qrGuaqxr/NQR+OeB6zc9vg5Yu4o5L1oL83NMTU5cNjY1OcHC/FxDFXVZVzXWVY11jZ86ztJ5ALghIl4NrAI3A+/ZMuc0cEtE3EP3zdqnM/OJGtbeEy69ETRuZwVYl3VZ1/jUNQqR2ffISrWdRBwBbqd7WubdmXlrRHwAIDPv7J2W+QngMN3TMt+fme2d9ttqtbLd3nGaJKknIs5kZqvftlrOw8/M+4H7t4zduel+Ah+qYy1J0tXxN20lqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgpxzTDfHBE/DfwdcAB4DPjNzPyfPvMeA34EXASey8zWMOteydLZVRaXV1jb6DA7PcXC/BzHDu7freUkqTa7nV/DvsL/CPCVzLwB+Erv8XZ+MTPfsNthf+LUOVY3OiSwutHhxKlzLJ1d3a0lJakWo8ivYQP/KPDp3v1PA8eG3N9QFpdX6Fy4eNlY58JFFpdXGqpIkgYzivwaNvB/NjOfAOh9ffk28xL4h4g4ExHHr7TDiDgeEe2IaK+vr1cqZm2jU2lcksbFKPJrx2P4EfFl4BV9Nv1RhXXemplrEfFy4EsR8Z3M/Gq/iZl5EjgJ0Gq1ssIazE5PsdqnObPTU1V2I0kjN4r82vEVfma+LTNv7HO7D3gyIq4F6H19apt9rPW+PgXcCxyq7U+wycL8HFOTE5eNTU1OsDA/txvLSVJtRpFfwx7SOQ28t3f/vcB9WydExE9FxMsu3Qd+BXhoyHX7OnZwPx971+vYPz1FAPunp/jYu17nWTqSxt4o8isyKx01ufybI34G+DzwSuC/gXdn5g8iYha4KzOPRMRr6L6qh+4hpM9m5q2D7L/VamW73b7q+iSpNBFxZruzIYc6Dz8zvw/8cp/xNeBI7/6jwOuHWUeSNDx/01aSCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRDXDPPNEfFu4E+AnwMOZWZ7m3mHgTuACeCuzLxtmHVVn6Wzqywur7C20WF2eoqF+TmOHdzfdFljy35VY7/Gy1CBDzwEvAv4y+0mRMQE8Eng7cB54IGIOJ2Z3x5ybQ1p6ewqJ06do3PhIgCrGx1OnDoH4A9lH/arGvs1foY6pJOZD2fmyg7TDgGPZOajmfkscA9wdJh1VY/F5ZXnfxgv6Vy4yOLyTn+lZbJf1div8TOKY/j7gcc3PT7fG+srIo5HRDsi2uvr67teXMnWNjqVxktnv6qxX+Nnx8CPiC9HxEN9boO+So8+Y7nd5Mw8mZmtzGzNzMwMuISuxuz0VKXx0tmvauzX+Nkx8DPzbZl5Y5/bfQOucR64ftPj64C1qylW9VqYn2NqcuKysanJCRbm5xqqaLzZr2rs1/gZ9k3bQTwA3BARrwZWgZuB94xgXe3g0htnnkUxGPtVjf0aP5G57dGVnb854teBjwMzwAbwjcycj4hZuqdfHunNOwLcTve0zLsz89ZB9t9qtbLd7numpySpj4g4k5mtftuGeoWfmfcC9/YZXwOObHp8P3D/MGtJkobjb9pKUiEMfEkqxFDH8HdbRKwD373Kb98HfK/GcupiXdVYVzXWVc2Lsa5XZWbfc9rHOvCHERHt7d64aJJ1VWNd1VhXNaXV5SEdSSqEgS9JhXgxB/7JpgvYhnVVY13VWFc1RdX1oj2GL0m63Iv5Fb4kaRMDX5IKsacDPyIOR8RKRDwSER/ps/2miHg6Ir7Ru/3xiOq6OyKeioiHttkeEfEXvbr/PSLeOCZ1NdWv6yPinyLi4Yj4VkT8Xp85I+/ZgHWNvGcR8ZMR8W8R8c1eXX/aZ04T/RqkrkaeY721JyLibER8oc+2Rn4mB6ir3n5l5p680f0gtv8EXgO8BPgm8PNb5twEfKGB2n4BeCPw0DbbjwBfpHutgDcDXx+Tuprq17XAG3v3Xwb8R5+/y5H3bMC6Rt6zXg9e2rs/CXwdePMY9GuQuhp5jvXW/jDw2X7rN/UzOUBdtfZrL7/CH9tLJ2bmV4EfXGHKUeAz2fU1YDoirh2DuhqRmU9k5oO9+z8CHuaFV0Ubec8GrGvkej34ce/hZO+29eyLJvo1SF2NiIjrgF8F7tpmSiM/kwPUVau9HPiDXjrxLb3/Yn4xIl47mtJ2VOmyjyPWaL8i4gBwkO6rw80a7dkV6oIGetY7DPAN4CngS5k5Fv0aoC5o5jl2O/AHwP9us72p59dOdUGN/drLgT/IpRMfpPu5Eq+n+7n9S7te1WAqXfZxhBrtV0S8FPh74Pcz84dbN/f5lpH0bIe6GulZZl7MzDfQvYLcoYi4ccuURvo1QF0j71dEvBN4KjPPXGlan7Fd7deAddXar70c+DteOjEzf3jpv5jZ/Uz+yYjYN7oStzWWl31ssl8RMUk3VP82M0/1mdJIz3aqq+nnWGZuAP8MHN6yqdHn2HZ1NdSvtwK/FhGP0T30+0sR8Tdb5jTRrx3rqrtfeznwn790YkS8hO6lE09vnhARr4iI6N0/RPfP+/2RV/pCp4Hf7Z0Z8Gbg6cx8oumimupXb82/Ah7OzD/fZtrIezZIXU30LCJmImK6d38KeBvwnS3TmujXjnU10a/MPJGZ12XmAbo58Y+Z+dtbpo28X4PUVXe/RnFN212Rmc9FxC3AMv9/6cRvRcQHetvvBH4D+GBEPAd0gJuz99b3boqIz9F9d31fRJwHPkr3DaxLdd1P96yAR4BngPfvdk0D1tVIv+i+0vkd4Fzv+C/AHwKv3FRbEz0bpK4menYt8OmImKAbAJ/PzC9see430a9B6mrqOfYCY9CvQeqqtV9+tIIkFWIvH9KRJFVg4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RC/B9ZDZG9WfmIsAAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "source": [
    "cents.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "array([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9],\n       [0, 0, 0, 0, 0, 1, 1, 1, 1, 1]])"
     },
     "metadata": {},
     "execution_count": 56
    }
   ],
   "source": [
    "orig.sindex.query_bulk(edges, predicate='within')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "axis = cents.drop(orig.sindex.query_bulk(cents, predicate='intersects')[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<matplotlib.axes._subplots.AxesSubplot at 0x128c12670>"
     },
     "metadata": {},
     "execution_count": 101
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAAuCAYAAAAlUx6oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAH10lEQVR4nO3dfYwVVx3G8e9T3ERCa7YKjbzU1mCF2CqUbihKYlAxIBqxpiZoim81qEkTjQZTjFHrX41EozVagrWJTdWGRERC2iDGauMfxe4C5UWKQYN2WcJLW4rKRlt8/GPm0svd2d3ZztzuzL2/T7Lh3plz955nf7CHO2fmjGwTQgihe1022R0IIYQwuWIgCCGELhcDQQghdLkYCEIIocvFQBBCCF0uBoIQQuhypQwEklZKOiLpqKQ7M/ZL0j3p/v2SFjXtOybpgKR9kvrL6E8IIYT8XlX0G0iaAvwQeC8wCDwhabvtPzc1ex9wXfp1M3Bv+mfDu2yfKdqXEEIIE1fGJ4LFwFHbf7P9X+AhYHVLm9XAA048DvRKmlnCe4cQQiiojIFgNvB00/PBdFveNgZ+I2lA0roS+hNCCGECCh8aApSxrXXdirHaLLU9JOkqYJekp2w/NuJNkkFiHcC0adNumj9/fpE+hxBC1xkYGDhje0br9jIGgkFgoaQjwBTgKPCHjDZfk/Qm4DzwGmAo3fc2SY+mrz1GcqhpxEBgezOwGaCvr8/9/TGvHEIIEyHp71nbyxgIBoAFwHvSx2dIJo+bDQEfBq4CPg183/YJSVcAPwKWA88Cx4EtJfRphG17j7Nx5xGGzg4zq3cq61fM40M3th7BqofIUj2dkgMiS1W1M0sZA8FNwH7gPpL/1T8G3CBpNoDtTcAsYC/Jp4XzwJl0svgdwHTgV2lf/gi8roQ+XWLb3uNs2HqA4RcuAHD87DAbth4AqN1fishSPZ2SAyJLVbU7S1mTxXttv9n2XOBBYLbtTekg0GjzLdtzbb8V+Gu6TcAW2wtsX994bQl9usTGnUcu/gAbhl+4wMadR8p+q7aLLNXTKTkgslRVu7OUMRAUmSzO89rkG0jrJPVL6j99+vSEOjh0dnhC26ssslRPp+SAyFJV7c5SaCCQ9FpgPfBRSbskXQnM4aWJ4GZbmq48brRZAdyWXlW8j2SeIeu12N5su89234wZIya9xzSrd+qEtldZZKmeTskBkaWq2p2l6CeCO4FtwElgD/BVYA2wvdEgvfJ4EfAX4C3A7cB/bJ8g+aX/b+AWkrOFbmx+bVnWr5jH1J4pl2yb2jOF9Svmlf1WbRdZqqdTckBkqap2Zyk6WbwaWAbsA34AXA3cZfuQpM+lbZ4kmUw+ChwGpgE7033/I5ko3kky0Xy/7UMF+zRCYzKlE84eiCzV0yk5ILJUVbuzqMg9iyWdtd3b9Pw521e2tLkVWGn7M+nztcDNtu+Q9E3gk8A5oB/4su3nxnvfuI4ghBAmTtKA7b7W7eMeGpL0W0kHM75a1xMa9VtkbGuMPvcCc4GFwAngO2P042VPFocQQhjduIeGbC8fbZ+kk5JmpheHzQROZTQbJDlk1HBxMtn2yabv9WNgxxj9uHhlsaTTo10hl8N0koveOkFkqZ5OyQGRpaqKZLkma2PROYLtwCeAu9M/f53R5gngOklvJLlyeA3wMYDGIJK2uwU4mOdNs9bKyEtSf9ZHozqKLNXTKTkgslRVO7IUHQjuJjkt9HbgH8BHACTNAu6zvcr2i5LuIHtC+NuSFpIcKjoGfLZgf0IIIUxQoYHA9jMk5/63bh8CVjU9fxh4OKPd2iLvH0IIobhuvGfx5snuQIkiS/V0Sg6ILFVVepZCp4+GEEKov278RBBCCKFJRw4Eku6XdEpS5llIStyTrn20X9KiV7qPeeXIskzS8431miR9/ZXuYx6Srpb0qKTDkg5J+kJGm1rUJWeWutTl1ZL+JOnJNMtdGW3qUpc8WWpRF0iW55G0V9KI0+pLr4ntjvsC3kmyvtHBUfavAh4hudhtCbB7svtcIMsyYMdk9zNHjpnAovTxFaRrT9WxLjmz1KUuAi5PH/cAu4ElNa1Lniy1qEva1y8BP8/qb9k16chPBE7uefzsGE1WAw848TjQm14QVzk5stSC7RO296SP/0my7lTrQim1qEvOLLWQ/qz/lT7tSb9aJw7rUpc8WWpB0hzg/SQ3/MpSak06ciDIYTbwdNPzQWr6Dzn19vTj8COSrp/szoxH0rUkK83ubtlVu7qMkQVqUpf0EMQ+kpUBdtmubV1yZIF61OV7wFdIFubMUmpNunUgyH1DnBrYA1xjewHJCrDbJrk/Y5J0OfBL4Iu2z7XuznhJZesyTpba1MX2BdsLSZZ/WSzphpYmtalLjiyVr4ukDwCnbA+M1Sxj28uuSbcOBKOuf1Q3ts81Pg47uXCvR9L0Se5WJkk9JL84f2Z7a0aT2tRlvCx1qkuD7bPA74GVLbtqU5eG0bLUpC5LgQ9KOgY8BLxb0oMtbUqtSbcOBNuBj6cz70uA5/3Smke1Iun1kpQ+XkxS02cmt1cjpX38CXDY9ndHaVaLuuTJUqO6zJDUmz6eCiwHnmppVpe6jJulDnWxvcH2HNvXkqzN9jvbt7U0K7UmRdcaqiRJvyA5O2C6pEHgGyQTR9jeRLLcxSqSm+WcBz41OT0dX44stwKfl/QiMAyscXpaQcUsBdYCB9JjuJDc0e4NULu65MlSl7rMBH6q5E6ClwFbbO9QemOpmtUlT5a61GWEdtYkriwOIYQu162HhkIIIaRiIAghhC4XA0EIIXS5GAhCCKHLxUAQQghdLgaCEELocjEQhBBCl4uBIIQQutz/AWh0RH7n18zRAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "source": [
    "axis.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "average = LineString([start] + [p for p in axis] + [end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "av = gpd.GeoSeries(average)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<matplotlib.axes._subplots.AxesSubplot at 0x12a4d6820>"
     },
     "metadata": {},
     "execution_count": 107
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAClCAYAAABfoUaZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAOHklEQVR4nO3db4wc9X3H8ffHf6JWNC2y7IJjn3M8sCK5UaFoZRohUGjdyHaj2kVNZaKWqIKerIJk0gfUDaJSG4HyJFXkiia+UtSStnGJUojTGBPcEtFIkPhMgdo4tBdD4Gq7NqElQX1ADd8+2DFZn3fv9m5m59/v85JWt7Mzt7/f2Ph9w9zujiICMzNrvyVVT8DMzMrh4JuZJcLBNzNLhINvZpYIB9/MLBEOvplZIpZVPYG5rFy5MsbHx6uehplZYxw5cuS1iFjVb12tgz8+Ps7U1FTV0zAzawxJ3x+0rpBTOpIekHRG0tEB6yVpj6RpSc9LurqIcc3MbHhFncP/K2DzHOu3AOuz2wTw+YLGNTOzIRUS/Ih4Enh9jk22AQ9G19PApZJWFzG2mZkNp6xz+GuAV3uWZ7LHTpU0frudPg07d3a/mrXN2Bjs3QsrVlQ9k8YrK/jq81jfT22TNEH3tA/r1q0b5Zza4fRpuOEGeOUVuO66qmdjVqwI+NrX4Hvfg0OHHP2cygr+DDDWs7wWONlvw4iYBCYBOp2OP8pzLr2xf/RRuP76qmdkVryDB2H7dti0ydHPqaw3Xu0Hbs5erfOLwBsR4dM5eTj2lorNm+GRR+CFF7rRf32uXxfaXIp6WeaXgKeAD0iakXSLpJ2SdmabHABOANPAXwC/V8S4yXLsLTWOfiEKOaUTETfNsz6A24oYK3mOvaXqfPR9emfR/Fk6TeLYW+p8pJ+Lg98Ujr1Zl6O/aA5+Ezj2Zhdy9BfFwa87x96sP0d/wRz8OnPszebm6C+Ig19Xjr3ZcBz9oTn4deTYmy2Moz8UB79uHHuzxXH05+Xg14ljb5aPoz8nB78uHHuzYjj6Azn4deDYmxXL0e/Lwa+aY282Go7+RRz8Kjn2ZqPl6F/Awa+KY29WDkf/XQ5+FRx7s3I5+oCDXz7H3qwajr6DXyrH3qxaiUffwS+LY29WDwlH38Evg2NvVi+JRt/BHzXH3qyeEoy+gz9Kjr1ZvSUWfQd/VBx7s2ZIKPqFBF/SZkkvSpqWtLvP+g9LekPSs9ntj4oYt7Yce7NmSST6uYMvaSlwH7AF2ADcJGlDn03/JSKuym5/knfc2nLszZopgegXcYS/EZiOiBMR8RawD9hWwPM2j2Nv1mwtj34RwV8DvNqzPJM9NtuHJD0n6VFJPzfoySRNSJqSNHX27NkCplcSx96sHVoc/SKCrz6PxazlZ4D3R8SVwJ8Bjwx6soiYjIhORHRWrVpVwPRK4NibtUtLo19E8GeAsZ7ltcDJ3g0i4ocR8WZ2/wCwXNLKAsaunmNv1k4tjH4RwT8MrJd0haT3ADuA/b0bSLpckrL7G7Nxf1DA2NVy7M3arWXRzx38iDgH3A48BhwHHoqIY5J2StqZbfYbwFFJzwF7gB0RMfu0T7M49mZpaFH0VefudjqdmJqaqnoaF3PszdJz8CBs3w4bNsChQ7BiRdUz6kvSkYjo9Fvnd9oulGNvlqYWHOk7+Avh2JulreHRd/CH5dibGTQ6+g7+MBx7M+vV0Og7+PNx7M2snwZG38Gfi2NvZnNpWPQd/EEcezMbRoOi7+D349ib2UI0JPoO/myOvZktRgOi7+D3cuzNLI+aR9/BP8+xN7Mi1Dj6Dj449mZWrJpG38F37M1sFGoY/bSD79ib2SjVLPrpBt+xN7My1Cj6aQbfsTezMtUk+ukF37E3syrUIPppBd+xN7MqVRz9dILv2JtZHVQY/TSC79ibWZ1UFP32B9+xN7M6qiD6hQRf0mZJL0qalrS7z3pJ2pOtf17S1UWMOy/H3szqrOTo5w6+pKXAfcAWYANwk6QNszbbAqzPbhPA5/OOOy/H3syaoMToF3GEvxGYjogTEfEWsA/YNmubbcCD0fU0cKmk1QWM3Z9jb2ZNUlL0lxXwHGuAV3uWZ4BrhthmDXCqgPEvNjnJ//3HNPfs+hzHjy+H40+NZJg6+fZL1X9Oh9moXHPFiqqnUIKf4crfvZc7//xOln3xi7BrV+EjFHGErz6PxSK26W4oTUiakjR19uzZxc1o2TKWv32Onz9+GKLvMGZm9RLBVceeZtk7b8PSpSMZoogj/BlgrGd5LXByEdsAEBGTwCRAp9NZXK1374aXXuLG++/nxqvXwqc/Der3M8fMrAYi4I474Ikvwyc/CbfdNpJhigj+YWC9pCuA/wR2AB+ftc1+4HZJ++ie7nkjIkZzOgdgyRLYu7d7/557ul8dfTOro/Ox37OnG/vPfnZkrcod/Ig4J+l24DFgKfBARByTtDNb/wXgALAVmAb+F/idvOPOy9E3s7orMfZQzBE+EXGAbtR7H/tCz/0ARvP/KHNx9M2srkqOPRQU/Fpz9M2sbiqIPaQQfHD0zaw+Koo9pBJ8cPTNrHoVxh5SCj44+mZWnYpjD6kFHxx9MytfDWIPKQYfHH0zK09NYg+pBh8cfTMbvRrFHlIOPjj6ZjY6NYs9pB58cPTNrHg1jD04+F2OvpkVpaaxBwf/xxx9M8urxrEHB/9Cjr6ZLVbNYw8O/sUcfTNbqAbEHhz8/hx9MxtWQ2IPDv5gjr6ZzadBsQcHf26OvpkN0rDYg4M/P0ffzGZrYOzBwR+Oo29m5zU09uDgD8/RN7MGxx4c/IVx9M3S1fDYg4O/cI6+WXpaEHtw8BfH0TdLR0tiDzmDL2kF8PfAOPAy8JsR8d99tnsZ+BHwNnAuIjp5xq0FR9+s/VoUe8h/hL8b+KeI+Iyk3dnyHwzY9oaIeC3nePXi6Ju1V8tiD/mDvw34cHb/r4FvMjj47eTom7VPC2MP+YN/WUScAoiIU5J+dsB2AXxDUgB7I2Jy0BNKmgAmANatW5dzeiVx9M3ao6WxhyGCL+kQcHmfVXctYJxrI+Jk9gPhcUnfjYgn+22Y/TCYBOh0OrGAMarl6Js1X4tjD0MEPyI2DVon6b8krc6O7lcDZwY8x8ns6xlJDwMbgb7BbzRH36y5Wh57gCU5v38/8Ins/ieAr87eQNIlkt57/j7wEeBoznHr63z0b721G/277+7+h2Rm9ZVA7CH/OfzPAA9JugV4BfgYgKT3AfdHxFbgMuBhdf/wlgF/FxEHc45bbz7SN2uORGIPOYMfET8AfrnP4yeBrdn9E8CVecZpJEffrP4Sij34nbaj5eib1VdisQcHf/QcfbP6STD24OCXw9E3q49EYw8OfnkcfbPqJRx7cPDL5eibVSfx2IODXz5H36x8jj3g4FfD0Tcrj2P/Lge/Ko6+2eg59hdw8Kvk6JuNjmN/EQe/ao6+WfEc+74c/Dpw9M2K49gP5ODXhaNvlp9jPycHv04cfbPFc+zn5eDXjaNvtnCO/VAc/Dpy9M2G59gPzcGvK0ffbH6O/YI4+HXm6JsN5tgvmINfd46+2cUc+0Vx8JvA0Tf7Mcd+0Rz8pnD0zRz7nBz8JnH0LWWOfW5L8nyzpI9JOibpHUmdObbbLOlFSdOSducZM3nno3/rrd3o33139x+CWZs59oXIe4R/FLgR2DtoA0lLgfuAXwFmgMOS9kfECznHTpeP9C0ljn1hcgU/Io4DaO4//I3AdEScyLbdB2wDHPw8Zkc/Au66q9o5mRUtAj71Kce+IGWcw18DvNqzPANcU8K47dcb/Xvv7d7M2sixL8S8wZd0CLi8z6q7IuKrQ4zR729o4ElnSRPABMC6deuGePrEnY/+ddfB6dNVz8aseGNjsGOHY1+AeYMfEZtyjjEDjPUsrwVOzjHeJDAJ0Ol0/NvIYSxZAjffXPUszKzmcr1KZ0iHgfWSrpD0HmAHsL+Ecc3MrEfel2X+uqQZ4EPA1yU9lj3+PkkHACLiHHA78BhwHHgoIo7lm7aZmS1U3lfpPAw83Ofxk8DWnuUDwIE8Y5mZWT5lnNIxM7MaUNT4XZqSzgLfX+S3rwReK3A6TeB9br/U9he8zwv1/ohY1W9FrYOfh6SpiBj4cQ9t5H1uv9T2F7zPRfIpHTOzRDj4ZmaJaHPwJ6ueQAW8z+2X2v6C97kwrT2Hb2ZmF2rzEb6ZmfVoXfBTvNiKpAcknZF0tOq5lEHSmKQnJB3PLsCzq+o5jZqkn5D0HUnPZfv8x1XPqSySlkr6V0n/WPVcyiDpZUn/JulZSVOFPnebTulkF1v5d3outgLc1PaLrUi6HngTeDAiPlj1fEZN0mpgdUQ8I+m9wBFge5v/ntW96MQlEfGmpOXAt4BdEfF0xVMbOUm/D3SAn46Ij1Y9n1GT9DLQiYjC33vQtiP8dy+2EhFvAecvttJqEfEk8HrV8yhLRJyKiGey+z+i+xlNa6qd1WhF15vZ4vLs1p6jtQEkrQV+Fbi/6rm0QduC3+9iK60OQeokjQO/AHy72pmMXnZq41ngDPB4RLR+n4HPAXcC71Q9kRIF8A1JR7LrgxSmbcFf0MVWrNkk/RTwFeCOiPhh1fMZtYh4OyKuontNiY2SWn36TtJHgTMRcaTquZTs2oi4GtgC3Jadsi1E24K/oIutWHNl57G/AvxtRPxD1fMpU0T8D/BNYHPFUxm1a4Ffy85p7wN+SdLfVDul0cs+bZiIOEP304g3FvXcbQu+L7aSgOwXmH8JHI+IP616PmWQtErSpdn9nwQ2Ad+tdlajFRF/GBFrI2Kc7r/lf46I36p4WiMl6ZLshQhIugT4CFDYq+9aFfxUL7Yi6UvAU8AHJM1IuqXqOY3YtcBv0z3ieza7bZ3vmxpuNfCEpOfpHtg8HhFJvEwxMZcB35L0HPAd4OsRcbCoJ2/VyzLNzGywVh3hm5nZYA6+mVkiHHwzs0Q4+GZmiXDwzcwS4eCbmSXCwTczS4SDb2aWiP8HoHrd9T7YVz8AAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "source": [
    "ax = av.plot()\n",
    "orig.plot(ax=ax, color='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## clean version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 1), (2, 1), (3, 1), (4, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "orig = gpd.GeoSeries([top, bottom])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "centroids = gpd.GeoSeries(triangulate(top.union(bottom), edges=True)).centroid\n",
    "axis = centroids.drop(orig.sindex.query_bulk(centroids, predicate='intersects')[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "average = LineString([start] + [p for p in axis] + [end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "'LINESTRING (0 0, 4 0, 3.5 0, 3 0, 2.5 0, 2 0, 1.5 0, 1 0, 5 0)'"
     },
     "metadata": {},
     "execution_count": 111
    }
   ],
   "source": [
    "average.wkt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "1.89 ms ± 43.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 1), (2, 1), (3, 1), (4, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "orig = gpd.GeoSeries([top, bottom])\n",
    "\n",
    "centroids = gpd.GeoSeries(triangulate(top.union(bottom), edges=True)).centroid\n",
    "axis = centroids.drop(orig.sindex.query_bulk(centroids, predicate='intersects')[0])\n",
    "\n",
    "average = LineString([start] + [p for p in axis] + [end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "316 µs ± 96.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 1),  (4, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "orig = gpd.GeoSeries([top, bottom])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "784 µs ± 71.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "centroids = gpd.GeoSeries(triangulate(top.union(bottom), edges=True)).centroid\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "646 µs ± 179 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "axis = centroids.drop(orig.sindex.query_bulk(centroids, predicate='intersects')[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "138 µs ± 12.7 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "average = LineString([start] + [p for p in axis] + [end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 1), (2, 1), (3, 1), (4, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "orig = gpd.GeoSeries([top, bottom])\n",
    "\n",
    "centroids = gpd.GeoSeries(triangulate(top.union(bottom), edges=True)).centroid\n",
    "axis = centroids.drop(orig.sindex.query_bulk(centroids, predicate='intersects')[0])\n",
    "\n",
    "average = LineString([start] + [p for p in axis] + [end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [],
   "source": [
    "union = top.union(bottom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "0      True\n1      True\n2      True\n3      True\n4      True\n5      True\n6      True\n7      True\n8      True\n9      True\n10    False\n11    False\n12    False\n13    False\n14    False\n15    False\n16    False\ndtype: bool"
     },
     "metadata": {},
     "execution_count": 123
    }
   ],
   "source": [
    "centroids.intersects(union)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 1), (2, 1), (3, 1), (4, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "orig = top.union(bottom)\n",
    "\n",
    "centroids = gpd.GeoSeries(triangulate(orig, edges=True)).centroid\n",
    "axis = centroids[~centroids.intersects(orig)]\n",
    "\n",
    "average = LineString([start] + [p for p in axis] + [end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "1.97 ms ± 24.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 1), (2, 1), (3, 1), (4, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "orig = top.union(bottom)\n",
    "\n",
    "centroids = gpd.GeoSeries(triangulate(orig, edges=True)).centroid\n",
    "axis = centroids[~centroids.intersects(orig)]\n",
    "\n",
    "average = LineString([start] + [p for p in axis] + [end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "996 µs ± 29.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 1), (2, 1), (3, 1), (4, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "orig = top.union(bottom)\n",
    "\n",
    "centroids = from_shapely(triangulate(orig, edges=True)).centroid\n",
    "axis = centroids[~centroids.intersects(orig)]\n",
    "\n",
    "average = LineString([start] + [p for p in axis] + [end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 1), (2, 1), (3, 1), (4, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "orig = top.union(bottom)\n",
    "\n",
    "centroids = from_shapely(triangulate(orig, edges=True)).centroid\n",
    "axis = centroids[~centroids.intersects(orig)]\n",
    "\n",
    "average = LineString([start] + [p for p in axis] + [end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [],
   "source": [
    "def consolidate(start, end, line1, line2):\n",
    "    orig = line1.union(line2)\n",
    "\n",
    "    centroids = from_shapely(triangulate(orig, edges=True)).centroid\n",
    "    axis = centroids[~centroids.intersects(orig)]\n",
    "\n",
    "    return LineString([start] + [p for p in axis] + [end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 263,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 3), (2, 2.5), (3, 1.2), (4.3, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<shapely.geometry.linestring.LineString at 0x12d371880>",
      "image/svg+xml": "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"-0.2 -0.7 5.4 2.75\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,1.3500000000000003)\"><polyline fill=\"none\" stroke=\"#66cc99\" stroke-width=\"0.10800000000000001\" points=\"0.0,0.0 4.65,0.5 3.15,1.7499999999999998 4.15,0.0 3.65,0.0 3.0,0.09999999999999998 3.6500000000000004,1.1 2.5,1.85 1.5000000000000002,0.6 1.0,1.25 2.5,0.09999999999999998 1.0,-0.5 5.0,0.0\" opacity=\"0.8\" /></g></svg>"
     },
     "metadata": {},
     "execution_count": 149
    }
   ],
   "source": [
    "consolidate(start, end, top, bottom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<shapely.geometry.linestring.LineString at 0x12d67f670>",
      "image/svg+xml": "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"-0.2 -0.2 5.4 1.4\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,1.0)\"><polyline fill=\"none\" stroke=\"#66cc99\" stroke-width=\"0.10800000000000001\" points=\"0.0,0.0 1.0,1.0 2.0,1.0 3.0,1.0 4.0,1.0 5.0,0.0\" opacity=\"0.8\" /></g></svg>"
     },
     "metadata": {},
     "execution_count": 262
    }
   ],
   "source": [
    "top"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<shapely.geometry.linestring.LineString at 0x129cf7a90>",
      "image/svg+xml": "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"-0.2 -1.2 5.4 1.4\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,-1.0)\"><polyline fill=\"none\" stroke=\"#66cc99\" stroke-width=\"0.10800000000000001\" points=\"0.0,0.0 1.0,-1.0 2.0,-1.0 3.0,-1.0 4.0,-1.0 5.0,0.0\" opacity=\"0.8\" /></g></svg>"
     },
     "metadata": {},
     "execution_count": 151
    }
   ],
   "source": [
    "bottom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [],
   "source": [
    "orig = top.union(bottom)\n",
    "\n",
    "centroids = from_shapely(triangulate(orig, edges=True)).centroid\n",
    "axis = centroids[~centroids.intersects(orig) & centroids.intersects(list(polygonize(orig))[0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<matplotlib.axes._subplots.AxesSubplot at 0x12d462580>"
     },
     "metadata": {},
     "execution_count": 172
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT0AAAD4CAYAAACeyTEuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO2deXiU5dX/Pyck7CA7soUQQBZFUSJCFRUFCrhg1VYp1bdWpVZpXd5WRfvWt7VUfesPl2q1iGuLWouAG+4bIKCyIxIREnaEhD2ENdy/P2YGwzDJPJl59jmf68qVmWeezH2PJl/Oec5zzleMMSiKomQKWV5vQFEUxU1U9BRFyShU9BRFyShU9BRFyShU9BRFySiyvd5AdbRo0cLk5eV5vQ1FUQLG/PnzS40xLRO95mvRy8vLY968eV5vQ1GUgCEia6p6TdNbRVEyChU9RVEyChU9RVEyChU9RVEyChU9RVEyirRFT0TqisgXIrJYRJaJyB8TnCMi8qiIrBSRJSJyWrrrKoqipIIdt6zsB84zxpSJSA4wS0TeNsbMrXTOMKBr9OsM4Inod0VRFFdJO9IzEcqiT3OiX/HzqkYAL0TPnQs0EZE26a4dRlZuKWP7ngNeb0NRQostNyeLSC1gPtAFeNwY83ncKe2AdZWer48e25TgvUYDowFyc3Pt2F5gGP/+Cp6fXczBQ4e55qxOXHdWPk0b1PZ6W4oSKmwpZBhjKowxvYH2QF8ROSnuFEn0Y1W81wRjTIExpqBly4RdJKHEGMO0hRvo0qoRA3u05u+frOKsBz7ir+8WauSnKDZia/XWGLMD+AQYGvfSeqBDpeftgY12rh10lm3cxdpt5fykoD2P//Q03rn5bM7t3krFT1Fsxo7qbUsRaRJ9XA8YBBTGnfY6cHW0itsP2GmMOSa1zWTeXLKJ7CxhSM/jAeh2fCMVP0VxADuu6bUBno9e18sCXjHGvCkiNwAYY54EpgPDgZVAOXCNDeuGBmMM05du4gddWhxzDS8mft98t5tHP/qWv3+yiuc+W83Pz8zTa36KkgLiZ2OggoICkwlTVr7asJML/zaLBy7rxRWnV1+8iYnf9KWbqJ9TS8VPURIgIvONMQWJXtOODB8Qn9pWh6a9ipIeKnoeU11qWx0qfoqSGip6HhOr2l7QK3mUlwgVP0WpGSp6HlOT1LY6VPwUxRoqeh6SampbHSp+ilI9Knoekm5qWx0qfoqSGBU9D7Erta0OFT9FORoVPY9wIrWtDhU/RYmgoucRTqa21aHip2Q6Knoe4UZqWx0qfkqmoqLnAW6nttWh4qdkGip6HuBValsdKn5KpqCi5wFep7bVoeKnhB0VPZfxU2pbHSp+SlhR0XMZP6a21aHip4QNFT2X8XNqWx0qfkpYsGNcfAcR+VhElkfNvm9OcM65IrJTRBZFv/6Q7rpBJCipbXWo+ClBx45x8YeA/zbGLBCRRsB8EXnfGPN13HkzjTEX2rBeYImltjcN7Oz1VtJGx9grQcUOs+9NxpgF0ce7geVEPG2VOIKa2lZHfOT3+MerOO/BT7hn2lccqjjs9fYU5RhsvaYnInnAqUC82TdAfxFZLCJvi8iJ1bzHaBGZJyLzSkpK7Nyep4Qhta2Obsc34rGRp3J+91Zs33uQ5+euYegjM/h64y6vt6YoR2Gb6IlIQ+BV4BZjTPxv+gKgozHmFOBvwLSq3iesZt9Bq9qmwjOfrebDwi388ux8jquXQ1HJHi7620z+33vfsP9QhdfbUxTAJtETkRwigjfJGDMl/nVjzC5jTFn08XQgR0Ra2LF2UAhjaluZ2atK+cv05Qzp2Zo7hnbnqasLEITWx9Xlbx+t5MJHZ7Fw7Xavt6kotlRvBXgaWG6MGV/FOcdHz0NE+kbX3Zru2kEh7Knt+u3ljHlxIZ1aNGD8Fb3JyhL6dmrGHcO6sXHHPkae3oGy/Ye49InZ3Pvm1+w9oFGf4h12VG/PBK4ClorIouixu4BcOGL2fTnwKxE5BOwFrjR+Nty1mTBVbePZd7CCG/41n4OHDvOPq/rQsM73v1LXD8hn3urt/Gf+ep75eQHvLtvM07OKef/rzdx/WS9+0Dmjgn3FJ6jZtwvc/3YhE2cW8eXdg0IV6RljuO2VxUxbtIGJVxdwfo/Wx5yzc+9BLn5sFvsOVvDWbwbw7eYy7pyyhDVbyxnZN5exw7vTuG6OB7tXwoyafXtImFPbZz5bzdSFG7h10AkJBQ/guHo5PDGqDzvKD3Lzywvp26kZ79x8NtcP6MS/v1zLkPEz+Khws8s7VzIZFT2HCWvVtnLhYszALtWe27NtY+695CQ+W7mVh95fQb3atbj7gp5MufFMGtfL5hfPzeOWlxeyTTs6FBdQ0XOYMFZtExUukvGTgg5cUdCBxz5eeSSy692hCW/8+ix+c35X3lyyicHjP+XNJRvx8yUXJfio6DlIGFPb6goXyfjjiBPp2aYxt/57Meu2lQNQJ7sWtw0+gTd+fRZtm9RjzIsL+eU/57Nl1z6nPoKS4ajoOUjYUltjDGOnLGXZxl08fGVvOrdsWKOfr5tTiyd/1ofDxnDjpAXsO/j9rSs92jRm6o0/YOyw7ny6ooRB4z/llXnrNOpTbEdFz0HCltpaKVwkI7d5fcb/pDdLN+zkT28ePZMiu1YWvzynM2/fPIDuxzfm9slLuPqZL45EhYpiByp6DhG21LYmhYtkDO7ZmhvO6cyLn69lyoL1x7ye37IhL4/ux70jTmTBmu388OEZPD97NYcPa9SnpI+KnkOEKbVNpXCRjN8OOYF++c24a+pSCr87dihBVpZwVf883r31bArymnHP68u4YsIcVpWUpb22ktmo6DlEWFLbdAoX1ZFdK4tHR55K47o5/OpfC9i972DC89o3rc/z15zOgz8+hW++282wR2byxCerdGyVkjIqeg4QltQ23cJFMlo1qstjPz2NtdvKuX3ykiqLFiLC5X3a88F/n8N53VrxwDuFXPL3z3RslZISKnoOEJbU1o7CRTL6dmrGHUO78fZX3/HMZ6urPbdVo7o8eVUf/j7qNL7buY+LH5ulY6uUGqOi5wBhSG3tLFwk4/oB+Qzp2Zr7pi9n3uptSc8f3qsN7996Dhef0lbHVik1RkXPZsKQ2jpRuKgOEeGvPz6Fdk3rcdOLCygt25/0Z5o2qM34K3rz7M9P17FVSo1Q0bOZoKe2ThUukhE/mKDC4u0pA7u34r1bz+anfXN5elYxP3x4BrNXlTq8WyXIqOjZTJBTW6cLF8mIH0xglUZ1cxj3o168dH0/ROCnT33O2ClL2VVFRVjJbFT0bCToqa0bhYtkJBpMYJX+nZvr2ColKW6ZfYuIPCoiK0VkiYiclu66fiTIqa2bhYtkJBpMYBUdW6Ukw45IL2b23QPoB9wkIj3jzhkGdI1+jQaesGFd3xHU1NbtwkUyqhtMYBUdW6VUhVtm3yOAF0yEuUATEWmT7tp+IqiprVeFi2RUN5jAKjq2SkmEW2bf7YB1lZ6v51hhjL1HIM2+g5jael24SEaywQRW0bFVSmXcMvtOlC8l/I0Lqtl3EFNbPxQukpFsMIFVdGyVEsMVs28ikV2HSs/bAxvtWNsPBDG19VPhojqsDiawio6tUlwx+wZeB66OVnH7ATuNMZvSXdsvBC219VvhIhlWBxNYRcdWZTZ2RHoxs+/zRGRR9Gu4iNwgIjdEz5kOFAErgaeAG21Y1zcEKbX1a+EiGZUHEzw9q9iW96w8tmrF5jIdW5UhpP0bb4yZReJrdpXPMcBN6a7lR4KU2lYuXEy8usB3hYtkXD8gn3mrt3P/24X07tCEgrxmab9nbGzV2Se04A/TlvHAO4VMX7qJBy47mZ5tG9uwa8VvaEdGmgQptQ1C4aI6UhlMYJXKY6s27dzLxY/NYryOrQolKnppEpTUNiiFi2SkOpjAKpXHVj2qY6tCiYpeGgQltQ1a4SIZqQ4msEr82KrLnpjNn3VsVWhQ0UuDIKS2QS1cJCOdwQRWiY2tGtk3l4mzihn6yAzmrNrqyFqKe6jopYHfU1u/d1ykSzqDCaxSeWwVwMin5nLX1KVp3y+oeIeKXooEIbUNeuEiGXYMJrBK5bFVL3+xliEPzeDjwi2Orac4h4peivg9tQ1L4SIZdgwmsErlsVWN6mZzzXNfcsvLC9muY6sChYpeivg5tQ1b4SIZdg0msErvDk1489cDuDk6tmqQjq0KFCp6KeDn1DashYtk2DWYwCq1s7O4dfAJvPmbs2jXVMdWBQkVvRTwa2ob9sJFddg9mMAq3Y9vzJRf6diqIKGilwJ+TW3DXrhIht2DCayiY6uChYpeDfFrapsphYtkODGYwCo6tioYqOjVED+mtplWuEjG9QPy+eGJrbn/7ULmrd7m6to6tsr/qOjVEL+ltplauKgOJwcTWEXHVvkXFb0a4LfUNpMLF8loXNfZwQRWiI2tev+2szmvWyseeKeQH/19Nl9vdL66rFSNXePinxGRLSLyVRWvnysiOysNGf2DHeu6jd9S20wvXCTD6cEEVtGxVf7CrkjvOWBoknNmGmN6R7/+ZNO6ruKn1FYLF9ZwYzCBVY6MreqtY6u8xBbRM8bMANy9YuwyfkpttXBRM9wYTGCVpg1qM/4nvXn2Gh1b5RVuXtPrLyKLReRtETmxqpP86nvrl9RWCxc1x83BBFYZ2C0ytuqnZ+jYKrdxS/QWAB2NMacAfwOmVXWiX31v/ZDaauEiddwcTGCVRnVz+PMlvXh5dD8EHVvlFq6InjFmlzGmLPp4OpAjIi3cWNsO/JLaauEiPdweTGCVfvnNefvmsxl9dr6OrXIBV0RPRI6P+uMiIn2j6wYmlvdDaquFC3twezCBVerVrsVdw3vo2CoXsOuWlZeAOUA3EVkvItfG+d5eDnwlIouBR4ErTYA6sr1ObbVwYR9eDSawio6tch7x83/MgoICM2/ePE/3YIzhnL9+Ql6LBrzwi76ur7/vYAWXPzmbNaXlTBtzpl7Hs4kvircx8qm5DOnZmr+POo1oIuIrCr/bxe2Tl7Bk/U6G9GzNny85iVaN63q9rUAgIvONMQWJXtOOjCR4mdpq4cI5vBxMYBUdW+UMKnpJ8DK11cKFs3g5mMAqVY2tWr9dx1aliopeNXhZtdXChfP4YTCBVeLHVg15SMdWpYqKXjV4ldpq4cI9/DCYwCpVja0q0rFVNUJFrxpcTW0nTYK8PPbVrssNdzzPwb37tePCJfwymMAq8WOrhurYqhqholcFrqa2kybB6NGYNWsYO+QmljVtz8OvPUDn915zdl3lCH4aTGCFqsZWLd/kn3sP/YqKXhW4mtrefTeUl/Nsn4uZetJ53DrrRc7/elbkuOIafhpMYJX4sVUX/U3HViVDRa8KXE1t165ldm4vxp13LUNWzGHM7H8fOa64hx8HE1hFx1ZZR0UvAW5Xbdd3782YEXfSadsGxr81niyiF9Nzcx1fWzma3Ob1eeCyk1m6YSc3/Gs+UxeuZ/G6Hb7r3EiEjq2yhl4lT0Astb1pYGfH19p3sIIbfnIPB3ceYMKUP9PwwN7IC/Xrw7hxjq+vfE9x6R5e+mIt/5m3DoBPvynhk2++H2/WqlEd8ls2oHPLhuS3bEh+ywZ0admQtk3qUctHFfbY2KoH3ilk4qxi3l++mfsvPZn+nZt7vTVfoKKXALdS2yMdF/uzmXjCPvIb58AOiUR448bBqFGOrq/A/kMVvLdsMy9+vpY5RVuplSUU5Dbl89Xb6NCsHs/8vC+rSsooKtlDUUkZq0rKeGvpJnaUfx/51c7OolPzBnRu1YD8Fg2PfM9v2YBGdXM8+VyxsVUXntyWO19dwsin5vLTM3IZO6y7Z3vyCyp6cbiZ2j4b7bi4bfAJnH9+V/jVFY6up3xPLKqbPH892/YcoH3Tevzuh934cZ/23DX1K7KzhFaN6tKlVUO6tDq2/W/bngNRMSxjVVQQCzft5t1lm4+61y8WHea3bBiNEN2NDmNjqx76YAUTZxbxceEW/vKjXgzs3srxtf2Kil4cbqW2s1eVMk47LlwlUVQ3uEdrRp6Ry4AuLcjKElaVlPHB8s20a1K3WlFq1qA2zRo04/S8ZkcdP3DoMGu3lR8THU73MDqMja0a3qsNd0xewjXPfcklvdtyz0Unem594AUqenG4kdpqx4W7VBfVxU8tmTiziNrZWbROcZpJ7ews30aHvTs04Y1fn8XjH6/k8Y9XMvPbUv444kQu6NXGl1NmnEJFrxJupLaVPS4maMeFY1iJ6uIp2b2fVxds4LLT2jvS2uWH6LB2dha3Dj6BYb2O5/bJSxjz4kJe77kxo8ZW6V9cJZxObSuPipp4dQH5OirKdmoS1cXzzzmrOVhxmOsGdOKuKUvd2TDORYedWzSkXdPE0WFsbNXTs4oZ//4KBo3/lN9f2JMf92kf+qjPFtETkWeAC4EtxpiTErwuwCPAcKAc+LkxZoEda9uJ06ntUYULHRVlG6lEdfGUHzjEC3PXMKhHa1/NLbQ7Ovz+lptIdPjLczozuGdr7nx1KbdPXsIbizdy36W9aN+0vtsf1TXsivSeAx4DXqji9WFA1+jXGcAT0e++wenUVgsX9pNOVBfP5Pnr2VF+kNFn5zu0W3uxOzoccmJruh3fkMnz1zN4/KfcOawHV/XrGMrrzbaInjFmhojkVXPKCOCFqC/GXBFpIiJtjDGb7FjfDpxMbbVwYR92RHXxVBw2TJxZzKm5TSjo2NSBXbtLutGhCNzz+jLeXLKRp39+Oo1Ddl+fW9f02gHrKj1fHz12jOiJyGhgNECui21YTqW2WriwBzujunjeXfYda7eVM3ZY91Bfz7IcHW4p47nZa/hy9XZ+98oi/nH16R7s1jnc+gtM9JuUcFqjMWYCMAEixkBObqrSmo6ktlq4SA8norp4jDH8Y0YRHZvXZ8iJ3ll8ek3l6PDlL9ZyoOIwIvDu11t4Y/FGLjqlrddbtA23RG890KHS8/bARpfWTopTqa0WLlLDyagunnlrtrN43Q7uHXGir/pnveK7nfsY99ZymtbPoWXDOtSrk83dU5dSkNeUNsfV83p7tuCW6L0OjBGRl4kUMHb66XqeE6mtFi5qhhtRXSL+8WkRTevncHmfDslPDjnGGH4/7SsOHj5MkwZ16NyqIbcP7c7wR2by2/8s5p+/OCMU16PtumXlJeBcoIWIrAfuAXIAjDFPAtOJ3K6yksgtK9fYsa4dOJHaauHCOm5GdfHEWs5+c35X6tWu5ehaQeDNJZv4YPlmxg7rzl/f/YaLT2lLpxYN+J8Le3LX1KU8O3s1157Vyettpo1d1duRSV43wE12rGU3dqe2WrhIjldRXTyxlrOr+3d0ZT0/s23PAf739WWc0qEJ5/dozX1vF9KpReQa9Mi+Hfhw+WYeeKeQs7q0oNvxjTzebXpk/F+knamtFi6qx8uoLp7KLWctGtZxdW0/8qc3lrFr30H+77KTWbN1DwD5LRsAET+O+y87maEPz+Dmlxfy2pgzqZMd3Mg4o0XP7tRWCxfH4peoLp7KLWeZzkeFm5m2aCO3DOpKt+MbMfPbyODU/BYNjpzTslEdHrjsZK57YR7j31vB2OE9vNpu2mS06NmZ2mrh4mj8FNXF49eWMy/Yve8gd0/9im6tG3HjuZHf26LSPTStn0OT+kcHAoN6tmZk31wmzCxiYPdW9MsP5iTmjBY9u1JbLVxE8GtUF0/QWs6c5P63C9m8ax9P/KwPtbMjljlFJWVVXpr5/QU9mLOqlP9+ZTHTbx7AcfWC162RsaJnV2qrhQt/R3XxhK3lLB3mrNrKpM/Xcv2ATvTu0OTI8eLSPQzo2jLhzzSok81DV/Tm8ifncM9rX/Hwlae6tV3byLy/0Ch2pLaZXLgISlQXT6a0nCVj74EKxk5ZQsfm9bltcLcjx8v2H2Lzrv1HihiJODW3KWMGduGRD7/l/B6tA9etkbGiZ0dqm4mFiyBFdfFoy9n3PPzBClZvLefF68846h7F1aXRym2LqkUPYMx5XfhkRUkguzUyUvTsSG0zqXAR1KguHm05i7B43Q6emlnEyL65/KBzi6NeK4qKXuwevarIqZXFw1f0DmS3RkaKXrqpbaYULoIc1SVCW84i46XueHUJrRrVZezw7se8XlRShgh0bJ58iGhQuzUyUvTSSW3DXrgIS1QXj7acRXjik1UUfrebp/+rIOGcvOLSPbRrUo+6Odb+GwWxWyNcf7EWSCe1DXPhImxRXTzacgYrNu/msY+/ZUTvtlVegy4q2VOj3+sgdmtkeb0Bt4mlthf0qnmUFytc3DoooIWLSZMgLw+ysiAvj/3/nMQbizcycsJcBj74CU/PKqZvXjOe/0VfZvxuIDcN7BIKwdOWs8itOrdPXkKjujn84cKeCc8xxlBcuidpESOeWLdG4Xe7Gf/eCju26ygZF+mlmtoGvnAxaRKMHg3l5RQ3bctLnc5j8rwsti1bGLqoLh5tOYNnPytm0bodPHJlb5pXIfwlZfsp23+ITjUUPQhWt0ZGiV6qqW0YChf7/+ce3svtw4u9hzKn4ynUOlzB4G/nMnLTQgbMnh7Iz2QFbTmDtVvLefC9bxjUoxUXV3NPXVHJ0YMGakpQujUySvRSqdoGvXBx5FrdhX9iW/3jaL/jO3736fP8eOkHtNqzPeICE1LBA205M8Zw55Ql5GRl8edLelV7Q3bxkdtVUhO9oHRrBOsvOE1qmtoGtXCRsAK7tYiRb05lQPFCsirbk7hovuQ22nIG//5yHbNXbeW+S3tx/HHVX7ooKimjTnYWbdO40TgI3Rp2TU4eSsTMuxYw0Rhzf9zr5wKvAcXRQ1OMMX+yY22rpJLaBq3jotoK7BvbYVohR/kx1a8P48Z5tl+nyfSWs5jfRf/85lx5evJ7E4tL99CpRYO0L3WMOa8Ln/q4WyNt0RORWsDjwGAiBkBfisjrxpiv406daYy5MN31UqWmqW1QCheW76sbNSry/e67Ye3aSIQ3btz3x0NGprecVfa7uP+y6tPaGEWle+jWOv377HJqZfGQj7s17Ij0+gIrjTFFAFHznxFAvOh5Sk1S2yAULlK6r27UqNCKXDyZ3nIW87v4/QU96Ng8+TW6gxWHWbu1nGEn2fMPhJ+7NewQvURG3mckOK+/iCwmYv34W2PMskRv5oTZd01SWz8XLsLaLeEEmdxyVtnv4pozrYnN+u17OXTYJO25rQl+7daw4y/aipH3AqCjMaZMRIYD04Cuid7MCbNvq6mtXwsXYe+WsJtMbzmr7HdhNcotKikDUq/cJsKv3Rp2iF5SI29jzK5Kj6eLyN9FpIUxptSG9ZNiNbX1U+FCo7rUyeSWs3i/C6vEblfpnOI9elXhR28NO0TvS6CriHQCNgBXAj+tfIKIHA9sNsYYEelLpP1tqw1rJ8VqauuXwoVGdemRyS1nu/cd5K4pR/tdWGVVSWJfDDvwW7dG2qJnjDkkImOAd4ncsvKMMWaZiNwQff1J4HLgVyJyCNgLXBn1wnUcK6mt14ULjersI5Nbzu5/u5Atu/fx5FXf+11Ypbi0al8MO/BTt4ZdZt/Tgelxx56s9Pgx4DE71qopyVJbLwsXGtXZSya3nFXld2GV6nwx7MBP3Rr+KU06QLLU1ovChUZ1zpGpLWdV+V1YJeaLYWcRIxF+6dYIteglS23dLFxoVOcsmdxy9lAVfhdWWe1QESMRfujWCLXoVZfaulG40KjOPTK15Wzxuh1MrMLvwiqrjtyu4nym44dujdCKXnWprdOFC43q3CVTW86S+V1Ypbh0j2VfDDuo3K3xzGfFXDfA3csRoRW9qlJbpwoXGtV5R6a2nCXzu7BKTX0x7CDWrfF/737DgK4tXe3WCK3oJUptnShcaFTnPZnYcmbF78IqRSV7HC9ixONlt0YoRa+q1NauwoVGdf4hE1vOrPhdWCXmi3F5n/Y27c46XnVrhFL0EqW2dhQuNKrzH5nYcmbF78IqJbtT98Wwg8rdGud2a0X/zs53a4RS9OJT23QKFxrV+ZdMbDmz6ndhlaLS9Hwx7OD7bo1FvH3L2Y53a4RO9OJT21QLFxrV+Z9Mazmrid+FVdL1xbADt7s1Qid6lVPbmhYuNKoLDpnYclYTvwur2OGLYQdudmuETvQqp7ZWCxca1QWISZPg7ruZ3Pwkdgz+FaP3rQQK7F/j8y2wbz/cN9IXY/Vr6ndhFbt8MezArW6NUIle5dR2+Xe7qi1caFQXQKKG5RV79zHxh7/n1A2FFPztD9Cgwj5Ripmij/ifyPM1ayLPwTPhi/hdLK2R34VVikr2+GaisVvdGjWbP+NzYqntD/KbVVm4KC7dw1+mL6f/fR/x65cWsm57Ob/7YTfm3HkeT17Vh3NOaKmC51fuvhvKy3n3hP6sbdqG0V9MQcrLI8dtXuMo7F6jhkT8Lrbw2yHdLPldWOVgxWHWbiv3tIgRT6xb47OVW3nms+LkP5ACoYr03lyyiVoCry3aeFThQqO6kLB2LQb4R99L6bh9I0O+nXvkuJ1r1Oi4w6Tid2EVJ3wx7MDpbg1bIj0RGSoi34jIShG5M8HrIiKPRl9fIiKn2bFuZYwxvLVkI80b1qFw824evrI3IqJRXZjIzWVeu54sbtuN676cRi1z+MhxO9eo0XEnmDQJ8vIgK4s/XTuOXeX7a+R3YXWNoouvAKDT6Ksia/qEWLdGozrZ3PzyQvYfqrD1/d3yvR1GxAioKxGntCdI7JiWMss27mLd9r0AXNCrDRNnFmtUFzbGjeMfr31L0/KdXL70w8gxuw3Lx437/hpeDDdN0WPXFMvL+Si/gGl5fbll7it0O2mn7dcti08cAkDnr+fB6C8ir/nEItTJbg07Ir0jvrfGmANAzPe2MiOAF0yEuUATEWljw9pHePTDbwHIqSW8tXSTRnUhZNWQEXyQfzpXFX1GvYoD0LEjTJhg7x/qqFGR96wbvdnZiTWqI3pNcXftetw1dAzdSlZz44wXHbluuapZO5qW76TJvjLPr1smItat8Y8ZRcxZZZ+ljlu+t4nOaQdsin+zVH1vG9eN9F22bFiHO4Z156KT26rIhYwjLWcvj4eGjzu30KhRUDYn8vjFY67WOEv02ufNLdgAAA0/SURBVOH9517DlgZNeXLqX6h9+JAj1y2Lm7Yjf9uGY477iav65TJ96SbK9x+y7T3tiPSs+N5aOSdy0JgJxpgCY0xBy5bWZ/ZfeHI7srOEjTv38bvJSxg3fTlbdu+z/POKv8mYlrPcXOZ06MWkU4dz7bzX6L1pxZHjdq4BUNysHZ22bzjmuF8oKinj589+SXaW0MHGWX92iF5S31uL56TFud1bMfXGM2lYuxbZWcJzn61mwAMfc++bX6v4hYBMaTnbe+84xg7/DR23b+S2mdHiggPXLcuOa8bmRs3ptG2jM2ukSVFJGVdOmEvFYcNLo/txQmv7Krh2iN4R31sRqU3E9/b1uHNeB66OVnH7ATuNMcektunSq/1xvDS6P9lZQstGtTmveyuem63iF3QyqeXsoeansbpJG+5b+Iqj1y1XPxi5PNB52wb3r1smwUnBA/d8b6cDw4GVQDlwTbrrVkWv9scx6bp+jJo4l6UbdvKva/vy6oINPDd7Nf+au4af9evIL8/Jp1UjbS8LCpnicnaU38X97zu61qqCAbByEZ0+fQd80pEBzgseuOd7a4Cb7FjLCpWF73eTl/DS9f0YM7ALj328UsUvYGSKy5ldfhdWcdsXwwpuCB6ErA2tMjHh27X3ICOfmkutLOHBH5/Ch7edw0WntNW0NyDEXM5GD8gPtctZzO/iz5eclJbfhVW88MWoDrcED0IsenCs8K3bVk5eiwYqfgEhU1zOYn4XF5/SlkE9nfVfjuGFL0ZVuCl4EHLRg8TCB6j4BYCYy9l1Z3UKrctZZb+Ley5Kz+/CKjFfDD8UhdwWPMgA0YOqhQ9U/PxMJricxfwu7rmoZ9p+F1bx2hcjhheCBxkielC98IGKn9+IuZxd1T8vtC5na7busdXvwip+8MXwSvAgg0QPkgsfqPj5hbC7nBljuPPVpeRkZXHvJSe5WqQpKvHWF8NLwYMMEz2wJnyg4uclmdBy9u8v1zGnaCtjh/dwbCx6VRSXeueL4bXgQQaKHlgXPlDx84Kwt5xV9rsY2df965Ve+WL4QfAgQ0UPaiZ8oOLnFmFvOXPS78IqXtyu4hfBgwwWPai58IGKn9OEveXsDYf8LqzihS+GnwQPMlz0IDXhAxU/Jwh7y5mTfhdWWbet3FVfDL8JHqjoAakLH6j42UnYW87+9MYydu87aL/fRQ0oLnWvcutHwQMVvSOkI3yg4pcuYW85+6hwM9MWbeSmgV089ZmNiV5nh9NbvwoeqOgdRbrCByp+qRLmlrPd+w5y15Sv6Na6ETeee6zxvJusKtlD0/o5NKlf27E1/Cx4oKJ3DHYIH6j41ZQwt5zd93YhW3bv44HLT6Z2trd/csWlZeQ7WBX3u+CBil5C7BI+UPGzQphbzuas2sqLn6/l2rM60btDE6+34+jtKkEQPEhT9ESkmYi8LyLfRr8nLLmJyGoRWSoii0RkXjpruoWdwgcWxK+SwTN5eb4yX3aasLac7T1QwdgpS+jYvD63De7m9XYo23+ILbv3OyJ6QRE8SD/SuxP40BjTFfgw+rwqBhpjehtjCtJc0zXsFj6oQvz+8gH3vjCLLaU7wRhYsyZi+JwBwhfmlrOHPljB6q3l3HdpL19EsKsdKmIESfAgfdEbATwfffw8cEma7+c7nBA+iBO/VXN57pRhDPjl09x73nVsadDEl+bLThDWlrOj/C46t/B6O0DkMgJg6z16QRM8SF/0WsdczaLfW1VxngHeE5H5UTPvKhGR0SIyT0TmlZSUpLk9e3BK+CAqfpP/wodP3cBFhTN4rs9F34vf1l22reNHwtpy5rbfhVXs9sUIouCBBdETkQ9E5KsEXyNqsM6ZxpjTgGHATSJydlUnpmr27TROCh+5ueTt2MSD0x85VvxCXPAIa8uZ234XVikqsc8XI6iCBxZEzxgzyBhzUoKv14DNItIGIPp9SxXvsTH6fQswFehr30dwD8eEb9y4iNkyfC9+/7qFi1qa0FZ7w9py5oXfhVVi01XSJciCB+mnt68D/xV9/F/Aa/EniEgDEWkUewwMAb5Kc13PcET4Ro2KmC137Egk/+hI3l/v5cE7Lg3trS5hbDnzwu/CKnb5YgRd8CB90bsfGCwi3wKDo88RkbYiEvPBbQ3MEpHFwBfAW8aYd9Jc11McE77Vq+Hw4cj3qNt8GO/zC2vLmRd+F1axwxcjDIIHaYqeMWarMeZ8Y0zX6Pdt0eMbjTHDo4+LjDGnRL9ONMaMs2PjXuPoNb4EhEn8wthy5pXfhVXS9cUIi+CBdmSkhdvCB+EQv7C1nHnpd2GVdHwxwiR4oKKXNl4IHwRX/MLYcual34VVUvXFCJvggYqeLXglfBA88Qtby5nXfhdWScUXI4yCByp6tuGl8EEwxC9sLWd+8LuwSk0HDYRV8EBFz1a8Fj7wt/iFreXMa78Lq8R8MayKXpgFD1T0bMcPwgf+E7+wtZz5we/CKjFfDCtz9MIueKCi5wh+ET7wj/iFreXsjz7wu7CKVV+MTBA8UNFzDD8JH3grfmFrOftw+WZe84HfhVVit6tUN1IqUwQPVPQcxW/CB96IX5haznbvO8jdU/3hd2GVotLqfTEySfBARc9x/Ch84J74ha3lzE9+F1YpLi2rMrXNNMEDFT1X8KvwgfPiF6aWM7/5XVilqGRPwiJGJgoeqOi5hp+FD5wTv7C0nO09UMGdPvK7sEpVvhiZKnigoucqfhc+sFf8wtRy9tAHK1jjI78LqxQnKGJksuCBip7rBEH4wB7xC0vLmR/9LqxSVHq0L0amCx6o6HlCUIQPUhe/sLSc+dXvwiqVfTFU8CKo6HlEkIQPai5+YWk586vfhVVivhgbd+xVwYuSrtn3j0VkmYgcFpEq/WxFZKiIfCMiK0WkOm/cjCJowgfWxC8sLWd7D1T41u/CKsWlezi+cV0VvEqkG+l9BVwKzKjqBBGpBTxOxAmtJzBSRPxlIOAhQRQ+qF78nv2sOBQtZ6tKy3zpd2EVYwwrt5SxbOMuFbxKiDEm/TcR+QT4rTFmXoLX+gP/a4z5YfT5WABjzH3J3regoMDMm3fMW4aSpet3MmriXOpkZ1lqDPcb+w5WsGHHXkrLDpAlUC8nixPbHhfYDoylG3ZSfqCCzi0bBPaa5KGKwyxYt4N6ObWYdtOZGSV4IjLfGJMw+8x2Yf12wLpKz9cDZ1R1ctQMfDRAbm6uszvzEb3aH8erv/oBd01d6vVWUqJuTi06t2xIuyYV7D1QQaO62YEVPIBe7Y5j/6EKatcK7mXv7FpZ9Mltyv2X9aJLq8wRvGQkFT0R+QBI1D90d9T7NulbJDhWZXhpjJkATIBIpGfh/UND19aN+M8NP/B6G4oSapKKnjFmUJprrAcq347fHtiY5nsqiqKkhBux+5dAVxHpJCK1gSuJmIQriqK4Trq3rPxIRNYD/YG3ROTd6PEjZt/GmEPAGOBdYDnwijFmWXrbVhRFSY20ChnGmKnA1ATHNwLDKz2fDkxPZy1FURQ7CG5pSlEUJQVU9BRFyShU9BRFyShU9BRFyShsaUNzChEpAdbU4EdaAKUObcctwvAZIByfIwyfAcLxOWr6GToaY1omesHXoldTRGReVf12QSEMnwHC8TnC8BkgHJ/Dzs+g6a2iKBmFip6iKBlF2ERvgtcbsIEwfAYIx+cIw2eAcHwO2z5DqK7pKYqiJCNskZ6iKEq1qOgpipJRhEL0wmA8JCLPiMgWEfnK672kioh0EJGPRWR51DDqZq/3lAoiUldEvhCRxdHP8Uev95QqIlJLRBaKyJte7yVVRGS1iCwVkUUikrZ/ROCv6UWNh1YAg4kMLP0SGGmM+drTjdUQETkbKANeMMac5PV+UkFE2gBtjDELRKQRMB+4JID/LwRoYIwpE5EcYBZwszFmrsdbqzEichtQADQ2xlzo9X5SQURWAwXGGFtusA5DpNcXWGmMKTLGHABeBkZ4vKcaY4yZAWzzeh/pYIzZZIxZEH28m8j8xHbe7qrmmAhl0ac50a/ARQci0h64AJjo9V78RBhEL5HxUOD+0MKGiOQBpwKfe7uT1IimhYuALcD7xpggfo6HgduBw15vJE0M8J6IzI8ah6VFGESvRsZDivOISEPgVeAWY8wur/eTCsaYCmNMbyKeLn1FJFCXHETkQmCLMWa+13uxgTONMacR8c6+KXopKGXCIHpqPOQjotfAXgUmGWOmeL2fdDHG7AA+AYZ6vJWaciZwcfR62MvAeSLyL2+3lBrRSewYY7YQmdTeN533C4PoqfGQT4gWAJ4Glhtjxnu9n1QRkZYi0iT6uB4wCCj0dlc1wxgz1hjT3hiTR+Rv4iNjzM883laNEZEG0aIYItIAGAKkdYdD4EUvLMZDIvISMAfoJiLrReRar/eUAmcCVxGJKhZFv4Yn+yEf0gb4WESWEPlH9X1jTGBv+Qg4rYFZIrIY+AJ4yxjzTjpvGPhbVhRFUWpC4CM9RVGUmqCipyhKRqGipyhKRqGipyhKRqGipyhKRqGipyhKRqGipyhKRvH/AWFCUj6F9mk1AAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "source": [
    "ax = gpd.GeoSeries(triangulate(orig, edges=True)).plot()\n",
    "gpd.GeoSeries(axis).plot(ax=ax, color='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<matplotlib.axes._subplots.AxesSubplot at 0x12bad0520>"
     },
     "metadata": {},
     "execution_count": 171
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 0 Axes>"
     },
     "metadata": {}
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<shapely.geometry.polygon.Polygon at 0x12a73f850>",
      "image/svg+xml": "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"-0.2 -1.2 5.4 4.4\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,2.0)\"><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"0.10800000000000001\" opacity=\"0.6\" d=\"M 0.0,0.0 L 1.0,3.0 L 2.0,2.5 L 3.0,1.2 L 4.3,1.0 L 5.0,0.0 L 4.0,-1.0 L 3.0,-1.0 L 2.0,-1.0 L 1.0,-1.0 L 0.0,0.0 z\" /></g></svg>"
     },
     "metadata": {},
     "execution_count": 167
    }
   ],
   "source": [
    "list(polygonize(orig))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## another approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 471,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance = 1\n",
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 3), (2, 2.5), (3, 1.2), (4.3, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "lines = [top, bottom]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 452,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = []\n",
    "ids = []\n",
    "\n",
    "extended_lines = [list(polygonize(lines))[0].buffer(distance).exterior] + lines\n",
    "\n",
    "for ix, line in enumerate(extended_lines):\n",
    "    current_dist = distance\n",
    "    line_length = line.length\n",
    "\n",
    "    while current_dist < line_length:\n",
    "        points.append(list(line.interpolate(current_dist).coords)[0])\n",
    "        ids.append(ix)\n",
    "        current_dist += distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 453,
   "metadata": {},
   "outputs": [],
   "source": [
    "voronoi_diagram = Voronoi(np.array(points))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate regions\n",
    "def _regions(voronoi_diagram, ids):\n",
    "    \"\"\"\n",
    "    Generate GeoDataFrame of Voronoi regions from scipy.spatial.Voronoi.\n",
    "    \"\"\"\n",
    "    import pandas as pd\n",
    "    from shapely.geometry import Polygon\n",
    "\n",
    "    # generate DataFrame of results\n",
    "    regions = pd.DataFrame({'ids': ids, 'region': voronoi_diagram.point_region})\n",
    "\n",
    "    # add vertices of each polygon\n",
    "    vertices = []\n",
    "    for region in regions.region:\n",
    "        vertices.append(voronoi_diagram.regions[region])\n",
    "    regions[\"vertices\"] = vertices\n",
    "\n",
    "    # convert vertices to Polygons\n",
    "    polygons = []\n",
    "    for region in regions.vertices:\n",
    "        if -1 not in region:\n",
    "            polygons.append(Polygon(voronoi_diagram.vertices[region]))\n",
    "        else:\n",
    "            polygons.append(None)\n",
    "    # save polygons as geometry column\n",
    "    regions[\"geometry\"] = polygons\n",
    "\n",
    "    # generate GeoDataFrame\n",
    "    regions_gdf = gpd.GeoDataFrame(regions.dropna(), geometry=\"geometry\")\n",
    "    regions_gdf = regions_gdf.loc[\n",
    "        regions_gdf['ids'] != 0\n",
    "    ]  # delete external cells\n",
    "    return regions_gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 455,
   "metadata": {},
   "outputs": [],
   "source": [
    "regions_gdf = _regions(voronoi_diagram, ids)\n",
    "tessellation = regions_gdf[[\"ids\", \"geometry\"]].dissolve(by=\"ids\", as_index=False)\n",
    "linestring = tessellation.geometry.exterior[0]\n",
    "clipped = linestring.intersection(list(polygonize(lines))[0])\n",
    "if hasattr(clipped, '__len__'):\n",
    "    lengths = []\n",
    "    for part in clipped:\n",
    "        lengths.append(part.length)\n",
    "    clipped = clipped[lengths.index(max(lengths))]\n",
    "coords = list(clipped.coords)\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 456,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<matplotlib.axes._subplots.AxesSubplot at 0x1303028e0>"
     },
     "metadata": {},
     "execution_count": 456
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT0AAAD4CAYAAACeyTEuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3de3RdZZ3/8fe3adqm16S00LRpkwqVQqu9hQJWXCBVEbmog4LOiCg/CiKjjvwQxJEBnVmyvKCAMzoFnZGRAWUYAbkJAnLxB5QktKWlLdSa9Eov5FJKL2mS7++PvROSNEmTnH3OPufsz2uts3JOzs7ez4H0k+e7936ex9wdEZGkGBJ3A0REMkmhJyKJotATkURR6IlIoij0RCRRhsbdgL5MmDDBKyoq4m6GiOSY6urqXe4+saf3sjr0KioqqKqqirsZIpJjzKyut/dU3opIoij0RCRRFHoikigKPRFJFIWeiCRKyqFnZiPMbJmZrTCz1WZ2Qw/bmJndYmbrzWylmc1P9bgiIoMRxS0rB4APuvseMysEnjOzR9z9hU7bfBSYET5OBH4WfhURyaiUe3oe2BO+LAwf3eerOhe4I9z2BaDYzEpTPXbeqa+Hn/0M3nwz7paI5K1IzumZWYGZLQd2AI+7+4vdNpkCbOr0enP4vZ72tcTMqsysaufOnVE0L3d8//tw+eVQUQHf+pbCTyQNIgk9d29197lAGbDQzGZ328R6+rFe9rXU3SvdvXLixB5HkeQnd7jnHjjxRPjYx+B731P4iaRBpFdv3b0R+BNwRre3NgNTO70uA7ZGeeyc9/LLsGEDXHIJ3H03vPKKwk8kDaK4ejvRzIrD50XAYmBtt80eAC4Mr+KeBDS5+7ZUj51XfvtbGDoUPv7x4PWsWQo/kTSIoqdXCjxlZiuBlwjO6T1oZpeZ2WXhNg8DG4D1wG3A5REcN3+0l7annw5HHNH1PYWfSKQsmxcGqqys9ETMslJTAwsWwO23w8UX973t6tXw3e8GPcNRo+ArX4Gvf/3QsBRJMDOrdvfKnt7TiIxs0L207Yt6fiIpUejFra/Sti8KP5FBUejFrf2q7ac+NbifV/iJDIhCL24DKW37ovAT6ReFXpwGW9r2ReEn0ieFXpxSLW37ovAT6ZFCL05RlbZ9UfiJdKHQi0s6Stu+KPxEAIVefNJZ2vZF4ScJp9CLSyZK274o/CShFHpxyHRp2xeFnySMQi8OcZW2fVH4SUIo9OIQd2nbF4Wf5DmFXqZlU2nbF4Wf5CmFXqZlY2nbF4Wf5BmFXqZlc2nbF4Wf5IkopoufamZPmdmacLHvr/awzalm1mRmy8PHdakeNyflSmnbF4Wf5LgoenotwJXufhxwEvBlMzu+h+2edfe54eM7ERw39+RaadsXhZ/kqCgW+97m7jXh87eANfSypm3i5Wpp25fewu+ee+JumUiPIj2nZ2YVwDyg+2LfACeb2Qoze8TMZvWxj/xc7DsfStu+dA6/2bPh858PnotkmchCz8xGA/cCX3P33d3ergHK3X0OcCtwX2/7ydvFvvOptO3LrFnwu99BcTH8zd/A7u6/CiLxiiT0zKyQIPDudPf/7f6+u+929z3h84eBQjObEMWxc0Y+lra9mTQJfvObIOS/+MWglyuSJaK4emvAL4A17n5TL9tMCrfDzBaGx03O2e58L217csopcOONcO+9cPPNcbdGpMPQCPaxCPgc8IqZLQ+/dy0wDcDdfw6cB3zJzFqAfcAFns0L7katvbS99tq4W5JZV14Jf/4zXHUVnHACLFoUd4tEUg89d38OsMNs81Pgp6keK2clqbTtzAz+4z+gshI+/ekg/I88Mu5WScJpREa6JbG07ay4OChx6+vhs5+F1ta4WyQJp9BLt6Rcte3LnDnwb/8GTzwB//RPcbdGEk6hl25JLW27+8IX4OKL4V/+BR56KO7WSIIp9NIp6aVtd7feCnPnwuc+B7W1cbdGEkqhl04qbbsqKgrO77W1wXnnwf79cbdIEkihl04qbQ/1rnfBHXdAdTV87Wtxt0YSSKGXLipte3fOOXD11fDv/w7/9V9xt0YSRqGXLipt+/bP/wynngqXXqqJCSSjFHrpotK2b0OHwl13aWICyTiFXjqotO0fTUwgMVDopYNK2/7TxASSYQq9dFBpOzBXXhn8t7rqqmCCApE0UuhFTaXtwLVPTFBeHkxMsGNH3C2SPKbQi5pK28HRxASSIQq9qKm0HTxNTCAZoNCLkkrb1GliAkmzTC32bWZ2i5mtN7OVZjY/1eNmJZW20dDEBJJGmVrs+6PAjPCxBPhZBMfNPipto6GJCSSNMrXY97nAHR54ASg2s9JUj51VVNpGSxMTSJpkarHvKcCmTq83c2gwtu8jNxf7VmkbPU1MIGmQqcW+e1o4qMcxRzm72LdK2/TQxAQSsYws9k3Qs5va6XUZsDWKY2cFlbbpo4kJJGIZWewbeAC4MLyKexLQ5O7bUj121lBpm16amEAiFEVPr32x7w+a2fLwcaaZXWZml4XbPAxsANYDtwGXR3Dc7KHSNv00MYFExDyL/2pWVlZ6VVVV3M3omzsccwzMmAGPPhp3a/KbO3zyk/Dgg/CnP8GiRXG3SLKUmVW7e2VP72lERqpU2maOJiaQCCj0UqXSNrM0MYGkSKGXCl21jYcmJpAUKPRSodI2PpqYQAZJoZcKlbbx0sQEMggKvcFSaRs/TUwgg6DQGyyVttlBExPIACn0BkulbfbQxAQyAAq9wVBpm300MYH0k0JvMFTaZh9NTCD9pNAbDJW22UkTE0g/KPQGSqVtdus8McFPfhJ3ayQLKfQGSqVt9rvySvjEJ+Ab34A//znu1kiWUegNlErb7KeJCaQPCr2BUGmbO8aN08QE0qOopov/pZntMLNVvbx/qpk1dZpk9LoojptxKm1ziyYmkB5E1dP7T+CMw2zzrLvPDR/fiei4maXSNvdoYgLpJpLQc/dngPoo9pW1VNrmLk1MIJ1k8pzeyWa2wsweMbNZvW2UteveqrTNXZqYQDrJVOjVAOXuPge4Fbivtw2zdt1blba5TRMTSCgjoefuu919T/j8YaDQzCZk4tiRUGmbHzQxgZCh0DOzSeH6uJjZwvC4b2bi2JFQaZs/NDFB4kV1y8pdwPPAsWa22cwu7rbu7XnAKjNbAdwCXODZvPZkdypt84cmJkg8rXt7OFrXNj89+yycdlrwh+yee4JRHJI3tO5tKlTa5idNTJBYCr3DUWmbvzQxQSIp9Pqiq7b5TRMTJJJCry8qbfOfJiZIHIVeX1TaJoMmJkgUhV5vVNomiyYmSAyFXm9U2iaPJiZIBIVeb1TaJo8mJkgEhV5PVNomV+eJCS64IBiju2wZNDXF3TKJyNC4G5CV2kvba6+NuyUSh3POge9+F66/Hu6//53vT5oEM2fCsccGj/bn5eVQUBBbc2VgFHo9UWkr//iPwU3LGzbA2rWwbl3wWLs2+P1oaHhn2+HDg2GKnYOw/TFuXHyfQXqk0OtOpa20GzYsCLGZM7t+3x127eoahOvWwcqVcN99Xe/1U+8w6yj0ulNpK4djBhMnBo/3v7/re83Ng+8dtn8dOzaznydhFHrdqbSVVETdO+zcS1TvMBIKvc5U2kq6DKR32P5VvcO0UOh1ptJW4tCf3mHnclm9w5REEnpm9kvgLGCHu8/u4X0DbgbOBPYCF7l7TRTHjpRKW8km6h2mRVQ9vf8Efgrc0cv7HwVmhI8TgZ+FX7OHSlvJJVH0DktLD72qnIDeYSSh5+7PmFlFH5ucC9wRrovxgpkVm1mpu2+L4viRUGkr+SDK3mFlJfz4x8F6InkkU+f0pgCbOr3eHH7vkNAzsyXAEoBp06ZlpHGASlvJfwPpHa5dGwzBO3gQfv3reNqbJpkKvZ5WXelxRSJ3XwoshWBhoHQ2qtNBVdpKcvXWO7z+erjhBjj7bDj//NiaF7VMTTiwGZja6XUZsDVDxz48TSMlcqhvfQtOPBEuuww2b467NZHJVOg9AFxogZOApqw6n6fSVuRQhYVBidvcDBddFEy5lQcytdj3w8AGYD1wG3B5FMeNhEpbkd7NmBFczHjiCbjllrhbE4mort5+5jDvO/DlKI4VOV21FenbJZfA738P11wDixfD7ENuxc0pmkRUpa1I38zg9tuDm5n/9m/hwIG4W5SSZIeeSluR/jnqKPjFL4KbnL/97bhbk5Jkh56u2or039lnw5Il8MMfwtNPx92aQUt26Km0FRmYH/0Ijj4aLrwQGhvjbs2gJDf0VNqKDNzo0cEIjS1b4Ior4m7NoCQ39FTaigzOiScGa4jceSf85jdxt2bAkht6Km1FBi+HR2skM/RU2oqkJodHayQz9FTaiqQuR0drJDP0VNqKROOSS+Css4LRGqtWxd2afkle6Km0FYlODo7WSF7oqbQViVaOjdZIXuiptBWJXg6N1khW6Km0FUmfHBmtkazQU2krkj45MlojWaGn0lYkvXJgtEZUMyefYWbrzGy9mV3Tw/unmlmTmS0PH9dFcdwBUWkrkhlZPloj5dAzswLgXwkW9D4e+IyZHd/Dps+6+9zw8Z1UjztgKm1FMiPLR2tE0dNbCKx39w3u3gzcTbC4d3ZRaTso/2/9Ln702Dqa9h2MuymSS7J4tEYUodfbQt7dnWxmK8zsETOb1dvOzGyJmVWZWdXOnTsjaB4qbQehad9Brv6flXz29he59cn1nPGTZ3ju9V1xN0tySZaO1ogi9PqzkHcNUO7uc4Bbgft625m7L3X3SnevnDhxYgTNQ6XtAD3+6nY+/OOn+Z+azXzp1KP57aUnM3JYAX/3ixe57v5V7G1uibuJkguydLRGFKF32IW83X23u+8Jnz8MFJrZhAiO3T8qbfvlzT0HuOK/a7jkjirGjxrOfZcv4uozZrJw+nge+sopfHHRdO54vo4zb36W6rr6uJsruSALR2tEEXovATPMbLqZDQMuIFjcu4OZTTIzC58vDI/7ZgTHPjyVtofl7ty/fAuLb3qax1Zv58oPvZsHrljEe8rGdWwzorCA684+nrsuOYmWNudTP3+eGx9Zy4GW1hhbLjkhy0ZrpBx67t4CXAH8AVgD/NbdV3db7Ps8YJWZrQBuAS4I18JNP5W2fdrWtI//86sqvnr3ciomjOKhr7yfvz99BoUFPf9qnHz0ETz6tQ9w/glT+fnTf+GcW//M6q1NGW615JwsGq1hmcqewaisrPSqqqrUdnLNNcF/8DfeUE+vk62N+/jD6je46bHXaGlz/u9HjuWi91VQMKSnU7Q9e2rtDq6+dyX1bzfz1dNn8KVTj2ZoL2EpwosvwqJFcMEFwciNNDKzanev7PG9vA49dzjmmODy+aOPRtewHNPW5qzfuYdlf62nqrael2ob2NK4D4D3HX0EN37yvUw7YuSg9t24t5nr7l/NAyu2MmdqMT/61ByOOXJ0lM2XfHL99XDDDXD33XD++Wk7THJDr6YGFiwIriBdfHF0DctyB1paWbWliZdqG3jpr/VU1TV03Gc3ccxwFlaMp7KihBMqxjNr8ljC060peWjlNv7xvlfY29zK1WfM5KL3VTBkAL1GSYiDB+GUU2DdOnjlFSgrS8th+gq9oWk5YrZIyFXb3fsPUlPXwEthL27FpkYOtAR3wb9r4ijOmDWJE6aP54SKEqaNHxlJyHX3sfeWcsL0Er557yt858FXeezVN/jBeXOYOn5wPUjJU+2jNebODUZrPPYYDMnsKZH87enlcWm7fff+IOD+GoTc2jd20+ZQMMSYPXksJ1SMpzLszU0YPTyjbXN37qnezHd+/yoA3z7rOD5dOTUtQSs5bOlSuPRSuOkm+Id/iHz3ySxv86S0dXf+svPtsBcXPDbVB+fjigoLmF9ezAkV4zmhYjzzphUzclh2dN43N+zlqntW8vyGN/ngzCO5/uxZgz5vKHnIHc45Bx5/HKqqYPbsSHefzNDL0au2B1vbWLWliaraoFytqmug/u1mAI4YNSzsxQXn446fPLbXW0uyQVubc8fztdz46Fr2H2yjdNwI5peXMH9aCQvKSzi+dCzDhmZv+yXNtm+H97wHSkth2TIYHl1VkrzQy6HS9u0DLdRsbOi46PDypgb2HwzOx5UfMTLsxQUhN33CqJwsEzfV7+WJNdup2dhIdd07V46HDx3Ce8vGMb+8hAXTSphfnvlyXGL2+98HPb6rroLvfz+y3SYv9LK4tN351gGqautZVltPVW0Dr27bTWubM8Tg+MljqSwfz8Lp46ksL+HIsSPibm5abN+9n5q6BqrrGqje2MCqLU0cbA1+D8uPGMmCaSXMC4Pw2EljBnTvoOSgSy+F226DJ5+EU0+NZJfJC70sKW3dndo393ZcdKiqa+Cvu94Ggl7OvGldz8eNGVEYW1vjtP9gK6u3NgUhWNdAdV0ju/YEg9NHDStg7rTijp7gvGkljCtK5n+nvLVnD8ybF0xIsHIlFBenvMtkhV6MpW1Laxtrtr0V9uKCK6vt/3iLRxaGvbgSKivGM3vyOJ3P6oW7s7lhH9V1DdRsDIJwzbbgCjXAjCNHsyA8Nzi/vISjJ+Zm2S+dRDxaI1n36bWPtb322rQfam9zC8s3NvJSbQNVdfXU1DXwdnMwAL+spIhTZkzoOCd39MTRulm3n8yMqeNHMnX8SD4+L5ia8e0DLazY3NhRFj+y6g3ufimYxrF4ZCHzphYHQVhewpyyYkYNz79f7e7cnTd272dfcyuTi4sYUVgQd5MGr31tjRtuCCYoSOdojbzr6aWxtK1/uzm4olpbz7LaBlZvaaKlzTGDY48aE5yLC0OudFxRpMeWrtranA273qamU2/w9R17ABhicFzp2I7e4ILyEspKinK+N7jnQAsrNzXy8qZGlm9qZMWmRna89c4cdRNGD2dKSRFlJUWUFRcxpaSIKZ2+Zv3pkwhHaySnvI2wtHV3NtXv63J/3F92Bufjhg0dwtyy4o5bR+aX6zxTNmjae5CXNzUEvcGNDSzf2NjR8544ZjjzpwW9wQXlJcyaPC7tPaP2Mr2qLrhoVV3XwMhhBR3ncSsrSigeOazHn21pbWPd9rdYsamJ5ZsaWL6pkdd37KH9n+v0CaOYUzaOuVOLGT2ikC0N+9jSuJctjfvY0rCPrY37aW7tujbFuKLCjhAsC4Mw+DqSKSVFlIwsjP8Pw+uvB6M1Tj45pdEayQm9FK7atrY5a9/YHYxyqGugqrae7buDv6JjRwztGOGwsGI8s6ek/x+MpK61zVn3xltUb2zg5TAI697cC8CwgiHMmjK24wLJgvISjkrxann7Od32kKuqe+d3aPTwocybVsze5lZWbm7suFr97qNGU1kxnoUV4xk2dAjLNzWyfGMjr2xpYt/BILBLRhYyd2oxc6YWMzd89BaW7dranF17DrCpYV9HEG5p3MuWhn1sDr+3t7nrXIgjhxV06RlO6QjGkZSVFDFx9PDMnKKJYLRGckJvAKXt/oOtrNjU2DFetaaugbcOBNOgl44b8c79cdPH8+4jx+h8XJ7Y+dYBXt4YBGBNXQMrNzd1jFOeUlwU3jNYzILy8cwsHdPnzd97DrTwcniPZXVdPS9vbOwIkinFRSwoL6GyIgjUmZPGdtx60/67V1XXwLK/1nf53Rs2dAizJo/tCLe5U4vTMl7a3Wnce5AtjUEIbm7Y2ykcg0fj3q6LQQ0rGEJp8YhDeojtryeNGxHNzfIRjNZIe+iZ2RnAzUABcLu739jtfQvfPxPYC1zk7jWH2++AQu8wpW3j3uZglEP4V7i3v7aVFSWUlWi4VFI0t7Tx6rbd71wprm3gjd37ARhROIQ5ZcUd5wanTxzF6q27qa7tOuZ5iMHMSWOprAiuzFeWlzC5uP/ndNurjNY2Z+ak7BmlsudAS9ceYhiK7T3FnW91XfNiiMGksSO69BC79xr7XSGlOFojraEXrnv7GvAhgvUyXgI+4+6vdtrmTODvCULvROBmdz/xcPseUOh1K223NO4LB+QHj9e2Bye5CwuM90wZF8w6Ut73eRVJpq2N+zoujtTUNbB6625a2t75dzJyWAHzwt5gZXlJYu+x3H+wlW1N+zuCcXPDvi7h+Mbu/bS2dc2XCaOHMaVkZJcLLWUlvVxsSWG0RrpvWelY9zY8WPu6t6922uZc4I5wivgXzKzYzErdfVsExwfgwH/fReGQAq5tfRfPfO8JtjYFf61HDx/K/PISzpkzmcqK8cydWqzzcdKnycVFTC4u4qz3TgZgX3Mrr2xponbX2xxXOpbjSsdohmiCdVOmTxjF9Amjeny/pbWN7W8dYHP9oaXzq9t28/ia7TS3dL3YMnbE0CAUS4qYUvwuzj/7Amb+4Acc/MhHGHb66ZG0O4rQ62nd2+69uN7Wxj0k9MxsCbAEYNq0af1uRMuUMupHFbPmtS3Mm388S8qDUuO40rEaxiQpKRpWwMLpwfBA6b+hBUOC0raXUr+tzdn19oF3wrBT6bzxzb08/5c3efGI93HXiAcZuWdPdO2KYB/9Wfe2P9sE33RfCiyFoLztbyNGfeD9FA1p5b57v41d8RRUVPT3R0UkBkOGGEeOGcGRY0Ywb1rJIe/7unX4qV+gbdQIhh5zTHTHjWAfh133tp/bpGbBAoY88UessRFOOw1qayPdvYhk0GuvYaedxpDWFoY+/SeYNSuyXWdk3dvw9YUWOAloivJ8XocFC+CPfwyWmFPwieSm114LZltpaYGnnoo08CBz694+DGwA1gO3AZenetxeKfhEcleaAw/y7ebkzqqrYfHiYJqap3SOTyTrRRh4fd2ykr/X3dXjE8kdGejhtcvf0AMFn0guyGDgQb6HHij4RLJZhgMPkhB6oOATyUYxBB4kJfRAwSeSTWIKPEhS6IGCTyQbxBh4kLTQAwWfSJxiDjxIYuiBgk8kDlkQeJDU0AMFn0gmZUngQZJDDxR8IpmQRYEHSQ89UPCJpFOWBR4o9AIKPpHoZWHggULvHQo+kehkaeCBQq8rBZ9I6rI48EChdygFn8jgZXnggUKvZwo+kYHLgcCDFEPPzMab2eNm9nr49dDVPYLtas3sFTNbbmaDnBU0wxR8Iv2XI4EHqff0rgGecPcZwBPh696c5u5ze5vNNCsp+EQOL4cCD1IPvXOBX4XPfwV8PMX9ZR8Fn0jvcizwIPXQO6p9VbPw65G9bOfAY2ZWHS7m3SszW2JmVWZWtXPnzhSbFxEFn8ihcjDwoB+hZ2Z/NLNVPTzOHcBxFrn7fOCjwJfN7AO9bejuS9290t0rJ06cOIBDpJmCT+QdORp40I/Qc/fF7j67h8f9wHYzKwUIv+7oZR9bw687gN8BC6P7CBmk4BPJ6cCD1MvbB4DPh88/D9zffQMzG2VmY9qfAx8GVqV43Pgo+CTJcjzwIPXQuxH4kJm9DnwofI2ZTTazh8NtjgKeM7MVwDLgIXd/NMXjxkvBJ0mUB4EH+bzYdyZoQXFJihwLvGQu9p0J6vFJEuRY4B2OQi9VCj7JZ3kWeKDQi4aCT/JRHgYeKPSio+CTfJKngQcKvWgp+CQf5HHggUIvego+yWV5Hnig0EsPBZ/kogQEHij00kfBJ7kkIYEHCr30UvBJLkhQ4IFCL/0UfJLNEhZ4oNDLDAWfZKMEBh4o9DJHwSfZJKGBBwq9zFLwSTZIcOCBQi/zFHwSp4QHHij04qHgkzgo8ACFXnwUfJJJCrwOqS72/SkzW21mbWbW63q2ZnaGma0zs/Vm1tfauMmi4JNMUOB1kWpPbxXwSeCZ3jYwswLgXwlWQjse+IyZHZ/icfOHgk/SSYF3iKGp/LC7rwEws742Wwisd/cN4bZ3EywS/moqx84r7cG3eDHMnw+lpXG3SPLFli0wbJgCr5OUQq+fpgCbOr3eDJzY28bhYuBLAKZNm5belmWTBQvgySfhhz+E5ua4WyP5Yv58+OY34XgVV+0OG3pm9kdgUg9vfStc+/awu+jhe72uRuTuS4GlECwM1I/954958+DOO+NuhUheO2zoufviFI+xGZja6XUZsDXFfYqIDEombll5CZhhZtPNbBhwAcEi4SIiGZfqLSufMLPNwMnAQ2b2h/D7HYt9u3sLcAXwB2AN8Ft3X51as0VEBifVq7e/A37Xw/e3Amd2ev0w8HAqxxIRiYJGZIhIoij0RCRRFHoikigKPRFJFHPP3vt/zWwnUDeAH5kA7EpTczIlHz4D5MfnyIfPAPnxOQb6GcrdfWJPb2R16A2UmVW5e6+zveSCfPgMkB+fIx8+A+TH54jyM6i8FZFEUeiJSKLkW+gtjbsBEciHzwD58Tny4TNAfnyOyD5DXp3TExE5nHzr6YmI9EmhJyKJkhehlw8LD5nZL81sh5mtirstg2VmU83sKTNbEy4Y9dW42zQYZjbCzJaZ2Yrwc9wQd5sGy8wKzOxlM3sw7rYMlpnVmtkrZrbczKpS3l+un9MLFx56DfgQwYSlLwGfcfecWoPDzD4A7AHucPfZcbdnMMysFCh19xozGwNUAx/Pwf8XBoxy9z1mVgg8B3zV3V+IuWkDZmZfByqBse5+VtztGQwzqwUq3T2SG6zzoafXsfCQuzcD7QsP5RR3fwaoj7sdqXD3be5eEz5/i2D+xCnxtmrgPLAnfFkYPnKud2BmZcDHgNvjbks2yYfQ62nhoZz7h5ZvzKwCmAe8GG9LBicsC5cDO4DH3T0XP8dPgG8AbXE3JEUOPGZm1eHCYSnJh9Ab0MJDkn5mNhq4F/iau++Ouz2D4e6t7j6XYE2XhWaWU6cczOwsYIe7V8fdlggscvf5BGtnfzk8FTRo+RB6Wngoi4TnwO4F7nT3/427Paly90bgT8AZMTdloBYB54Tnw+4GPmhmv463SYMTzsSOu+8gmKl9YSr7y4fQ08JDWSK8APALYI273xR3ewbLzCaaWXH4vAhYDKyNt1UD4+7fdPcyd68g+DfxpLv/XczNGjAzGxVeFMPMRgEfBlK6wyHnQy9fFh4ys7uA54FjzWyzmV0cd5sGYRHwOYJexfLwcebhfigLlQJPmdlKgj+qj7t7zt7ykeOOAp4zsxXAMuAhd380lR3m/C0rIiIDkfM9PRGRgVDoiUiiKPREJFEUek9hmM4AAAAYSURBVCKSKAo9EUkUhZ6IJIpCT0QS5f8DKbK04Jo5OOsAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "source": [
    "ax = gpd.GeoSeries(average).plot()\n",
    "gpd.GeoSeries(orig).plot(ax=ax, color='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 470,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<matplotlib.axes._subplots.AxesSubplot at 0x1332d0760>"
     },
     "metadata": {},
     "execution_count": 470
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAClCAYAAABfoUaZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAN+UlEQVR4nO3db6xkd13H8fdntyWaijZNKyztLsuDhqQaweZmkTQloJVsK7HFiLZGQUKzbGyTgg9koalGCQQfYExNA12RaFFpSbBQZGmhCkESkO5ii1sWzNIUum7rLhSBxgdN4euDmYXb7dy/c2bOmXPer2Ry58w5c36/2T/ve3Lu3DmpKiRJ/bel7QlIkubD4EvSQBh8SRoIgy9JA2HwJWkgDL4kDcQZbU9gNeeee27t3Lmz7WlI0sI4dOjQt6rqvEnrOh38nTt3cvDgwbanIUkLI8k3VlrXyCmdJO9PciLJ4RXWJ8nNSY4m+XKSi5sYV5K0fk2dw/9bYPcq6y8HLhzf9gDvaWhcSdI6NRL8qvos8Pgqm1wJ3FYjXwDOTrKtibElSeszr3P45wOPLFs+Nn7s0TmN32+PPQZ7946+Sn2zfTvceiucc07bM1l48wp+Jjw28VPbkuxhdNqHHTt2zHJO/fDYY/CKV8A3vwmXXtr2bKRmVcHHPgZf/zrce6/Rn9K8gn8M2L5s+QLg+KQNq2o/sB9gaWnJj/JczfLYf+IT8LKXtT0jqXl33w1XXQWXXWb0pzSvX7y6C3jt+N06vwR8t6o8nTMNY6+h2L0bPvIR+MpXRtF/fLUfF2o1Tb0t84PA54EXJjmW5A1J9ibZO97kAPAQcBT4a+APmhh3sIy9hsboN6KRUzpVdc0a6wu4romxBs/Ya6hORd/TO5vmZ+ksEmOvofNIfyoGf1EYe2nE6G+awV8Exl56OqO/KQa/64y9NJnR3zCD32XGXlqd0d8Qg99Vxl5aH6O/bga/i4y9tDFGf10MftcYe2lzjP6aDH6XGHtpOkZ/VQa/K4y91AyjvyKD3wXGXmqW0Z/I4LfN2EuzYfSfweC3ydhLs2X0n8bgt8XYS/Nh9H/E4LfB2EvzZfQBgz9/xl5qh9E3+HNl7KV2DTz6Bn9ejL3UDQOOvsGfB2MvdctAo2/wZ83YS900wOgb/Fky9lK3DSz6Bn9WjL20GAYU/UaCn2R3kq8lOZpk34T1L0/y3ST3j29/3MS4nWXspcUykOhPHfwkW4FbgMuBi4Brklw0YdN/q6oXj29/Nu24nWXspcU0gOg3cYS/CzhaVQ9V1ZPA7cCVDex38Rh7abH1PPpNBP984JFly8fGj53upUkeSPKJJD+30s6S7ElyMMnBkydPNjC9OTH2Uj/0OPpNBD8THqvTlr8EPL+qXgT8FfCRlXZWVfuraqmqls4777wGpjcHxl7ql55Gv4ngHwO2L1u+ADi+fIOq+l5VPTG+fwA4M8m5DYzdPmMv9VMPo99E8O8DLkzygiTPAq4G7lq+QZLnJsn4/q7xuN9uYOx2GXup33oW/amDX1VPAdcD9wBHgA9V1YNJ9ibZO97sN4HDSR4AbgaurqrTT/ssFmMvDUOPop8ud3dpaakOHjzY9jSeydhLw3P33XDVVXDRRXDvvXDOOW3PaKIkh6pqadI6f9N2o4y9NEw9ONI3+Bth7KVhW/DoG/z1MvaSYKGjb/DXw9hLWm5Bo2/w12LsJU2ygNE3+Ksx9pJWs2DRN/grMfaS1mOBom/wJzH2kjZiQaJv8E9n7CVtxgJE3+AvZ+wlTaPj0Tf4pxh7SU3ocPQNPhh7Sc3qaPQNvrGXNAsdjP6wg2/sJc1Sx6I/3OAbe0nz0KHoDzP4xl7SPHUk+sMLvrGX1IYORH9YwTf2ktrUcvSHE3xjL6kLWoz+MIJv7CV1SUvR73/wjb2kLmoh+o0EP8nuJF9LcjTJvgnrk+Tm8fovJ7m4iXHXZOwlddmcoz918JNsBW4BLgcuAq5JctFpm10OXDi+7QHeM+24azL2khbBHKPfxBH+LuBoVT1UVU8CtwNXnrbNlcBtNfIF4Owk2xoYezJjL2mRzCn6ZzSwj/OBR5YtHwNeso5tzgcebWD8Z9q/n9+++Pfh9b8AR86EI5+fyTCS1JyfgT/5MBw+zB0f+ADccEPjIzRxhJ8Jj9UmthltmOxJcjDJwZMnT25uRmecAVXw+Hc293xJasPjj4/atXXrTHbfxBH+MWD7suULgOOb2AaAqtoP7AdYWlqa+E1hTfv2cccb3wh//la48UZ4+9shk77nSFIHVMGb3gQ33wxvfjNc95aZDNNE8O8DLkzyAuC/gauB3zltm7uA65Pczuh0z3erajancwC2bIFbbx3df8c7Rl+NvqQuOj327373zFo1dfCr6qkk1wP3AFuB91fVg0n2jte/FzgAXAEcBf4PeP20467J6EvqujnGHpo5wqeqDjCK+vLH3rvsfgHXNTHWhhh9SV0159hDQ8HvNKMvqWtaiD0MIfhg9CV1R0uxh6EEH4y+pPa1GHsYUvDB6EtqT8uxh6EFH4y+pPnrQOxhiMEHoy9pfjoSexhq8MHoS5q9DsUehhx8MPqSZqdjsYehBx+MvqTmdTD2YPBHjL6kpnQ09mDwf8zoS5pWh2MPBv/pjL6kzep47MHgP5PRl7RRCxB7MPiTGX1J67UgsQeDvzKjL2ktCxR7MPirM/qSVrJgsQeDvzajL+l0Cxh7MPjrY/QlnbKgsQeDv35GX9ICxx4M/sYYfWm4Fjz2YPA3zuhLw9OD2IPB3xyjLw1HT2IPUwY/yTnAHcBO4GHgt6rqOxO2exj4PvAD4KmqWppm3E4w+lL/9Sj2MP0R/j7gX6rqXUn2jZffssK2r6iqb005XrcYfam/ehZ7mD74VwIvH9//O+AzrBz8fjL6Uv/0MPYwffCfU1WPAlTVo0l+doXtCvhkkgJurar9K+0wyR5gD8COHTumnN6cGH2pP3oae1hH8JPcCzx3wqobNzDOJVV1fPwN4VNJvlpVn5204fibwX6ApaWl2sAY7TL60uLrcexhHcGvqstWWpfkf5JsGx/dbwNOrLCP4+OvJ5LcCewCJgZ/oRl9aXH1PPYAW6Z8/l3A68b3Xwd89PQNkpyV5Nmn7gOvBA5POW53nYr+tdeOon/TTaN/SJK6awCxh+nP4b8L+FCSNwDfBF4DkOR5wPuq6grgOcCdGf3hnQH8Y1XdPeW43eaRvrQ4BhJ7mDL4VfVt4FcmPH4cuGJ8/yHgRdOMs5CMvtR9A4o9+Ju2s2X0pe4aWOzB4M+e0Ze6Z4CxB4M/H0Zf6o6Bxh4M/vwYfal9A449GPz5MvpSewYeezD482f0pfkz9oDBb4fRl+bH2P+IwW+L0Zdmz9g/jcFvk9GXZsfYP4PBb5vRl5pn7Ccy+F1g9KXmGPsVGfyuMPrS9Iz9qgx+lxh9afOM/ZoMftcYfWnjjP26GPwuMvrS+hn7dTP4XWX0pbUZ+w0x+F1m9KWVGfsNM/hdZ/SlZzL2m2LwF4HRl37M2G+awV8URl8y9lMy+IvE6GvIjP3Utkzz5CSvSfJgkh8mWVplu91JvpbkaJJ904w5eKeif+21o+jfdNPoP4LUZ8a+EdMe4R8GfgO4daUNkmwFbgF+FTgG3Jfkrqr6ypRjD5dH+hoSY9+YqYJfVUcAsvof/i7gaFU9NN72duBKwOBP4/ToV8GNN7Y7J6lpVfC2txn7hszjHP75wCPLlo8BL5nDuP23PPrvfOfoJvWRsW/EmsFPci/w3Amrbqyqj65jjEl/QyuedE6yB9gDsGPHjnXsfuBORf/SS+Gxx9qejdS87dvh6quNfQPWDH5VXTblGMeA7cuWLwCOrzLefmA/wNLSkj+NXI8tW+C1r217FpI6bqp36azTfcCFSV6Q5FnA1cBdcxhXkrTMtG/LfHWSY8BLgY8nuWf8+POSHACoqqeA64F7gCPAh6rqwemmLUnaqGnfpXMncOeEx48DVyxbPgAcmGYsSdJ05nFKR5LUAakO/5ZmkpPANzb59HOBbzU4nUXga+6/ob1e8DVv1POr6rxJKzod/GkkOVhVK37cQx/5mvtvaK8XfM1N8pSOJA2EwZekgehz8Pe3PYEW+Jr7b2ivF3zNjentOXxJ0tP1+QhfkrRM74I/xIutJHl/khNJDrc9l3lIsj3Jp5McGV+A54a25zRrSX4iyReTPDB+zX/a9pzmJcnWJP+R5J/bnss8JHk4yX8muT/JwUb33adTOuOLrfwXyy62AlzT94utJHkZ8ARwW1X9fNvzmbUk24BtVfWlJM8GDgFX9fnvOaOLTpxVVU8kORP4HHBDVX2h5anNXJI/BJaAn66qV7U9n1lL8jCwVFWN/+5B347wf3Sxlap6Ejh1sZVeq6rPAo+3PY95qapHq+pL4/vfZ/QZTee3O6vZqpEnxotnjm/9OVpbQZILgF8D3tf2XPqgb8GfdLGVXodg6JLsBH4R+Pd2ZzJ741Mb9wMngE9VVe9fM/CXwB8BP2x7InNUwCeTHBpfH6QxfQv+hi62osWW5KeADwNvqqrvtT2fWauqH1TVixldU2JXkl6fvkvyKuBEVR1qey5zdklVXQxcDlw3PmXbiL4Ff0MXW9HiGp/H/jDwD1X1T23PZ56q6n+BzwC7W57KrF0C/Pr4nPbtwC8n+ft2pzR7408bpqpOMPo04l1N7btvwfdiKwMw/gHm3wBHquov2p7PPCQ5L8nZ4/s/CVwGfLXdWc1WVb21qi6oqp2M/i//a1X9bsvTmqkkZ43fiECSs4BXAo29+65XwR/qxVaSfBD4PPDCJMeSvKHtOc3YJcDvMTriu398u2KtJy24bcCnk3yZ0YHNp6pqEG9THJjnAJ9L8gDwReDjVXV3Uzvv1dsyJUkr69URviRpZQZfkgbC4EvSQBh8SRoIgy9JA2HwJWkgDL4kDYTBl6SB+H9Kf9OpmY1SDgAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "source": [
    "distance = 1\n",
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 1), (2, 1), (3, 1), (4, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "points = []\n",
    "ids = []\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "for ix, line in enumerate(extended_lines):\n",
    "    current_dist = distance\n",
    "    line_length = line.length\n",
    "\n",
    "    while current_dist < line_length:\n",
    "        points.append(list(line.interpolate(current_dist).coords)[0])\n",
    "        ids.append(ix)\n",
    "        current_dist += distance\n",
    "\n",
    "voronoi_diagram = Voronoi(np.array(points))\n",
    "\n",
    "regions_gdf = _regions(voronoi_diagram, ids)\n",
    "tessellation = regions_gdf[[\"ids\", \"geometry\"]].dissolve(by=\"ids\", as_index=False)\n",
    "linestring = tessellation.geometry.exterior[0]\n",
    "clipped = linestring.intersection(poly)\n",
    "if hasattr(clipped, '__len__'):\n",
    "    lengths = []\n",
    "    for part in clipped:\n",
    "        lengths.append(part.length)\n",
    "    clipped = clipped[lengths.index(max(lengths))]\n",
    "coords = list(clipped.coords)\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)\n",
    "\n",
    "ax = gpd.GeoSeries(lines).plot(color='r')\n",
    "gpd.GeoSeries(average).plot(ax=ax)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 472,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "22.3 ms ± 337 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "points = []\n",
    "ids = []\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "for ix, line in enumerate(extended_lines):\n",
    "    current_dist = distance\n",
    "    line_length = line.length\n",
    "\n",
    "    while current_dist < line_length:\n",
    "        points.append(list(line.interpolate(current_dist).coords)[0])\n",
    "        ids.append(ix)\n",
    "        current_dist += distance\n",
    "\n",
    "voronoi_diagram = Voronoi(np.array(points))\n",
    "\n",
    "regions_gdf = _regions(voronoi_diagram, ids)\n",
    "tessellation = regions_gdf[[\"ids\", \"geometry\"]].dissolve(by=\"ids\", as_index=False)\n",
    "linestring = tessellation.geometry.exterior[0]\n",
    "clipped = linestring.intersection(poly)\n",
    "if hasattr(clipped, '__len__'):\n",
    "    lengths = []\n",
    "    for part in clipped:\n",
    "        lengths.append(part.length)\n",
    "    clipped = clipped[lengths.index(max(lengths))]\n",
    "coords = list(clipped.coords)\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 473,
   "metadata": {},
   "outputs": [],
   "source": [
    "polygons = []\n",
    "for point_region in voronoi_diagram.point_region:\n",
    "    region = voronoi_diagram.regions[point_region]\n",
    "    if -1 not in region:\n",
    "        polygons.append(Polygon(voronoi_diagram.vertices[region]))\n",
    "    else:\n",
    "        polygons.append(None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 483,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<shapely.geometry.polygon.LinearRing at 0x134827700>",
      "image/svg+xml": "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"-0.24919796515140274 -0.20636609208842577 5.5718844863874954 2.0272556666524197\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,1.614523482475568)\"><polyline fill=\"none\" stroke=\"#66cc99\" stroke-width=\"0.11143768972774991\" points=\"4.237499933207152,1.2678930852277572 4.331442318208224,1.280486018754476 5.064185473657881,0.7280881851489126 5.116320429147667,0.030344170978790275 5.10013238199036,0.0 3.6035533905932735,0.0 3.085786437626905,0.0 2.0857864376269046,0.0 1.4309644062711508,0.0 0.0004402470253619395,1.1102230246251565e-16 -0.04283187306297698,0.14747595138299852 0.21720035761886924,1.0973618550409823 0.9149214894707098,1.3753753291912947 0.9814893563858217,1.3484251496559871 1.7305380699108404,1.614523482475568 2.085786437626905,1.3854765175244317 2.7305380699108404,1.614523482475568 3.085786437626905,1.385476517524432 3.719062072322486,1.610446653064392 4.237499933207152,1.2678930852277572\" opacity=\"0.8\" /></g></svg>"
     },
     "metadata": {},
     "execution_count": 483
    }
   ],
   "source": [
    "from_shapely(polygons)[np.array(ids) == 1].unary_union().exterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 508,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "4.26 ms ± 70.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "points = []\n",
    "ids = []\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "for ix, line in enumerate(extended_lines):\n",
    "    current_dist = distance\n",
    "    line_length = line.length\n",
    "\n",
    "    while current_dist < line_length:\n",
    "        points.append(list(line.interpolate(current_dist).coords)[0])\n",
    "        ids.append(ix)\n",
    "        current_dist += distance\n",
    "\n",
    "voronoi_diagram = Voronoi(np.array(points))\n",
    "\n",
    "polygons = []\n",
    "for point_region in voronoi_diagram.point_region:\n",
    "    region = voronoi_diagram.regions[point_region]\n",
    "    if -1 not in region:\n",
    "        polygons.append(Polygon(voronoi_diagram.vertices[region]))\n",
    "    else:\n",
    "        polygons.append(None)\n",
    "\n",
    "linestring = from_shapely(polygons)[np.array(ids) == 1].unary_union().exterior\n",
    "\n",
    "clipped = linestring.intersection(poly)\n",
    "if hasattr(clipped, '__len__'):\n",
    "    lengths = [part.length for part in clipped]\n",
    "    clipped = clipped[lengths.index(max(lengths))]\n",
    "coords = list(clipped.coords)\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 502,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "116 µs ± 38.8 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "poly = list(polygonize(lines))[0]\n",
    "points = []\n",
    "ids = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 503,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "110 µs ± 67.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 504,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "1.29 ms ± 87 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "for ix, line in enumerate(extended_lines):\n",
    "    current_dist = distance\n",
    "    line_length = line.length\n",
    "\n",
    "    while current_dist < line_length:\n",
    "        points.append(list(line.interpolate(current_dist).coords)[0])\n",
    "        ids.append(ix)\n",
    "        current_dist += distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 511,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "693 µs ± 26 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "voronoi_diagram = Voronoi(np.array(points))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 513,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "406 µs ± 35.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "polygons = []\n",
    "for point_region in voronoi_diagram.point_region:\n",
    "    region = voronoi_diagram.regions[point_region]\n",
    "    if -1 not in region:\n",
    "        polygons.append(Polygon(voronoi_diagram.vertices[region]))\n",
    "    else:\n",
    "        polygons.append(None)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 515,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "1.69 ms ± 458 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "linestring = from_shapely(polygons)[np.array(ids) == 1].unary_union().exterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 517,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "150 µs ± 11.5 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "clipped = linestring.intersection(poly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 522,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "143 µs ± 3.87 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "clipped = linestring.intersection(poly)\n",
    "if hasattr(clipped, '__len__'):\n",
    "    lengths = [part.length for part in clipped]\n",
    "    clipped = clipped[lengths.index(max(lengths))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 523,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "29.2 µs ± 8.06 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "coords = list(clipped.coords)\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 527,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygeos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 528,
   "metadata": {},
   "outputs": [],
   "source": [
    "py_line = pygeos.io.from_shapely(line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 529,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "array([<pygeos.Geometry POINT (0.707 -0.707)>,\n       <pygeos.Geometry POINT (1.59 -1)>,\n       <pygeos.Geometry POINT (2.59 -1)>,\n       <pygeos.Geometry POINT (3.59 -1)>,\n       <pygeos.Geometry POINT (4.41 -0.586)>], dtype=object)"
     },
     "metadata": {},
     "execution_count": 529
    }
   ],
   "source": [
    "pygeos.linear.line_interpolate_point(py_line, [1, 2, 3, 4, 5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 530,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "[<shapely.geometry.polygon.LinearRing at 0x134b96c70>,\n <shapely.geometry.linestring.LineString at 0x1332e36a0>,\n <shapely.geometry.linestring.LineString at 0x1332e3730>]"
     },
     "metadata": {},
     "execution_count": 530
    }
   ],
   "source": [
    "extended_lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 537,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "array([0.        , 1.45710678, 2.91421356, 4.37132034, 5.82842712])"
     },
     "metadata": {},
     "execution_count": 537
    }
   ],
   "source": [
    "np.linspace(0, line.length, num=int(line.length//1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 535,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3"
     },
     "metadata": {},
     "execution_count": 535
    }
   ],
   "source": [
    "pygeos.linear.line_interpolate_point(py_line, [1, 2, 3, 4, 5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 597,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "1.25 ms ± 298 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "points = []\n",
    "ids = []\n",
    "for ix, line in enumerate(extended_lines):\n",
    "    current_dist = distance\n",
    "    line_length = line.length\n",
    "\n",
    "    while current_dist < line_length:\n",
    "        points.append(list(line.interpolate(current_dist).coords)[0])\n",
    "        ids.append(ix)\n",
    "        current_dist += distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 609,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "333 µs ± 88 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "points = np.empty((0, 2))\n",
    "ids = []\n",
    "pygeos_lines = pygeos.io.from_shapely(extended_lines)\n",
    "lengths = pygeos.measurement.length(pygeos_lines)\n",
    "for ix, (line, length) in enumerate(zip(pygeos_lines, lengths)):\n",
    "    pts = pygeos.linear.line_interpolate_point(line, \n",
    "        np.linspace(0, length, num=int(length//distance))\n",
    "        )\n",
    "    points = np.append(points, pygeos.coordinates.get_coordinates(pts), axis=0)\n",
    "    ids += [ix] * len(pts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "2.89 ms ± 464 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "\n",
    "distance = 1\n",
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 1), (2, 1), (3, 1), (4, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "points = np.empty((0, 2))\n",
    "ids = []\n",
    "\n",
    "pygeos_lines = pygeos.io.from_shapely(extended_lines)\n",
    "lengths = pygeos.measurement.length(pygeos_lines)\n",
    "for ix, (line, length) in enumerate(zip(pygeos_lines, lengths)):\n",
    "    pts = pygeos.linear.line_interpolate_point(line, \n",
    "        np.linspace(0, length, num=int(length//distance))\n",
    "        )\n",
    "    points = np.append(points, pygeos.coordinates.get_coordinates(pts), axis=0)\n",
    "    ids += [ix] * len(pts)\n",
    "\n",
    "voronoi_diagram = Voronoi(np.array(points))\n",
    "\n",
    "polygons = []\n",
    "for point_region in voronoi_diagram.point_region:\n",
    "    region = voronoi_diagram.regions[point_region]\n",
    "    if -1 not in region:\n",
    "        polygons.append(Polygon(voronoi_diagram.vertices[region]))\n",
    "    else:\n",
    "        polygons.append(None)\n",
    "\n",
    "linestring = from_shapely(polygons)[np.array(ids) == 1].unary_union().exterior\n",
    "\n",
    "clipped = linestring.intersection(poly)\n",
    "if hasattr(clipped, '__len__'):\n",
    "    lengths = [part.length for part in clipped]\n",
    "    clipped = clipped[lengths.index(max(lengths))]\n",
    "coords = list(clipped.coords)\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance = 1\n",
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 1), (2, 1), (3, 1), (4, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "points = np.empty((0, 2))\n",
    "ids = []\n",
    "\n",
    "pygeos_lines = pygeos.io.from_shapely(extended_lines)\n",
    "lengths = pygeos.measurement.length(pygeos_lines)\n",
    "for ix, (line, length) in enumerate(zip(pygeos_lines, lengths)):\n",
    "    pts = pygeos.linear.line_interpolate_point(line, \n",
    "        np.linspace(0, length, num=int(length//distance))\n",
    "        )\n",
    "    points = np.append(points, pygeos.coordinates.get_coordinates(pts), axis=0)\n",
    "    ids += [ix] * len(pts)\n",
    "\n",
    "voronoi_diagram = Voronoi(np.array(points))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "[[-1, 0],\n [0, 4],\n [-1, 4],\n [-1, 2],\n [2, 6],\n [-1, 6],\n [1, 5],\n [-1, 1],\n [-1, 5],\n [8, 11],\n [8, 10],\n [9, 12],\n [9, 10],\n [11, 13],\n [12, 13],\n [-1, 10],\n [2, 8],\n [-1, 9],\n [-1, 12],\n [-1, 7],\n [7, 13],\n [-1, 15],\n [-1, 14],\n [14, 16],\n [15, 16],\n [7, 15],\n [0, 3],\n [3, 14],\n [17, 20],\n [17, 18],\n [20, 21],\n [19, 21],\n [18, 19],\n [1, 17],\n [-1, 18],\n [3, 20],\n [16, 21],\n [6, 19],\n [11, 21],\n [22, 25],\n [22, 23],\n [23, 27],\n [24, 26],\n [24, 25],\n [26, 27],\n [4, 22],\n [-1, 25],\n [20, 23],\n [5, 27],\n [-1, 24],\n [-1, 26]]"
     },
     "metadata": {},
     "execution_count": 88
    }
   ],
   "source": [
    "voronoi_diagram.ridge_vertices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "array([[11, 12],\n       [11, 25],\n       [11, 10],\n       [ 3,  2],\n       [ 3, 18],\n       [ 3,  4],\n       [ 6, 20],\n       [ 6,  5],\n       [ 6,  7],\n       [22, 18],\n       [22,  2],\n       [22,  0],\n       [22,  1],\n       [22, 23],\n       [22, 15],\n       [ 2,  1],\n       [ 2, 18],\n       [ 0,  1],\n       [ 0, 15],\n       [15, 14],\n       [15, 23],\n       [13, 14],\n       [13, 12],\n       [13, 24],\n       [13, 23],\n       [14, 23],\n       [12, 25],\n       [12, 24],\n       [19, 20],\n       [19,  5],\n       [19, 24],\n       [19, 18],\n       [19,  4],\n       [ 5, 20],\n       [ 5,  4],\n       [24, 25],\n       [24, 23],\n       [18,  4],\n       [18, 23],\n       [26, 10],\n       [26, 25],\n       [26, 20],\n       [26,  8],\n       [26,  9],\n       [26,  7],\n       [10, 25],\n       [10,  9],\n       [25, 20],\n       [20,  7],\n       [ 8,  9],\n       [ 8,  7]], dtype=int32)"
     },
     "metadata": {},
     "execution_count": 62
    }
   ],
   "source": [
    "pts = voronoi_diagram.ridge_points\n",
    "pts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(row):\n",
    "    if (ids.index(1) <= row[0] <  ids.index(2)) and (row[1] >=  ids.index(2)):\n",
    "        return True\n",
    "    if (ids.index(1) <= row[1] <  ids.index(2)) and (row[0] >=  ids.index(2)):\n",
    "        return True\n",
    "    return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test(row):\n",
    "    a = ids[row[0]]\n",
    "    b = ids[row[1]]\n",
    "    if a in [1, 2] and b in [1, 2] and a != b:\n",
    "        return True\n",
    "    return False    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "array([False, False, False, False, False, False, False, False, False,\n        True, False, False, False, False, False, False, False, False,\n       False, False, False, False, False, False, False, False, False,\n       False, False, False,  True, False, False, False, False, False,\n       False, False,  True, False, False,  True, False, False, False,\n       False, False,  True, False, False, False])"
     },
     "metadata": {},
     "execution_count": 139
    }
   ],
   "source": [
    "mask = np.apply_along_axis(test, 1, pts)\n",
    "mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [],
   "source": [
    "rigde_vertices = np.array(voronoi_diagram.ridge_vertices)\n",
    "verts = rigde_vertices[mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [],
   "source": [
    "lss = []\n",
    "for f, l in verts:\n",
    "    lss.append(LineString([voronoi_diagram.vertices[f], voronoi_diagram.vertices[l]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<shapely.geometry.linestring.LineString at 0x121ecbd30>",
      "image/svg+xml": "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"0.17267903255066605 -0.17429518882466202 4.705970098265869 1.0840762415072622\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,0.7354858638579381)\"><polyline fill=\"none\" stroke=\"#66cc99\" stroke-width=\"0.09411940196531737\" points=\"0.34697422137532785,0.6819560703481348 1.0008820789065405,0.0 1.771446609406726,-2.220446049250313e-16 3.228553390593273,0.0 3.9991179210934593,0.0 4.704353941991873,0.7354858638579383\" opacity=\"0.8\" /></g></svg>"
     },
     "metadata": {},
     "execution_count": 143
    }
   ],
   "source": [
    "linemerge(lss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = list(linemerge(lss).coords)\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<shapely.geometry.linestring.LineString at 0x121e55e80>",
      "image/svg+xml": "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"-0.2 -0.20000000000000023 5.4 0.40000000000000024\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,-2.220446049250313e-16)\"><polyline fill=\"none\" stroke=\"#66cc99\" stroke-width=\"0.10800000000000001\" points=\"5.0,0.0 1.0008820789065405,0.0 1.771446609406726,-2.220446049250313e-16 3.228553390593273,0.0 3.9991179210934593,0.0 0.0,0.0\" opacity=\"0.8\" /></g></svg>"
     },
     "metadata": {},
     "execution_count": 145
    }
   ],
   "source": [
    "average"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [],
   "source": [
    "def edge(row): \n",
    "    a = ids[row[0]]\n",
    "    b = ids[row[1]]\n",
    "    if a in [1, 2] and b in [1, 2] and a != b:\n",
    "        return True\n",
    "    return False  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "2.46 ms ± 262 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "\n",
    "distance = 1\n",
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 3), (2, 2.5), (3, 1.2), (4.3, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "points = np.empty((0, 2))\n",
    "ids = []\n",
    "\n",
    "pygeos_lines = pygeos.io.from_shapely(extended_lines)\n",
    "lengths = pygeos.measurement.length(pygeos_lines)\n",
    "for ix, (line, length) in enumerate(zip(pygeos_lines, lengths)):\n",
    "    pts = pygeos.linear.line_interpolate_point(line, \n",
    "        np.linspace(0, length, num=int(length//distance))\n",
    "        )\n",
    "    points = np.append(points, pygeos.coordinates.get_coordinates(pts), axis=0)\n",
    "    ids += [ix] * len(pts)\n",
    "\n",
    "voronoi_diagram = Voronoi(np.array(points))\n",
    "\n",
    "pts = voronoi_diagram.ridge_points\n",
    "\n",
    "mask = np.apply_along_axis(edge, axis=1, arr=pts)\n",
    "rigde_vertices = np.array(voronoi_diagram.ridge_vertices)\n",
    "verts = rigde_vertices[mask]\n",
    "\n",
    "edgeline = linemerge([LineString([voronoi_diagram.vertices[f], voronoi_diagram.vertices[l]]) for f, l in verts])\n",
    "if hasattr(edgeline, '__len__'):\n",
    "    lengths = [part.length for part in edgeline]\n",
    "    edgeline = edgeline[lengths.index(max(lengths))]\n",
    "coords = list(edgeline.coords)\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "<shapely.geometry.linestring.LineString at 0x118c09a90>",
      "image/svg+xml": "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"-0.2 -0.22797463279663638 5.4 0.8236731551167618\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,0.36772388952348906)\"><polyline fill=\"none\" stroke=\"#66cc99\" stroke-width=\"0.10800000000000001\" points=\"5.0,0.0 3.97806558865857,0.02195533473655309 3.553022986316233,0.09674918196636484 3.2285533905932735,-0.027974632796636367 1.771446609406726,0.38505711494923256 1.7602605914494638,0.3956985223201255 1.2211026809627277,0.22966657252741365 0.0,0.0\" opacity=\"0.8\" /></g></svg>"
     },
     "metadata": {},
     "execution_count": 187
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "2.34 ms ± 340 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "\n",
    "distance = 1\n",
    "start = (0, 0)\n",
    "end = (5, 0)\n",
    "top = LineString([start, (1, 3), (2, 2.5), (3, 1.2), (4.3, 1), end])\n",
    "bottom = LineString([start, (1, -1), (2, -1), (3, -1), (4, -1), end])\n",
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "points = np.empty((0, 2))\n",
    "ids = []\n",
    "\n",
    "pygeos_lines = pygeos.from_shapely(extended_lines)\n",
    "lengths = pygeos.measurement.length(pygeos_lines)\n",
    "for ix, (line, length) in enumerate(zip(pygeos_lines, lengths)):\n",
    "    pts = pygeos.linear.line_interpolate_point(line, \n",
    "        np.linspace(0, length, num=int(length//distance))\n",
    "        )\n",
    "    points = np.append(points, pygeos.get_coordinates(pts), axis=0)\n",
    "    ids += [ix] * len(pts)\n",
    "\n",
    "voronoi_diagram = Voronoi(np.array(points))\n",
    "\n",
    "pts = voronoi_diagram.ridge_points\n",
    "\n",
    "mask = np.apply_along_axis(edge, axis=1, arr=pts)\n",
    "rigde_vertices = np.array(voronoi_diagram.ridge_vertices)\n",
    "verts = rigde_vertices[mask]\n",
    "\n",
    "strings = pygeos.linestrings([[voronoi_diagram.vertices[f], voronoi_diagram.vertices[l]] for f, l in verts])\n",
    "edgeline = pygeos.line_merge(pygeos.union_all(strings))\n",
    "\n",
    "\n",
    "coords = list(pygeos.get_coordinates(edgeline))\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "metadata": {},
   "outputs": [],
   "source": [
    "realdf = gpd.read_file('/Users/martin/Dropbox/Python/Jupyter/edges.geojson').to_crs(5514)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {},
   "outputs": [],
   "source": [
    "top = realdf.loc[0, 'geometry']\n",
    "bottom = realdf.loc[1, 'geometry']\n",
    "start = top.coords[0]\n",
    "end = top.coords[-1]\n",
    "distance = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 223,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "308 ms ± 106 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "# FULL PYGEOS\n",
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "points = np.empty((0, 2))\n",
    "ids = []\n",
    "\n",
    "pygeos_lines = pygeos.from_shapely(extended_lines)\n",
    "lengths = pygeos.measurement.length(pygeos_lines)\n",
    "for ix, (line, length) in enumerate(zip(pygeos_lines, lengths)):\n",
    "    pts = pygeos.linear.line_interpolate_point(line, \n",
    "        np.linspace(0, length, num=int(length//distance))\n",
    "        )\n",
    "    points = np.append(points, pygeos.get_coordinates(pts), axis=0)\n",
    "    ids += [ix] * len(pts)\n",
    "\n",
    "voronoi_diagram = Voronoi(np.array(points))\n",
    "\n",
    "pts = voronoi_diagram.ridge_points\n",
    "\n",
    "mask = np.apply_along_axis(edge, axis=1, arr=pts)\n",
    "rigde_vertices = np.array(voronoi_diagram.ridge_vertices)\n",
    "verts = rigde_vertices[mask]\n",
    "\n",
    "strings = pygeos.linestrings([[voronoi_diagram.vertices[f], voronoi_diagram.vertices[l]] for f, l in verts])\n",
    "edgeline = pygeos.line_merge(pygeos.union_all(strings))\n",
    "\n",
    "\n",
    "coords = list(pygeos.get_coordinates(edgeline))\n",
    "# resolve start/end join\n",
    "coords[0] = start\n",
    "coords[-1] = end\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "295 ms ± 16.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "# PYGEOS INTERPOLATE + SHAPELY LINEMERGE\n",
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "points = np.empty((0, 2))\n",
    "ids = []\n",
    "\n",
    "pygeos_lines = pygeos.io.from_shapely(extended_lines)\n",
    "lengths = pygeos.measurement.length(pygeos_lines)\n",
    "for ix, (line, length) in enumerate(zip(pygeos_lines, lengths)):\n",
    "    pts = pygeos.linear.line_interpolate_point(line, \n",
    "        np.linspace(0, length, num=int(length//distance))\n",
    "        )\n",
    "    points = np.append(points, pygeos.coordinates.get_coordinates(pts), axis=0)\n",
    "    ids += [ix] * len(pts)\n",
    "\n",
    "voronoi_diagram = Voronoi(np.array(points))\n",
    "\n",
    "pts = voronoi_diagram.ridge_points\n",
    "\n",
    "mask = np.apply_along_axis(edge, axis=1, arr=pts)\n",
    "rigde_vertices = np.array(voronoi_diagram.ridge_vertices)\n",
    "verts = rigde_vertices[mask]\n",
    "\n",
    "edgeline = linemerge([LineString([voronoi_diagram.vertices[f], voronoi_diagram.vertices[l]]) for f, l in verts])\n",
    "if hasattr(edgeline, '__len__'):\n",
    "    lengths = [part.length for part in edgeline]\n",
    "    edgeline = edgeline[lengths.index(max(lengths))]\n",
    "coords = list(edgeline.coords)\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 227,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "533 ms ± 54.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "# PYGEOS INTERPOLATE + SHAPELY INTERSECTION\n",
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "points = np.empty((0, 2))\n",
    "ids = []\n",
    "\n",
    "pygeos_lines = pygeos.io.from_shapely(extended_lines)\n",
    "lengths = pygeos.measurement.length(pygeos_lines)\n",
    "for ix, (line, length) in enumerate(zip(pygeos_lines, lengths)):\n",
    "    pts = pygeos.linear.line_interpolate_point(line, \n",
    "        np.linspace(0, length, num=int(length//distance))\n",
    "        )\n",
    "    points = np.append(points, pygeos.coordinates.get_coordinates(pts), axis=0)\n",
    "    ids += [ix] * len(pts)\n",
    "\n",
    "voronoi_diagram = Voronoi(np.array(points))\n",
    "\n",
    "polygons = []\n",
    "for point_region in voronoi_diagram.point_region:\n",
    "    region = voronoi_diagram.regions[point_region]\n",
    "    if -1 not in region:\n",
    "        polygons.append(Polygon(voronoi_diagram.vertices[region]))\n",
    "    else:\n",
    "        polygons.append(None)\n",
    "\n",
    "linestring = from_shapely(polygons)[np.array(ids) == 1].unary_union().exterior\n",
    "\n",
    "clipped = linestring.intersection(poly)\n",
    "if hasattr(clipped, '__len__'):\n",
    "    lengths = [part.length for part in clipped]\n",
    "    clipped = clipped[lengths.index(max(lengths))]\n",
    "coords = list(clipped.coords)\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "716 ms ± 44.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "# SHAPELY INTERPOLATE + SHAPELY INTERSECTION\n",
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "points = []\n",
    "ids = []\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "for ix, line in enumerate(extended_lines):\n",
    "    current_dist = distance\n",
    "    line_length = line.length\n",
    "\n",
    "    while current_dist < line_length:\n",
    "        points.append(list(line.interpolate(current_dist).coords)[0])\n",
    "        ids.append(ix)\n",
    "        current_dist += distance\n",
    "\n",
    "voronoi_diagram = Voronoi(np.array(points))\n",
    "\n",
    "polygons = []\n",
    "for point_region in voronoi_diagram.point_region:\n",
    "    region = voronoi_diagram.regions[point_region]\n",
    "    if -1 not in region:\n",
    "        polygons.append(Polygon(voronoi_diagram.vertices[region]))\n",
    "    else:\n",
    "        polygons.append(None)\n",
    "\n",
    "linestring = from_shapely(polygons)[np.array(ids) == 1].unary_union().exterior\n",
    "\n",
    "clipped = linestring.intersection(poly)\n",
    "if hasattr(clipped, '__len__'):\n",
    "    lengths = [part.length for part in clipped]\n",
    "    clipped = clipped[lengths.index(max(lengths))]\n",
    "coords = list(clipped.coords)\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 230,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "1.16 s ± 148 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "# SHAPELY INTERPOLATE + SHAPELY INTERSECTION + GEOPANDAS\n",
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "points = []\n",
    "ids = []\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "for ix, line in enumerate(extended_lines):\n",
    "    current_dist = distance\n",
    "    line_length = line.length\n",
    "\n",
    "    while current_dist < line_length:\n",
    "        points.append(list(line.interpolate(current_dist).coords)[0])\n",
    "        ids.append(ix)\n",
    "        current_dist += distance\n",
    "\n",
    "voronoi_diagram = Voronoi(np.array(points))\n",
    "\n",
    "regions_gdf = _regions(voronoi_diagram, ids)\n",
    "tessellation = regions_gdf[[\"ids\", \"geometry\"]].dissolve(by=\"ids\", as_index=False)\n",
    "linestring = tessellation.geometry.exterior[0]\n",
    "clipped = linestring.intersection(poly)\n",
    "if hasattr(clipped, '__len__'):\n",
    "    lengths = []\n",
    "    for part in clipped:\n",
    "        lengths.append(part.length)\n",
    "    clipped = clipped[lengths.index(max(lengths))]\n",
    "coords = list(clipped.coords)\n",
    "coords[0] = end\n",
    "coords[-1] = start\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 232,
   "metadata": {},
   "outputs": [],
   "source": [
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 233,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "33.4 ms ± 11.6 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "points = np.empty((0, 2))\n",
    "ids = []\n",
    "\n",
    "pygeos_lines = pygeos.io.from_shapely(extended_lines)\n",
    "lengths = pygeos.measurement.length(pygeos_lines)\n",
    "for ix, (line, length) in enumerate(zip(pygeos_lines, lengths)):\n",
    "    pts = pygeos.linear.line_interpolate_point(line, \n",
    "        np.linspace(0, length, num=int(length//distance))\n",
    "        )\n",
    "    points = np.append(points, pygeos.coordinates.get_coordinates(pts), axis=0)\n",
    "    ids += [ix] * len(pts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "397 ms ± 200 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "points = []\n",
    "ids = []\n",
    "\n",
    "for ix, line in enumerate(extended_lines):\n",
    "    current_dist = distance\n",
    "    line_length = line.length\n",
    "\n",
    "    while current_dist < line_length:\n",
    "        points.append(list(line.interpolate(current_dist).coords)[0])\n",
    "        ids.append(ix)\n",
    "        current_dist += distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 371,
   "metadata": {},
   "outputs": [],
   "source": [
    "# helper functions\n",
    "def get_ids(x):\n",
    "    return ids[x]\n",
    "mp = np.vectorize(get_ids)\n",
    "\n",
    "def dist(p1, p2):\n",
    "    return np.sqrt(((p1[0]-p2[0])**2)+((p1[1]-p2[1])**2))\n",
    "\n",
    "def get_verts(x):\n",
    "    return voronoi_diagram.vertices[x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 277,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "274 ms ± 27.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "# PYGEOS INTERPOLATE + SHAPELY LINEMERGE = OPTIMUM\n",
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "points = []\n",
    "ids = []\n",
    "\n",
    "pygeos_lines = pygeos.from_shapely(extended_lines)\n",
    "lengths = pygeos.length(pygeos_lines)\n",
    "for ix, (line, length) in enumerate(zip(pygeos_lines, lengths)):\n",
    "    pts = pygeos.line_interpolate_point(\n",
    "        line, np.linspace(0, length, num=int(length // distance))\n",
    "    )\n",
    "    points.append(pygeos.get_coordinates(pts))\n",
    "    ids += [ix] * len(pts)\n",
    "\n",
    "voronoi_diagram = Voronoi(np.concatenate(points))\n",
    "\n",
    "pts = voronoi_diagram.ridge_points\n",
    "\n",
    "mask = np.apply_along_axis(edge, axis=1, arr=pts)\n",
    "rigde_vertices = np.array(voronoi_diagram.ridge_vertices)\n",
    "verts = rigde_vertices[mask]\n",
    "\n",
    "edgeline = linemerge(\n",
    "    [\n",
    "        LineString([voronoi_diagram.vertices[f], voronoi_diagram.vertices[l]])\n",
    "        for f, l in verts\n",
    "    ]\n",
    ")\n",
    "if hasattr(edgeline, \"__len__\"):\n",
    "    lengths = [part.length for part in edgeline]\n",
    "    edgeline = edgeline[lengths.index(max(lengths))]\n",
    "coords = list(edgeline.coords)\n",
    "\n",
    "if dist(coords[0], start) > dist(coords[0], end):\n",
    "    coords[0] = end\n",
    "    coords[-1] = start\n",
    "else:\n",
    "    coords[0] = start\n",
    "    coords[-1] = end\n",
    "average = LineString(coords)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 305,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "213 ms ± 37.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "edgeline = linemerge(\n",
    "    [\n",
    "        LineString([voronoi_diagram.vertices[f], voronoi_diagram.vertices[l]])\n",
    "        for f, l in verts\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 301,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "77.9 ms ± 21 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "strings = pygeos.linestrings([[voronoi_diagram.vertices[f], voronoi_diagram.vertices[l]] for f, l in verts])\n",
    "edgeline = pygeos.line_merge(pygeos.union_all(strings))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 296,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_verts(x):\n",
    "    return voronoi_diagram.vertices[x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 304,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "48.7 ms ± 7.48 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "strings = pygeos.linestrings(get_verts(verts))\n",
    "edgeline = pygeos.line_merge(pygeos.union_all(strings))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# helper functions\n",
    "def get_ids(x):\n",
    "    return ids[x]\n",
    "mp = np.vectorize(get_ids)\n",
    "\n",
    "def dist(p1, p2):\n",
    "    return np.sqrt(((p1[0]-p2[0])**2)+((p1[1]-p2[1])**2))\n",
    "\n",
    "def get_verts(x):\n",
    "    return voronoi_diagram.vertices[x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 401,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "131 ms ± 31.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
    }
   ],
   "source": [
    "%%timeit\n",
    "# FULL PYGEOS - SUPER PERFORMANCE\n",
    "lines = [top, bottom]\n",
    "\n",
    "poly = list(polygonize(lines))[0]\n",
    "\n",
    "extended_lines = [poly.buffer(distance).exterior] + lines\n",
    "\n",
    "points = np.empty((0, 2))\n",
    "ids = []\n",
    "\n",
    "pygeos_lines = pygeos.from_shapely(extended_lines)\n",
    "lengths = pygeos.length(pygeos_lines)\n",
    "for ix, (line, length) in enumerate(zip(pygeos_lines, lengths)):\n",
    "    pts = pygeos.line_interpolate_point(line, \n",
    "        np.linspace(0, length, num=int(length//distance))\n",
    "        )\n",
    "    points = np.append(points, pygeos.get_coordinates(pts), axis=0)\n",
    "    ids += [ix] * len(pts)\n",
    "\n",
    "voronoi_diagram = Voronoi(points)\n",
    "\n",
    "pts = voronoi_diagram.ridge_points\n",
    "\n",
    "mapped = mp(pts)\n",
    "mask = np.isin(mapped[:, 0], [1, 2]) & np.isin(mapped[:, 1], [1, 2]) & (mapped[:, 0] != mapped[:, 1])\n",
    "rigde_vertices = np.array(voronoi_diagram.ridge_vertices)\n",
    "verts = rigde_vertices[mask]\n",
    "\n",
    "edgeline = pygeos.line_merge(pygeos.multilinestrings(get_verts(verts)))\n",
    "coords = pygeos.get_coordinates(edgeline)\n",
    "if dist(coords[0], start) > dist(coords[0], end):\n",
    "    coords[0] = end\n",
    "    coords[-1] = start\n",
    "else:\n",
    "    coords[0] = start\n",
    "    coords[-1] = end\n",
    "average = LineString(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "[(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)]"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "source": [
    "list(itertools.combinations(range(4), 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.2-final"
  },
  "orig_nbformat": 2,
  "kernelspec": {
   "name": "geo_dev",
   "display_name": "geo_dev"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}