polygonlength.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "polygonlength.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyO61y7OHSB8XdCwyndP4Qh1",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/felluksch/length_of_a_polygon/blob/main/polygonlength.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pdAjIt-_j3Db"
      },
      "source": [
        "In this notebook, we're trying to estimate the \"length of a polygon\". This is not a well defined metric; one might describe length in terms of the bounding box of a polygon or approximate a curve that follows the primary axis of the polygon. Depending on the shape of the polygon, humans might intuitively give preference to one of the definitions; e.g. if the polygon approximates a curve, we would describe the \"length\" in terms of that curve. If the shape is more complex, however, most people would probably simplify the geometry mentally and simply draw a straight line along the principal axis of the polygon.\n",
        "\n",
        "Here, we try a method more akin to the \"curve fitting\" approach. We compute the skeleton of the polygon and take the longest path in the skeleton graph as the length of the polygon.\n",
        "\n",
        "But first we need to load the required packages."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "uAI8OvaJl6z9",
        "outputId": "88690336-0fbd-4068-9df1-5cfd483e34d7"
      },
      "source": [
        "!pip install geopandas\n",
        "!pip install alphashape\n",
        "!pip install skan\n",
        "#!pip install https://github.com/scikit-image/scikit-image/archive/refs/tags/v0.18.3.tar.gz"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Requirement already satisfied: geopandas in /usr/local/lib/python3.7/dist-packages (0.10.2)\n",
            "Requirement already satisfied: fiona>=1.8 in /usr/local/lib/python3.7/dist-packages (from geopandas) (1.8.20)\n",
            "Requirement already satisfied: pyproj>=2.2.0 in /usr/local/lib/python3.7/dist-packages (from geopandas) (3.2.1)\n",
            "Requirement already satisfied: pandas>=0.25.0 in /usr/local/lib/python3.7/dist-packages (from geopandas) (1.1.5)\n",
            "Requirement already satisfied: shapely>=1.6 in /usr/local/lib/python3.7/dist-packages (from geopandas) (1.7.1)\n",
            "Requirement already satisfied: certifi in /usr/local/lib/python3.7/dist-packages (from fiona>=1.8->geopandas) (2021.5.30)\n",
            "Requirement already satisfied: six>=1.7 in /usr/local/lib/python3.7/dist-packages (from fiona>=1.8->geopandas) (1.15.0)\n",
            "Requirement already satisfied: munch in /usr/local/lib/python3.7/dist-packages (from fiona>=1.8->geopandas) (2.5.0)\n",
            "Requirement already satisfied: click-plugins>=1.0 in /usr/local/lib/python3.7/dist-packages (from fiona>=1.8->geopandas) (1.1.1)\n",
            "Requirement already satisfied: attrs>=17 in /usr/local/lib/python3.7/dist-packages (from fiona>=1.8->geopandas) (21.2.0)\n",
            "Requirement already satisfied: setuptools in /usr/local/lib/python3.7/dist-packages (from fiona>=1.8->geopandas) (57.4.0)\n",
            "Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.7/dist-packages (from fiona>=1.8->geopandas) (7.1.2)\n",
            "Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.7/dist-packages (from fiona>=1.8->geopandas) (0.7.2)\n",
            "Requirement already satisfied: numpy>=1.15.4 in /usr/local/lib/python3.7/dist-packages (from pandas>=0.25.0->geopandas) (1.19.5)\n",
            "Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.7/dist-packages (from pandas>=0.25.0->geopandas) (2.8.2)\n",
            "Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.7/dist-packages (from pandas>=0.25.0->geopandas) (2018.9)\n",
            "Requirement already satisfied: alphashape in /usr/local/lib/python3.7/dist-packages (1.3.1)\n",
            "Requirement already satisfied: click-log>=0.3.2 in /usr/local/lib/python3.7/dist-packages (from alphashape) (0.3.2)\n",
            "Requirement already satisfied: shapely>=1.4.0 in /usr/local/lib/python3.7/dist-packages (from alphashape) (1.7.1)\n",
            "Requirement already satisfied: networkx>=2.5 in /usr/local/lib/python3.7/dist-packages (from alphashape) (2.6.3)\n",
            "Requirement already satisfied: Click>=6.0 in /usr/local/lib/python3.7/dist-packages (from alphashape) (7.1.2)\n",
            "Requirement already satisfied: trimesh>=3.9.8 in /usr/local/lib/python3.7/dist-packages (from alphashape) (3.9.35)\n",
            "Requirement already satisfied: scipy>=1.0.0 in /usr/local/lib/python3.7/dist-packages (from alphashape) (1.4.1)\n",
            "Requirement already satisfied: rtree>=0.9.7 in /usr/local/lib/python3.7/dist-packages (from alphashape) (0.9.7)\n",
            "Requirement already satisfied: numpy>=1.8.0 in /usr/local/lib/python3.7/dist-packages (from alphashape) (1.19.5)\n",
            "Requirement already satisfied: setuptools in /usr/local/lib/python3.7/dist-packages (from trimesh>=3.9.8->alphashape) (57.4.0)\n",
            "Requirement already satisfied: skan in /usr/local/lib/python3.7/dist-packages (0.9)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fDhIRH8HOaoc"
      },
      "source": [
        "import numpy as np\n",
        "from shapely.geometry import Polygon, MultiPoint, LineString, MultiLineString, Point\n",
        "import alphashape\n",
        "import matplotlib.pyplot as plt\n",
        "import geopandas as gpd\n",
        "import pandas as pd\n",
        "from skimage.morphology import medial_axis, skeletonize\n",
        "import scipy.sparse as sparse\n",
        "import skan\n",
        "import gdal\n",
        "import ogr"
      ],
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1mwasI_olm_G"
      },
      "source": [
        "We write two classes that we use to generate synthetic data: `GenericPolygon` is a parent class that takes a polygon geometry. In the following, we work with the child class `RandomPolygon` that inherits from `MyPolygon` and generates a polygon (concave or convex) consisting of random coordinates."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dDoKekws1fcr"
      },
      "source": [
        "class GenericPolygon:\n",
        "  \"\"\"\n",
        "  Class that represents a generic polygon.\n",
        "  Args:\n",
        "    - geometry: Polygon geometry (GeoSeries)\n",
        "  \"\"\"\n",
        "  def __init__(self, geometry):\n",
        "    self.geometry = geometry\n",
        "    self.bbox = self.geometry.bounds\n",
        "\n",
        "class RandomPolygon(GenericPolygon):\n",
        "  def __init__(self, extent, crs, polytype):\n",
        "    \"\"\"\n",
        "    Class that represents a randomly generated polygon.\n",
        "    Args:\n",
        "      - extent: range for randomly generated coordinates as 4-uple of floats: (xmin, xmax, ymin, ymax), format must be aligned with crs\n",
        "      - crs: coordinate reference systems as string, e.g. \"epsg:4326\" for WGS84\n",
        "      - polytype: two options \"concave\" or \"convex\"\n",
        "    \"\"\"\n",
        "    self.polytype = polytype # concave or convex\n",
        "    geom = self.create_geom(extent, crs)\n",
        "    GenericPolygon.__init__(self, geom)\n",
        "\n",
        "  def generate_coords(self, extent):\n",
        "    lon = np.random.uniform(low=extent[0], high=extent[1], size=10) \n",
        "    lat = np.random.uniform(low=extent[2], high=extent[3], size=10)\n",
        "    return list(zip(lon,lat))\n",
        "\n",
        "  def create_geom(self, extent, crs):\n",
        "    self.coords = self.generate_coords(extent)\n",
        "    if self.polytype == \"convex\":\n",
        "      poly = Polygon(self.coords) # to polygon\n",
        "      geom = gpd.GeoSeries(poly.convex_hull) # take convex hull\n",
        "    elif self.polytype == \"concave\":\n",
        "      geom = gpd.GeoSeries(alphashape.alphashape(self.coords))\n",
        "    geom.set_crs(crs)\n",
        "    return geom"
      ],
      "execution_count": 12,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "W_YU75DzlQwf",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 295
        },
        "outputId": "91a8cc16-89a5-4a17-fa20-616a0c6b1046"
      },
      "source": [
        "random_poly = RandomPolygon((600000.0, 620000.0, 5000000.0, 5020000.0), \"epsg:32630\", \"concave\")\n",
        "random_poly.geometry.plot()"
      ],
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.axes._subplots.AxesSubplot at 0x7f3edd661990>"
            ]
          },
          "metadata": {},
          "execution_count": 13
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ8AAAEDCAYAAAAryNZsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAddUlEQVR4nO2deZRU1b3vP7+uqp6gaaBp5qGhG5wgInbAAWVwxQEnjENwBB8EZ9+KL3mJ9yZxXc27K8l76yV6vYmXaIwap8iLCTcqPpwi+gQEnGdEEQjKDCJj07/3x9mFh6ar61R11dlVp/ZnrVp9ap+z9/7t7m/vfc4+e/9+oqo4HDYos22Ao3Rx4nNYw4nPYQ0nPoc1nPgc1nDic1gjMuITkd+LyHoReTvg9ReJyLsi8o6IPJRv+xyHIlGZ5xORk4EdwP2qOjLNtcOBPwGTVXWLiPRW1fVh2On4msj0fKr6IrDZnyYijSIyX0SWichCETncnPou8O+qusXkdcKzQGTEl4I5wA2qeizwfeA3Jn0EMEJEXhaRRSJyujULS5i4bQPyhYh0BU4AHhORZHKF+RkHhgMTgYHAiyIySlW3hm1nKRNZ8eH16ltVdXQ759YAi1V1H/CJiHyIJ8ZXwzSw1InssKuq2/GEdSGAeBxtTv8Fr9dDRHrhDcMrbdhZykRGfCLyMPAKcJiIrBGRmcClwEwReQN4BzjXXP40sElE3gWeB36gqpts2F3KRGaqxVF8RKbncxQfkXjg6NWrlzY0NNg2w9EOy5Yt26iq9e2di4T4GhoaWLp0qW0zHO0gIqtSnXPDrsMaTnwOazjxOazhxOewhhOfwxpOfA5rOPE5rOHE57BGJCaZc8mmHXu48ZHXqErEqCqPU5UooyoRo7I8RlUiRrX5WZmIUVXufa9MxMz1bX4mYsRj7v87FU58baitSrBo5Wb2t+ZmwUUiJlSmEO2B70nB+kTbnpArffn8/xCJIhV4IPGJyKfAl8B+oEVVm9ucF+B2YAqwE5ihqsvNufnAccBLqnqWL89CoMZ87Q0sUdWpIjIR+CvwiTn3Z1W9NavWZUE8Vkb/7pWs3rwrJ+Xt26/s29/Cl7tbclJeeyQFnhRreayMikQZ5bEyyuNllMdjVMS94wqTlvxeHi+jPBZrc713viLNea+82IG0WJmkN9ZHJj3fJFXdmOLcGXgrgYcD44Dfmp8A/xOoBq7yZ1DVk5LHIvJ/8ASXZKFfqGEzqEd1zsQXBmEIPB29ayp49r9NoKYyEThPrvrrc/G2LKqqLgK6i0g/AFV9Fq/XbBcR6QZMxltdXBAM7llt24Si48dnHZmR8CC4+BT4v2YL4ux2zg8AVvu+rzFpQZgKPGuWvSc5XkTeEJGnROSogOXkjEFOfBlxYlMdZ3+jX8b5gg6741V1rYj0BhaIyPtmn2wuuBi42/d9OTBEVXeIyBS8HnF420zmn2A2wODBg3NkisfAHlU5LS/KJGLCreeOxLdDMDCBej5VXWt+rgceB8a2uWQtMMj3faBJ6xCzeWcs8ISvru2qusMcPwkkzHVtbZqjqs2q2lxf3+5axaxxPV9wrjq5kcb6rlnlTSs+EekiIjXJY+BUoK0/lHnAFWaH2HHANlVdF6D+C4C/qepuX319zdMzIjLW2Bjq5h53zxeMgT2quG5SU9b5gwy7fYDHjR7iwEOqOl9ErgZQ1buAJ/GmWVbgTbVcmcxsplQOB7qKyBpgpqo+bU5PA37epr4LgGtEpAXYBUzTkHc51XUppyoRY9e+/WFWW3T8yzlHUVUeyzp/WvGp6krg6HbS7/IdK3BdivwntZduzk1sJ+1O4M50duUTEWFQzyo+/GKHTTMKmlOP7MMpR/TpVBnFOTUeAoN6uKE3FVWJGD89+8hOl+PElwL30JGaG08ZzsAc/HM68aXAia99mnp3Zeb4oTkpy4kvBYPcXF+73HbuSMrjuZGNE18KXM93KOcdM4DjG+tyVp4TXwqc+A6mpjLOzVMOT39hBjjxpaBrRZyeXcptm1Ew/OC0w+hdU5nTMp34OsDd93mMHNCNS8cNyXm5TnwdMNANvYjAz6aOynihaBCc+DrAveOFS8YOZvSg7nkp24mvA0r9LUddl3L++2m5fcjw48TXAYN6lvY93z9NOYLa6sxWJ2eCE18HlHLPN7ahJ98eE3QxenY48XVA/+5V5OE+u+CJlwm3Tc1udXImOPF1QHm8jH61pTf0zhw/lMP61qS/sJM48aWh1PZz9Kut5MZTDtkykxec+NJQaq/Zbjn7SLpUhOPIwokvDaU01zfxsHpOO6pvaPU58aWhVKZbyuNl/Ms5R+X9IcNPIPGJyKci8paIvC4ih8QcMLvW7hCRFSLypoiM8Z2bLyJbReRvbfL8QUQ+MWW+LiKj05Vlg1KZbrluYhND6rqEWqc1Xy2GH6jq3AzKCp1SuOdrqKvmqgnDQq/Xuq+WTMqyQX3XCipytHK3ULn13JFUJrLfApkttn21/A8ztP5KRJKBmAOVJSKzRWSpiCzdsGFDsFZkQVmZRHq65cxR/Th5RG49PgQlqPjGq+oYvCHxOhE5OQd134y3mfybQE/gh5lkzqe7jLZEdejtUh7jJ2d1fgtktljz1aKq68zQuge411dmVn5f8klUHzq+960R9K3N7erkTLDmqyV5H2f8skz1lZmt35e8EcW5vsP71jDjhAarNtj01fKgiNQDArwOXG2ypCzLFlGc6/vZ1JHWnZVb89WiqpNTpKcsyxa52J1fSFzUPJDmhp62zXBvOIIQpQeO7tUJfnTGEbbNAJz4AlFblaC2Kn8resPkh6cfXjBbQp34AhKF+77Rg7rzneZB6S8MCSe+gBT7dEuZeA8ZZQW0NNuJLyDFft93xfENjBxQa9uMg3DiC0gxi6++poKbTh1h24xDcOILSDG7zvjxmUfQLcMALWHgxBeQYu35Tmyq45yj+9s2o12c+AIyoHvx9XydCdASBk58AalMxOjbzd5L+GzoTICWMHDiy4BimuvrbICWMHDiy4BimuvrbICWMHDiy4Bi8df3rRwEaAkDJ74MKIZ1fZWJMm7JQYCWMHDiy4BimOvLVYCWMHDiy4BCn+trrO/CrPHhb4HMFie+DOjTrZJErDDnzABum5q7AC1hUDyWFgCxMinYIW3q6P6c0HhITOyCxokvQwpxD29NRZx/OrMwVidngk1fLQ+KyAci8raI/F5EEiZ9oohs8/lw+WlnG5lLCvG+7/t5CNASBjZ9tTwIXGaOHwJmmXwAC1X1rAxsC41Cm2geOaAblx2X+wAtYWDNV4uqPmmuV2AJ3ubwgqeQ5vryGaAlDGz7asEMt5cD833Jx4vIGyLylIgclSJfKL5a2lJI73fzGaAlDIIOu+NVda2I9AYWiMj7qvpijmz4DfCiqi4035cDQ1R1h4hMAf6CN5wfhKrOAeYANDc3a45sSUuhDLv5DtASBtZ8tQCIyC1APXCTr67tqrrDHD8JJESkYOYQulcn6BqSz+KOuDnPAVrCwKavllnAacDFqtrqS+9r/LcgImONjZsyaFNeERHrT7xjG3pyfp4DtISBTV8tdwGrgFdM2X9W1VuBC4BrRKQF2AVMMw8lBcOgHlW8t267lbpjIQVoCQObvlrarVtV7wTuTGeXTWz2fGEFaAkD94YjC2ytbunbrZL/GlKAljBw4suCwXV2er4wA7SEgRNfFtiYbpkwop7TR4YXoCUMnPiyIOyVLTYCtISBE18WVJXH6NW1Iv2FOeLaiY009Ao3QEsYOPFlyeCQXrM11FVz9YTGUOoKGye+LAlruuXmKUdYCdASBk58WRLWQ8fqzTtDqccGTnxZEtbqlgcXf0Zra0G94MkZTnxZEtaw+8nGr3j541RreIsbJ74sCXOu74FXVoVWV5g48WVJv9rK0FYQP/PeF6zduiuUusLEiS9L4rEy+ncPZ9NOq8LDiz8Lpa4wceLrBGHu53jk1c/Y29Ka/sIiwomvE4R537dxx17mv/N5aPWFgRNfJwh7Xd8fI/bg4cTXCcL2XrDk0828/7mdFdT5wImvE9jYw/vHRdHp/Zz4OoGN5fSPL1/Ll7v3hV5vPrDpq2WoiCw2eR4VkXKTXmG+rzDnGzrXxPxR16WcqpBf+n+1dz+Pv5Z2V2pRkEnPN0lVR6tqczvn/L5aZvO1zxXwfLVc3k6eXwC/UtUmYAsw06TPBLaY9F+Z6woSbxtl+Ps5HnhlFQW2oS8rrPhqMftyJwNzTdJ9wFRfWfeZ47nAKVLAS3ht3Pd9tH4Hiz/ZHHq9ucaWr5Y6YKuqtrRz/YGyzPlt5vqDsOWrpS22nEU+EIEHj6DiG6+qY/CG1+tE5OQ82hQIVZ2jqs2q2lxfX2/NDlt7eJ9++3PWb99tpe5cYctXyya8oTnezvUHyjLnaykgdxltsbWHt6VVeeTV1ekvLGCs+GoxHg6ex3ONATAd+KuvrOnm+ALguUJzl+HH1h5egIcWf0bL/uJ93xuk5+sDvCQib+A5cXwi6asl6a8Fz1fLSjxfLb8Drk1mNr5aHsN7cFgjIqeZUz8EbhKRFXj3dPeY9HuAOpN+E/CjTrUwz9h0mfb59t08894X1urvLDZ9tazk0OEbVd0NXJjOrkKhujxGbVWCbbvsTPw+sGgVp4/sZ6XuzuLecHSCbbv2ce2Dy60JD+DlFZtYsX6Htfo7gxNflryxeitn/dtCnnrb/jKnBxcX57SLE1+GqCr3vPQJF9z1/1i9uTCWts9dtoade1vSX1hgOPFlwNade5n9wDJu+9u77NtfOA/gX+5uYd7r/7BtRsY48QVk+WdbOPOOl1jwbmE+Xd5fhO97nfjS0NqqzHnxYy6665WC3kH27rrtLP9sq20zMiI6ngbzwOav9vL9x97guffX2zYlEH9ctIpjh/SwbUZgXM+Xglc/3cyZdywsGuEBPPHmOjbt2GPbjMA48bWhtVX5zQsrmDZnEeu2FdeL+737W/nT0jW2zQiME5+PjTv2MOMPr/LL+R+wv0id8zy4eFXR2O7EZ1i0chNTbl/Iix/aWxuYC9Zs2cXfPyyOW4WSF9/+VuWOZz/ikt8tYv2XxXO/1BHF4liopJ9213+5m+89+jovryjY5YJZ8cKHG/hs006ry72CULI938srNjLl9pciJzwA1eJ431ty4tvfqvzvBR9y2T2L2VhE0xKZ8ujS1ezet9+2GR1SUuL7YvtuLr17EXc8+xFF9iYqY7bu3McTb3YY+NM6JSO+v3+4gSm3L2TRyuLfchiUQt/hFnnxtexv5Zfz32f675ew6au9ts0JlddXb+WtNdtsm5GSSItv3bZdXPy7RfzmhY9tm2KNQnYsFIavluki8pH5TDdpNaas5GejiPzanJshIht852Zl07Dn31/PlNsX8uqnW7LJHhn++sZatu0sTMdCmczzTVLVVD75/b5axuH5ahknIj2BW4BmPK8Hy0RknqpuAUYnM4vIMuDPvvIeVdXrM7DtILbt3MeVf3g12+yRYve+Vh5btppZJw2zbcoh5NtXy2nAAlXdbAS3ADjdn1FERgC9gYU5soVuVXFqqxK5Kq7oKdRAMvn21RLEh8s0vJ7O/9s53wzfc0VkEO3Qka8WEaGpd9dADSsFCjWQTCH4apkGPOz7/p9Ag6p+A6+nvK+9TOl8tTTVO/H5KcT3vfn21dKhDxcRORqIq+oyX12bVDX56uFu4NhALWmD6/kO5pn3vuAfBbYNIN++Wp4GThWRHiLSw+R92pfvYg7u9Uj69TOcA7yXYZsAJ762tCo8vKSwAskEedrtAzxu/DPGgYeSvlrggNuMJ4EpeL5adgJXmnObReQ2IPnoeauq+l8xXGTy+blRRM4BWoDNwIws2kWjG3YP4eElq7lh8nDK44UxvSvFtt2uPZqbm3Xp0oOnH/e3Kkf+dD57Iha1p7PccfExnHN0/9DqE5FlKVwpR/cNR6xMGOZ6v0MopEAykRUfuPu+9iikQDLRFp/r+dqlUN73Rlp8jb272DahICmUQDKRFl//7nb8JRc6X+3dz18KIJBMpMUXDykSeDHywCL7joUiLb6ywo0dY50Pv9jBEsuBZCItPqe9jrnf8oNHpMXner6OsR1IxomvhLEdSCbS4otFunW5wWYgmUj/eQo4WGXB4AWSseNYKNLic8NuMGy98Yi4+GxbUBy8tGIjH28IP5BMxMXn1BcUG71fpMXntBccG4FkIi0+1/MFx0YgGSc+xwHCDiQTcfHZtqC4eHfddl5bHV4gGSu+Wkz6CyLygc8nS2+TXiEij5qyFotIQ7aNc/N8mRPmMvtMer5Jqjo6xWYQv6+W2Xi+WvD5ahmHt9f3FrOFMsmlpszRZk8wwExgi6o2Ab8CfpFRi3y4ni9z/vbmOjaH5ErOuq+WFGUlvRTMBU6RLLuwmFNfxuzd38qjIb3vte2r5V4z5P7EJ7ADeVS1BdgG1LWtsCNfLb5rAjXOcTBhBZKx6avlUlUdBZxkPpdnkjmdrxZww262hBVIxpqvFl+ZXwIP+co8kEdE4kAtkFW8AjfVkj1hOBay4qtFROIi0suUmQDO8pU5D0g+FV8APKdZTj458WVPMpBMPrHiq8WI+GkjvBjwDPA7c809wAMisgLPV8u0bBvntJc9qvDgklXcfMYReasjsr5aAPa2tDLix09ZsCga9KhO8MrNp1CZiGVdRkn6agH3wNFZtuzcx5Nv5S+QTMTF59TXWe7P44NHpMXntNd58hlIJuLiEzf05oB8LTSNtPjADb25IF+BZJz4HGnZva+VucvX5LzcyIvPaa/zfOvIPhw/7JDX650m8mHuXc+XPacf1ZcbTmniqP61eSm/BMRn24LiY8qovtwweThH9OuW13pKQHxOfUEQgSmj+nHj5OEc1rcmlDojLz6nvY4RgbO+0Z8bJjcxok84oksSefGVuXG3XcoEzj7aE11T73BFlyT64nNd30GUCZw7egDXT26yHqWpBMRn24LCIFYmnDu6P9dPaiqY4DglIL7SVl+sTDjvmAFcP6mJhl6FFRrCiS+ixMuE88cM5NpJjQypKyzRJSkB8dm2IFziZcKFzQO5dmITg3pW2zanQyIvvlLZPpmICRc2D+LaiY0M7FHYoktixV2GiFSLyBMi8r6IvCMiP/ddP0NENvjcaMzqVAMj/va6PFbGZccN5oUfTOJfzxtVNMKDzHq+Saq6McU5v7uMcXjuMsb53GU04208XyYi84A9wP9S1edFpBx4VkTOUNXkhotHVfX6LNpzCFG95yuPlTFt7CCuntBYtGG+cjXsHnCXASwSkaS7jIkYdxkAIrIAOF1VHwaeB1DVvSKyHG9Pb86JmvjK42VcMnYwV09opG9tpW1zOkVQ8SXdZSjwH6o6p835bN1lICLdgbOB233J5xuvCB8C31PVrJ2HREV7FfEyLhnnia5Pt+IWXZKg4huvqmuNG7MFIvK+qr7Y2cqNR4KHgTtUdaVJ/k/gYVXdIyJX4TkNmtxO3tl4HrEYPHhwyjqKveerTJRx2bghzD55GL0jIrokgcTnd5chIkl3GX7xdeQuY2Kb9Bd83+cAH6nqr311+V1j3A38MoVNc0x+mpubU24+LtaplqpEjMuPH8J3TxpGfU2FbXPyQlrxGe8CZar6pc9dxq1tLpsHXC8ij+A9cGxT1XUi8jTwrz6ffKcCN5tyf4bnh+Wgp1kR6WdcbQCcA7yXXdM8iq3nqy7/WnS9ukZTdElsucsYCPwz8D6w3JR9p6reDdwoIucALXjuMmZ0poHFIr4u5TGuOKGBWeOHUhdx0SWJtLsMgLP+bSFvr90eskXB6VoRZ/oJQ5g5fhg9u5TbNifndOQuI/JvOAq156upiDPjxAZmjh9K9+roiS4IkRdfob1eq6mMc+WJQ5l54lBqqxO2zbFK5MVXKE+73Srj/JfxQ7nyxKHUVpW26JKUgPjsqq+2KsHM8UOZcWID3Sqd6PyUgPjs1Nu9OsGs8UOZfkIDNU507RJ58YV9z9ejOsGsk4Yx/YQGulZE/tfbKSL/2wmr5+vZpZzvnjSMK44fQhcnukBE/reU73u+Xl3LmX3yMC47bgjV5ZH/deaUyP+28hWFqFfXCq6eMIxLxw2hqjx7n8WlTOTFl+t7vvqaCq6e0MglYwc70XWSyIsvVx1f75oKrpnYyMVjB3fKO7vja0pAfJ1TX99ulVwzsZHvfHOQE12OKQHxZZevX20l105s5KJvDqIi7kSXDyIvvkzv+QZ0r+KaiY1c2DzQiS7PRF58QXu+gT2quG5SE+ePGUh5POL7LQuEEhBfx+ob1LOK6yc18e0xA0nEnOjCpGTFN6SumusmNXHeMQOc6CwRefG11V5DXTXXTx7O1NH9iTvRWSXy4kv2fMN6deH6yU2cc7QTXaFgxVeLST/WlLnC5BWT3lNEFpjrF/h2vmXF8N5d+fV3RrPgpgl8e8xAJ7xCQlXTfoBPgV4dnJ8CPAUIcByw2KT3BFaanz3McQ9zbom5VkzeM0z6L4EfmeMfAb9IZ9+xxx6rjsIEWKop/m656gYO+GpR1UVA0lfLaRhfLaq6BVgAnG7OdVPVRcbA+4GpvrLuM8f3+dIdESOo+JK+WpYZNxVtydRXywBz3DYdoI9+vWn8c7x9w4cgIrNFZKmILN2wYUPAZjgKiaDiG6+qY/BcoV1nnPjkHdMrtruxWFXnqGqzqjbX19eHYY4jxwQSn/p8tQBJXy1+OvLVkip9YDvpAF+YYRnzc30QGx3FR1rxiUgXEalJHuP5W3m7zWXzgCvMU+9xGF8twNPAqSLSwzy1ngo8bc5tF5HjzFPuFcBffWUln4qn+9IdEcOKrxZzfC3wB6AK72k36ZX058CfRGQmsAq4qDMNdBQukffV4rBLR75a3IyrwxpOfA5rRGLYFZENePeHYdELSOWZv9jry3VdQ1S13bmwSIgvbERkaar7mGKvL8y63LDrsIYTn8MaTnzZ0TYOSZTqC60ud8/nsIbr+RzWcOJzWKPkxGeCEs414VbfE5HjReRCE3q1VUSafdd+y6xhfMv8nOw794KIfOALzdrbpFeIyKNme8BSEXnSV9dvRWSdiOxpp64GEdnlK+8u37mgWw6G+Nq2WUTWiMgqEfm4nfou9dX1ujk/OoO2LRaRBl95N5v0D0TktEB/jFRLnKP6wVsdPcsclwPdgSOAw/BCczX7rj0G6G+ORwJrfecOutaXfi1wlzn+O7DEHB8NvAV8Ay+W3C5grC9fA/B2CpuDbjl4Cy+i0xRgvmnbNOCNVPaavKOAjzNs2zS80LQAR5o6KoChwMdALN3foqR6PhGpBU4G7gEv3KqqblXV91T1g7bXq+prqvoP8/UdoEpE0oUHOhe4z9Q1GBhmeqopeCuC3lTV5/DEd1QAm4NuOfgz3j/QPcl007ZHgEq8f7RUXAw8ks6WNvXNBU4xbTsXeERV96jqJ3irm9qu+TyEkhIf3n/lBuBeEXlNRO42axSDcD6wXFX3+NLuNcPST5JDIV9vHUjWVY7XK1zGwQtj9wBtXzsNNXb9XURO8pUXZMtBNd7f8168XukSX9vW4PVKqfgOXvRPPx21DVVtAbYBdQQIbdsepSa+ODAG+K2qHgN8hTdcdYiIHAX8ArjKl3ypqo4CTjKfy1PUtR1vmG0BzuqgmnXAYGPXTcBDItItSKN89cXworwvBHYTrG3jgJ2q6l8gnK5tOaHUxLcGWKOqi833uXgCSYl4QQofB65Q1Y+T6fr11oIvgYf4ephJbh1YYz6VwCZgEd69ZZIKvJ4xWd4eNeFeVXUZ3n3TCAJuOQD2Ai2mbWvxbhPG+PL4e2w/02jT6wVoWzJWcq1pW6rtEh1SUuJT1c+B1SJymEk6BXg31fXiRUF/Au+m/mVfelxEepnjBF6Pluw55gHTTV27gFfNvVoL0NM8MQ7FW8H9jq/MehGJmeNhwHBgpQbfcjAF+Ny0bR5wKfCumG0NeOJs274yvJXij/jS0rbNHF8APGfaNg+Y5mvbcLyHpI6x/fQZ9gcYDSwF3gT+greZ/Ty8XmoP8AXePhOAH+MNza/7Pr2BLsAyU8Y7wO2Ypzu8nu4xvJvut801ybpuw9sOus+IwV/X+aas14HlwNk+m5tNWR8Dd/L1m6k64FngI+AZvIepZNtWAp/gLTX7om3bTP6JwKI2v5+gbVsCDPPl+2dj3weYp/F0H/d6zWGNkhp2HYWFE5/DGk58Dms48Tms4cTnsIYTn8MaTnwOa/x/qeeaEirFrWwAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "12zGP1mDnkrk"
      },
      "source": [
        "We save the random poly to a file and rasterise in the next step."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YfDmGHighz1d"
      },
      "source": [
        "random_poly.geometry.to_file(driver='ESRI Shapefile', filename=\"polygon1.shp\")"
      ],
      "execution_count": 14,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sDpXHcP6CTtF",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 36
        },
        "outputId": "1faf8767-f912-4271-c027-d27a913f5811"
      },
      "source": [
        "def rasterise(shp_path, output_path):\n",
        "  \"\"\"\n",
        "  Save raster representation of ESRI Shapefile vector layer to GeoTIFF. \n",
        "  Args:\n",
        "    - shp_path: Path to input layer (ESRI Shapefile) as string\n",
        "    - output_path: Path to output layer as string\n",
        "  \"\"\"\n",
        "  shp = ogr.Open(shp_path)\n",
        "\n",
        "  shp_layer = shp.GetLayer(0)\n",
        "  xmin, xmax, ymin, ymax = shp_layer.GetExtent()\n",
        "  pixel_size = 100\n",
        "\n",
        "  ds = gdal.Rasterize(output_path, shp_path, xRes=pixel_size, yRes=pixel_size, \n",
        "                    burnValues=1, outputBounds=[xmin, ymin, xmax, ymax], \n",
        "                    outputType=gdal.GDT_Byte)\n",
        "  ds = None\n",
        "  return output_path\n",
        "\n",
        "rasterise(\"polygon1.shp\", \"raster1.geotiff\")"
      ],
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            },
            "text/plain": [
              "'raster1.geotiff'"
            ]
          },
          "metadata": {},
          "execution_count": 15
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "O2fYoJSlnQ5J"
      },
      "source": [
        "Now we open the newly created raster file."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mVIZB96Ijsmm",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 324
        },
        "outputId": "7391d720-7e66-492f-8742-334a1f310bb2"
      },
      "source": [
        "raster = gdal.Open(\"raster1.geotiff\")\n",
        "gt = raster.GetGeoTransform()\n",
        "print(gt)\n",
        "\n",
        "raster_arr = raster.ReadAsArray()\n",
        "print(np.unique(raster_arr))\n",
        "plt.imshow(raster_arr)"
      ],
      "execution_count": 16,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "(610439.8540557557, 100.0, 0.0, 5017977.990950754, 0.0, -100.0)\n",
            "[0 1]\n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.image.AxesImage at 0x7f3edd567d50>"
            ]
          },
          "metadata": {},
          "execution_count": 16
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI8AAAD8CAYAAABO1ff8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAXCElEQVR4nO2deXxU5bnHv89kIwmGVZBVYgiKoFZlERcEFUFq1dp7RaWKiFJEXCp17b3X++lyq6UfW1uVVioVlaoUd0EBFcXKIqCgLAIBIoQtrJEkkGXmuX/MCUZMyMw5M3POmXm/nw+fzJzt/QV+vOc973mf5xFVxWCwQ8BtAQb/YsxjsI0xj8E2xjwG2xjzGGxjzGOwTdzMIyLDRGSdiBSJyAPxasfgHhKPeR4RSQPWA0OAEmApcJ2qrol5YwbXiFfP0w8oUtVNqloNvARcGae2DC6RHqfrdgK21vteAvRv7OBMydJm5DpqsF3vw7QMhKI+r1yV7V86azuZOcj+Pap6fEP74mWeJhGRscBYgGbk0F8udnS9214v4qrc8qjPq9IaJh8odNS2XSa/ehnd/nuRK21Hyns68+vG9sXLPNuALvW+d7a2HUFVnwaeBsiT1q69YMuSDO5uVexK22VXfMS0LgNcaTtiRs1sdFe8zLMUKBSRfMKmuRa4Pk5t+ZaHj1/Dw5d68xlic005t2+6hka7HeJkHlWtFZEJwBwgDZiqqqvj0ZYhPsz45kyCg7cf85i4jXlUdTYwO17XN7iPmWE2fI/u80ezYMhJTR7n2tOWwZt0f3EcJ808RO3OXU0ea3oew3do9ynIopURHZs05plR2peNNdHP8xi+5ZmyE8gsD0Z8fNKYZ+95+/nx57e6LcPXvDr4DLJmLY34+KQxjyHxGPMYWF5VzSXX30xwz96ozjNPWynOpH0FPDt9KJ0/XEi074hMz5PivF5yBp1/t9DWucY8KUxpsIKyymzb5xvzpDAX//leOl1t/5WjMU8q43AhjDFPinLqU+Pp/O4+R9dIKvPoJ60YveUCt2V4miqtoecnN5A/bQuhVV85ulZSmafjpIWseP40t2V4mrJQNV1HrKF2a4njayWVeQxNUx3DUCvb5hGRLiIyX0TWiMhqEbnL2t5aROaJyAbrZ6uYqTU44v5dP+DWXpdBKPKXn8fCSc9TC0xU1VOBc4DbReRU4AHgfVUtBN63vhtcpucnN/D5+DMIHTwYs2vaNo+q7lDVz6zPB4G1hOO1rgSmWYdNA65yKtLgnKodORGv04mUmLzbEpFuwJnAEqC9qu6wdu0E2jdyznfitgzx49G9heRuSYv5dR2bR0SaA68Ad6vqNyJyZJ+qqog0OELzStxWKvDOA4PoOMve+6tj4ehpS0QyCBtnuqq+am3eJSIdrP0dgFJnEg1excnTlgDPAGtV9bF6u94ERlmfRwFv2JcXPR2mr2bQGLOiECCoIYZedQPZ8+MTMuek5zkPuAG4SERWWH+GA48AQ0RkA3CJ9T1hBA+UkVN8IJFNepLlVdX0+d0E5PN1hCor49KG7TGPqv4bkEZ2O8taYHBMcU1b2j0R/QKvaDAzzElISW0588t6xr0dY54kZOjSn7Ghb1Xc2zHmMdjGLIBPMgpeHkfh9PK4jnXqMD1PEnHS3DEUzDyMLluVkPZMz5NEnPJ/Bwiu35iw9pKz5wmGKA1WuK0ioZQGKyDB5a+S0jzB9Ru5Mf9CNqdI4oONNeXccOJAghs2JbTdpDQPgNbWui0hscRogVc0JK15UoW/HujET/5wnyttmwGzj5m0r4C/zRlCwV9iv9wiEkzP42OmvHkpBRMXu9a+MY/BNua25VMGjh9LwdyVRF9tI3aYnsenZJbVxm2dTqQ4No+IpInI5yLytvU9X0SWWEXaXhaRTOcyDXUENUT+7FvILHF/wVssep67CIfd1PEo8EdV7Q7sB8bEoA1bPLJrCFtqk2eicH+wkv8u/QGnTFiV8AnBhnC6AL4z8EPg79Z3AS4C6kqluBq3VdzvEA9vv8yt5mPOywcLWX5mgNDhw25LAZz3PH8C7oMj47Y2wAFVrZveLSEcCGhIQpxET1wOlKrqcpvnjxWRZSKyrIb4r3rzO2ctG8Eb13orfYyTR/XzgCusiIlmQB7wONBSRNKt3ud7RdrqMEF/kdN9/mg6zsgg9MWnbkv5Dk5i1R9U1c6q2o1wMbYPVHUkMB/4D+uwhMdtJSNt5jQj+w1vGQfiM89zP3CPiBQRHgM9E4c2UoYVVVWkVXuzY47JDLOqfgh8aH3eRLg0tiEG/HLgTzhuq3vvr46FmWE22MaYx6NsrCnnnHvHEdq9x20pjZL05lk8+zTu3N7XbRlR8Up5HkPenkiL6Ys9MyHYEElvnq6/Wsg77/VxW0ZUPL11IIW3L3FbRpMkvXn8Ro0GqQ7FPotXPDDm8RinP30HWcN3NH2gBzDm8RhSC1pT7baMiDDm8RD5b4yl04JDbsuIGLMM1QMENcSITZfS80/7CK4rcltOxBjzeIByraJ88AG01rtzOg1hblsG26SEeQofWcOAiePcltEgj+07ievPH+HL8OiUME/wQBnZu2vcltEgB4PNqC3e4rYMW6SEebzKo3sLeX7uhW7LsI0xj0vMrczgmbcuoeDeRW5LsY3T6ImWIjJTRL4SkbUiMsDU24qMX91/M/kP+dc44LzneRx4V1VPAc4gHL9l6m2lCE6iJ1oAA7GWmapqtaoewNTbapILJvyMvA83uC3DMU56nnxgN/APK9z47yKSS4T1tlKRjTXl9Jh2G83f/ZLgXmdlqb2AE/OkA2cBk1X1TKCCo25Rqqo0Uvo9FeO2Nta0Iv/BRa4nKIgVTsxTApSoat2qpZmEzRRRvS1VfVpV+6hqnwyyHMiIjEBQXU1wWRmqprjmeNfajwdO4rZ2AltF5GRr08XAGlyut9UYafM/Y/ypQ6lSdyYLL1zxU17p2c6VtuOF0xejdwDTrTQqm4DRhA05Q0TGAF8D1zhsw+BRHJlHVVcADS0QNvW26lHwwWi6TE++BQxmhjnOnPfF1XR5MZ2sd5a6LSXmJN9/B4+R94tMQquSzzhgep74k8R/w0n8q3mDabOeoXTCuW7LiAvGPHGmXVouoSQdHBjzJIBWw7ez99YBbsuIOallnmCQO7cNZE+Ca3F92Pt15Mq9CW0zEaSUeUKHD1Pc7xBvlBe4LSUpSCnzGGKLMU+CeO30qfRcnlwjZ2OeBNE5vTnDW6x0W0ZMMeZJIJ3Sv2H3uAFIVvyXoCQCY54E0iszm8/+ZzKBvDy3pcQEYx6DbYx5DLZJSfO8cPflDF17uWvtd3q7kuCgs1xrP1Y4Dfr7uYisFpFVIvKiiDTzQ7G2zDnL2LirrWvtT+nyCYfbZLjWfqxwErfVCbgT6KOqvYE0wjUoPFOszcvs75FGWg9/z3Q7vW2lA9kikg7kADvwULE2L7P6jqdYe08bt2U4wkn0xDbgD8AWwqYpA5ZjirWlDE5uW60IhxbnAx2BXGBYFOenXNDf0bw97HGq5nZzW4ZtnNy2LgE2q+puVa0BXiVcwK2ldRuDJoq1JTLoz4v0ysxmcPv1bsuwjRPzbAHOEZEcqzBtXdCfKdaWIjgZ8ywhPDD+DPjSutbTmGJtUXFc2mHST+zitgxbOA36exh4+KjNplhbFNzTehNjPlnNiG4X+C6pZUrOMBtiQ0qZR7KySJvfkeyP2vNEvxfdlnOE5pJFzgetSDu5u9tSoiK5lrYdg/RuXSm6pRNrejxJmnjr/0yaBHi1+zyG5d3gtpSoSAnzpHfpzPbhnVl381OkWGcbV5LaPHUr9tbd2YWikU+5rKZpQplpSCANQkG3pUREUpvn9+s+ont6gAxZSPi9rbd5a8YUznriLjr/bqHbUiIiKfvw9PwTuWz1AXpmZJATyCRDvG8cgCzJQH30L+IjqZFRe/HZbJ50HHe3KvaNaerT/0dfsu0BfyRGSJrbVuXV/anJFnYNDLL53BfclmObf3T9mKFDWsAjbitpGv+bR4T0jh345e+fZVhOcrydz82ooqrDCdTu2Om2lGPie/OktWzJrKWz3ZYRU17tPo/lC6t5KN/bb3l8PeapHtaXX382x20ZKYsve571U/rS/aSd9Mr7grOzPLe+PmXwl3lE2HbfAP446DmuynUvm3siOD6tmpKHzqXr4ysJVSQ2n1Ck+Oa2JVlZ6Dmn8/GEPyS9cQC6pjdn9YSnCLT2brmyJs0jIlNFpFREVtXb1mBBNgnzZytm6wsRcR7ZFkhD0tPh9B7MfWUardJyHF/SV2Skg4jbKhokkp7nWb6/sL2xgmyXAYXWn7HAZKcCy97O5+Xij5n52hSnl/Il/1zwImUj+7sto0GaNI+qLgCOLg7VWEG2K4HnNMxiwovhO9gV1+LfbXiu1zRaBLJpHmhm9zK+pkUgGw14s+exO2BurCBbJ2BrvePq4rZ2EAWB447jq0k9Wdj1MTqkN7cp0RBvHD9tqaqKSIMF2Y6FiIwlfGujGeFxjJzdi/L85lTlBdh8xWTAGAdg72nQqk9vdNmqpg9OIHbNs0tEOqjqjqMKsm0D6ocCHDNui3C0BXnSWtPatGb9xAyKBv3NpqTkpWjkZApOGE13jy00tPuo3lhBtjeBG62nrnOAsnq3t2NyyUfFFA161qYcgxtE8qj+IrAIOFlESqwibI8AQ0RkA+HI0bp3wLMJh94UAVOA8ZGIOKF3JeNafmVDfurw7wv/QufF3rqNN3nbUtXrGtn1vYJsVkHa26MVkReAnIB5zXAsOqQ3p2/eZkrwTp1S38wwG+CUrB3s/Pm5EPDGIjdPmGdvKJ2y0CG3ZXieQdkhPv3F4wSaeSMxhDfMsyqTZw70cluGIUo8YR6DPzHm8RlZksHpn1QifXq7LcWYx4882n4Ftce5/3RqzONTdp3dzPVsqsY8PuXLe56i6KZ2rmow5jHYxpjHYBvPmOflR4dyddEQt2X4ipE//IhNj7hXNdkz5mn5/CI+39jVbRm+4uHj19D3AvdeKHvGPAZ7tM6sdO2py5jH5zzRaQl/mTet6QPjgDGPwTbGPAbb2A36myQiX1mBfa+JSMt6+x60gv7WicjQeAk3fEvrQIBNvx9AWpvWCW3XbtDfPKC3qp4OrAceBBCRUwkXbOtlnfOUiA/Tc/mMVmk5bPjpZGiT2NBkW0F/qjq3Xk2txYSjJCAc9PeSqlap6mbCa5m9nWQmidCcrHBodoKIxZjnZuAd63NjQX+GBPD2rBcovaVvwtpzZFMR+SVQC0y3ce73gv4Mzkh0Znsnlf5uAi4HRlpRExBl0F+qF2uLB71HrWbHxMRkU7VlHhEZBtwHXKGqlfV2vQlcKyJZIpJPOFvGp85lGiLluRMXEDy3LCFtNXnbsoL+BgFtRaSEcH2tB4EsYF64yB+LVXWcqq4WkRmEK/7VArerqj9y4Ruixm7QX6PV+1T1t8BvnYgy+AMzw5yELOn/DJesOhj3dox5kpDmgWaclLk77u0Y8yQpp2Tuovg3A46UjYoHxjxJSq/MbNaMfpJAdvzS8RnzGGxjzGOwjTFPEpMmAa5atIHgYOfpsBvCmCfJGdtiO7U58VkVY8yTAmw/Pz0uiRGMeVKA9aMms/7G3Jhf15jHYBtjnhRh8Y8fI/uj9k0fGAXGPClCu7RcCprH9pWFMY/BNsY8KUS7zIPQ77SYXc9W3Fa9fRNFREWkrfU99sXaDDHj/jYbeP6Vv8Ysj7PduC1EpAtwKbCl3uaYF2szeBe7xdoA/kh4HXP9ckkxLdZmiD1tAtlctPKbmGTWsLsA/kpgm6quPGqXidvyOGkS4P42G9Bs59lUo47bEpEc4CHCtyzbmLgtd6jRINdtGkrg4CFCDq9lJ+ivAMgHVlqRE52Bz0SkHw6KtdnQYYiCPcEKdgUD7Aw25+CF+yG0x/E1ozaPqn4JHMnhKiLFQB9V3SMibwITROQloD9RFGszxJd+r91D4Z1LrG+xiYayFbelqo2F3swGhhNOcFAJjI6JSoMtghri8suuJ3DwEKccWB8jy3yLk2Jtdfu71ftsq1ibIbbc+PVA1k4JVxFqs2YptbW1TZxhj8Tl4zDEnTOXXktFZRbZn+ZywtSFwHfnUWKNMY/PKQ1WMLsiH4COd1VSW5y41LrGPD6lxkoBMOHrKyg7f6+1dUvjJ8QBYx6fMuSWcWQvWAPBCtc0GPP4iFfK83hy/H8CkPPpOoIV7hkHjHl8Qe/FI6len0fWXqHje+GBsBfy1hjzeJTy0GHGbw2/AWo7JYesdxa5rOj7GPN4jD3BCspCyvKqTuwa8A0AWSx1WVXDGPN4jIF/u5cuv1notoyIMObxAGWhQ1x36SioDdKtdLUnxjORYMzjIkPXXs7uf3VBQtB27WJQfy0uMOZxgbOWjaCqJp2seXkc/1fvDYQjxZgnQZTUlrPwcHhR5QmjdxPc29DKXn9hzBNnghper3flyptp+6P11lb/GweMeeLO8J/cRNqaYtrVljhe9uk1jHniwKR9Bcy67yIAsr9YTbCysokz/EkkKwmnEq4xUaqqvettv4Pwwq8gMEtV77O2PwiMsbbfqapz4iHci5z88Y3Ubs8htyRAh9nhuZpk623qE0nP8yzwBPBc3QYRGUw4RusMVa0SkXbW9vrF2joC74lIj2QuIbAnWMH/7gz3Mif+SZBFi11WlDgiWYa6QES6HbX5NuARVa2yjim1th8p1gZsFpG6Ym3+fR5thP3BSg5riBe/OZ0NfasAEI4OY0tu7CY66AFcICJLROQjEamrEBZx0J+IjBWRZSKyrIYqmzLc45JfT+SmEy9gzmkt3JbiGnYHzOlAa+AcoC8wQ0ROiuYCfozb2lhTzoQf3gJA+5I1BH02Ixxr7JqnBHjVipb4VERCQFuiCPrzE2cvv4ba99sSqIZ2q/zx0jIR2DXP68BgYL6I9AAygT2Ei7X9U0QeIzxg9m2xtqCGOG/lNQRDATJebk3bF4xpjsZusbapwFQrZ081MMrqhXxfrK2ktpxV1W2oCGXR8sdb0Sr/jccShZOgv582cryvi7UN/OAuCm9abn0zxjkWZobZYviQEbBrDycf+iqpJ/ZiScqaJ3/WrXSe/e1MRe76z9GaahcV+Y+UMk+V1nDKW7dDCPJfD5Exd8mRfan90G2PlDDPgsMw95vTOFCTQ4/xyyHkqzG8Z0lq8+wPVhJEuemtiRTetRgzAI4tSW2ekX2uIrR3H4XBJU0fbIiapDPPkLU/Iu2+PAB091fmFhVHksY8hc/dRtYBoeWGILnLTU+TCHxtnhVVVdy9YQQAPZ7cSu3WEpcVpRa+NM8X1Yc5GMrkvzZeQ9alxUD4XYghsfjSPHffcjsZ7y0nk6/dlpLS+MY8Lx1sxfOD+gOQuXulmdTzAJ43T8GMcbRfDBkVIZrt8OXqjqTFk+bZH6zkrLl3AlD4r8PIJytcVmRoCE+ZJ2NnJpP2FVBU2Y4eY5a5LcfQBKIeWIebJ621v1zstgxDA7ynM5erap+G9nnCPCKyG6ggvJQ1lWiL93/nE1X1+IZ2eMI8ACKyrDGHJyt+/51NgVqDbYx5DLbxknmedluAC/j6d/bMmMfgP7zU8xh8huvmEZFhIrJORIpE5AG39cQLESkWkS9FZIWILLO2tRaReSKywfrZym2d0eCqeUQkDXgSuAw4FbjOyvGTrAxW1R/Uezx/AHhfVQuB963vvsHtnqcfUKSqm1S1GniJcI6fVOFKYJr1eRpwlYtaosZt80SczycJUGCuiCy3asoDtK9X/Xkn0N4dafbw1IvRJOd8Vd1mpeCbJyLfqeeoqioivnr0dbvnScp8Pg2hqtusn6XAa4Rv2btEpAOA9bO08St4D7fNsxQoFJF8EckknAzzTZc1xRwRyRWR4+o+A5cCqwj/rqOsw0YBb7ij0B6u3rZUtVZEJgBzgDRgqqqudlNTnGgPvCYiEP47/6eqvisiSwmn5BsDfA1c46LGqDEzzAbbuH3bMvgYYx6DbYx5DLYx5jHYxpjHYBtjHoNtjHkMtjHmMdjm/wGHbkwyasCbEgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Fp12cqRrF6uP"
      },
      "source": [
        "We write a class that estimates the length based on the skeleton of the polygon."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ITR5WXzgF6QT"
      },
      "source": [
        "class LengthEstimator:\n",
        "  def __init__(self, raster_arr, gt):\n",
        "    \"\"\"\n",
        "    Initialise an estimator for the length of a polygon by calculating the longest path between endpoints of its skeleton.\n",
        "    Run method \"estimate\" to obtain results.\n",
        "    Args:\n",
        "      - raster_arr: 2-dimensional numpy array representation of a raster layer, where 0 = null and 1 = polygon shape\n",
        "      - gt: GDAL GeoTransform object of the raster layer\n",
        "    \"\"\"\n",
        "    self.raster = raster_arr\n",
        "    self.gt = gt\n",
        "    self.skel = skeletonize(raster_arr)\n",
        "\n",
        "  def plot_skeleton(self):\n",
        "    nonzeros = np.nonzero(self.skel.astype(int))\n",
        "    for x,y in zip(nonzeros[0], nonzeros[1]):\n",
        "      self.raster[x,y] = 0.5\n",
        "\n",
        "    plt.imshow(self.raster)\n",
        "\n",
        "  def estimate(self):\n",
        "    \"\"\"\n",
        "    This method uses the skan package to generate a DataFrame where each row corresponds to an edge in a graph representation of the skeleton\n",
        "    (nodes = junctions and endpoints, edges = path between them). We rename the nodes and create a sparse distance matrix that stores the weights of \n",
        "    the edges (negative distance between nodes). We run a shortest path algorithm on the matrix, the size of the smallest value denotes the longest path,\n",
        "    since shortest path with negative weights = longest path.\n",
        "    \"\"\"\n",
        "    branches = skan.summarize(skan.Skeleton(self.skel)) # generate df of all edges in skeleton\n",
        "\n",
        "    all_nodes = sorted(np.unique(np.concatenate([branches[\"node-id-src\"], branches[\"node-id-dst\"]]))) # all unique nodes in the skeleton\n",
        "    branches[\"endpoints\"] = [(all_nodes.index(a),all_nodes.index(b)) for a,b in zip(branches[\"node-id-src\"], branches[\"node-id-dst\"])] # to tuple indices\n",
        "\n",
        "    matrix = np.zeros((len(all_nodes), len(all_nodes))) # empty matrix\n",
        "\n",
        "    for i in range(len(branches)):\n",
        "      matrix[branches[\"endpoints\"][i]] = -(branches[\"euclidean-distance\"][i] * self.gt[1]) # negative distance between all nodes as weights\n",
        "\n",
        "    graph = sparse.csr_matrix(matrix)\n",
        "\n",
        "    paths = sparse.csgraph.shortest_path(graph) # distance matrix with longest paths (because of negative weights)\n",
        "\n",
        "    return round(-np.amin(paths), 5) # longest path in the graph == longest path between endpoints"
      ],
      "execution_count": 17,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 305
        },
        "id": "1aeOP9IiY0R_",
        "outputId": "68cf4b31-517a-41b9-f7c8-f03e3ef936dd"
      },
      "source": [
        "estim_skel = LengthEstimator(raster_arr, gt)\n",
        "estim_skel.plot_skeleton()\n",
        "print(\"longest path in graph = length of the polygon\")\n",
        "\n",
        "length = estim_skel.estimate()\n",
        "\n",
        "print(\"length:\", length)"
      ],
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "longest path in graph = length of the polygon\n",
            "length: 21090.27057\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI8AAAD8CAYAAABO1ff8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dd3hUZRaH3zMzSUgh9B5KgCBNerMAIlVURFQsqDSliVhQLFtc17XiuuKuoCIgNkRBsdBFEJFi6FVIQgKEXgPpmZlv/5gJREzI9Htnct/nyZOZ286ZyS9fP98RpRQGBp5g0toBg+DFEI+BxxjiMfAYQzwGHmOIx8BjDPEYeIzfxCMi/URkr4gki8iz/rJjoB3ij3EeETED+4DeQDqQCNyrlNrtc2MGmuGvkqcTkKyU2q+Uyge+AG7zky0DjbD46bl1gENF3qcDnUu6OFwiVDmivTJYvWUuFU12t+/LVIojO7yzHcpc4OwppVS14s75SzylIiKjgFEA5Yiis/T06nljFyQzMDrT7fvyVAHTziV4ZdtTpn19Ew3+tk4T267yo5p3oKRz/hLPYaBukfdxzmMXUUp9AHwAECuVNZtgi5AwHq+UpontjAE/M7vuNZrYLokm7+ajEne4dK2/xJMIJIhIPA7R3APc5ydbJTIouTf/qPc9rcLLBdq0S7xQbTcv9NG+D3HWls2Q5DsBONCrAeWbdCHivJ1y3/92xfv80mBWSlmB8cBSYA/wpVJqlz9sXYncsZV5Pu12jlrdr87KAgetmWzNy+PzC1ch94PcD+0H7GT95Pfo/8rKUu/3W5tHKbUIWOSv57vCouVz6TZuFN26PE3Sg9O0dEWX3PzOJGq/uRZLwwYsTFzg9v0hP8L8/ttvY4uxc+0TY7R2RRecsmVx87UDuLnzLeS0zeGl1EReXjH3D9c0Xjmc1b0blvoszXpbgaJZeBQIlDtdoLUrmtFk9lgq7XG8tlsg9sPDmER4q95cOkWE/eHaxnPG0HBeDtZjx0t9bsiLB6BVyzT2nGlIwqph/N59JmYJ+QKX1IJMeq1+FACLDTLjBAB7GKxvtqDE76D6byDrtrlkI2TE8+WJjlwd9z2NwmL+dG5BwlLeqtKQxWO6MbV1PGFipXtUkqNUCjFW5JhJyqtJ4vl4Et6xAtD5w9W8WK1of6V44czIqEl4ps1lW36Z23KXWKmsvB0kBDi6oBnbO80p8fxZWzb3JdyIslpJ/rgle7rPIEzMXtvVA3nKUS13efkxanzwG9IigcWLS/4uiuPmdn3/VF39qOZtUkp1KO76kCl5XKGSOYovkh1d0Gvf7USXVePZ9I/g74WlFmQyrlkfsNth7hnmPv8L8AsQ6Ve7ZUo8ABVMji/0jZEzmbj5Lq6fMJo177yvsVfu0/blcVTZlQuANcpMj3W/AjC4wkwqmNyrjjfl5fPM8DFYTm13674yJ55Cbo7K5XSrxbxwZhBXzRoLwMRB3zKqwhGNPSuZ9bk2hs4ZD0BBcxvnrnb++SJsvFCtcKTaPeFMPtOIjz7rS9yqtbjbgCmz4gF4MPYUHfu9zcNPPQHA5Pg+mNsuYWSFYxp79kdeOXUV2y/UYf+5KlTb4lg58NC/vmFY7Amvn70gvTVxr6716N4yLR5wjAMVVlsdXhjLm3sHceOIN4gvptcWSGzKToo1B4Cv37mR6t8mE9O2Aj995Lsq9oQti4zsSDz9pKE/4OEGG1+cRl58HqPveURrV9iab+XR+tfxaP3rqHRvOou2Leenjz70qY2e7zxNnUGeTzmW+ZLnclbfOIWXmvemf8+7+HL5J8SYAjsjn/DJWBKmH8NWtTyTUj4DoF34Wtxty7iEl6M0hnguI84Sw6Qay7n3zeF0f+kJxA6V7klnebPv/WJvfmYsr04ecvG9rY2ds/81ERV2jp6RhQN2vhdO86njqL/kDO6vvbxESIlH/VqJ4TW7MqveL149p1FYDL+2+YK2qx/FZIWUnXUYaOnLgoSlPvIUBiT1Y9/xauRnhxMTKxePd2mzj8/jS18O4Sl5qoA2a0fQcPZBrIfSvXpWSImn9uS1bM2+Fv7inXgAwsTMzsemAtDwm9Ec+qwh85+K5Y6Y8x4/M08VMC+zJgBHPoknPvEcZ1rHsOG1qV776yoZ9nzq3b0bq931aYiSCCnx+Iv9t7/Pna16MaNrF+7Ysszt+23Kjh3FypwYPm7WAICma/bw6T9X+dZRF8j34XSUx+IRkbrAx0ANHE2vD5RSU0SkMjAXaACkAYOVUme9dzV4aTNlPHHTdkB8HT45MAOASqZyQGDn1Z453obt3cqD/YJPnudNV90KTFRKNQe6AI+ISHPgWWCFUioBWOF8H/S8Xm8BVRfk0mP4Q2zKyy/1+lO2LHqMeJgewx8iq3kezVdl0eWT7VQ3R1PdHB3wCdlmvz7AlnGtsV/wjXDAi5JHKXUUOOp8fUFE9uCI17oNuMF52WxgFfCMV17qgEZhMbwVt5hONzfj7vkTQAlxrY+yquWl5ZtNfh6KPd3ZMxKFvb+jihjd/ieeqZKkhdsXyTsahazb4NNn+qTNIyINgLbABqCGU1gAx3BUa8Xd84e4rWCgqjma/Xe+T8e/jCXqpI3jWbUZXfFS6Ez0mmhiDzjW0FijTKz+zzRdLDx7/XQC0Qd9X9J5LR4RiQHmA48rpc6LXOp2KqWUiBTbQtNL3JYnJL7sWMbRdM0DpHXOxVyxIgC9Vqxjcs0tRa7UXjgAi5+9gdoLPZu/uhJeiUdEwnAI5zOl1NfOw8dFpJZS6qiI1AK8n73TMaZWTfl+kWMkWA+lTCDx+NOKo4iZAexRSr1V5NR3wFDn66HAt5675z61PtvFDSMfDqRJzGLSpXBsyk7fgQ8QudI/IXPefOLrgAeAG0Vkq/OnP/Aa0FtEkoBezvcBw3Yug6i0c4E0qUs25eXT4dXxyJa92LOz/WLDm97WGkBKOO39gmQDr0grqEr1/7m/wMsd9FfWGnhNujWTlRnN/G7HEI+3SEmFr3b0TRxNUsc8v9sx5rY8pM1r4wgrB9O/ewM8XosX3Bji8RBLjsIeJsRZ9CWcRnPHkPBZpl/bOoUY1VYI0XDZSBrNy0Vt3BkQe0bJE0I0feUctn0pAbMXmiWPzc4JW5bWXgSUE7YsCHDoeEiKx7YvhQfju5NaUDZ2BEspyOSB+t2wJe0PqN2QFA+Aslq1diGw+GBZqbuErHjKCu+dq8Mdb07SxLbRYPaAvntuITNOKGiUo6kfk8804v2lvWn0X98vt3AFo+TxAPPoMCp1Ok7KjbM09WP6d31oNHG9ZvYN8Rh4jFFtBSndxo2i0bJtXkV8eotR8gQp4RlWv63TcRWvxSMiZhHZIiI/ON/Hi8gGZ5K2uSIS7r2bBoXYlJ34RQ8Rnq79gjdflDyP4UgRUMjrwH+UUo2Bs8BIH9jwiNeO9+ZgCKUOOGvL5m8n2tB0/M6ADwgWh1fiEZE44GbgQ+d7AW4E5jkvmQ0M9MaGN6R1yuGFIzdpZd7nzL2QwKa2Juy5uVq7Anhf8rwNTIKL7bYqwDln4hJwJGmr46UNA53iTfTELcAJpdQmD+8fJSIbRWRjAf5f9eYLUgsy6TvwARrOOcLnzWcH1Ha7jXfz7T1dA2qzNLzpql8HDHBGTJQDYoEpQEURsThLnz8laSskGIP+cpUJftvB7ZV3BHTPwsYrh1P7yzDs26+c/yrQeFzyKKWeU0rFKaUa4EjG9pNSagiwErjTeVnA47ZCkSpLyxH5rb6EA/4Z53kGeFJEknG0gWb4wUaZYWteHuZ8fRbMPhlhVkqtwrEbBkqp/ThSYxv4gL90u4Pyh7Sbv7oSxvSEiww/2JVjD9XhldSPaRVuJtAbM+kRY3rCRS4URCAnz9A+IjwgGzOlFGTS5ekx2E+e8rstTwl58axfdDUTjnTU2g23mJ8ZS+8fJlLhs/W6GRAsjpAXT71/rmXxj8Wmi9ItHxzqRsIjvt3Fyx+EvHiCjQJlI98eHO0pQzw6o9UHjxLR/2jpF+oAQzw6Q6ygCkrfbVUPGOLREfHfjqLOam0X1buDMc7jAjMyarI9vQ6xt/hn11absnP3/j40e/sMtr3JfrHhDwzxuMAb826n2g7Fr2/7J5ltpsojs8c5lFW/YzrFYVRbBh5TJsST8Npurpk4Rms3iuWtMw257/q7gzI8ukyIx3Yug8iTBVq7USwXbOWwph3U2g2PKBPi0Suvn07gk2XdtXbDYwzxaMSy7DBmfN+LRk+v09oVj/E2eqKiiMwTkd9FZI+IXCMilUVkuYgkOX9X8pWzocQ/nxlB/PPBKxzwvuSZAixRSjUFWuOI3wrJfFsGf8ab6IkKQDecy0yVUvlKqXM48m0VhhZoGrelV7qOH03sKm3zb/kCb0qeeOAkMMsZbvyhiETjYr6tskhKQSZNZo8lZskObKfPaO2O13gjHgvQDpimlGoLZHFZFaWUUpSQ+j0Y47a8JaWgEvHPrdN8gwJf4Y140oF0pVThqqV5OMR03Jlniyvl21JKfaCU6qCU6hBGhBduuIbJpjTd4DLbnk9aQTXN7PsDb+K2jgGHROQq56GewG40zrdVEuaVmxnXvC95SpvBwu5b72d+s+qa2PYX3k6MPgp85txGZT8wHIcgvxSRkcABYLCXNgx0ilfiUUptBYpbIGzk2ypCo5+GU/ez0FvAYIww+5nrtg+i7hwLEYsTtXbF5xjicYGCSnbyKgjfZbm/GCz2qXAiFoaecMAQj0vsH/Q+DR5M4r3rr3f/5hD+hkP4o+mD2QtncGL8tVq74RcM8fiZ6uZo7KHXVgYM8QSESv2PcPrha7R2w+eULfHYbEw43I1TAc7FtarlAuS20wG1GQjKlHjsubmkdcrh28xGWrsSEpQp8Rj4FkM8AeKbVjNptim0Ws6GeAJEnCWG/hW2ae2GTzHEE0DqWM5zcsw1SIT/l6AEAkM8AaRFeCSb/z4NU2ys1q74BEM8Bh5jiMfAY8qkeD59/Bb67rlFM/t1fsjGdkM7zez7Cm+D/p4QkV0islNE5ohIuWBI1ha+dCMpx6tqZn963V/JrRKmmX1f4U3cVh1gAtBBKdUSx67W96CjZG165mwTM+YmwT3S7W21ZQEiRcQCRAFH0VGyNl8SbrIhkeV8Ni+269Gp7Hmyik+epRXeRE8cBt4EDuIQTQawiRBN1vZ5/EpeXvkV98d3Z1VOmWwq/glvqq1KOEKL44HaQDTQz437gy7oL0JsKKsVG+KT5/3Qbwp5yxr45Fla4M2/UC8gVSl1UilVAHyNI4FbRWc1BqUkawtk0J8eaREeSY8a+7R2w2O8Ec9BoIuIRDkT0xYG/RnJ2soI3rR5NuBoGG8Gdjif9QFGsja3KG/OxVK/rtZueIS3QX8vAC9cdthI1uYGT1bez8hfd3F3g65Bt6ml0W0w8JgyJR6JiMC8sjaRP9fgf53muH1/fYuFyJ9r8Myro3j40HU+8ytGIoj6qRLmqxr77JmBoMyIx9KgHql/a8f3TX5gQcJS+kW5PzwQZQpnQcJSqq86xs4zNX3mm1lMfN14OfbYSJ89MxCE1rrIErDUjeNI/zj2jphKGfp/8TshLZ7CFXt7J9QlechUjb0pHXu4GTGZwW7T2hWXCGnxvLH3ZxpbTITJWoIhG/H3X06n3f8eI+7VtVq74hIhWYZb4utz065zNAsLI8oUmGzEviBCwlBB9BcJIlddw9qzPamTy/N4pbSgEU1ROt+6g8PPBsfGCCEjnuxBnckY0oX9g03svvZTv9o60r8Wx45VZPIZ36/HmVXvF+r0Do5EJsEvHhEsdWrzlzc+Yv3k90i9dbrfTW59dipV1obzyay+fnl+dFgellq+GwrwF0HfYDZXrMjCxEVau+FTvm68nE1r83k+Xt+zPEFd8uT368hLm5dq7UaZJShLnn3TO9K44TFaxG6nfYTu1teXGYJLPCIcnnQN/7nhYwZGa7ebeyCoZs4n/flrqTdlG/aswO4n5CpBU21JRASqSyt+Gf9myAsHoJ4lhl3jp2KqrN90ZaWKR0RmisgJEdlZ5FixCdnEwTvOmK3tIuJ9ZJvJjFgs0KoJy+bPppLZP7nN3cVuBiX4Px1BmAXEN2umfY0rJc9H/Hlhe0kJ2W4CEpw/owCvE5Fn/BDP3LRfmPeN/7vg7rDmr1PI65DJgNtH+NXO56vnkDGks19teEqp4lFKrQYuTw5VUkK224CPlYP1OBbD1/LUuQprqvBxi9lUMEUSYyrn6WP8QpQpHJNJYcrz7+q/CqZIlEmfJY+nDeaSErLVAQ4Vua4wbusobmAqX57fJzdjbb23qGWJ8dBFA3/jdW9LKaVEpNiEbFdCREbhqNooh6MdI+1bkBkfQ16sidQB0wBDOACnr4ZKHVqiNu4s/eIA4ql4jotILaXU0csSsh0GioYCXDFuC0e0BbFSWZmrVGbfxDCSb3jfQ5cCy1FrJjarGWtF/492JA+ZRqOaw2n8gN9NuYWnXfWSErJ9Bzzo7HV1ATKKVG9XpNfPaSTf8JGH7gSeW/71NJGbolg+Z5bWrmiGK131OcA64CoRSXcmYXsN6C0iSTgiR19zXr4IR+hNMjAdGOeKEzVbZjOm4u8euK8d7lfU3rGm+3+JW6+varzUMlcpdW8Jp/6UkM2ZkPYRd52INTl6LwYlU8sSQ8fYVNLRT57SoBlhNoCmEUc59sS1YNLHIjddiOe03UKGPUdrN3TPDZF2fntqCqZy+tgYQh/i2RnOjHMttHbDZWzKjhLH9ERZJrhm1XXCLTfdBy+fZl3bjwF9jXwHEl2UPMGGZOcRbrFpMmUSIWG0+jUb6dAy4LYvxxBPEPJ6ja1Yy2vfOzXEE6Qcb19O891UDfG4yUFrJiq6HBFmbUOCdzw5leRh1TX1wWgwu8Ge/Gwej+/KxKR59Iny8yKwIMAoedxFKcxi19oLXaAb8cx9vS+Dkntr7UZQMeTmn9n/mnZZk3UjnoqfrGNLSj2t3QgqXqi2m45dtZtQ1o149M4pWxY/ZjXD3CyBcqKf9k7l8GzNel1Gg9lFxqTdRvadJhZt+Qo9/c/9r84GUpavYFz96wNuWz/fgkHQYYjHwGM8DfqbLCK/OwP7vhGRikXOPecM+tsrIv7Zg8TgD1Q2mdj/xjWYq1QOqF1Pg/6WAy2VUq2AfcBzACLSHEfCthbOe6aKBOH2XJfxXVYUv5+sQUa3eK1dKZZK5iiS7p8GVQIbmuxR0J9SalmRnFrrcURJgCPo7wulVJ5SKhXHWmZ9bzLjAk9/MZQKX5Tn1yn6juxQURGO0OwA4Ys2zwhgsfN1SUF/BgHgh4WfcuKhjgGz55VMReQvgBX4zIN7/xT0Z+AdZgls/8ebTH/DgFuAIc6oCXAz6K+sJ2vzBy2H7uLoxMDspuqReESkHzAJGKCUyi5y6jvgHhGJEJF4HLtl/Oa9m9pxZ0ovKiZp7YXrfFx/NbZrMwJiy9Ogv/8B5YHlIrJVRN4DUErtAr7EkfFvCfCIUio49sIvhlU5Jk78uyFVvtuttSu6xNOgvxKz9ymlXgZe9sYpPZBtz+f1Ft24OXElMxf0ouqOAIeIBgHG3NZlJKwaRsJjh5HwMF7c8y1tw02MHLaVIfvupn/rPy4ZSRudwO5x+kuIsqHzDKbtbMGPLcv71Y4hHicd/zqWiPN2pL2J2AU2TFhpH27GLCaqmqOZ3Ggeryzo/4d7cnYXcP2E0SiT8MO/39LNlncxpnI0DD+Jo2XhP8q8eI5aM7l20ZOENQZRJup0OMIX8T85z15qErYKL1fkuINnos6w4IRjMVa7hY9fPF6pdgabO8z1u+9Xomn4cdL+NYT4lzaj8vyTt77Mi+eANZImY37j9dQNtIlwb8jg9RpbeX34VmzKTs+HRmPJcfQNzjapzIv1mwPweJVNVDAFPoNfi/BIdg9/l1vf7IHNEI/vKVA2ztmjMEVFYfJizxSzmFg149KGmxOOdGR951gA5m2J587yqRfPaSEkf1GmxdNn9yCiRgsLkr4iQnwX/flO7UTyUhyrDQfePJSvdu4DwJQQz6IVX/nMjtaUafEoJWCzESFhPn924TPvn7uUCzaHMLdnWeh3m2NvuBbv7eLftTb73G5RzGJi4LokvhzXD/NK39sqs+Jpt/FuMndVxvLvC361M6T86Yuv02P20n1CMwBSl3RmSVYXAOwW2D76v35JLjeqwhE+jzL7JUlmmV1JmL+2ChX24ffEbkWJs8SQ0nMWKT1nYQ+DmHTl+Dmk6LN7EL333Mr8zFif2z1yvcUvGyOUyZJne34udjPkV9Bug52k+y9tjp9hz+G+bvdCgeKpFwdTt8elc4VjTd6wb+g0GkaPJmGjV4/5E2VSPJPuGEn+xGx262T31QqmSBauWQBA83fH8bcRl9bkvJSaSKcIfVYQZVI8embZ6Dc49dClBvyzd47AfOwsyePqs3e456k81t/+FqNaDyKn+3FfuAmUMfFk2nPp+dwT5P/9LG81/0Zrd4olzhJDXJG/SuyUY2QWRJC/10aXp8dcPP7o37/6Q2O8NKqbo2kUcxJf7iFfpsSTq2xUmpPIXc+nMyA6u/QbdMCXDVcAMDr6Gn453Pbi8b+uGsTfIuzUq32aVS0XaOJbmRHPKVsW751tj71zSyqa92ntjtu8H7cOxq67+L79P8ZSZUc2p1rX4u06DS4evzVmJ43Cit/su3r4BejUBX7b4ROf5NIK0hIuEJmJY7npCaVUy8vOTQTeBKoppU6JiABTgP5ANjBMKVXq6FSsVFadpSf7ZrUntW+JS4W84pnjbdh+fTQL9q3yy6CgVtyfdgOnul8qRU/Mj2dNu0+A4jdGP2HL4oH63cDu2hq9H9W8TUqpDsWdc6Xk+QjHysGPix4UkbpAH6BoBvmiydo640jWps9MYyHC7Po/kZOaf/F911c6MOiurtC8MUsWuh2X4BaurCRcLSINijn1HxzrmL8tcuxisjZgvYhULMyO4wtnPeWqXx4kdkk0gzcsC6lSBxxTEDFF5uVefOIj9oytw5bzR+k9eNjF4ymjTOzvNZMqpkhu3Haen+9ohW1file2PWrziMhtwGGl1Db5Y/5LnyRr8zUFpyIpfyifpyt792UFAwOisxkQnURK7BZ6Dbu0xsh02kLT6eOwhyv2PPguqyKLrYncwm3xiEgU8DyOKstjAhW39fzxVmAVjnUpW+E9jcJiSO3/4cX3bRPvIfLzitgtcEfyTZgu5ODt5nielDyNgHigsNSJAzaLSCe8SNbmgR8usebFLpi7iC7XGgeS5e1mcLy1iWO2GN5IaONyg/lKuC0epdQO4OIeriKSBnRw9ra+A8aLyBc4GsouJ2sz8C+dvnmShAkbnO98Ew1VqniccVs3AFVFJB14QSlVUn96EY5uejKOrvpwn3jpITf1v4+otw6zovFcymK+Upuyc8tN92G6kEPTc/t8JJlLeJOsrfB8gyKvPUrW5i/k0DFqRhVQr4xlSH7wQDf2THdkEaqyOxGr1T/pu0NyhDndmsmNvz5C2Kho7qjwtdbuBIy2ifeQlR1B5G/R1Jy5FgB/hiqGpHg251Un/t5tvHdgDfElDNWHCidsWSzKcmw6VfuxbKxpgdtaNyTFU6AsAd3kSAsKnFsAjD8wgIzrC2fXD5Z8gx8IuW+44+bBVP2rhU9Tp1PVHLqlTu+HxhC5ejfYsjTzIeTEY7WZkLwCqpqjtXbF58zPjOXdcXcBEPXbXmxZ2gkHQkw8nbbcxYWkipx/3j8RklrRcv0Q8vfFEnFaqP2joyGsh31rQko8+UuqUd4OW+6bpbUrXpNpz2XcIccMUNXpUUQsXlfKHYEnZMRz0JqJsoAtyD/RKVsWGXbFprw6HL/mPAARJGrsVfEE+Vd9iTH9RqBeOcuOTnO0dsUrur3/NHX/tVZrN1wiZMQTzGTYc7i3z1Cw2mhwYpcu2jOuEPTiOWXLos8rT5HzYiavttBnRERJ9N1zCye/qovYoeqe9VDKkmC9EfTiyVaKatPWMfbJZAZGZ2rtjku023g3eQUWIpbHUu09/TWEXSWoxXPKlsW3F1pgatmUKNmltTtXJN2aydpcx2b4NYefxHb6TCl36J+gFs8zh/ty9KYwFu/6QmtXSsSmHOv1bts2gqq3Fob8BL9wIMjFEwz0v2MY5t1pVLeme73sU28ErXgSVg0jel0UPX9cr7Urf2LymUYsnHQjAJHbd2HLDo7oVHdxZSVhsUF/IvIojoVfNmChUmqS8/hzwEjn8QlKqaX+cJz0SGIPWv2+u5Y7XPXLg1iPRBGdbqLWIsdYTaiVNkXxKOhPRHrgiNFqrZTKE5HqzuNFk7XVBn4UkSa+TiHw+ukElElxuoX2BecpWxb/OOYoZeq/Lcg6/ZWE/sLToL+xwGtKqTznNSecxy8mawNSRaQwWZtP+6NLn+iG/TbFrvHaRUSctWWTq+zMOd+KpI6OiVhhm2b+aIGnuwY1AbqKyAYR+VlECncjcjlZm4iMEpGNIrKxgOCbBe/10kSG1e/K0qsraO2KZnha7luAykAXoCPwpYg0dOcBnsRtFSgbtw4cRsybh1nSYD4Q2DU7KQWZjL/5IQBqpO/GFmQjwr7GU/GkA187oyV+ExE7UBU3gv48RXbvp1lsFk3CAiec9psGY11RFVM+VN8ZHJOWgcBT8SwAegArRaQJEA6cwpGs7XMReQtHg9lnydoOWjO5fdsI7MOqMjTK/zPnNmXnum2DsdlNhM2tTNVPDdFcjkdBf8BMYKYz13o+MNRZCu0SkcJkbVZ8mKztp+yGVL0thc8OzPLrEtN0ayY786uQZY+g4u2H/Jb0IxTwJujv/hKu91uyNjH5f+vbbj89RsKwTc53hnCuhPYDJS7QYt0Qar8dznsp//VbRET/3nfD8VNclfN7SA/s+ZKgEE9BvgVzZn6Je+15QvzCh4lbdGmkInrfFlRB/hXuMLgc3Yvn+u2DsB2PZO9Y75tOeaqApt8/AnaIX2AnbNmGi+fKdqfbM3Qvnry5NYiqJux4/D2Pn7E6F5adv5pzBVE0GbfJJ3vTGOhcPGdt2SizIyuMp/fbUAz7fiIJj63HaAD7Fl2LZ0iPIUS9e4wNLefjyUzKkA4DsZ8+Q4JtQ3PA31YAAAK9SURBVOkXG7iNrsUjVhsmMbmV9aX3nlsxT3KkHVInfzeqKD+iS/EctWbSe+oksv+Zy+QG81y6J+HjsUScEyom2YjeZJQ0gUCX4jlpt1DntbX8ff9mritXcqmzNS+Px5PuBqDJu4ewHkoPlIsG6FA8GfYcEnPisTSoR5gkUlxbZ3t+Lhfs4fw1ZTARfdIAx1yIQWDRnXiGJA9C3atYuOk7oPjd2h9/6BHCftxEOAcC65zBH9CdeEriiwuV+OQGRxqL8JPbjEE9HaAr8Vw1NZ/UXvG0nrfn4rFGX46hxnoIy7JT7qhPVncY+AhdiSe7TiT5V2fzbr2FxC+eAEDCV7nIr1s19sygOHQlnmOdzYxq9QvPHu1Jk5E+TsVr4HNKTdYWCAqTtRnojysla9OFeETkJJCFYylrWaIq+v/M9ZVS1Yo7oQvxAIjIxpIUHqoE+2fWZ7Z3g6DAEI+Bx+hJPB9o7YAGBPVn1k2bxyD40FPJYxBkaC4eEeknIntFJFlEntXaH38hImkiskNEtorIRuexyiKyXESSnL8rae2nO2gqHhExA+8CNwHNgXude/yEKj2UUm2KdM+fBVYopRKAFc73QYPWJU8nIFkptV8plQ98gWOPn7LCbcBs5+vZwEANfXEbrcXj8n4+IYAClonIJmdOeYAaRbI/HwNqaOOaZ+hqYjTEuV4pddi5Bd9yEflDPkellBKRoOr6al3y+H0/H72glDrs/H0C+AZHlX1cRGoBOH+fKPkJ+kNr8SQCCSISLyLhODbD/E5jn3yOiESLSPnC10AfYCeOzzrUedlQ4FttPPQMTastpZRVRMYDSwEzMFMppe88AJ5RA/hGRMDxnX+ulFoiIok4tuQbCRwABmvoo9sYI8wGHqN1tWUQxBjiMfAYQzwGHmOIx8BjDPEYeIwhHgOPMcRj4DGGeAw85v/q6GvoaejG7AAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}