{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "# Point in Polygon & Intersect\n", "\n", "Finding out if a certain point is located inside or outside of an area,\n", "or finding out if a line intersects with another line or polygon are\n", "fundamental geospatial operations that are often used e.g. to select\n", "data based on location. Such spatial queries are one of the typical\n", "first steps of the workflow when doing spatial analysis. Performing a\n", "spatial join (will be introduced later) between two spatial datasets is\n", "one of the most typical applications where Point in Polygon (PIP) query\n", "is used. \n", "\n", "For further reading about PIP and other geometric operations, \n", "see Chapter 4.2 in Smith, Goodchild & Longley: [Geospatial Analysis - 6th edition](https://www.spatialanalysisonline.com/HTML/index.html).\n", "\n", "## How to check if point is inside a polygon?\n", "\n", "Computationally, detecting if a point is inside a polygon is most commonly done using a specific formula called [Ray Casting algorithm](https://en.wikipedia.org/wiki/Point_in_polygon#Ray_casting_algorithm).\n", "Luckily, we do not need to create such a function ourselves for\n", "conducting the Point in Polygon (PIP) query. Instead, we can take\n", "advantage of [Shapely's binary predicates](https://shapely.readthedocs.io/en/stable/manual.html#binary-predicates)\n", "that can evaluate the topolocical relationships between geographical\n", "objects, such as the PIP as we're interested here.\n", "\n", "There are basically two ways of conducting PIP in Shapely:\n", "\n", "1. using a function called\n", " [within()](https://shapely.readthedocs.io/en/stable/manual.html#object.within)\n", " that checks if a point is within a polygon\n", "2. using a function called\n", " [contains()](https://shapely.readthedocs.io/en/stable/manual.html#object.contains)\n", " that checks if a polygon contains a point\n", "\n", "Notice: even though we are talking here about **Point** in Polygon\n", "operation, it is also possible to check if a LineString or Polygon is\n", "inside another Polygon.\n", "\n", "Let's import shapely functionalities and create some points:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from shapely.geometry import Point, Polygon\n", "\n", "# Create Point objects\n", "p1 = Point(24.952242, 60.1696017)\n", "p2 = Point(24.976567, 60.1612500)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's also create a polygon using a list of coordinate-tuples:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Create a Polygon\n", "coords = [(24.950899, 60.169158), (24.953492, 60.169158), (24.953510, 60.170104), (24.950958, 60.169990)]\n", "poly = Polygon(coords)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POINT (24.952242 60.1696017)\n", "POINT (24.976567 60.16125)\n", "POLYGON ((24.950899 60.169158, 24.953492 60.169158, 24.95351 60.170104, 24.950958 60.16999, 24.950899 60.169158))\n" ] } ], "source": [ "# Let's check what we have\n", "print(p1)\n", "print(p2)\n", "print(poly)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "- Let's check if those points are ``within`` the polygon:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Check if p1 is within the polygon using the within function\n", "p1.within(poly)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Check if p2 is within the polygon\n", "p2.within(poly)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Okey, so we can see that the first point seems to be inside that polygon\n", "and the other one isn't.\n", "\n", "-In fact, the first point is quite close to close to the center of the polygon as we\n", "can see if we compare the point location to the polygon centroid:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POINT (24.952242 60.1696017)\n", "POINT (24.95224242849236 60.16960179038188)\n" ] } ], "source": [ "# Our point\n", "print(p1)\n", "\n", "# The centroid\n", "print(poly.centroid)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "It is also possible to do PIP other way around, i.e. to check if\n", "polygon contains a point:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Does polygon contain p1?\n", "poly.contains(p1)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Does polygon contain p2?\n", "poly.contains(p2)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Thus, both ways of checking the spatial relationship are identical; [contains()](https://shapely.readthedocs.io/en/stable/manual.html#object.contains) is inverse to [within()](https://shapely.readthedocs.io/en/stable/manual.html#object.within) and vice versa.\n", "\n", "Which one should you use then? Well, it depends:\n", "\n", "- if you have **many points and just one polygon** and you try to find out\n", " which one of them is inside the polygon: You might need to iterate over the points and check one at a time if it\n", " is **within()** the polygon.\n", "\n", "- if you have **many polygons and just one point** and you want to find out\n", " which polygon contains the point: You might need to iterate over the polygons until you find a polygon that **contains()** the point specified (assuming there are no overlapping polygons)\n", " \n", "## Intersect\n", "\n", "\n", "Another typical geospatial operation is to see if a geometry intersects\n", "or touches another one. Again, there are binary operations in Shapely for checking these spatial relationships:\n", "\n", "- [intersects():](https://shapely.readthedocs.io/en/stable/manual.html#object.intersects) Two objects intersect if the boundary or interior of one object intersect in any way with the boundary or interior of the other object.\n", "\n", "- [touches():](https://shapely.readthedocs.io/en/stable/manual.html#object.touches) Two objects touch if the objects have at least one point in common and their interiors do not intersect with any part of the other object.\n", " \n", "\n", "Let's try these out.\n", "\n", "Let's create two LineStrings:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "deletable": true, "editable": true }, "outputs": [], "source": [ "from shapely.geometry import LineString, MultiLineString\n", "\n", "# Create two lines\n", "line_a = LineString([(0, 0), (1, 1)])\n", "line_b = LineString([(1, 1), (0, 2)])" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Let's see if they intersect" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "line_a.intersects(line_b)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Do they also touch?" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "line_a.touches(line_b)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Indeed, they do and we can see this by plotting the features together" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Create a MultiLineString from line_a and line_b\n", "multi_line = MultiLineString([line_a, line_b])\n", "multi_line" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Thus, the ``line_b`` continues from the same node ( (1,1) ) where ``line_a`` ends.\n", "\n", "However, if the lines overlap fully, they don't touch due to the spatial relationship rule, as we can see:\n", "\n", "Check if `line_a` touches itself:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Does the line touch with itself?\n", "line_a.touches(line_a)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "It does not. However, it does intersect:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Does the line intersect with itself?\n", "line_a.intersects(line_a)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Point in Polygon using Geopandas\n", "\n", "Next we will do a practical example where we check which of the addresses from [the geocoding tutorial](geocoding_in_geopandas.ipynb) are located in Southern district of Helsinki. Let's start by reading a KML-file ``PKS_suuralue.kml`` that has the Polygons for districts of Helsinki Region (data openly available from [Helsinki Region Infoshare](http://www.hri.fi/fi/dataset/paakaupunkiseudun-aluejakokartat).\n", "\n", "Let's start by reading the addresses from the Shapefile that we saved earlier." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "deletable": true, "editable": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/aagesenh/.conda/envs/python-gis/lib/python3.8/site-packages/geopandas/_compat.py:111: UserWarning: The Shapely GEOS version (3.9.1dev-CAPI-1.14.1) is incompatible with the GEOS version PyGEOS was compiled with (3.9.1-CAPI-1.14.2). Conversions between both will be slow.\n", " warnings.warn(\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
addressidaddrgeometry
0Ruoholahti, 14, Itämerenkatu, Ruoholahti, Läns...1000Itämerenkatu 14, 00101 Helsinki, FinlandPOINT (24.91556 60.16320)
1Kamppi, 1, Kampinkuja, Kamppi, Eteläinen suurp...1001Kampinkuja 1, 00100 Helsinki, FinlandPOINT (24.93169 60.16902)
2Kauppakeskus Citycenter, 8, Kaivokatu, Keskust...1002Kaivokatu 8, 00101 Helsinki, FinlandPOINT (24.94179 60.16989)
3Hermannin rantatie, Verkkosaari, Kalasatama, S...1003Hermannin rantatie 1, 00580 Helsinki, FinlandPOINT (24.97783 60.18892)
4Hesburger, 9, Tyynenmerenkatu, Jätkäsaari, Län...1005Tyynenmerenkatu 9, 00220 Helsinki, FinlandPOINT (24.92160 60.15665)
\n", "
" ], "text/plain": [ " address id \\\n", "0 Ruoholahti, 14, Itämerenkatu, Ruoholahti, Läns... 1000 \n", "1 Kamppi, 1, Kampinkuja, Kamppi, Eteläinen suurp... 1001 \n", "2 Kauppakeskus Citycenter, 8, Kaivokatu, Keskust... 1002 \n", "3 Hermannin rantatie, Verkkosaari, Kalasatama, S... 1003 \n", "4 Hesburger, 9, Tyynenmerenkatu, Jätkäsaari, Län... 1005 \n", "\n", " addr geometry \n", "0 Itämerenkatu 14, 00101 Helsinki, Finland POINT (24.91556 60.16320) \n", "1 Kampinkuja 1, 00100 Helsinki, Finland POINT (24.93169 60.16902) \n", "2 Kaivokatu 8, 00101 Helsinki, Finland POINT (24.94179 60.16989) \n", "3 Hermannin rantatie 1, 00580 Helsinki, Finland POINT (24.97783 60.18892) \n", "4 Tyynenmerenkatu 9, 00220 Helsinki, Finland POINT (24.92160 60.15665) " ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import geopandas as gpd\n", "\n", "fp = \"data/addresses.shp\"\n", "data = gpd.read_file(fp)\n", "\n", "data.head()" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "\n", "### Reading KML-files in Geopandas\n", "\n", "It is possible to read the data from KML-files with GeoPandas in a similar manner as Shapefiles. However, we need to first, enable the KML-driver which is not enabled by default (because KML-files can contain unsupported data structures, nested folders etc., hence be careful when reading KML-files). Supported drivers are managed with [`fiona.supported_drivers`](https://github.com/Toblerity/Fiona/blob/master/fiona/drvsupport.py), which is integrated in geopandas. Let's first check which formats are currently supported:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'ARCGEN': 'r',\n", " 'DXF': 'rw',\n", " 'CSV': 'raw',\n", " 'OpenFileGDB': 'r',\n", " 'ESRIJSON': 'r',\n", " 'ESRI Shapefile': 'raw',\n", " 'FlatGeobuf': 'rw',\n", " 'GeoJSON': 'raw',\n", " 'GeoJSONSeq': 'rw',\n", " 'GPKG': 'raw',\n", " 'GML': 'rw',\n", " 'OGR_GMT': 'rw',\n", " 'GPX': 'rw',\n", " 'GPSTrackMaker': 'rw',\n", " 'Idrisi': 'r',\n", " 'MapInfo File': 'raw',\n", " 'DGN': 'raw',\n", " 'PCIDSK': 'rw',\n", " 'OGR_PDS': 'r',\n", " 'S57': 'r',\n", " 'SQLite': 'raw',\n", " 'TopoJSON': 'r'}" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import geopandas as gpd\n", "gpd.io.file.fiona.drvsupport.supported_drivers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Let's enable the read and write functionalities for KML-driver by passing ``'rw'`` to whitelist of fiona's supported drivers:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "deletable": true, "editable": true }, "outputs": [], "source": [ "gpd.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's check again the supported drivers:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'ARCGEN': 'r',\n", " 'DXF': 'rw',\n", " 'CSV': 'raw',\n", " 'OpenFileGDB': 'r',\n", " 'ESRIJSON': 'r',\n", " 'ESRI Shapefile': 'raw',\n", " 'FlatGeobuf': 'rw',\n", " 'GeoJSON': 'raw',\n", " 'GeoJSONSeq': 'rw',\n", " 'GPKG': 'raw',\n", " 'GML': 'rw',\n", " 'OGR_GMT': 'rw',\n", " 'GPX': 'rw',\n", " 'GPSTrackMaker': 'rw',\n", " 'Idrisi': 'r',\n", " 'MapInfo File': 'raw',\n", " 'DGN': 'raw',\n", " 'PCIDSK': 'rw',\n", " 'OGR_PDS': 'r',\n", " 'S57': 'r',\n", " 'SQLite': 'raw',\n", " 'TopoJSON': 'r',\n", " 'KML': 'rw'}" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gpd.io.file.fiona.drvsupport.supported_drivers" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Now we should be able to read a KML file using the geopandas [read_file()](http://geopandas.org/reference/geopandas.read_file.html#geopandas.read_file) function.\n", "\n", "- Let's read district polygons from a KML -file that is located in the data-folder:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "deletable": true, "editable": true }, "outputs": [], "source": [ "# Filepath to KML file\n", "fp = \"data/PKS_suuralue.kml\"\n", "polys = gpd.read_file(fp, driver='KML')" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of rows: 23\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
NameDescriptiongeometry
0Suur-EspoonlahtiPOLYGON Z ((24.77506 60.10906 0.00000, 24.7766...
1Suur-KauklahtiPOLYGON Z ((24.61578 60.17257 0.00000, 24.6155...
2Vanha-EspooPOLYGON Z ((24.67576 60.21201 0.00000, 24.6752...
3Pohjois-EspooPOLYGON Z ((24.76792 60.26920 0.00000, 24.7699...
4Suur-MatinkyläPOLYGON Z ((24.75361 60.16631 0.00000, 24.7537...
5KauniainenPOLYGON Z ((24.69075 60.21958 0.00000, 24.6924...
6Suur-LeppävaaraPOLYGON Z ((24.79747 60.20827 0.00000, 24.7954...
7Suur-TapiolaPOLYGON Z ((24.84436 60.16598 0.00000, 24.8443...
8MyyrmäkiPOLYGON Z ((24.82459 60.29025 0.00000, 24.8243...
9KivistöPOLYGON Z ((24.94309 60.33845 0.00000, 24.9421...
10EteläinenPOLYGON Z ((24.78277 60.09997 0.00000, 24.8197...
\n", "
" ], "text/plain": [ " Name Description \\\n", "0 Suur-Espoonlahti \n", "1 Suur-Kauklahti \n", "2 Vanha-Espoo \n", "3 Pohjois-Espoo \n", "4 Suur-Matinkylä \n", "5 Kauniainen \n", "6 Suur-Leppävaara \n", "7 Suur-Tapiola \n", "8 Myyrmäki \n", "9 Kivistö \n", "10 Eteläinen \n", "\n", " geometry \n", "0 POLYGON Z ((24.77506 60.10906 0.00000, 24.7766... \n", "1 POLYGON Z ((24.61578 60.17257 0.00000, 24.6155... \n", "2 POLYGON Z ((24.67576 60.21201 0.00000, 24.6752... \n", "3 POLYGON Z ((24.76792 60.26920 0.00000, 24.7699... \n", "4 POLYGON Z ((24.75361 60.16631 0.00000, 24.7537... \n", "5 POLYGON Z ((24.69075 60.21958 0.00000, 24.6924... \n", "6 POLYGON Z ((24.79747 60.20827 0.00000, 24.7954... \n", "7 POLYGON Z ((24.84436 60.16598 0.00000, 24.8443... \n", "8 POLYGON Z ((24.82459 60.29025 0.00000, 24.8243... \n", "9 POLYGON Z ((24.94309 60.33845 0.00000, 24.9421... \n", "10 POLYGON Z ((24.78277 60.09997 0.00000, 24.8197... " ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Check the data\n", "print(f\"Number of rows: {len(polys)}\")\n", "polys.head(11)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nice, now we can see that we have 23 districts in our area. \n", "Let's quickly plot the geometries to see how the layer looks like: " ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM8AAAD4CAYAAABL/rJKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3HUlEQVR4nO2deXhjV3nwf0e7ZUuW5H23ZzJLZl88noRAWMKSBpoAAQoBmrI0LXu3r1/a0gW+j3yFUlqgNDSEJZQQCAkEAmkSCIQtIbPv+4w93rfxLttaz/eHZMdjS/LV1b1a7Pt7Hj+W5HvufSXrveecdxVSSgwMDNLHlGsBDAwKFUN5DAxUYiiPgYFKDOUxMFCJoTwGBiqx5FqARJSXl8vm5uZci2FgwMGDB4ellBWJ/paXytPc3MyBAwdyLYaBAUKIy8n+ZizbDAxUYiiPgYFKDOUxMFCJoTwGBioxlMfAQCWG8hgYqESR8gghPEKIR4QQZ4QQp4UQ1wshfEKInwohzsd/e1OMNwshDgshfqyd6AYGuUXpzPN54Ekp5UZgO3AauBt4Rkq5Dngm/jwZH4uPMTBYMSyrPEIIN3Aj8FUAKWVQSjkG3AY8ED/sAeCNScbXA68H7s9cXINkzIYiPHmij/t/fSnXoqwalEQYrAGGgK8LIbYDB4nNJFVSyj4AKWWfEKIyyfh/B/4acKW6iBDiLuAugMbGRkXCr3YiUckLl67w2JEe/udEP5OzYT5z+7Zci7VqUKI8FmAX8BEp5QtCiM+Teok2jxDiDcCglPKgEOIVqY6VUt4H3AfQ2tpqpLem4GTvOD841MPjx3oZmAhc9bfT/RM5kmr1oWTP0w10SylfiD9/hJgyDQghagDivwcTjL0BuFUI0QF8B3iVEOJbGUu9ivnthWF+/4u/4f7ftC9RHICH93exv2MkB5KtPpZVHillP9AlhNgQf+km4BTwI+DO+Gt3Aj9MMPZvpJT1Uspm4O3Az6WU79JC8NXIwMQsH/vOYaIp5mV/MMJ7v7GfC4NT2RNslaI0qvojwINCCBtwCXgPMcV7WAjxPqATeCuAEKIWuF9KeYsO8q5aQpEoH/72IYangsseOzkb5lDnKNdUlmRBstWLIuWRUh4BWhP86aYEx/YCSxRHSvks8Gxa0hnM89mnzrK/Y1Tx8Q6rWUdpDMCIMCgInj7Zz3/9Kj0TdHmxTSdpDOYwlCfP6bwyzV9+72ja4z7y0GF6xmZ0kMhgDkN58pjZUIQPfvsgk7PhtMde8Qf50/8+yGwoooNkBmAoT17zf358ihM96v02x3vG+bsfnCAfq8IGw9Fci5AxhvLkKY8d7uHBFzozPs+jh7r55vNJ0/A1Z7GiPn/xCn/x8BGmgy/Onvc+e5HbvvRbxmdCWZNLD/KyAMhq5/zAJH/z/eOane+TPz5FU5mTV2xIFkGlDWf6J/jrR47x5p11NJUXs6nGzSd/fIrTfROc7Z/kTTvrePxYH0e7xjCbBM9fvMLNW6p1lUlPDOXJM/yBMB948BAzGu5VIlHJBx88xMN/cj1b6ko1O+9CRvxBbv/P5/AHIxzrHsdlt1Bst9A/MQvAyd4JTva+uASNRCWPHe7htZuqMJmELjLpjcjH9XBra6tcjaWnpJT8+XeP8NiRXl3OX15i50cfvoFaT5Eu5//SLy7wL0+dTWtMsc2Mu8hKsd1Csc1Msd3CJ27dzLqqlHHEWUMIcVBKmcjHaex5lBKNSg53jnL3o8f4zJNndLlG58i0booDMDwV4G++f1w3A8IHXr6Wl60rT2uMPxihb3yWC4NTHO0e57mLVzjcNaaLfFpjLNsSMBuK4LCaCYaj/O7SFZ4+1c/TJwcYnIwFYpoEbKv38LrNVQihzZLjRM843/qd/hv7X54b4vuHerh9d73m5zaZBP/6tu3c8vlfKwojSsbjR3s50TNOldvB5lo3G6pdhMKSUDRKeYmd0iKrhlKrx1i2LWJoMsArP/ss3mIrY9OhlD6WOk8RlW47L1lbxntuaKG8xK7qmhOzIW794m/ouDKtVuy0KC2y8tO/uJFKl0OX8//y3BB3fm2f5uddU1HMDz90Ay5H9pTHWLalwReeOc9UIEzXyMyyzsmesRkOd47xpV9c5PZ7n2MmmP4mX0rJ3Y8ey5riAIzPhPinH53U7fwvX1/BXTeu0fScJXYL9727NauKsxyG8izg0tAU396nzrdy+co0t9/7HF945jzj08r9F994roMnjverumYmPHG8nydP9Ol2/r967Qa21Wtn2fvc27bnXZS4oTwL+OLPLxBJlSyzDKf6JvjcT8/x0s/8nH99+iyj/tTr/n3tI3zqJ7mri/Lxx04yNq1+b5IKm8XEF96+k2Jb5tHdH71pHa/dnH/+IEN54swEIzx1UpsZYHI2zBd/foGXfeYX/PP/nGFocmnGZ9/4DB988CDhDJQ1U4anAnzy8VO6nb/SbafYnrlNqrTIyrHuscwF0hjD2hbn52cGmVaxZ0nFVCDMl395ka//tp079jbyJzeupbrUQSAc4U+/pSyxTW++f7iHV26s5Pe312p63uGpAO/7xn7GZ0KsqyyhZ2xG9ef7f358iiKrGatZcOuOWv7vG7dqKqtaDOWJ8+Nj+vlXAuEoX/9tBw/+rpPbd9cTlZKjeeTLuPvRY1xb49ZsT3G0a4x7f3kRfzDClrpSDl6OJfFVuuw0lTkVJfVtry/laPf4/POZUISZEFSU6GMhVIOhPMRCYn5xNlH9Em0JRqI8tK+TWk/+fAEg5qj8wLcO8oMP3UBJBsus9mE/n33qLD85ntgQMTgZwOtUZi37xG1b+MqvL/GTY1efq86rT3SEGgzlIbZkmw1lL0R+SkV+Tjb4s+8c5j/u2JV2Cvfg5CxfeOY839nXtewezuNcPsO1qczJ2opiqt0v3mRsFhN337yR23fVpSWbnqx65ZFS8vCBrqxdTwATeag8ZpPgZ6cH+eCDh/jyu3ZjsyxvS5qcDfGVX13iK79uVxzIemFwirYWHz1jM/SMxjJdi2xmLCbBp964hQqXg+pSB9FozFDwj7+/CZfDSp2niOvXlmX0HrVm1SvPs2eH+PX54axdz+O0MpqGH0gPrGZBKPLiDGESsR+IzcJ3fm0fzeXF+ANhpoNhghFJNCqJRCVRKWnwFnFhyM/Foam0s1yv+INcaY/Vlat02WkuL+Zo1yizYbj/N+187Y/2zEdqfPSmddq8YZ1Y1coTiUo+rVOQZzJcjuwrT0t5MVJKHFYzUsLIdBCnzUyFy04gHse3cBP//KUrPH/pStLzRaWPIxoYPAYnA/PxggDHusd5838+xwPv2UNLRX45RBNRcMoTikSxmrVxT/3gcA9n+ic1OZdSnBo4DdOl84ofm9W8JHzosoqQoJbyYi4M6veZdY5M84nHT/HhV11Da7NPt+toQcE5ST/2ncPc/egxDnWOZhRaf2UqwOeeTi/3RAscVhNmk2Bno4dXrC+nvET/ElERCRurtcmPqXDZGZ0OUVZsm49uFgLqNbSCjc+GuOP+F/jRUf3cB1pQcMpT5yniO/u7ePN/Psfr/v1XfP237WmFmESikm8+38ErP/ssveOzOkqamMHJAJ4iK4c7x+ifCBAIR9muYQxYMqwmbf7Vczesem8R4zMhzAJsZsGoP0iRRoUWL1+ZJhiO8tGHDvMfPz+flwVMoACVZ0O1e/7xuYEpPvH4KdrueYaPPHSY5y4ME01hKj3UOcqt//Eb/uGHJ3Nm8eodm+VKPOata2SaydkwR7vH2dPsxa7AwqWWUEQbU/xc/pI9rigRCYGwxB+M4C6y0FzmpLXJi9o0p52NHkYWxAR+9ulzvOcb+zndl3/dHwpuz7O+aulGMhiO8vjRXh4/2kujz8nbWut5a2sDVXE/wfh0iE8/dYaH9nWSTzcxfzBCvbeI7tEZ9neM0ugrwiSE5ukJ2+tLOaJRbJjNLNjV6OFI59IogVBEMjAxTceVaYptZjbVurkwOKXYQGIxQV+C1cCzZ4f45bkh3rijjr94zXoafM6M34cWFFwy3EwwwqZ/fHJZJTCbBK/aWMmN68v5/M/O50UcWSJ2NXo41Dk2/9xqFuxo8KRVlzoVe1t8vNCuXcuRTTUuTvUpNxg0+IoYngomzHUyCa7q+NDW7GPfMu1RrGbBO/c28aFXXkOFS13yYTqsqGS4IpuZJgV3nkhU8tNTA/zDYydpKS/OgmTqsCyyHIYikv0do2yqcVNTmlkYz7U1Lk0VB+BU3yTX1ig3PgxNBAgkcaBurnWztyVmUXM5LIoac4Uikm8818HL/+UXfO7psxmlkGRKwSkPkFZlFQns7xilrdlLPlY4Spbzc6pvglF/kLZmn+r9g1mnN5xO/NtsOEpr04sm57YWL1vr3HidVo73TLC/Y4RGn5Nra9xpOVyngxG+8PMLfPuF7BV0XIyureSFEA4hxD4hxFEhxEkhxCe0EDrRvmc59nWMsrnWTVEO/Cyp6Ljix5LkSz4bjrKvY4Q9KvwdLeXFGZXqTUW6UQVHukbn/2e9o7Mc75mY3wdFJVSU2OYjr9PlM0+dTZgvlQ30biUfAF4lpdwO7ABuFkJcl6nQ61XW9DreM0GVy54V34pSQhFJU1nqZWi6Xw6zQJMMzmQk2tSnIhiRXBqaYk+zF8nSZZYE1cuvydkw/+9/cpONq2sreRljrr+fNf6T8SJ1QwYOv44r0wRCESpddpzW2NuvdtvZ2VDKnmYvbc2+2E9L9rzb3mV66bQP+xXPtlaTYHNtKSd69Zl1Wpu8qmpMh6Ox5XPP2FLFyzRi5PuHenghRTiRXujeSl4IYY6PuQb40oLGwIuPU9xKvqW8GLNJqL9bBSJMBl7cxPZPBOhf1Bx3TrGygZKdiZJ9hs0sWFfl4ljP+LLHqsEsYstMrZkNZ57B+/c/PMFPPvoyzUK3lKDkSnOt5O+VUu4E/ChsJQ8gpYxIKXcA9UCbEGJLkuPuk1K2SilbKyoqUp7TbjHTvMxSJ1OKbNlzgSkpnHi8exzPMolkZSX2q+pBa82ORq8uJn81/YcWc25gim/8tiPlMd98voNv/LadUCTK4c5RPv+z8xlVENK7lfw88aXes8DNaoVdiNp9j1K0CjVRQlSBry0UlVyTJNLYahZsrHYRjUrWlBezu9HLmgrtzfN6VdpJp1RXKv7tZ+foG0/cDW9oMsC9z17knx4/xZZ/fIo33/scPWPT1HvV34R1bSUvhKgQQnjij4uAVwMZ5wBEolL3yvr2LC7blHZvW6hkRTYzm2vd7G70UFtaxJn+SQYmAwxNBTjYOcqlIT/XVJawQaObTKXLzsUh7ZdsgliKhBZMByN870D3ktcvDE7x0YcOzxs6AuEoUkJTWXFGXSP0biVfAzwQ3/eYgIellD9WLS1wbmCSv37kmCb5JKmwZXHtPDGjbNnSNfLiXXVjtYvDCyIT5li4BLowOEVb8xIPQtrUlDpS+ozWlBdTVmLjRM9E2q1RvE6bZsoDL65IhiYDPH60l8eO9HCsO/Ee8LNPn+V1m6tVFz7RtZW8lPIYsFOVZEm4MDiFr9iGxSR0rXmWzY2nUlP08FSA8hIbaytKFEcOaBEOOjYdXLIHbC5z4nHamJoN0TU6zaVhP3WeIsamg2yuK2V8OsTZgeXDeNwOi6bKs6MhNpM8ebKfT/44eU06k4BP374to4pBBRcYesvWGm7ZWsOIP8jee352VTqxlpjN2QlHcBdZFM88kpixJJ2Qm97RzDtiz4SibKkrvira2VdsuyomD2K1ux0WEyd7xtlY40YJJQ7tvoJVbjs/ONzLzkYP//7TcymPXV/lyrhWXcEpzxxep5Uqt4NuDb4cichWKE9FiV2x8gBptYf3OK2a5SwttAjubvImjQiYjTfqPdY1tqRWQiKcGlo1az1FfOapM7gd1mV9UWf6JzndN8HORvXL2oKMbQN44LkO3RQHslceSss772IWlm5SQrHNTG2pA2uCWdcfCLO7ycO6yhJFoTShqGRHg2dZq2Uwwzwju8WEw2piW30po/4gUqLIiVvtdmTc56cgledQ5yifekLfkAx3kTUrcXB2i37XUFI+CmJKs6fZi8Us6B2fTWiBCoajHLw8xvnBqQRnSMyBjlHcDgtrK4oTRoXUe4o4lUGSm90SU5pAKMqx7vG08qAGJmf54s8v4A+ov0kW5LLtI98+rNteZ465fJpqt50Kl50im4VoVHJuYFLbLFQd30YqC5klXkdhcGKWwcnAVflDl4aWKohFxR5QAgOTAQbiBpHSIgtmYWImFGZzbSmRqKQ7jWXoHL5iK+sqXZzpn1Cd9yQlvHDpiiIfWzIKUnnmOixn51pXh+40lzk1VZ5Mly2pSFS3oNHnxOu0cm5gMukXzx8IU+12zH/ODquJsxpUGRpfsLc7cHmU1qb09huNPieVLjtHu8ZU5yltrnXz7uuaaCxzcv2asozaYhak8uQy+3VMRVBkyvNNByktWn6Dq4aTvbGQnjpPEcV2CxMzQc70T9G5zPcuHI21B5lTHrfDij8Qxq9xF4kDl0cVGRUgFpB64PIonSOZpajfvquet7eljp1USkHueXKJ1l/ymWA4bceiUvzBCLPBMIFwlEhUcqZf+X7lWPc4uxs9VLjsWMyCcFT7GTJL3oCr0LJbnaE8aaL1/3ttpYtgWJ+lm89po87r5MLgFAcvxzbv6XCwc4xgOEK9p4hAWPvZvsHnVLx31SorVssbVUEqTy5Llpg0ah0PsKnGzXMXr7BJoUMxHRp9RVgt4qp4NDXBtLF9ij5TRFmx8gIelzNcrgG8bF05u9PcZ6WiMJUnh9qjVUCq1SwYm4l57KMy/ZmnrcVHQ4oqnRaTiYFFOUqTs+qWnKMz+kRTByPKZ4HBiVm21ZfiyKC2XccVv6bfnYJTnlyXytIq8mBXo5feeFblmf4p1qaRQlBaZGFf+wi947O0tfhw2a/2Fa0pL+bS8NII6EvDflXLn4uDUwkdp5mSTnRBVMb2YU1l6lMtukZm+MhDhzmvIOZOCQWnPEIIvvCOnYodgFrT6HPS1uxjV6OHjdWupHLYLKakX7i2Zu8SU2uZwroKFS47O+o9QCw1Y1/7CCaTibYWL9vqS7mmsjipczcUkaqKp0QlrKvUPn9qVEVA6NmBSfZkECn+8zOD/PfvtKm4U5Cm6lu311JT6uCubx7IeruOcwNLLVYldjPhiMRmNWEzx5TGbDJR4bIvSRvY1ehhXwL/Ssdw6jW9y2Hh2ho3hztHmVmUtjw+E2JfuzJnocuefkhKkdXE8FRmFWrcRRbKiu2M+AOMz4RxWk2q84MyjaZ/9uxQRuPnKLiZZ449zT4e+9ANrMmDgoZTgQiz4SgTM2GGp4L0jQfoHp1haCJwVVmpbYua1C5kcDLAugSzgs1iYm+LDyljs0woIpkNqrfOdY2m/4XdWu+5qo+OGjZUuWgf9jM+E8ZuMVHhctDoU9dZYSRJrTulDE5q42QvWOWBWCbg9z/4kqymTKdD99jMvHVnc62L030TKYuWeIteXLoJEXMMljqsvNA+wtSCgiVDGcwCfeMBGtL80k5nEP8F4HNaaV+wBwuEo1wemWbEH0yrFkW128HORg+TaUShL8TlsLCn2ctbdzeoGr+Ygly2LcTjtFFkM+vmaMyUQ5dHua7Fx6GusWV9GnNfsG31pYxNBzmQJHq5b3wWm8Wk2j9UUWK/Kit1OTKJqqgpdWA1mxJGBozPhJmcDbOtrhSzSXCqb4JAkvc0V3NbbWjW9gYP373rurSbFaei4JUHspd7o4Zra9wc7BxV5AwcmgrQ5HMmTRteSG2pQ3U3BUuavXrSbS8viBlAqksdXB6eZjLFzBWVzJfKKraZ2drkZSoQvqpjn9thQRBLwMtkyaZ1C5eCXrbNoVdN5kzwOK20tXg51TeRVgS4UmfgcoUSk2Ezi7Qjic8NTLKjwUOdR9lyb21lMcNTQU70TKRUnMX4gxEOXB7lTP8kNaV29rb48DitbKxx87v2kYz2t1JKzVcnK0N5NPT6a0GdpwiH1cy+9lHd6iyk+44FsLOhFF+JPelyMBlRCUe6xvAHw4oS7KSMFXT3FatPNusbD/BC+whj0yEOx+U90zfBjgZP2ud63eYqvnvX9ZpmrcIKUR69y1Clw/b6UgYmZunXuWVjsuLwyXDazFwY9Gck19h0iEYFG/yLQ372tY9it5gyztaEWFYqwFQwQtfINBurXTQqbHB1141ruPedu3VJbCx45ZkNRbimsoSdjR62ZVCDSwuafE6Odo/rWtVHLf5ghPXVmbdnT0dl+8YDrMugOk0irviDnOmfpHNkmi11brbVlyZswVLpsnPPm7byt7dcq9vNtaANBuFIlA9/+9BVTq+2Zm9CJ6TeVLvtWLMY9XBp2I/bYUkrMe9I5xjlJTbVJXN9TquiclIL6VWRKaqUuRYqlS47zeXFdAz7kUCjtwi71cQde7XJ20lGQc88n3nqLD87fXWV330do1ntcACxfkH9EwEuDE5hMcVK35ap3NArZXgqSGV8/6F0yxeRsDZJyV4lNJY5GUsjoqPR58xKx/HByQD72kcY9QcZ8wc52DlGtVu71vbJKFjl6Ruf4RvPdST82772EU1Dz5djxB/E47Ris5jYVOPmTP8kZpOg1uNgQ7WL7fWl7Gz0aH5dl93C3hYfrjRMyefjCp4MIWKp5nuavexp9rKtrhSTiC3XhibTm7GqM2wLmS6hqJzfHynZm2VKwS7b/uuXl1I6CQ93jrK1rpTjOrXbWMjwVJAdDR76xmfmfRYvhrPE7ryt8WDGrXWlDE7OLkkXUEP/xAy1nvRqKoz4g+xs8HAuvvxamFpdbDOztrJkSSWa8hIbW+tK+UWaMWFaFXBXw3INw7SgIJWn88o031omMjYq4fzAJNvrS7k45GcqwxCTZLiLLKyvdC1r/j3QMUpLeTEne8fZ0eDRRHn6xgPUedL/kkhiStPa5GUmFKHIasYkBO3D/oQO2uGpYNolmjbXurmSww7k6STaqaXglKd3bIa//N4RRRat2XB0PhCzymWnwefEZBIIIBSJ/S2Tbsp7W3wc7xlX7DfpGvGzo8GzpExtJqgpnTSXKpGOv2duL7lPQdWaBl8Rl4f9rKksyWqlo4V890AXN65P3ecpUwpOeTpHplXV6lpYP2yOLXVu1U1vK132tMsfXVNZoqniQKytRjpYTOKq0Jf0SKyoVrNgZ6OXkz3jWMwmwpEoU8EIDh0LOi7HT4718ZFXTbCxWvsU9zkKzmCwvd7DH17fpEk821zFSTVd5oqsJnwpLGqJ5HNr4DBcTLrFCDfXpteyfSGneydwJnA27mr0sq99BLvFhNdppW88dpPKdeDHvy1T7D1TCm7mKbKZ+eRtW3hbawMff+xERn16Dl6OjVXT4/TyyAxOayzX5lDnKC6HlUhUMj4Toq3FR9/4DIFQ9Ko8mOAyFWiKrCaaypzYLGbsFhNCCKJRSVRKpgJhpoMRfE4bDpuZEX+AC4P+tIM8M3EYSsmSqOc1FcXs74jNwCPTIUYWGAlO9U3Q2uSd/wyGpgLMaFz7LRVPnRzgePc4WzUsN7UQoaQmQLy72/3AFmJz93uBs8B3gWagA3iblHJ00bgG4JtANbFWMfdJKT+/3PVaW1vlgQMHlpUrGpV872AXY9Mh7BYT9zxxRnUFTqXr+US47GYmAxFsZsHmulIOd47hspuJSvCV2KhyOZgJRTjZO4G7yMKGKheDkwGkjLVvtFtNSClpH55Oy7DR4CuiprQobblbm70cULH0vaayhAsLalVbTIJ6b5Hi6O71VSUJM3H15JUbKvj6e9pUjxdCHJRSJupNpVh5HgB+LaW8P94dzgn8LTAipfxnIcTdgFdK+b8XjasBaqSUh4QQLmJdsd8opUzedQjlyrOYP/7mAX56aiDtcQB2i8BmMWvSXBZiLVDWV7noHZ+ha2SGBm8RjT4nz1+6gpbRO7sbvRzsTF8R1EYazOXVLH6shF2NHi4O+XWpjpqKRz/wEtV+v1TKs+ycL4RwAzcCXwWQUgbjzXlvAx6IH/YA8MbFY6WUfVLKQ/HHk8BpoC79t6CMTKKrA2Gpaf200ekQUkrWlJdwTWUJXaMz/PbiFc2dt+n061lIuu1H5nihfYTtDaW0lBdzoCO9Ge9Q55iqAiRq2FpXOr8/02vvo2TBvAYYAr4uhDgshLhfCFEMVEkp+yCmJEBlqpMIIZqJtVh8Icnf7xJCHBBCHBgaUlegIdOSsKd6x9MqAbUc+zpG6RufYcT/4r5Hy2jrKpddtSk4k+pDZ/smEUKiplHF8FQw7eS6dNnZ6OF4zziBUIRra1wZpUakQsknaCHWOv5eKeVOwA/cnc5FhBAlwKPAn0kpE9qGpZT3SSlbpZStFRXK7POBcITLV/wc6hzl/Q8cUF05f47JQITByfRz/FPhddrm80hcdjNdGjbkqlWYnJaIYCSKVaXxYG1lCZeG1GWxtg/7ubZG+zJWc1jNgr74bByRcLpvkrISfRymSpSnG+iWUs7NGI8QU6aB+J5mbm8zmGiwEMJKTHEelFJ+P3ORX8RuMdNUVsyuRi9v2V1HIJR5zefJ2TCBUJRKV/of+JzzceF3cngqMN/BbjoY0TSvJBNT8ImeCbY1lOJyWNK+M2fqJugenaam1J5R/bVk7Gr0XtUSBmBvS5nm1wEFyiOl7Ae6hBAb4i/dBJwCfgTcGX/tTuCHi8eKWPOTrwKnpZSf00TiJNy8pYYH3tuWVpBkMuYsYZsU3CHNIuY7afQVUVZsjz92sqaimEqX/araZBEJW2u121dFMqie6iu2YreYmQ1F0ipoWO8t4rhKx/IcfeMB+sYD8wG0WuFxWjmRIJZRryh7pd+0jwAPxi1tl4D3EFO8h4UQ7wM6gbcCCCFqgfullLcANwDvBo4LIY7Ez/W3UsontHsLL3L92jK+94HrufNr+zKOHRuaCmCzLP3H1nmKqPU4mAlGcFjNhMJR7DYzp/smiMrlG28tLCGVKVGVZrvdjR4Odo7x3MUrcZmUWxhrSrVrolznKcog2mEp6ypLlkSfbKhypXRmZ4Ii5ZFSHgESmetuSnBsL3BL/PFv0KvEfhI2Vrt59AMv4c6v7VNdkXIOsUj0Pc1ejnWPJbRwtTX72H95ZNlC4pevZCbTQtJ1afmKbayrLFmyN1SasOYrtnFcQWUfpXSPTLOpxp1RX9I5Gn3OhL6rvWv0y+0quPAcJdR7nfznO3dnnALcPTbD9WvLaGv2UeWys79jNGmfmn0dI2yqcS/bA6e8xKZZqax0IiKKrGY8RdaERhWln1NLWfF8q3gtmApGONM/MZ+ukQluhyVh5J1e+x1YocoDsKHaxdN/fiP3/2GrqpJFjT4nzWVOnr94hX0dI0uCShPhtJmXLTN1eWSG1iZt7oZKTfMmEVOQRJ0TIFbUcH1VScpiHRuqXaqcscsRlbF0jbZmr2oDyJZaNyd6E89eemYVr1jlgVhHhVdvquLd1zelNa7R50RKmVZRwbYWH/s7RhXVBjtweUSTdXhEoaNld5N3PkkvEecGpjg3MMWGKhdCxCLGaz1XO1H1aDGykH0do2ypdScMPE2FScBEkr5DayuKqVBhNVV8bd3OnEek8w/ZVl9K9+h02v4YpV9kiN1tF3851RBSMPPsjSu1EvZ1jOB12hiaDNA7NovPaWVno4c9zT46ksxaWnK8ZyJtH9DuJi+dSUoH712j35INCjCqWg1v2d3AA89dTrox3dXoieX1C8H+9hFVbRsPdipPFoPMK/0DBMNRttWXYjHFOkovTjnf0eBJ23G8UK6R6RAj8fyjmlKHppbCZKRjui62ma8KVF3MXp0LwawK5RmaDHBxaOmH7LCa2FpXqiq5LhHHu8fwOq3L9gwqtpk1CVGJxLulQWz5Mheo6XVaWVflYn+GERcLqXTZ6dO5Eo5ZQP+4chfDlrrSlDeH64yZJ3OePTu4JA+lyGbGahaaKQ7ATCiKr9iOxWxiKIWBwR+MaDLzLEzBjko42jWG1SwYnQ6pTq9IxvBUYH5JNWfC9wfDXFZZbD4RG6rdigNdq90ODqVII28uc1KlMvhVKatCef5gTwNjMyE+/eQZNlbHNqWhcDTlJlotPWMzXFNZzKg/mLLOgs2c+XZzcWduLc3Ii+kZm6Vn7OqZZ0eDR1PlmVtWr68qweu0zX9+8+9SvBg5PxuKpHRI62minmNVKI8Qgj99+Vru2NuI22Hl3MAEr/23X+t2vQuDfvY0e1POag0aFASsKXVoMoOpRYsbQCK0SJjT0zk6x6qwts3hdsT8GOur3Lz/pS3cuL5Ct94+qZZtAAcvj7Atw/RgvUP7lyV/6usvQW9LG6ySmScRH3/DJgBO903wDz88oeneB6DjyjSba92cTOK8C0djdeV2NHhU12EIhHPbDe/w5VFqSu24HFZcdqsuTlQ11HuLFPcSyoRVNfMk4toaNw//yfV88rbNmld7We50M6EoR7rG2KWyFG+yFoTZIhSV9I0H4k7WzOPTtCIb+x0wlAeI7Yn+8PpmPv/2nWn3vUnFid4JRefzq/SfZJCRsKLJxn4HDOW5ilu31/Lld+3WVIGUUOpMP03Yaha6tu8oZPR2js6xavc8yXj1pipeubFSdRWexexo8Cxb1rZ9yM/6qhJMQiBEzOEphAmrWXBhYCphX8/yEv2dloVIlduuuGtcphjKk4A3bKvRTHkGJ5f/gg9NBRiaSmyd21hdwqXh6SUdIdwOa34pT67Lg8bZWleKyJIsxrItAa++tgqHVZuPpnNkhtYMyk2d6Z9iY7VriUm9yJq7OtD5zFqN2zimwlCeBBTbLbxqY8pKWmlxqHM0I5/Ose5xdjV6KS+x4Su24XFacdjy6183NRtmY7V+VXGUkk49hkzJr/9AHvG+l67RrARVVMb8SZl8uQ5cHmV4KsiIP8jYdCjvLG2SWGhSrmdErRsIp8JQniTsbvJy/x/u0ex8oYike2SaFhVZrYnIpK+QXkzOhilZJg1db4xlW57w9Ml+Tc83FYwwNh1rwZgpagva6409ix3BF1Nb6shqyJKhPCn4g7YGzc85Oh3iSNdYRrn1Lrs5b308Vp2CRZVwTVV291yG8qSg0uXQbRN8sGNEVX69xQR1XqeqDgfZwKZzrYNUZHO/A4byLMvLN+jT1zIiUVXVZ1u9R9NCgVpjyeHMYyhPnvGya/RrCnuydzwt61SsC92YbvJoQbZDmxZyjaE8+UVrszejdhypmApE2FqnrHb17qb0i3nkgmx59xNhKE+e4bCaaWvWL9Cwa3Rm2ciWTTVujnZpnzKuBW0tPrbHmyIX2cwZ9YjNhAqXHY9Tn5rUyTCURwFaRhsspm98lq11yaMPGnxFdI1Mp6yHkCvKim3sax/haPc4HVems9qsdzHZ3u+AoTyKePW1VbqeP5nD0+e0EgxHE0ZV5wNrNOyilynZXrKBoTyKaCxz6npnO9k7sSSM3m4R+ErsGbdK0RMlpYWzRd7OPEIIjxDiESHEGSHEaSHE9UIInxDip0KI8/HfCUOHhRBfE0IMCiFOaCt6drlJ59mnyn21z2djtTtlNcxcU2Izc6Yvf0zm12QxIHQOpTPP54EnpZQbge3EulrfDTwjpVwHPEPyPqXfAG7OUM6cU6xhO8REHO8enw8taWv2cVTDPjh6sL7alVf7sHVZ6rK9kGUDgRa0kv8jiLWSB4JCiNuAV8QPewB4Fvjfi8dLKX8V74Rd0CRrz6EVs+EobfUeIFZwPd/Jp6huj9NKmU7d31KRtVbyy6FFK3m9mA1F+PX5Yd2v0zs2kzflm1JhNQnODuTPkm1dZUlO/EtZaSWvBDWt5LPFd/d3MZwkTVpLusdm2Kxhw1+92FjjYjqHZunF5GK/A1loJV/oBMIRvvzLi1m73uJaBfmII89SwHNhpgadW8mvBPa1j2S10MaZ/knW5pH/JBGJ2rXkkppSfbshJEOptW2ulfwxYAdwD/DPwGuEEOeB18SfI4SoFULMt4oXQjwEPA9sEEJ0x1vPFwx7W8qWmJH1xp2iN2iuWV9Vwog/df+hbOPKUfaqrq3k48/foVa4fMBmMfG+l7ZwzxNnsnbNo11jVLrsDCpoIpxtsh0/poS5Av7ZxogwUMA72hpxZTG9NyqhuSw/l275mMGaq5nHUB4FuBxW7tjbmNVrnuwdT7sztN7UeRx0p9noOBu4jJknv3nPDS26t1NfiD8YYUttZv17tKbOk50ytulizDx5TnWpg9t21GX1mh1XpnRrvqWGyUB+GQogtifNlencUJ40+OOXrcnq9QYng2yPh+zkmhK7WZN2h1qTK2MBGMqTFhuqXbxCp4IgyZiYzY+7/boqV14WWnTnsMiioTxp8ic3rs3q9S4O+fOiBnQerR6vIlf7HTCUJ22uW+PLuBFvuuhVgEQpJgHn83DJBrl1KBvKkyZCCO66Mbt7n2Pd41lpUJuM9VWuvE0FN2aeAuPmzdVZ6z42R60nN/FbkNtN+XIYBoMCw2I28f6XtWT1mke7xijN0RKlJw+jCuYwZp4C5K27G/CqaMSrlmBE5sRwUO8pynPlMWaegqPIZubd1zdn9ZrnByazGuUAUOvN3V5LCYapukC58/omPFmcfUamQ1lxmlrNgj3NXmo9DvbleYlfY+YpUMpK7Hz69m1ZvaaS7tpqcTks7G3xUWy3sL9jlN6xPOq2nQTDVF3AvG5zNe+6LnsR150j2tc5qHLZ2dPsJRyRvNA+wth0fkQ1KMEwGBQ4H3/9JtZnsW6YRJswmTXlxexs9DA0FWB/x2heVQBVimGqLnAcVjNffMeurPXjPNU7SXOZej/Tllo3m2rcXBr2c7hzjDwMWVOMMfOsADZUu/j4GzZl7XrpFvkzCdjV6KGl3MmJ3glO9U3oJFl2MWaeFcK79jby2k361rSe40jXmCIFKrGZ2dPspcJl51DnGO3D01mQLnvksnW9oTwaIoTg07dvo9qtfyhNRMLaiuT7rA1VLnY3eQlLyf6O0bzutqCWErsFcw6zBQ3l0RhvsY3P/cH2Zbu9acHpvqt7mnqcVtqavTT6nJwdmOTg5VFmQ/lfRFEtuXSQgqE8uvCSteV86BXX6H6dyXhP0211pexoKGVqNsy+jlE6R1bW0iwZuXSQgsK6bQbp87FXr+O3F4c5rHP36s4RPwOTwbzqWpAtcmlpA2Pm0Q2r2cQX3r5T93pv/RNBtqXoabqSMZRnBdPgc/J/37RF9+us5H1NKnK9bDOUR2du21HH7bvqdb3G2YHJnPTkzDXGzLMK+MRtmzOKCFBCSRbLAecLxsyzCiixW/jCO3bqmotzpGssK/6lfMKYeVYJ2+o9/K/XbVj+QJVIoNGX34lrWmP4eVYR73/pGm64pky385/oGV9Vy7eCWLYJITxCiEeEEGeEEKeFENcLIXxCiJ8KIc7Hf3uTjL1ZCHFWCHFBCKF5L9NCwmQSfPiV63Q7/3QoyqYC6GmqFbm+USideT4PPCml3AhsB04Ta+r7jJRyHfAMCZr8CiHMwJeA3wM2Ae8QQmQv9DgPuW6NjzU6tk3sGPbnNN4rm+T9nkcI4QZuBL4KIKUMSinHgNuAB+KHPQC8McHwNuCClPKSlDIIfCc+btUihOCONv0yTwcnA2zPckXTXFEIy7Y1wBDwdSHEYSHE/UKIYqBKStkHEP9dmWBsHdC14Hl3/LUlCCHuEkIcEEIcGBoaSutNFBpv3lWPzazfdrOQ0qgzIe9nHmLxb7uAe6WUOwE/CZZoSUi0fkgYhSWlvE9K2SqlbK2oyG4ngmzjK7bxe1urdTv/pWE/19bkvji83uS6kqkS5ekGuqWUL8SfP0JMmQaEEDUA8d+DScY2LHheD/SqF3floOfSDWKxdSudXCbCgQLlkVL2A11CiDknxU3AKeBHwJ3x1+4Efphg+H5gnRCiRQhhA94eH7fqaWvxsVZHw8Gx7nEa8rxgYSYU28w5N4wovT19BHhQCHEM2AHcA/wz8BohxHngNfHnCCFqhRBPAEgpw8CHgaeIWegellKe1PQdFChCCN6h8+xTtYIjDnJtLACF+TxSyiNAa4I/3ZTg2F7glgXPnwCeUCnfiuYtu+v5zFNnCYb1iYo+2j2G12lldAUaEHJtLAAjwiCneJw2Xr+1RrfzhyKS9VUr03BgKI8Bd+zVd+l2pn8ia/Xkskk+LNtW3qdaYLQ2eXXNxRmfCbN1BTpNjZnHICuGg/4CKNieLsbMYwDA7bvqdV1adY/NrLiQnVynI4ChPHlBqdPK67fpZzgACOhk0csVuY6oBkN58oZ36m44mNTVKZttjD2PwTy7Gr26tynJdSyYlhh7HoN59E5VgJjTtNJl1/Ua2cKYeQyu4k276nFY9fuXRCW0lK+MpZsx8xhcRWmRlTdsq9X1Gid7xnHazMsfmOcYM4/BEvT2+UwFI2xZAeV582H/ZihPnrGr0cPGan3j0Tqv+Cn0MgfGzGOwBCGE7vFu/RMBttd7dL2G3uQ6EQ4M5clL3rizTlfDAYA/GNb1/HpSZDXnRaZs7iUwWILbYeX3dTYcnBuYKtji8PmwZANDefIWvZduACX2wrO6Vbjs3Lg+PwrE5IcKGyxhR4OHa2vcnNax5fuRrnGq3Q76J/Iv6tokYg2LN9W6ubbGzaaa2O+KPHLyGsqTp8wZDv7+sRO6XSNWHN6Zc+UpsVu4tsY1ryCbat2sr3LhsOb3zGgoTx5z245a7vnJaWZCEd2ucaJnjGKbGX9Qv2sspM5TNK8gm2pcbKoppd5bhKkAbeeG8uQxboeVW7fX8t0DXcsfrJLpUJS9LT5eaB/R9LxWs2BdpSuuJC8uvUqduXduaoWhPHnOHXsbdVUegPbhmNM0qrKjttdpnVeOuT3K2ooSbCuwdsJCDOXJc7bVl7K51s3JXv0MB4OTAXY2epZtey8ENJcVz+9P5hSl2u1AiMJbdmWKoTx5zlyNg4/raDgAmJy5urabw2piY7V7wf7EzcZqF8V5kMGZLxifRAFw245a7nniNNM6beorXXbqfU7esL123jzcXFac83K2+Y6hPAWAy2Hlth21PLRPu72Px2nlLbvqecfeRtZWFGakQa4xlKdAuKOtSRPl2d3k5Z17G7lla03e+1HyHUN5CoSt9aVsqXNzoid9w0GJ3cKbdtZxx95Grq1ZPT1L9cZQngLijrYm/vYHxxUfv7nWzbuua+LW7bXGRl8HjE+0gLh1Ry2f+smplNEADquJW7fX8s69TWyrL12VJuRsoUh5hBAdwCQQAcJSylYhxHbgy0AJ0AG8U0q5ZE0hhPgY8MfEWix+RUr575pIvgopsVu4bWcd336hc8nf1lWW8M69jbxpVz2lRSvHi5/PpDPzvFJKObzg+f3AX0kpfymEeC/wv4C/XzhACLGFmOK0AUHgSSHET6SU5zOUe9VyR1vjvPJYzYLf21LDu65rYk+z15hlskwmy7YNwK/ij39KrPvb3y865lrgd1LKaQAhxC+BNwGfyeC6q5otdaW8fmsNW+tLecvuespL8idEf7WhVHkk8LQQQgL/JaW8DzgB3EqsF+lbubpx7xwngE8JIcqAGWId4w4kuoAQ4i7gLoDGRv0TwQqZL71zV65FMEB5JukNUspdwO8BHxJC3Ai8N/74IOAitiy7CinlaeDTxGamJ4GjQMLk+dXUSt5gZaBIeeJ9RpFSDgI/ANqklGeklK+VUu4GHgIuJhn7VSnlLinljcAIYOx3DFYEyyqPEKJYCOGaewy8FjghhKiMv2YCPk7M8pZo/NxxjcCbiSmagUHBo2TmqQJ+I4Q4CuwDfiKlfBJ4hxDiHHAG6AW+Dle3ko/zqBDiFPA48CEp5aim78DAIEcIKVVmQOlIa2urPHAgoV3BwCCrCCEOSilbE/1tZaf6GRjoiKE8BgYqMZTHwEAlhvIYGKgkLw0GQogh4HKWLlcODC97VG7IZ9kgv+XTSrYmKWVCr31eKk82EUIcSGZNyTX5LBvkt3zZkM1YthkYqMRQHgMDlRjKA/flWoAU5LNskN/y6S7bqt/zGBioxZh5DAxUYiiPgYFKVqTyCCEahBC/EEKcFkKcjBchWfj3vxJCSCFEeZLxHiHEI0KIM/FzXJ9n8v15fNwJIcRDQgiH3rIJIf5JCNEjhDgS/7klyfibhRBnhRAXhBB3ayVXprIt95mrQkq54n6AGmBX/LELOAdsij9vIFZv4TJQnmT8A8D7449tgCdf5APqgHagKP78YeCP9JYN+CdiBV9SjTUTS4pcE//cjs69rzyQLelnrvZnRc48Uso+KeWh+ONJ4DSxLx3AvwF/TawuwxKEEG7gRuCr8fFBKeVYvsgXxwIUCSEsgJNYPlU2ZFuONuCClPKSlDIIfAe4LR9ky/B9JWRFKs9ChBDNwE7gBSHErUCPlPJoiiFrgCHg60KIw0KI++MZtHkhn5SyB/gs0An0AeNSyqf1li3+0oeFEMeEEF8TQngTDKkDFhbU7ibDL6iGsqUaq4oVrTxCiBLgUeDPiBUe+TvgH5YZZgF2AfdKKXcCfkDTtXsm8sW/GLcBLUAtUCyEeJeesslYMct7gbXADmJK+6+JhiV4TXNfiErZko1VzYpVHiGEldiH9KCU8vvEPtwW4Gi8Amo9cEgIUb1oaDfQLaWcuys9QkyZ8kW+VwPtUsohKWUI+D7wEp1lQ0o5IKWMSCmjwFeILdEW083VJcjq0XBJmaFsCcdmwopUHhErnflV4LSU8nMAUsrjUspKKWWzlLKZ2D96l5Syf+HY+PMuIcSG+Es3AafyRT5iy7XrhBDO+HluIrZ+1022+Os1Cw57E7GafIvZD6wTQrQIIWzA24Ef5YNsycZmhFaWkHz6AV5KbLlwDDgS/7ll0TEdxK1ZxJY/Tyz42w5ixRmPAY8B3jyT7xPECq+cAP4bsOstW/w6x+Ov/wioSSLbLcQsWReBv8vG56ZENiWfebo/RniOgYFKVuSyzcAgGxjKY2CgEkN5DAxUYiiPgYFKDOUxMFCJoTwGBioxlMfAQCX/H5/YlfYKffjvAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "polys.plot()" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "We are interested in an area that is called ``Eteläinen`` (*'Southern'* in English).\n", "\n", "Let's select the ``Eteläinen`` district and see where it is located on a map:\n" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "# Select data \n", "southern = polys.loc[polys['Name']=='Eteläinen']" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "# Reset index for the selection\n", "southern.reset_index(drop=True, inplace=True)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
NameDescriptiongeometry
0EteläinenPOLYGON Z ((24.78277 60.09997 0.00000, 24.8197...
\n", "
" ], "text/plain": [ " Name Description geometry\n", "0 Eteläinen POLYGON Z ((24.78277 60.09997 0.00000, 24.8197..." ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Check the selction\n", "southern.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Let's create a map which shows the location of the selected district, and let's also plot the geocoded address points on top of the map:" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "deletable": true, "editable": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOYAAAEYCAYAAABIhYpmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+wklEQVR4nO29aXRb15Xn+9uYSRAcRBLgAJKSLMmyZVmiRDuyZct2FDu2Kim7otix064oncrLWpXqrHRS/Tp5Ve+91f2ha1W/rkpVVq0kFcfuRHEGl60odjzEkuNYlgfJluRB1mjNnOdBJMAJwHkfSKA4ACCGiwuAvL+1uEgA995zQOB/9zn77LO3KKUwMDDILUzZ7oCBgcF8DGEaGOQghjANDHIQQ5gGBjmIIUwDgxzEku0ORKOiokItX748290wMMg4x44d61VKVc59PieFuXz5co4ePZrtbhgYZBwRuRLteWMoa2CQgxjCNDDIQQxhGhjkIIYwDQxyEEOYBgY5iCFMA4McxBCmgUEOkpAwRaRURPaIyBkROS0it4jIMhF5RUTOTf8ui3O+WUTeF5EXtOu6gcHiJVGL+X3gZaXUWmADcBr4LvCqUmo18Or041h8c/ocAwODBFhQmCJSDGwDngBQSk0opQaB+4Hd04ftBh6Icb4X+BPg8fS7a6AlgUAAv9+f7W4YRCERi7kS6AF+Oj0cfVxEnIBHKdUBMP3bHeP8fwb+KxCK14iIfE1EjorI0Z6enoTfgEFyhEIhLl26xHPPPce//Mu/4PP5st0lgygkEitrATYB31BKvSMi3yf+sDWCiHwG6FZKHRORO+Mdq5R6DHgMoKmpych3oiFKKTo7O/noo484ceIEw8PDANjtdoqKirLcO4NoJCLMVqBVKfXO9OM9TAmzS0SqlVIdIlINdEc5dyvwpyKyA3AAxSLyC6XUo1p03iA+IyMjvPfee3z00Uf09vbOe318fBy/309BQUEWemcQjwWHskqpTqBFRK6dfmo7cAr4HbBr+rldwHNRzv2/lFJepdRy4GHgj4Yo9SEQCPCrX/2K1157Laoowzz99NOMj4/r2DODREjUK/sN4JcichzYCPwd8PfA3SJyDrh7+jEiUiMiL2WgrwZJ8PLLL9PR0bHgcd3d3fziF78gGAzq0CuDREloP6ZS6gOgKcpL26Mc2w7siPL8AeBAUr0zSImPPvqIY8eOJXx8a2srXV1d1NTUZLBXBslgRP4sMnp6enj++eeTPs+wmLmFIcxFxMTEBE8//TSTk5NJnWe323G7Y612GWQDQ5iLBKUUL7zwQlxHTyzGx8c5efJkBnplkCqGMBcJ4WWRVHnxxRdpaWnRsEcG6WAIcxHQ0dHB73//+7SuEQqFeOaZZxgZGdGoVwbpYAgzzxkbG+OZZ57RxHkzPDys2bX0IBSKG+WZ1xjCzGOUUjz33HMMDAxods3m5mb27dun2fW04OzZs7zxxhvMrEw3NjbGL37xC954440s9ixz5GReWYPEOHz4MGfOnNH8ukeOHMHtdtPUFG3pOnOEQiHeeOMNrrnmGkKhEPX19fj9fvbt28fAwADd3d3ce++9vPfee5w8eZKuri7Gxsa49dZbMZvNuvY100gu1sdsampSRsLn+LS0tPCzn/0sY8M5EeELX/gC11577cIHa4DP52PPnj1cvnwZk8lEKBTi7rvv5ujRo7NGBGazed5Qe8uWLXz605/WpZ9aIyLHlFLz7oCGxcxD/H4/e/bsyegcSynFb37zG3bt2kVtbW3G2gljs9kiW9DC7+uVV16Zd1y0+W9vby/vvvsuDocDm82G3W7HZrPhdDopLS3NaL8zhWExNcLv9+NwODCZMj9t/+Uvf8n58+cz3g6A0+nk61//OoWFhRlvq6enh5/85CdJB0jEwuVy8e1vf1uTa2WKWBbTcP6kQSgU4sqVK+zbt48f/vCH/PrXv854m36/n4sXL2a8nTA+n083Z1BlZSX33XefZtcbHh4mEAhodj09MYayCaCUYnR0lMLCQgKBABcvXuTMmTOcPXt2VmqO8+fPc+jQIW6++eaMOCMuXbrE+++/H5mD6cXx48e54YYbWL16dcbb2rhxIxcvXuTEiRNpX0tEeOqppzCbzRQXF+PxeKivrycUChEMBlFKUVtbi4ho0HNtMYayCXDy5En27NlDaWkpfr+fiYmJuMcXFRVhs9moqanhhhtu0MSBMjIywo9//OOsBQAUFxfz9a9/HbvdnvG2xsbG+PGPf8zg4GBG27njjju48847M9rGQhhD2RQJBoO8+uqrAAwODi4oSpgSUX9/PydOnOCpp57irbfeSqsPoVCIvXv3ZjUq5+rVq1GdMZnA4XCwc+fOjM7X16xZwx133JGx66eLIcwFmOuuT4XXXnuNw4cPpxRgDnDw4EEuXbqUVh+04NixY1y+fFmXtrxeL3fddVdGrl1eXs6f/dmf5eQQNowhzDiMj49z8ODBtK8TDAbZt28fP/jBD/jNb35DMlkAz507x+uvv552H7Ti+eef18xruhBbt25lxYoVml7TZrPxhS98AYfDoel1tcYQZhxOnz6ted7VEydO8MMf/pBnnnmG7u5o+cv+nf7+fvbu3atp++nS39/PgQMHdGlLRLj33ns1veZ9991HZeW8yuo5h+GVjYMWnsFYnDp1ilOnTrF27Vpuv/32eWk9JiYm+Ld/+zfGxsYy1odUOXToEGvXrqWuri6j7SilNF+vPXToEKdPn2bbtm1MTk7i8XhyMkugIcwY6LVeeObMGc6cOcOqVavYtm0bdXV1KKX43e9+t6BFzRZKKfbu3cvXvva1jH2plVLs37+fw4cPU1ZWhsvlYmxsjKGhobSy+nV3d9Pd3c3HH38MwCOPPMKaNWu06rZmGMKMwenTp9FzKen8+fOcP3+e5cuXU19fn/MZBQYHB3n22Wd5+OGHNXeiXL16lddff53Ozk5qamoYHR2lubk58rrL5cJisSTklBMRVq5cyYULF6K+XlFRoVm/tcQQZgxOnTqVlXYvX76sm3MlXT7++GPefPNNbr/9dk2uNzY2xltvvcXhw4fjRuwMDw9TVhazuNwszGYzjzzyCE8++SRXrlyZ9ZrFYsFms6XV50xhCDMKPp8vq8sT+RRGdvbsWUpKSrjxxhtTvkYgEODIkSO88cYbjI6OLnh8QUFBwktYtbW1mM1mGhoaZgnT5XLx2c9+NmdLRBjCjMKpU6d0HcbOJV8spsfjoa2tjWeffRar1cp1112X1PlKKY4fP85rr73G0NBQwueJCF6vl87Ozqg3seuvv57CwkKsVisbNmygu7ubM2fO4PV6sdvtOBwO7rrrLsrLy5Pqr54YwpzDxMRE1nfF50vJgr6+PhwOB2NjY+zZs4eHHnqI6upqJiYmGB8fZ3JyklAohFIq8iMiBAIBRkZGOHr0KF1dXUm36/f78fv9WK3WiGe4v78/IjSbzcY999yD1WqNnPOXf/mX2rxpnTCEOYdDhw5FqmFli1xYInG73QQCASYnJzGbzbhcLgKBwKyyC9XV1ZHMeqFQiKeeemrB69bX189y5KTD5ORkpP2CgoLIdZubm+np6eGLX/yiLtvVMoEhzBn4fD7efvvtrPbBarXmxFA2EAjQ398feTw4OIjFYsHtdkccJrmU7nLu3LStrY0nnniCRx99NGFHUS6Rl5E/mQrmPnjwYEJB6pkkV0LFojlFAoEA3d3dtLa20tramvQ1PR4P7e3tWnQvIfr7+3n66adT6mu2yUthHjhwgJ/85CccOXJEs2FfR0cHubDVLLytqqKiAq/Xy4oVK7ISbJ2JEvAOhyPirJm5fmg2mzNi1axWK11dXezevTvn14XnkpdDWbfbzbFjx2hvb2f//v2sW7eOzZs34/V6k/4Sj4+Pc+DAAd59992cyFM6Pj5ORUVFZCdKdXU1FRUVjI2N6Tr37e3tjZr4Kh1merrDyx02m41QKMTY2JjmG8ArKiro6OggEAiwZ88eBgYG2Lp1a07vKgmTt8IMEwgE+PDDD/nwww+prKyksbGRDRs2LDjpV0px4sQJ9u/fn1PZx4eHh2cJ0Gaz0dHRgc1mo7a2lra2Nl36UV5eTl9fX0aubbPZIlOG8O9AIIDT6cTn81FSUsL4+Hjao6G5+zlfffVVurq62L59e84n6cp7Yc6kp6eH/fv384c//IHrrruOxsZGVq5cOe8O2dfXx4svvpgTexwXIuwImpiYoK2tjbq6Ojo7OzPuIMpEDKyIUFVVhdVqjeo4WrZsGVarlcHBQSoqKnC73fT19UWy5yWD2+2OehM7ceIEp0+f5qabbuL222/PWa9tXgqzsLAwcneNRigU4uTJk5w8eZLS0lIaGxtpbGykoKCAN998kzfffDNvygDMjXBpaWmhpKQEi8WSEYtmNpuprq7OmMOks7Mz5mszxRoeyptMJjweT9LrnfGGxMFgkMOHD/P+++9z6623smXLlpwLzcvbnD8///nPk7J44WDmS5cu5cRcMhmKiormDbdFJLImqNVnKCKUl5ennGlhIcxmM3a7PWnHUkFBQdxQvZqaGrq6uiI3W6/Xm9SNxel0sm3bNjZv3qx7RvdFl/Mn2UKrSikuXLhAZWWlLgmltKSkpGTec0oprly5Qnl5uWYeTa/XmzFRwpSlSmU3R6xhtcPhwOPx0NnZicfjAaYsbDLhfTC1fv373/+eH/zgB0lll8gkS0aYYbq6unA4HBQXF2vco8wR7y7e29vL4OAgDQ0NaXsb9fD6phJu2N/fT0NDQ+RxTU1NZH9mV1cXoVCI9vZ2qqqq8Hq9Kb+PgYEBXnzxxazGSYdJSJgiUioie0TkjIicFpFbRGSZiLwiIuemf8+7bYuIQ0TeFZEPReSkiPx3rTqeTnqIoaEhJicn8yLFBLBg0EPYeqaTUaCmpibj6SKBlOdyLS0tEatosViiii8YDKYdwHDlypW0CgBrRaIW8/vAy0qptcAG4DTwXeBVpdRq4NXpx3MZBz6plNoAbATuFZEtafea9IQJUyFcfX19eL1eLbqTUWaGxsUjkS1T0TCZTLpFPKW6BBIKhejq6or7eZnNZk22zO3fvz/r8coLemVFpBjYBnwZQCk1AUyIyP3AndOH7QYOAN+Zea6aGhOEvRbW6R9Nxgnh4ejVq1dTvkYoFKK1tZXa2lqUUpGgbYfDgdVqjTgTwkMbpdSsIG69mJiYoKSkZMG5U09PT9Lrj2azmcrKyrjeUq3QwrEUz6ljsWizyODz+Xjttdc0LdeQLIm8k5VAD/BTEdkAHAO+CXiUUh0ASqkOEYk66RMR8/Q5q4AfKKXe0aTnTM0z0xFmmEQX7bPpUne5XAk5NQoLCxMWptVqpaysTBdRwpRwMjl/09KjeuTIETZu3Eh1dbVm10yGRIayFmAT8COlVCPgI/qwNSpKqaBSaiPgBW4WkRuiHSciXxORoyJyNFHPmN5zxJn7+/Qm0azk7e3tCd9AKisrdUv4VV5entLey2TQcjiulErYEdTd3T1vXTzdWONELGYr0DrD0u1hSphdIlI9bS2rgbifsFJqUEQOAPcC8/JCKqUeAx6DqXXMRDqfqmc2VbIpzEQjfYLBILW1tXH3PJpMJt2XjfRIEan1BvO2tjbef/99Nm3aFPMYpRTPP/88ExMT3HbbbXR2dnL+/Hl6enr49re/nXLqkgVvw0qpTqBFRMKVcbYDp4DfAbumn9sFPDf3XBGpFJHS6b8LgE8BmtUm19tiajWHSYVknBHRRDyz78XFxXR1dTExMUF9fT21tbXU19dn7P0VFBTost0rVedXPP7whz/EtX6HDh2itbWV7u5u9u7dy9tvv43H4+Gb3/xmWvmEEv0kvgH8UkRswEXgPzIl6qdF5C+AZuBBABGpAR5XSu0AqoHd0/NME/C0UuqFlHs7Bz3c+zPROypkJskE2vf392OxWCIeSrfbzfDwMCUlJQwMDET+b3Pn1jabTdMMA2HKy8szvidSRDIizNHRUa5cuTIvn5Hf7+fDDz/kj3/847xzLl26lPZcOiFhKqU+AOaFDTFlPece2w7smP77ONCYRv+i4vP5ePnllzOaKT0a2RKm3W5Papg2c+tYQUEBIyMjjI6OLvjFnZiY0Dw4Phx0Hw+Px4PdbqetrS3lGObCwsKUgt0TIVzqfnJyko8//pjjx49z/vz5mKGdIyMj7Nu3j4ceeijloI+8DGK/ePGi7tYSEnfAaI3T6Ux6/lRQUBDJ1ZOMg0frebRSap7YwzG5DoeDycnJiFOorKyM4eFhioqKKCkpmZcHNh4OhyMjwnS5XJEosYsXL7Jnz54Fz/F4PPzJn/xJWpFYeSnM9evXs379egYGBti3bx9nz57Vpd1sCTOVdCPDw8NUVFQk7QnV2mK2trbO2n8JUxuYo3neBwYGKC4uxm6309zcnFT+o0w55srKyhgcHMRut3PkyJEFjzeZTDzyyCNp56vN21hZmPqnJZvLNB9JxSkzODiY0vJEurVAozE3yD7elGBw8CoXLowQCqmk1hAztcYcLhf/5JNPxiyzMJNQKMQ776S/VJ+XFnMmehbeybftYsmSbiRVLPx+P8XFxRHxxJpzhkKwe/cuWlrqqKtr4W/+5g8Jt5GJwAWTyYTP50s6WkmLsMC8FubZs2eznm5SD/Ta1O1wOJIWpslkiuRa6u/vjxpc7nA4EtpO5fc7aWmpIxQy09JSx7lzg6xb14Df71/w/ETjiRci/H7C10xWlCKiSXnCvBXmwMAAzz77rK5tiggej4ehoSFdg5z1ysyezDzNZrNRVVXF0NBQZHklPD+c299EgwucTh91dS0Ri+l0+rhyxYfL5aK+vh6YGqKPjY3hdDrx+/2Mj49rtsRTVFREWVlZWtdSSvHb3/4Wn8/Hli2p79fIW2G+/PLLuu8AmJn6oqCgIDI8ExHN1/5mkqllgLksZJktFksk+/rExMS89zwyMjJvSFlSUpLwGqYI7Nq1G7/fidPpI+zUnJugDKaWdsLe23SnMxUVFTgcDlpbWzVJzGaz2dKe8+atMLO9LWfmumAmw81MJlNGFs6jEctiVlRUYDab6e7ujpt9PRQKzdvdkmjwfRiTCYqKErsRjYyMMDIykrLFrKmpiSSxThe3283NN9+M3W7nuuuuS3vNO2+FmQu7zMM4HI6MiScUCmUkGicaHR0dmM1mRAS73c6yZcvw+/1pbdXSY4kp2VQiFosFp9OpaZjgjTfeyObNmzW7Xl4vl+QKmQ4G18tiBgIBSkpKqK6uZnR0lNbW1qQz8fX19VFfXx+xGMmsi4ZCMDLiJNl7biq7SrT2Pq9cuVLT6xnC1IBMZ/bWMz+Rw+Ggo6ODUCiUcnaH5uZm7HY7tbW1CS8xhZdKvve9b/Gzn+0imZWpaMnK4pGJG6nWIwNDmDlOTU1NRtYWo1FfX097e3tkHS6dKCC/34/P50s4yGHuUonf70y4rWSF5vf7NS3eFB7+a0neCjOX5piZ6ovVauXq1av09PRkfItbQ0PDvHlsZ2dnWqFlg4ODCafWDC+VmEzByFJJoiT7/1dKpZRGMxZ333235nuD89b5k0tkSpjV1dURsaTq5XO73ZHiQLG8j06nM6ZzadmyZWktIbhcroTC/GItlWQKn8+XdGLoWLz55pts2LBBUyuctxZzsePxeGaJpbOzM6VCOOE5Y3d3N16vN6oFrKioiHlzSTcMMZnhcHipJFlRppqrtrW1NRK4kA4jIyP86le/0rQWTt4K8zOf+UzOJG3u7e2lurqa+vp6qqurKSsrW9AZEG8BetmyZVGtTLJOjurq6llLCa2trYyNjVFfX8+yZctwu91UVFTE9bymWx9loTXCVD2xYcxmc1rZ01tbWzUZhra0tHD48OG0rxMmb4eyHo+Hr371q/z617/OSkrJmQSDwXl9sFqtEQsnIphMpsiPiGC1WqPuNywuLmZsbCzqEkCiTqBwQEBHR8c8x0ggEJhliReqgTk6OkplZWXKX36n0xmz33OD1nft2k2yzk23253W5x8KhTTz0l66dIlAIKBJipa8FSZMzV++/OUvs3fvXt32ZCbK5ORk3MBqk8k0r1iQ0+lEKRUzx8zAwABlZWUx52wlJSUUFRXNShkSLWQtFCIyl0skQD6dyKbCwsKYwozmiY0X9VNVVUVvby9OpxOHw4HNZtNEVFqVhpicnGRoaIjy8vK0r5W3Q9kwNpuNhx56iJtvvjnbXUmKUCg0a2gaTjK90Jck2vC9oKCA+vp6rl69Oi+Pz9wvbirrhclG1oSx2WyzlkvmDltnemK93haUIu6QNpzHaGhoiK6uLlpaWjh//nxa80S73a7pGqRWWwPz2mKGMZlM3Hrrrbz77rvZ7kpStLW1UVNTQ39/P4WFhQltXZopEqvVSk1NDW1tbTG9qnPX15K1UuE2CwsLk86V6nK5InPUWMPWXbt24/M52bPn8/zTP30r5pA2Xlhic3Mz9fX19Pf3J+xBrq+vRylFV1dXWlvG7HY7t912G263G7fbnbQfIBaLQpiQ+eibTDE+Pp6wKGFqbbC8vByn00lPT8+CeXHmekWjba2C2cPbaP/KoqKipIU5Myt8rBuCyTS1VBLvZpFIrHD4dY/HE7HUsUL1ampqNIs93rFjBzfeeKMm15rJohFmtvLxpENhYSHBYDDpxGIDAwMJe0vnbhmLtl6YiBMmlZw6Q0NDuFwuhoeHY94QIPrNorKyEqUUFoslKRH9+9DZhMu1CqfTR1dX56zlIJPJpNnGgIsXLxrCjEe+WcyGhgb6+vpSyvaXzDzm6tWr8zyvc7dWLTS8LSwsTGnuFG67traWgYGBmAEE0W4Wk5OTC/5vYln5qRvNn0eE/pd/+Qxu99SyUFFREa2trRQWFib9fqKh1XXmsmiEmU8Ws7q6OqnUjOliMpniel9jWbOwqDo7O1NekggGg7S1teFyuXA6CzCZos9nZ94sTCYTLpcLq9Uac5kmnpWfe6Pp7RXGxqasY3h46/f7qa6u5urVqyltRBcR7rvvPm666aakz02ERSPMfLCYy5YtY3R0VPd114XSQIrAl760m97eSioreyLWZ6EaKLGIZsmGh4cTHj6GQqHIhuxY58Sz8vGGzTP/Dx0dHTQ0NDA2NpZUXiWbzcbnP/95Vq9enfA5ybKohCkiORXcPhO3243P59Ntb+VMFnLchELw85/Ptz49PT0UFBQk1ed4liyVz6a5uTlq9r544ksm7ra1tZVgMIjL5aKsrIyOjo64N7Hi4mIeeeQRqqqqkn4vybBohHnw4EGUUlitViwWC5WVlbrs+k8Uu92ua6rNmSwUXB3L+oyOjibtJIlnyVLJ9hcrtcpC4ks0RUm4T+G8Qlarlfr6eoaHh2cFcpSWllJSUsLOnTtxuVxJv49kyZ+JWRzefvvtSBrLyclJRkdHaW5upqamJss9mwqPCyewyhbRbggzF/vjbbnq7u5OamdLrGs1NDSklMrD4/FEtWALLe+kyuTkJM3NzQwMDFBSUoLX66Wuro7BwUFGRkZ0ESUsAovZ3d3NK6+8EvW1zs5O3G531iwVTN3xw3PKgoKCyL5KPa352NgYVqsVpRQ2m42REX9kuOn1tvDgg3tiWp+xsbEFt0fNFUmsTHepEG2ZJpHlHS2EOzQ0NCugI9G9pVqQ98J8/fXXY74WCoUYGhqKG1+aaWw2WyQutrS0NCLIurq6SMBzZ2cny5Yty2jF5ZKSEmw2G729vbOGm83NDXzve9+ivn7qCx7tSxxvW9VMkVxzTRff+c4+xsdHcTgcjI876e7uxuVypVwEKlrgxULLO1oEx0dDT2Hm9VC2p6eHU6dOxT1mfHw8koM0G7S1tUXmLDOtZEtLCx0dHbS0tDA5ORnZBubxeDJSKXtsbAyz2czExMSs4SaAUrPTecyNae3p6aG8vByz2Twv9namSC5c8HDmTB89PT20tLTQ3d1NaWkpZWVlKa2DlpeXRw2xWyjbQTppSuJhWMwEiTWEnYvP56O0tBSn06lb8mSY+iADgQCXL19e8NjW1laqqqro7u6O1GPUkpGRkcgXKzzc9PmcPPPM52lt/XfPZixrU1RURF9fH7W1tfT09OB0OrHZbJjNFurr22huro0qknClrFQoKyuLGuG0kOMnnsc2HfTc/5uXwlRKceLECc6dO5fwOeGhlNPpjGRQ9/l82O12QqEQJpNp3s6MVCkuLqa4uDiptBVKKTo7O6mtrc3okDaMyQQul48vf3n2F9zniz5MDHsvw8m6ZnpKv/Sln8adzw0NDeHxeJJ6X/X19XH3gMbzumYqTcnp06dZt26dNhdbgLwU5tWrV9m7d29K5/p8vqhW02w2J1WPMR7BYDClXDJlZWWa3RyiEW3ZYe4XPJq1sVqtkQpd0SpZLbQ0MTY2hsfjidu3maOZmYENJpMppWFwMhndE+XkyZPcfvvtC74XLcjLOabL5dJ8zhgMBjX5hxcXF6c8XM7WPHgmYWvz7W//E1/+8pQzyOPxpF1aLt4aZl1dHWNjYzQ0NODxeGaV6YsXPJ9uWpJUOHDggC7t5KUwTSYTf/VXf8VNN92kaSje+Pg41dXVaZVRC4eexcrpIyIxC7Jq8V4KCgooLy+PKvJE1yPnJsXSogxguLL0XIqKiujs7CQYDHLlyhW6urpmtRcrAVmqCaLTFfOZM2d0CanMy6EsTEWz7Nixg8bGRl566SVN0hDOnNMkOycKo5SiubmZgoIC6urqGB4exuVy4ff76evri2xsbmhomBXIXlBQMGutz2QyYbVa5y1VlJWVRRa5Q6EQExMTkRxB4eCK8JA1XA1LRCJ1SVIh3QI5MBUnHG3pY256lbn09fXNulH6fD76+/tT2vCt1TLKgQMHeOSRR5I/MQkSEqaIlAKPAzcACvgKcBb4N2A5cBl4SCk1MOe8OuDnQBUQAh5TSn1fm65PUV1dzVe+8hWOHz8e2T2+d+/etPO4hJcXUrUWo6OjkWifsOOpvLw8MocMD3fDaTHa2toYHR2lqKiI8vJyBgYGIvGh4fw2wWCQgYGBhNdkw9WwYCokMFWBtba2UldXl1b0UrQhaSLXDAQCs44JZypIxfOaipij8fHHH9Pa2ppyCYlESNRifh94WSn1eRGxAYXA3wCvKqX+XkS+C3wX+M6c8wLAXyul3hMRF3BMRF5RSsVffEwSEWHDhg2Rx+vWrUs7leDQ0JDmVbZmuv6dTidOpzNiNcvLyykuLubSpUvzLMjY2FjaZQfHx8fTykfT0tKCzWZLqYAPTG1gnjlKcDqdKY1IRCSy+TpZz6uWyygHDhzg0UcfTfn8hVhQmCJSDGwDvgyglJoAJkTkfuDO6cN2AweYI0ylVAfQMf33sIicBmoBTYUZpc+aXCeT0UKTk5ORnfThUnfhSlmZCtdL9/2UlZVFxJRKyNuVK1ci4X3FxcUpzdWuXLlCVVUVw8PDSXte011GCd8QAC5cuBDJNZQJEhlhrwR6gJ+KyPsi8riIOAHPtPDCAowbriIiy4FG4J0Yr39NRI6KyNF0EvhO9yet88OEHTmZoL29nZ6eHnw+36z6k1o4WqJRWlqadM6euYTzpaZTmau9vZ1rrrkmLQdKV1dXypkDUs32Xl5ePp0ixUldXR21tbWa+DVikYgwLcAm4EdKqUbAx9SwNWFEpAj4DfCflVJRk4wqpR5TSjUppZqSLaATDAa5evUq3d3dnD17VtMF+kzeFUOh0KzkwCaTSdNiqjPRYldEuK/phLyFQqG4WQkS8ZgqpXTb5REm/N59Ph8tLS20tbWl5b1fiESE2Qq0KqXClm4PU0LtEpFqgOnfUbdwiIiVKVH+UimVWlRAHJRSBAIBzGYzDoeD119/XdMaEkDGtpC53e5ZNxGz2Zyxjd5apF5pbm6mqqoqrcpcscoKJmuFu7q6qKioyKgDJkxNTc28m304dWimWPDTUkp1Ai0icu30U9uZmiP+Dtg1/dwu4Lm558rUZO8J4LRS6nua9Hh+G9jt9kio3a5du1ixYoXm7XR2dqa1a93r9WKxWCKpJysrK+el0p+cnNQki3emUEphNploGhxgz8ZvzwpCCDPX6s19HB6qz30+FSvc29sb8RhnCpPJFDVgxOv1arKMFItEvbLfAH457ZG9CPxHpkT9tIj8BdAMPAggIjXA40qpHcBW4M+Bj0Tkg+lr/Y1S6iXt3sJs7HY7X/ziF3n22Wc5efKkZtcNhUJ0dnYm5ZypqanBYrGglIp4NYeHhwkGgzGjg2bmYtWSdC2xKxRiXW8va/fupeH4cYJmM3/42/WMy7/f2+euE37pS7tnpSz51reep6urK+p6Yjoe03SjkuLh9Xqjft4NDQ0ZaxMSFKZS6gOgKcpL26Mc2w7smP77TUD3LFkWi4WdO3dSVFTEO+9E9TWlTKxY1traWgKBAH6/n9LSUkSE3t5eHA5HZGE9kaWG/v7+nMtdVBIMsutf/5WyGXNDczDI2pERPpyx42J+drrKWY9DoXIgdnBAqh7TVHLeJoLdbo/pr1i+fHlG2gyTlyF5iSAifPrTn+ZTn/qU5tedic1mw+v10tbWRldXF8PDw7S0tNDc3Izf72dkZCSp/ZU+ny8juUpTFfoav58/efXVWaIMUz3Hszp37llZ2RN5vHx5O8PDF6IeF7aOqXpMw8WWtKaqqirqJvFwWs9MkrcheYkgImzduhWfz8ehQ4c0uWYgEMDr9RIIBLBarZF5TiwmJibo6elJeAhsNpszkkkvFWGuGxnh/n/+Z6xRhopBs5l3V62KPA6va37pS7sZHZ2fZmTNmlLa26c8Olpvywon0aqurtYsjrWkpCTm5xX2F2SSRS3MMPfccw/r16/n7bff5sSJEylfx2azUV1dTVdXV1KROEopxsfHcTgcC56nlIo5r0mHZIVZOzHBff/6r1FFCVNDWfv0NaPNGcNiM5ngmmuKaG+fPQXQelvW5OQkHR0dmgVoOJ3OmFXOMj2/hEU8lJ1LdXU1O3fu5LOf/WzK16isrOTKlStJh8fV19cnLOZQKERzc7Pm3tlkhFkcCvG5xx/HuUDlrD//x3/kGz/8IaVDtrge1VTD+FIhvO6cTvSX2+2Ou55sCDMDJFqmbSbhsuipbmKeub8wUZxObfLUhEk0TtYSCvHQU0+xLIHMggWjoyzr7ubzL/+cdZ4LmCVAU+GH7Gw5xhq/H0IhiouLdU+E1tzcHKn6lQrxoq9MJlNGl2fCLImh7EyuvfZaDh48mHDoW0lJCYFAIK0ailVVVUkPr7T2yiYkTKV46JVXqP3448Svi2A5O8DLbMOEwjPSjfwbbAFGCwo4/8lPsjdD9T3i0dnZSVlZWdLWeqFUnbW1tRnzAs9kyVlMj8fDDTfcEPeY+vp63G43dXV1+P3+mKXKE6W5uTnp4Y/WtVgqL17kb/7H/+C//K//xa7f/paCuUJVigeOHmV1gk6yACaOs45tvEYN7Xhp5ws8jZqxOlYwOkq9Rk63VEj2c7NYLAtadz2GsbAELabP5+P06dNRX3M4HHg8noxU4urv709r21S6hEwmrJOTWCcncX74IX994gT//Ld/SwConpzktv37WXnsWELXCmCigj6GCFdPFhTwNrfSQyWeGdGZJf39OEMhfDpXY0vFCVRbW7vgZ28IM0OcPXs2qjhKSkoYGhrKWHm84eFhysrKsFgsae/ySIXgHAtsDgb5+j/8AwVJ9iWE8Ca3TYtSmNo3D6Bo4gjuKCHT/+EXvyAwJ4Vl99q1vDBjD62WmEympK2l0+lc0IcgIrrML2EJCrOxsZGhoSEOHjwITKW2mJycTDn3aTIMDAxQWVmZFWFGm7MmI8oQQhduHuYp3mYrZgIEsWAmQAi4mXd5i9uihnlVX7w47zlzIAAZEmYoFGJwcDCS+8jv98/KfhieJlgsFqxWK4FAALvdvmA2herqal2+J7AEhSki3HXXXdTV1VFUVITH4+H555/n/fff16X9RIMNtN6XWZ1iiQKYEuVd/JG32EoIMwoTFiZ5jTu5jTfooxI33UnFXhbo4Knt6+vTNO5Yr2EsLEHnT5hVq1ZRVVWFiNDY2Khr283NzQu68tva2jQN+6pMY9G9h0reZitBrCjAzCS38hZ3cBALCk+SogRwDg2ha95JDch0fOxMlqwwZ+L1evnqV7/Krl27WLlypS5txkrLOJO2tjbN6piURRlOJoqbbm7hLabyqcFNHOGPfDKt3Qm2yUnubG1lw+AgDePjFKWRj0gvMrVhPhpLbigbDRGJWKeGhgZOnTrF/v37014miUeiC/4DAwOahJmVpZHhXYBf8whe2gAT77CFbtxUk16miDueeCLy94jLxT/+9V+ndb1MUlVVtWABYC0xLOYcRIR169bx9a9/PaN3yN7eXioqKhY8LlxINZFj42FOs/SDmX+/kSjCCyRLBz3nl2AIMyZ2u51HH32Ua665JmNtJBp2Z7Va07bekuZ8zkM3t3MQM5PczsFZa5VLAUOYOYTVauXhhx9m7dq1Gbl+omF3k5OTkWTWqbBqdBRrmhZTgAN8kja8vM5d+u9+zzJ6C9OYYy6AxWLh/vvv59y5c5ovYSSz79JisUSG1uEMB+H1uEAgEHdx3NudvNc0GqZpD+xSo7KyMiOb1+NhCDMBHA4Hq1at4uzZs5pet6enB7fbTXcCOzkW2gBcW1sbU5zONMtFLHUyna0gGsZQNkEyVbBUq53wHR0dMcsIFsTY8JtLKI2D9rVEqyWrZDCEmSBr1qzJSDqJ9vZ2Tby/oVCI/v7+qBusHXkgTFMwSFIp3XUk2QTkWmAIM0HsdjurV6/OyLWbm5s1SVw8OTmJ3++nrq4u8uP1enGkEY6nF06fjwfeey/b3YhKuktVqWDMMZPgzjvvZGxsTPNM7zBV6q6mpibtEgkzy/+FyYehLMCGF17gxLp1nC8oyHZXIlit1rQ84qliWMwkcLvdPPjggxm7fldXV0bmM4V5IkwAV5b2q8aisrJS803riWAIM0kyuQslGAwyODioqRfQGgziSLO2pp6ku96qNdkYxoIhzKRZv369JgV6YjExMRHZWaLFnXpFHokSiJkuM1tkw/EDhjCTxuVyUV1dnfF2tNj2VRIIcO/PfqZNh3TClmMW0xBmHpEp7+xcYhUeSgRLMMhDTz0VtbRBLmPJsTmmMZTNI/TaszkwMJBa6T+l+ML+/dScP699pzJMurtgtMRsNmekJkoiGMJMgdraWt1yv6Qyn73/vfdYpXGVM72w5pDFrKioyKg/IR6GMFPAZDLplmaivb09oWwHYe66coWNzz+fuQ5lkIDZTDCDxWCTJVvDWDACDFJm5cqVmge1x8LlcjGYQPRO48AAW3fvznyHNKKvqoqT993HQFkZnYWFdJpMU9WGcoRsOX7AsJgps2bNGt3aam9vX3DovHJsjHt+/GPMORpvGo3LW7fyWkMDHxQX02mx5JQowRBmXlJaWqrbroNgMBjXCVQeCPCnjz2WV4EEAKd1vLmlQjaHsoYw00BPqxkuAT8XezDIg08+SUkaRY+ywUBlJRdSrMalByKieSnEZDCEmQZ6CnN4eJiamppZz0koxMPPP48nQ2UdMklrUxNpl5LOIOXl5Ziz6IhKSJgiUioie0TkjIicFpFbRGSZiLwiIuemf0dd8BGR/y0i3SKSeinnHEWvkmxhAnPC1T53+DDLP/hAt/a15Nz112e7C3HJ5jAWEreY3wdeVkqtBTYAp4HvAq8qpVYDr04/jsbPgHvT7GdOEggEZtXEyDRdXV2RpZN7zp3jhv37dWtbS0aKijihcWFercmm4wcSEKaIFAPbgCcAlFITSqlB4H4g7JvfDTwQ7Xyl1EEgvyZACZJOMdtUKS4u5hM9PXziV7/SvW2taNu0CZVjHti55LwwgZVAD/BTEXlfRB4XESfgUUp1AEz/TstFKSJfE5GjInK0J0/iO8+cOaN7m+1tbdz25JOY8qzux0wu3nhjtruwIPkwlLUAm4AfKaUaAR+xh60po5R6TCnVpJRqyvbdKhHGx8d5Jwthb4FgkMtbtujerlZM2Gx8sGxZtruxIPkgzFagVSkV/hbuYUqoXSJSDTD9e0klHH333XcZy9K64dtNTYRyfCgYi7b165nI8b6XlZXp6tSLxoL/IaVUJ9AiItdOP7UdOAX8Dtg1/dwu4LmM9DAHmZiY4PDhw1lrv8Nm4/LmzVlrPx2u6FzyMBWytaNkJoneur4B/FJEjgMbgb8D/h64W0TOAXdPP0ZEakTkpfCJIvJr4BBwrYi0ishfaNj/rGA2m7O6xgVw7I47stp+KgRNJt6fsxabixTkQDKwhILYlVIfAE1RXtoe5dh2YMeMx4+k2rlcxWw2s2XLFl555ZWs9eFUUREd11xD9YULWetDsnSsWcPVHB/GArpt6YtH7v+XcpTNmzdn/QP86J57stp+srTmyfBbzzqYsTCEmSJ2u53NWf6iHa6sZDDL3sNkOK5zxaxUyfYNFwxhpsWWLVuytsMdQJlMnP70p7PWfjL01tTQkcNB6zMxLGae43K5uDHLi+UHV65kNAecFQvRtmlTtruQMIYwFwG33HJLVtsfM5v5ePs8H1zOcS5DxX8zgSHMRYDb7WbVqlVZ7cNbGzbkVK6cuYwUFXFS58Kv6WDMMRcJt956a1bb77FaufiJT2S1D/Fo37gx59KGxMOwmIuE5cuX65KdPR5Hbr89q+3H4/L69dnuQlIYwlwkiEjWrea5ggJac3AeF7BY+CAPNiXMxBjKLiKuv/76pPK/ZoLjd9+d1faj0b52LaN5NIwFQ5iLCpPJxJYsb8c6smwZfamUVMggrRs3ZrsLSWG327NSD3MuhjA1pLGxMbsB0CKcuje3sric0CljvVbkwvwSDGFqis1mo6kpWqy/frxRV4evqCirfQjT4/XSYcmvZP+GMBcpN998c1a3hE2azZzNkblmPkX7hMmF+SUYwtScoqIiNmzYkNU+vLluHZM5YKnOXXvtwgflGIbFXMRke+lkwGLhwtatWWt/0mrl1PbtnMqDGN655Iows39bXYSUl5ezceNGPshiMuZ3br2Vta+/rmubfqeTs3ffzRvr1zOQwyGC8ciVoawhzAxx33330dzcnJXcswCX7Xau3HgjDcePZ7ytwfJyTu3YwRsrVjCWZ2uWczEs5iLHZrOxc+dOnnjiCUJZKo334fbtGRVmd10dH+3YwVseT84ncE4Uw2IuAWpqati+fXvWcgO9X1zMlro63C0tml63ed06Prj7bt4vKcnpwkCpYFjMJcItt9zCxYsXuZCNpFkinPr0p3E//njalwqJcPETn+DIHXfwcR46dRLFEOYSQUR44IEH+NGPfoTf79e9/Teqq9lcUoJraCil8ydsNs7dcQdv3XRT3qQGSYdcGcoujolBjlNUVMQDDzyQlbZDZjNnUsim53M6ee9zn+OH3/kOe7ZuXRKiBMNiLjlWr17NJz7xiazUO3lj7VputNmwT0wseGy318u5T36SN5YvZ3yROHSSwRDmEuRTn/oUly9fpqurS9d2h81mzm/bxro//CHq6+M2G5duvZUPb7qJMzletzLT5MpQ1hCmjlgsFnbu3Mljjz02rzp0pjl0881c9+qrs8r3da5Ywflt2zjU0IB/CVrHaBgWc4lSWVnJvffeywsvvKBru202G5c3bcJ99iyXbr+d9zdu5FKOWIdcwWQyYcmBGGMwhJkVNm3axIULFzh9+rSu7f7xwQfp9PsJGtYxKg6HIyc2SYPhlc0KIsJnP/tZiouLdW23bWyM8hzLcJBL5Mr8EgxhZo2CggI+97nP6d5urgzVcpFcmV+CIcys0tDQwO06p51sb2/X3VLnC4bFNIhw55134vV6dW0z29n8chXDYhpEMJlM7Ny5U9e7dXt7e05Zh1whl/4nhjBzgNLSUj7zmc/o1l4gEKDKcALNwxCmwTxuuOEGXXMF9fb2ZrW2Zy5iDGUNonLfffexbNkyXdry+XzU1NTo0la+YFhMg6jY7Xbu1jH15OjoqG5t5QN5J0wRKRWRPSJyRkROi8gtIrJMRF4RkXPTv8tinHuviJwVkfMi8l1tu7/4WLNmDS6XS5e2+vr68Hg8urSVD+TjUPb7wMtKqbXABuA08F3gVaXUauDV6cezEBEz8APgPuB64BERuV6Lji9WTCYTjY2NurZnMEVeWUwRKQa2AU8AKKUmlFKDwP3A7unDdgMPRDn9ZuC8UuqiUmoCeGr6PIM4bNIxg3lHR4exrjlNvlnMlUAP8FMReV9EHhcRJ+BRSnUATP92Rzm3FpiZCap1+rl5iMjXROSoiBzt6elJ6k0sNkpKSli9erVu7ek1dM518spiMrUDZRPwI6VUI+AjyrA1BtFC9VWU51BKPaaUalJKNVXmWaHTTLB582bd2mpra8tulbIcId8sZivQqpQK58TYw5RQu0SkGmD6d3eMc+tmPPYC7al3d+mwevVq3SxZKBTC7Y424Fla5JXFVEp1Ai0iEq4Qsx04BfwO2DX93C7guSinHwFWi8gKEbEBD0+fZ7AAejuBuru7s1qlLNuICFarNdvdiJCoS+4bwC9F5DiwEfg74O+Bu0XkHHD39GNEpEZEXgJQSgWA/wTsY8qT+7RS6qSm72ARs2nTJt027o6Oji7pgINc2iQNCWYwUEp9AESryLo9yrHtwI4Zj18CXkqxf0uakpISVq1axblz53Rp7+rVq7q0k4vk0jAWjMifnEdPJ9DQ0NCStZqGMA2SQk8nEEAwGNStrVwilzyyYAgz5zGZTLoGHHR1dVFeXq5be7mCYTENkqaxsVFXx8RSXNM0LKZB0ugdCdTa2kpRUZFu7eUChsU0SAk9h7OAbvtCcwVDmAYpsXr1al2z23V2dmJbIhW+wBjKGqSI3pFAExMTSyovkGExDVJGz0gggIGBgZyKhskkhsU0SJni4mJdnUDDw8PU1kbdpbfoMCymQVroGQkEMD4+rmt72cIQpkFarFq1SlcnUE9Pz5LYEmYMZQ3SQm8nEJBT26EyhWExDdJGbydQW1vboi9EZFhMg7TR2wkEi78QUa6t2RrCzFP0dgJ1dHTk3JdXK+x2e84tCxnCzFP0dgJNTk4u2oCDXBvGgiHMvEXv7WAA/f39OWdZtCDXHD9gCDOv0Xs72MjIyKIMODAspoGmFBcXs2bNGl3bHBsb07U9PcjFnTQJJeMyyF02b97M2bNndWuvt7cXt9tNd3e0NMK5jdVqxePx4PF4qKqqoqqqCrfbnZNOLUOYec4111xDSUkJQ0NDurWZDwEHLpeLqqqqWSIsKyvLmyJKhjDznHAk0IEDB3RrMxxwkAvpLkWEysrKWQL0eDw4nc5sdy0tDGEuAhobG3n99ddRKmpZmIxQWlqquzDtdvs8AbrdbiyWxfc1XnzvaAkSdgLpOddsb2/HZrMxMTGRkeuXlJTMG4qWlpYuyuWaaBjCXCTo7QQKBAI0NDRw5cqVtK5jNpuprKycJUKPx7MkM/XNxBDmIiEbTqC+vj5EJOEhdEFBwTwrWFFRsaSLGcXCEOYiIRwJ9Nprr+nW5sjICF6vl9bW1nmvLVu2bJ4IXS7XkhmKposhzEVE2DurpxMoEAhQW1s7b20wF8Pc8glDmIsIl8vFtddey5kzZzLWhslkYu3ataxdu5bq6mqWLVuWN2uD+YQhzEXGpk2bMiLMkpISNm3aRGNjo65FjpYqhjAXGVo7gVavXk1TUxOrVq0yLKOOGMJcZGjhBHI6nTQ2NrJ58+ZFn7kgVzGEuQhJ1Qm0fPlympqaWLt2rbGEkWUMYS5CknECORwONm7cyObNm6moqNChdwaJYAhzkbJ58+a4wvR6vWzevJl169blxW6RpUZCwhSRy8AwEAQCSqkmEdkA/CtQBFwG/oNSal5Us4h8E/g/AAF+opT6Z016bhCXaE4gq9XKjTfeSFNT06LN37NYSMZi3qWU6p3x+HHgvyilXheRrwD/J/D/zDxBRG5gSpQ3AxPAyyLyolLqXJr9NlgAEYk4gTweD01NTaxfv95Y+M8T0hnKXgscnP77FWAfc4QJXAccVkr5AUTkdeDPgP8vjXYNEmTz5s2sWLECr9drhMLlGYkuTClgv4gcE5GvTT93AvjT6b8fBOqinHcC2CYi5SJSCOyIcRwi8jUROSoiR3t6ehJ/BwYxcTqd1NXVGaLMQxIV5lal1CbgPuCvRGQb8JXpv48BLqaGqrNQSp0G/idTFvVl4EMgEK0BpdRjSqkmpVRTZWVl8u/EwGARkZAwlVLt07+7gd8CNyulziil7lFKbQZ+DVyIce4TSqlNSqltQD9gzC8NDBZgQWGKiFNEXOG/gXuAEyLinn7OBPzfTHloo50fPq4e+BxTIjYwMIhDIhbTA7wpIh8C7wIvKqVeBh4RkY+BM0A78FMAEakRkZdmnP8bETkFPA/8lVJqQNN3YGCwCBE99+4lSlNTkzp69Gi2u2FgkHFE5JhSqmnu88Z2AQODHMQQpoFBDmII08AgBzGEaWCQgxjCNDDIQXLSKysiPUB6mYRTpwLoXfCo3MLosz5kos8NSql5oW45KcxsIiJHo7mvcxmjz/qgZ5+NoayBQQ5iCNPAIAcxhDmfx7LdgRQw+qwPuvXZmGMaGOQghsU0MMhBDGEaGOQgS0KYIlInIq+JyGkROTmduW/m6/9FRJSIRE2sKiKlIrJHRM5MX+OWPOjzt6bPOyEivxYRR7b6LCL/TUTaROSD6Z8dMc6/V0TOish5EflupvurRb8X+pxSRim16H+AamDT9N8u4GPg+unHdUwlErsCVMQ4fzfw1em/bUBpLvcZqAUuAQXTj58GvpytPgP/jamMivHONTOVBWPl9P/4w/D7zfF+x/yc0vlZEhZTKdWhlHpv+u9h4DRTX16AfwL+K1MJx+YhIsXANuCJ6fMnlFKDudznaSxAgYhYgEKmNrNnlAX6vBA3A+eVUheVUhPAU8D9menpbNLpd5rvOSZLQpgzEZHlQCPwjoj8KdCmlPowzikrgR7gpyLyvog8Pp1iRTeS7bNSqg34B6AZ6ACGlFL79ehrmJl9nn7qP4nIcRH53yJSFuWUWqBlxuNWNPiCJ0sK/Y53bsosKWGKSBHwG+A/M5Wt72+B/3eB0yzAJuBHSqlGwAfoOf9Jus/TX6D7gRVADeAUkUcz29NZ7Uf6rKay8/8IuAbYyNSN4h+jnRblOV3X8lLsd6xz02LJCFNErEz9436plNrL1D98BfDhdAkIL/CeiMytHdAKtCqlwnfBPUwJNZf7/CngklKqRyk1CewFbs1Sn1FKdSmlgkqpEPATpoatc2llds5hLzoMv8Ok0e+o56bLkhCmTGU8fgI4rZT6HoBS6iOllFsptVwptZypL8YmpVTnzHOnH7eIyLXTT20HTuVyn5kawm4RkcLp62xnau6je5+nn6+ecdifMZUIfC5HgNUiskJEbMDDwO8y2d8Z/Uu537HOTRs9vF7Z/gFuY2pYdBz4YPpnx5xjLjPt4WRq+PfSjNc2Akenz38WKMuDPv93pjIYngCeBOzZ6vN0+x9NP/87oDpGn3cw5dW8APxttr8fifQ7kc8plR8jJM/AIAdZEkNZA4N8wxCmgUEOYgjTwCAHMYRpYJCDGMI0MMhBDGEaGOQghjANDHKQ/x+zQ3hj3hXKsAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "# Create a figure with one subplot\n", "fig, ax = plt.subplots()\n", "\n", "# Plot polygons\n", "polys.plot(ax=ax, facecolor='gray')\n", "southern.plot(ax=ax, facecolor='red')\n", "\n", "# Plot points\n", "data.plot(ax=ax, color='blue', markersize=5)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Okey, so we can see that, indeed, certain points are within the selected red Polygon.\n", "\n", "Let's find out which one of them are located within the Polygon. Hence, we are conducting a **Point in Polygon query**.\n", "\n", "First, let's check that we have `shapely.speedups` enabled. This module makes some of the spatial queries running faster (starting from Shapely version 1.6.0 Shapely speedups are enabled by default):" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#import shapely.speedups\n", "from shapely import speedups\n", "speedups.enabled\n", "\n", "# If false, run this line:\n", "#shapely.speedups.enable()" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "- Let's check which Points are within the ``southern`` Polygon. Notice, that here we check if the Points are ``within`` the **geometry**\n", " of the ``southern`` GeoDataFrame. \n", "- We use the ``.at[0, 'geometry']`` to parse the actual Polygon geometry object from the GeoDataFrame." ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 True\n", "1 True\n", "2 True\n", "3 False\n", "4 True\n", "5 False\n", "6 False\n", "7 False\n", "8 False\n", "9 False\n", "10 True\n", "11 False\n", "12 False\n", "13 False\n", "14 False\n", "15 False\n", "16 False\n", "17 False\n", "18 False\n", "19 False\n", "20 False\n", "21 False\n", "22 False\n", "23 False\n", "24 False\n", "25 False\n", "26 False\n", "27 False\n", "28 False\n", "29 False\n", "30 True\n", "31 True\n", "32 True\n", "33 True\n", "dtype: bool\n" ] } ], "source": [ "pip_mask = data.within(southern.at[0, 'geometry'])\n", "print(pip_mask)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "As we can see, we now have an array of boolean values for each row, where the result is ``True``\n", "if Point was inside the Polygon, and ``False`` if it was not.\n", "\n", "We can now use this mask array to select the Points that are inside the Polygon. Selecting data with this kind of mask array (of boolean values) is easy by passing the array inside the ``loc`` indexer:\n" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "deletable": true, "editable": true }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
addressidaddrgeometry
0Ruoholahti, 14, Itämerenkatu, Ruoholahti, Läns...1000Itämerenkatu 14, 00101 Helsinki, FinlandPOINT (24.91556 60.16320)
1Kamppi, 1, Kampinkuja, Kamppi, Eteläinen suurp...1001Kampinkuja 1, 00100 Helsinki, FinlandPOINT (24.93169 60.16902)
2Kauppakeskus Citycenter, 8, Kaivokatu, Keskust...1002Kaivokatu 8, 00101 Helsinki, FinlandPOINT (24.94179 60.16989)
4Hesburger, 9, Tyynenmerenkatu, Jätkäsaari, Län...1005Tyynenmerenkatu 9, 00220 Helsinki, FinlandPOINT (24.92160 60.15665)
10Rautatientori, Kaisaniemi, Kluuvi, Eteläinen s...1011Rautatientori 1, 00100 Helsinki, FinlandPOINT (24.94410 60.17133)
30Kampin keskus, 1, Urho Kekkosen katu, Kamppi, ...1031Urho Kekkosen katu 1, 00100 Helsinki, FinlandPOINT (24.93312 60.16909)
31Ruoholahdenkatu, Hietalahti, Kamppi, Eteläinen...1032Ruoholahdenkatu 17, 00101 Helsinki, FinlandPOINT (24.92496 60.16497)
32Easy Cycles, 3, Tyynenmerenkatu, Jätkäsaari, L...1033Tyynenmerenkatu 3, 00220 Helsinki, FinlandPOINT (24.92119 60.15890)
33Oluthuone Kaisla, 4, Vilhonkatu, Kaisaniemi, K...1034Vilhonkatu 4, 00101 Helsinki, FinlandPOINT (24.94709 60.17191)
\n", "
" ], "text/plain": [ " address id \\\n", "0 Ruoholahti, 14, Itämerenkatu, Ruoholahti, Läns... 1000 \n", "1 Kamppi, 1, Kampinkuja, Kamppi, Eteläinen suurp... 1001 \n", "2 Kauppakeskus Citycenter, 8, Kaivokatu, Keskust... 1002 \n", "4 Hesburger, 9, Tyynenmerenkatu, Jätkäsaari, Län... 1005 \n", "10 Rautatientori, Kaisaniemi, Kluuvi, Eteläinen s... 1011 \n", "30 Kampin keskus, 1, Urho Kekkosen katu, Kamppi, ... 1031 \n", "31 Ruoholahdenkatu, Hietalahti, Kamppi, Eteläinen... 1032 \n", "32 Easy Cycles, 3, Tyynenmerenkatu, Jätkäsaari, L... 1033 \n", "33 Oluthuone Kaisla, 4, Vilhonkatu, Kaisaniemi, K... 1034 \n", "\n", " addr geometry \n", "0 Itämerenkatu 14, 00101 Helsinki, Finland POINT (24.91556 60.16320) \n", "1 Kampinkuja 1, 00100 Helsinki, Finland POINT (24.93169 60.16902) \n", "2 Kaivokatu 8, 00101 Helsinki, Finland POINT (24.94179 60.16989) \n", "4 Tyynenmerenkatu 9, 00220 Helsinki, Finland POINT (24.92160 60.15665) \n", "10 Rautatientori 1, 00100 Helsinki, Finland POINT (24.94410 60.17133) \n", "30 Urho Kekkosen katu 1, 00100 Helsinki, Finland POINT (24.93312 60.16909) \n", "31 Ruoholahdenkatu 17, 00101 Helsinki, Finland POINT (24.92496 60.16497) \n", "32 Tyynenmerenkatu 3, 00220 Helsinki, Finland POINT (24.92119 60.15890) \n", "33 Vilhonkatu 4, 00101 Helsinki, Finland POINT (24.94709 60.17191) " ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pip_data = data.loc[pip_mask]\n", "pip_data" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Let's finally confirm that our Point in Polygon query worked as it should by plotting the points that are within the southern district:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "deletable": true, "editable": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAAEYCAYAAACnR+FYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7pUlEQVR4nO29eXBb55Xg+ztYSYI7wQVcJdmOZFnWRsqyIu+KPY6yWI7jztKZOEl3u1PpSXdPvZ6Z9Ey/VzM1U13trvfyxtWdTmIncdwvjmNbjjuyrUiyvEiOJCuWZVmStVjWQhLcKRHcRBIE8L0/CCCkCJBY7r0AyPurYoEA7r3fAYlzz/ed7yyilMLExCT7sGRaABMTk9iYymlikqWYymlikqWYymlikqWYymlikqXYMi1ALNxut1qyZEmmxTAx0Z333nuvXylVGeu9rFTOJUuWcOTIkUyLYWKiOyLSGu89c1prYpKlmMppYpKlmMppYpKlmMppYpKlmMppYpKlmMppYpKlmMppYpKlJKScIlIqIttF5IyInBaRTSJSLiKvici58GPZHOdbReR9EXlFO9FNTBY2iVrOx4FdSqkVwBrgNPA94HWl1A3A6+Hn8fir8DkmJiYJMq9yikgxcAfwUwCllF8p5QMeAJ4OH/Y0sC3O+fXAZ4CfpC+uiZaEQiGGh4czLYZJHBIJ31sG9AFPicga4D2mLGG1UqoLQCnVJSJVcc7/38B/BormGkREHgUeBWhsbExIeJPkUUrR2dnJiRMnOHXqFA888ABFRXP+a0wyRCLKaQPWA99VSh0WkceZewobRUQ+C/Qqpd4TkbvmOlYp9QTwBEBLS4tZO0VjLl++zIkTJzhx4gRXrlwBQERMxcxiElFOL+BVSh0OP9/OlHL2iIgnbDU9QG+MczcDnxeRrUAeUCwiv1BKfU0L4U3mZmJigvfff58TJ07Q2dk5632lFD6fj6qqeJMek0wy75pTKdUNtIvI8vBLW4BTwA7gkfBrjwC/iXHu3yql6pVSS4AvA2+YimkMSileeukldu/eHVMxI7zyyisMDAwYKJlJoiTqrf0u8IyIHAfWAn8P/ANwr4icA+4NP0dEakVkpw6ymiTBwYMHOXv27LzHDQ8P86//+q+Mjo4aIJVJMiSUz6mUOga0xHhrS4xjO4GtMV5/C3grKelMUqK1tZXXX3894eN9Ph/nz59n9erVOkplkixmhNACY2RkhO3bt5NsPeJgMKiTRCapYirnAiIUCvHrX/+akZGRpM81t6+yD1M5FxBvvfUWFy9eTOnco0ePaiyNSbqYyrlA+Pjjj3n77bdTPv/gwYOcOHFCQ4lM0sVUzgXA4OAgv/71r9O+zo4dO+ju7tZAIhMtMJUzxwkGg7zwwguMjY2lfa1AIMBzzz3H1atXNZBMf0KhUNKOr1zCVM4c57XXXqOjo0Oz6/l8Pl588UVCoZBm10yXvr4+Xn31Vfx+f/Q1pRS7d+9m+/btBAKBDEqnH1lZt9YkMU6dOsXhw4fnPzBJLly4wJ49e7j//vs1v/Z8nDx5kkAgQEVFBXV1dSil2Lt3Lx999BEdHR18+ctf5vTp01y8eJGzZ89it9sZHx+nsLDQcFn1RrJxWtDS0qLMotJzc/nyZZ544okZ1kRr7rvvPjZt2qTb9acTDAbZuXMnR48exWKxEAqF2LhxI729vTM80FarddaebH19Pd/61rcQEUNk1RIReU8pFSvAx7Scucjk5CQvvPCCrooJsGfPHoqLi7npppt0HQfAYrEwMTEBEJ1Sx5oVxAqW8Pv9HDp0iIKCApxOJw6HA4fDgdPpzOmgflM5NWJsbAyHw4HVatV9rF27dtHT06P7OAAvvfQSbreb6upqXccRET73uc/R2dmZdCB+b28vr732Wsz3/vIv/5KysrgVdLIa0yGUBkopvF4ve/fu5Uc/+hFPPfWU7t7Dq1evcunSJV3HmE4wGGTHjh2GOIicTicPP/wwFot2X8tczrgxLWeCXL16lYKCAoLBIJcuXeL06dOcPXt2Rqjc0NAQb775Jrfffjt2u11zGfr6+ti1a5fha6vOzk4OHTrE5s2bdR/L4/HwqU99ij179mhyvTfffJNDhw5RUlJCQUEBN998M8FgkFAoRDAYpLKykry8PE3G0hrTIZQA58+f5xe/+AUlJSWMj49H10bxcLlc2O12PB4PN910kyZrNr/fz5NPPkl/f3/a10oFm83Gt7/9bSoqKnQfSynFL3/5Sz7++GNdx1m6dClf+9rXNLXUyTKXQ8ic1s5DKBSKrmcGBwfnVUyA0dFRfD4fp0+fZvv27WlbAaUUr776asYUE6YCFF5++WVDNv1FhG3btum6PVJSUsIXv/jFjCrmfGSvZFnCiRMn0na+HDp0iP3796ccGvf+++9z/PjxtGTQgtbWVsP6prpcLh588EFdrm2z2fjSl75EQUGBLtfXClM55yAQCPDGG29ocq0333yTH//4x/zqV7+as2zItXR2dvLb3/5WExm0YO/evQwODhoy1rJly3RZ537uc5/D4/Fofl2tMZVzDs6dO8fQ0JCm1zx79ixPPvkkv/zlL+cNuxsdHeX555/PqvA0v9/PK6+8YlhM6913362pw2bt2rU5U/HB9NbOwcmTJ3W79rlz5zh37hzLli3j9ttvp6mpaYYXNhQK8eKLLxpmpZLh448/5tixY6xbt073sS5evKjpzenChQv8/Oc/Z/PmzTgcDoqKiigvL9fs+lpiKmcc/H4/H330ke7jXLhwgQsXLtDY2Mjtt9/Oddddh4iwd+/elBOnjeC3v/0t9fX1VFZW6jbGBx98wG9+8xscDge1tbX4/X6Gh4fTKkY2NDTE0NAQra2tAGzatIn77rtPK5E1xVTOOHz00UeGTifb2tp45plnqK2tZeXKlRw6dMiwsVNhcnKS559/nj/7sz/D4XBoem2/38+BAwe4dOkS1dXVWCwW2traou87HA4qKysTzsZZunQpbW1tMUP/3G63ZnJrjamccTh16lRGxu3s7MRmy41/S39/Py+//DJf+MIXNAmMCAaDHD16lH379s1pHf1+f1Jhklu2bOH06dMcOHBg1nvZ7LE1HUIx8Pv9nDt3LmPjZ1Mu5Xx0dnayb9++tK6hlOLDDz/kX/7lX9i5c2dC01afz5fQtV0uF6WlpSxbtmzG61arlS1btrB8+fI4Z2ae3LhFG8yZM2cy6iHNlTKVLpeL4eFh9u3bh9Vq5fbbb0/6GhcvXmTv3r1JbS8BlJaWEggEYlZtaGpqoqKiAofDwZIlS3A6nbz11lvU1taSl5eH0+lkzZo1Wa2YYCrnLAKBAG+99VbGZcgFRkdHcbvd9Pf388Ybb2Cz2Vi1ahV+vz/6EwqFouVEItsvwWCQiYkJTp48mXKIXmQNWltbi81mo7+/n8rKSkKhEE6nk9tuu21GNsq3vvWt9D+wwZjKeQ1HjhzJeCZDIiGCeuN2u1FKMTk5icViweVyYbFYaG9vjx7T2Ng4w1GzZ8+eeUMVKysr6evr00zOiMXNz8+PemAjr3/1q1/NiWCDeJjKOY3x8XH279+faTEYHx/PtAg4nc4Z3lCfz4eI4Ha7cTqdWK1WvF5v0tfVK9/12gJnIyMj/PznP+dLX/rSrPVmrpCTyjkyMqJLUPSBAwc0qWKXDiKie4WDRIgVEK6USiv4vrCw0NBO2n6/n3379jEwMEBzc7Nh42pFTnpr9+/fzxNPPMG7776rmTL19fXxzjvvaHKtdMjPzwegrKyMxsZGrrvuuoxsrfT392tu5YqLi6Oe2IqKihnbL3qlol29epVXXnmFvXv35lwZzZy0nFVVVbz77rt0dXWxZ88eVq5cSXNzMw0NDUnvt/n9ft5++20OHTqUFV7SUChEdXU1PT09DAwM0NjYSFFRESIS7UhtBGNjY5SVlWm6/p5uja1WK0opbDYbFouFy5cvU1pamvAWSSI4HA4uX74MTM2KBgYG2LZtmy6J8HqQs8oZIRAIcPz4cY4fP47b7WbdunWsWbMGl8s15zWUUpw+fZrdu3drHtyeDuPj4zPWnEopBgYGsFgssxwwemKz2TRVFJi5RRSZuk/3TBcUFODz+SgoKMBms6X9f6mpqZnx9zp16hQDAwN8+tOfpqGhIa1rG0FOVkIYHx/nsccei/u+xWJh+fLlrF+/nmXLls1aPw0MDLBz507dM+21oLy8fIbFrK2tZWBgQPe1cVNT0wzvpxZ4PB4CgQAFBQUxr11bW0sgEKC3t5eSkhJKSkoYHh5OyXo7nU4gvud7+fLl3HPPPRmvzjdXJYScVE6A73//+wk5F0pKSli7di3r16/H5XJx8OBB9u/fnzN7iTBlxa61MMXFxbr1NdFDMWFqG6W/vz/ptV8qM4ZEPoOIsGbNGu68805KS0uTur5WLEjl/MUvfsH58+cTvqaIsHTpUlpbW7NibZkM8fYGGxsb6ejo0PTz1NfXp7RFkihVVVX09vYmfV5xcXHcaW51dTXDw8PRaKGioiJGR0cTDoO0Wq1s2LCB22+/3fBY2wVZQyjZ6YhSigsXLlBRURH1iOYK8eRta2ujsLBQs6lZZWWlroqZKiISc+/XbrdTVlbG+Pg4Vqs16tUuKSlJKj45GAzyzjvv8Pjjj+uaw5ssi0Y5I/T29mK32zM2jUmFub5og4OD9Pb20tjYmLYXUuvUr1iMjo4mPY5SasZWS3FxMSUlJcCU/2BwcJDh4WHq6upwu90p32D8fj+vvvpqWvmiWpKQcopIqYhsF5EzInJaRDaJSLmIvCYi58KPs8pqi0ieiPxeRD4QkQ9F5H9oJXg6Sb5DQ0OMjY3pXsVcKxJpI9/W1paWBS0pKdG0W1k8KioqUgqy6OrqorGxEZiatg4ODjI5OTnjmNbW1rR9CePj4+zduzeta2hFopbzcWCXUmoFsAY4DXwPeF0pdQPwevj5tUwA9yil1gBrgftF5Na0pSY95YQpL15vb29OuNR9Pl9C1iYSYpcKudClq62tjcbGxrjlLC0WiybbP8eOHTNsy2ou5t3nFJFi4A7gGwBKKT/gF5EHgLvChz0NvAX8l+nnqilvU+S2bw//aOKBcjgcaW9aK6Vob2+ntrYWmNpzm5ycJC8vD4fDQTAYnJFNASSd2qQV5eXl83pnR0dHqaurS9oCNjQ0zAho1wun05l20PtcShNJYdOCV199lT//8z/PaF3bRIIQlgF9wFMisgZ4D/groFop1QWglOoSkZhzKhGxhs+5HviBUkqzhpJVVVWa3CkzpXDJkOg6LVnPbX19vSGKCVNeVT0tUlFRkWbK2dvby+HDhw1rgRiLRG4LNmA98EOl1DpglNhT2JgopYJKqbVAPXCLiKyKdZyIPCoiR0TkSKJ3Vz2LS8UiF8qHdHd3J+zsampqMsw763A46Orq0nUMrf8/b731VkLKPjQ0NMsv4Pf7Z62JkyUR5fQC3mkWbztTytojIh6A8OOcm1dKKR9TU9+Y7ZKVUk8opVqUUi2JKp3R0R25EpNZXFw87zHV1dWGtCuMUFNTk/aX1Wj8fn9CrTSOHj3KD37wA9555x0OHTrE008/zWOPPZZ2Hap5lVMp1Q20i0ikpsMW4BSwA3gk/NojwG+uPVdEKkWkNPx7PvAp4ExaEk/DaMtpxFZDPJL5YsfarJ9uVTweDz09Pfh8Purq6qirq6OpqUnXbltax+nGQo/gkpMnT3LhwoW473d1dfH2228zPj7O7t272bNnD6Ojo3zzm99Mu3h1ovOA7wLPiIgDuAB8kynFfl5E/gRoAx4GEJFa4CdKqa2AB3g6vO60AM8rpV5JS+JpGPEPn04mp7XJxNL6fD5cLld0v664uBi/3095eTmXL1+OTi+vzXKxWCy6hO5pnd0SD70qSJw+fXpWwnYgEODcuXO89tprs/ah+/r6GB4eTrsiYULfNqXUMSBWiNGWGMd2AlvDvx8HNC8LfvXqVXbv3m14cx8jp4HXkshe53TKysoYHR2NRs4MDQ3NW2EhFAppnpZWXV09r+zl5eW4XC6uXLmSVgCAXhUk6urqgD949z/44ANOnTo153g7d+6kqakprXDA7PdwxODixYsZqfOTKeXMy8tL+osXkdXj8STl9NFyOwKm1m2xFM7tdiMiWCwW+vv7uXLlCvn5+ZSVlTEyMkJNTU3SXmS9lXNoaIjnnnsuZsW/6TgcDh566KG043RzUjkjDWl9Ph979uzh9OnThoybKeV0uVxJf/F6e3upq6tL2hur9bp6YGBgVrC73W7n8uXLs7JTIlP3qqoq2tvbkwqSvzZzRysie7Nut5v33ntvXsUE+PSnP82SJUvSHjtnY2thqnbpypUrDRvP6HbvESK5ickwNjaWUjieHnmi1zqaamtr46aNTZc7Ge94JNZWa8rLy/m3f/s3XnzxRd5+++2Ezjl48KAmY+ek5ZxOuo1tkyFTOaBGbuHoUQolEAjgcrlwOBxxE61jMTo6mrBFLCgoiJYk0RK73c7k5CQffvhhwudo9T3JaeU8d+4cv/vd7wwbL1OW06ic2/z8/JQsZ319PRaLhatXr8asznf58mVCoRCjo6NJ+Qp8Ph8NDQ0Eg8F5o7i0bGFRX1+PiDA2NpZSRFNTUxNKKWO8tdmIz+fjpZdeMnRMi8VCTU0Ng4ODhpbQNMpiJ6OcVqsVj8fDxMREdF2bn59PRUXFLAtWUVGRcohke3s7IhLNSBkdHWV4eJi8vDxCoRAjIyNUVVVpEoJpt9vxeDy0t7endUM8duwYIyMjfPWrX01LQXNWOXfv3m14jdnpzpW8vDyKi4txOp2ICF6vV7cGREZVgC8sLJxzWisiNDQ0RKsvXOtsGhsbm7U+tlqtaTcAVkrNsmCRTmMNDQ0EAoG0lMnlclFRUUF3d7cmsb+RCvmL1nJmuir6tVXy9CQRD6GeFBcXU1hYSG9v77xf3sLCwhnBIYWFhbp15w4Gg7S3t1NfX5/S+ZGqGB0dHWknWIsIW7duRSnFqlWrNKm2kbPKmU21j6ZH4+hBZWWlIfmFbW1t0cLPDoeD8vJygsEg3d3dCZep7OnpmeHEcblcuilnBK/Xm/R6Wes0uaamJlpaYpYCSpmc3krJFvSMSYWpvUKj9lgtFguVlZWMjY0xMjKSdIW/yclJqqqqopYj4nHVm2T3Z7VeKujRj8VUTg3Q24tbXV1tWMVAq9UatZKp7h12dnYyOTlJY2MjJSUluju0bDZb0tZZ6wwZPZKyTeXUAD2n2BUVFfT39xuSkV9bW8vg4GB0jZvOnmcgEKCtrc2QRPbpfTgTResSmHrMDnJWObNpzaknSil8Pp/ufSYbGhro7OycYeUiMa7pkO75iZBKBJWWN7sVK1Zovt6EHFbObEKvG0VjY2PUeqXqcCouLqampiYavB0LEYnb2i/d6CQjbqKpTFH7+vpobGzUxF9w5swZXUq9mMqZpZSVlc2YEvp8vpSsUKRtQ0dHB9XV1ZSXl886pra2Nq6nM10vtBE5t6ko2Pj4OG1tbfM2vEqUl156id///veaXCtCzirnZz7zGd2CnZMlojhNTU3U1tZSXl4+r3d1Lu9iUVERfr9/liMlWSvmdrtnWJVIW8GGhgYqKiqoqqqitLR0Ts9lJJUrVUZHR3Uv4J3ODeTy5cvR6KN0GBoaYu/evZo6v3J2n7Oqqoo//dM/5dlnn8149bxQKDRry8Fms1FWVobFYkEphcViwWq1YrFYEBGsVmvMqVB+fj4WiyVmTmVfXx8iMu9UsaSkhMLCQjo6OnC73TPeiyQMR7BarfN6gt1ud8rTNqvVqmun7uk9OFNFq8iuyclJWltbue666zS5Xs5aTpiKPvnGN77BjTfemGlRZhEIBBgYGODy5ctcuXKF/v5+enp66OrqorOzk/b29lndnJ1OJ/n5+XG3BcbHx+ec2hYUFNDQ0MDg4GA07Wq+KV8iWzTpbOPU1dVpFuFkt9uxWCwUFRVRVVVFfX19NEg9Hbxer2bT23Tr8k4np5UTpv5hDz/8MLfeqkkheUOZPvW12WwUFxfPu30Ry2XvcDhoamrC7/fPsnBaBC90d3en7BjSMjrI4/EQCoUYHh6mt7cXr9fLhQsXqKmpScv72tDQoFmEl5bx1Tk7rZ2OiLB582beeeedTIuSFJF2EF1dXVRUVCSUm9rb24vFYiEUCmGxWKivr6evry9ujqQWa6BQKERlZWXSubMioplyVlVVxU0e7+zsxOPxMDk5GdfrfC0ejydaWyldT+uGDRtobGyksrJy1jIiHRaEckLmci3T5cqVKwkrJkyFndXW1mKxWBgcHJw35jbZwmDxSKV8SVlZmSbJ25E92Lmm15GKgqWlpdEZSLzPXlBQQE9PjyZWbvXq1WzdujXt68RiwShnJntapEqkMkCyFqm7uzvhL9bQ0FDU0qZDKje/oaGhuI1/56OsrAyr1UpBQUFSQf8+ny+6fVNdXY3D4aC7u3uG17q6upqrV69qUkXj0qVLBAIBXSKEFoxy5prlbGxsZGhoKKUvbjKKppSipKQkrf1Gu92eUjBBIBCgr68vWu09Ga9qYWFh2tPNiPLZbDbq6+sZHx/H5/PR1tZGMBhMuxEWEI0dNsP35iCXLGd9fT1tbW2GFcVO58ZVX1+Pw+FIS1EiKWfJrMcsFku0+1u6BAIBvF4v/f39iEh0eqyUSqtrwM0338zXv/513bKSTMtpIMXFxYRCIcNbu6fqadUy5zHSWjFRIg6u+vp6Tf9e06e3g4ODuN3u6F50MrODO++8kzvvvFPX713umJt5iBQozlbKysqiNW+MJpXAcJjygmrpfUylrIzX650zLjhd+vv7o1XxGxsb520ibLFY2LZtG3fddZfuBmHBWM79+/cTCoWw2+3YbDbDqgckSnFxseY9SPQmGAymrNixKCwsTCmaR88II/iDNY18XyJ1dbu7u6PW1Ol0Ul1dzSc/+UmWL18e91pakr2mJgkOHz4cLfg7OTkZLWmYam0ZLSktLaWuri6jitne3p6yknV0dFBUVJS2DLW1tSlNkfPz8zWNukmEzs5Ourq6osEdkdjbtrY2Q1LgIuS8cvb397Nr166Y70UyMTKJy+WaUcG8sbFRk0DrVOQQkXmnbbGIlcmSLFarNaXtHC2n1ckyMTFBa2srbW1tTExMREMHjSLnp7X79u2L+55SKrrJr0c18EQIBoNYrdZok55IQnNjY2PUBd/d3U1VVZWujqK8vDzq6upSalI0PSopHjabDbfbjdPpZHx8PBrg39XVhYik3HhKr3KjqVBaWmqoXyOnlfPy5cucPHlyzmMmJye5evUqxcXFCVeQ05Lu7m6WLFlCe3v7DGfQtevhyD+9oqICp9OpeabN0NBQypXwxsbGqK+vp6Ojg7y8vFmOHYfDgdvtjilzfn4+dXV1fPzxx0mPa7FYEm5kZASplENJh5ye1r722msJHTc2NpbylC4dioqKKCsr49KlS/NmdkTWM8PDw7pYi5GRkbSuGync7Ha7cTgcVFRU4PF4aGxsxOFwxL2ZjI2NpezQaWxszKpW9cXFxYaOl5OWUynFqVOnOHv2bMLnRCyGy+WiuLgYu93O6OhotKy/1WrVbFrpcrkoLy/H6/UmtXfW3d1NTU2NbhY+naJWEYdSV1cXgUAgqWWC1+tNes/U4/Fw6dIlnE6nYRXv5+PcuXNMTk4a1lgqJ5VzZGSE7du3p3Tu6Oho3PSgVBv5XIvT6UzJM+l0OpOuE5sM6eRlRm5uqWS5hEKheW9S0wtzl5eXR7NLSkpKsmZqOzIywpEjR9i0aZMh4+XktLawsFAXr1k6oVwRnE5nypkYWnhF9aCioiLtUMORkZG4m/Zut5urV69Gy6dMTk5Gp7PZFlhy4MAB3fddI2TXJ08QEeE73/kOGzdu1DRKY3R0lNra2rT2RycmJuat6havzKUWXaWdTicVFRUx10ep/q20qBLg8/li3vwsFgvBYDBaPuXy5cszPMpaBkFowejoqOaFvOKRk9NamNoauP/++1m7di07d+7UJAZ0+jqqtrY2ZY9pW1sbTqeTxsZGfD4fJSUljI+PR8sxtrW1RR8jOByOGVNqEcHhcMxab0WcTJEA7kjQRcTaTExMRM8pKCigtLQUq9VKb29vyso5MTERbSKbKhaLJeZ2SkNDw5wBGq2trdTV1WG1WlFKEQgEormbmeLgwYNs2LBB9xuHJOKwEJFS4CfAKkAB3wLOAs8BS4BLwB8ppQauOa8B+FegBggBTyilHp9vvJaWFnXkyJGEP4RSihMnTlBUVITFYuGll15KOwPf5XIxOTmp6RRm+n5rRPkbGxsREbq7u5mYmCA/P5/KykpGRkai02OHwzGjH2Uq2Gw2GhoauHjxYkrnV1dX09/fn/K61eVy4ff7Zyh4KvvPWjcgSpW77rqLO++8M+3riMh7SqmYFakTndY+DuxSSq0A1gCnge8BryulbgBeDz+/lgDwfyilbgRuBf5CRFYm+wHmQ0RYvXo1S5cupampiZtvvjnta46OjmoeqjX9i2i1WqPWs7W1lfz8fJqamqKhh9PXrX6/n6GhobSC5gOBQFoOoZ6enrTWxKOjozOmtakmgCuldC+1mQiHDh3SvT/svMopIsXAHcBPAZRSfqWUD3gAeDp82NPAtmvPVUp1KaWOhn8fZkqp9Usx+IPMmlynr69Pt+5efr+fYDBIU1MTVVVV+Hw+WltbdQ3tS7efaLp5i52dnTQ1NQFTFjCVqCGv16tru8W5cLlcUQfVxMQEhw4d0nW8RCznMqAPeEpE3heRn4iIC6hWSnXBlBICVXNdRESWAOuAw3Hef1REjojIkXQDnbVqATA2NqZbulKkTOb4+LghWwVWqzXh4lfx0MJz2trayrJly9LKGJqcnNQsETtRrFZrtN5wXV0dDQ0NugfkJ/LXtgHrgR8qpdYBo8SewsZFRAqBF4G/VkrF3GFXSj2hlGpRSrUku6UxvVziRx99pKnDoK2tLXq315rITWS6RdLrH15ZWZl25JFWbQgnJibSvoH29/cb0vczQn19PSMjI0xOTtLR0UF7e7vuSRWJKKcX8CqlIhZvO1PK2iMiHoDwY8zbv4jYmVLMZ5RSv05f5JkopZicnMRisZCfn8/+/fs5f/68pmO0trbqkn4WKTc5fbqpVy9LLaJavF5v2tPuoqKiuCUuk8Hv91NTUxMto6InBQUFMT33et20I8yrnEqpbqBdRCIZpluAU8AO4JHwa48Av7n2XJla/P0UOK2U+r4mEs8eA6fTicvloqioiK9//eualcOfTroZ+fX19TidTsrKyigoKKCqqirmGk6LQIhYaBWvOzw8zPLRUR46eBBSuKaWweNerxev1xsNx9SLSGDEdCI1g/Uk0XnBd4FnRMQBXAC+yZRiPy8ifwK0AQ8DiEgt8BOl1FZgM/DvgRMicix8rf+qlNqp3UeYicPh4Ctf+Qo7duzg+PHjml67o6Nj1v7kXHg8Hux2e7RukNVqZXx8HL/fH7dFgV4pUulOI22hEOsGBli6bx83vvkmAMfXr+dcEk4iu92uS3jiyMiIbp2/y8vLY27d1NXV6R5jm5ByKqWOAbH2YrbEOLYT2Br+/XeA4ZW3rFYr27Zto6ioiAMHDmh67ba2tpjbAJGon5GREUpKSqIb/0VFRVGHTzAYnNf93t/fn/aGfyzSUU5rKMQ3nnuOumsSDZZ6vZy7/vqEr+PxeHQpHeNyudL2RMcj3pRZ7ykt5HCE0HyICJ/61KcoKiqKWykhnWtHsNls1NbWzvjSTQ8/8/v9SUUbBQIB3G532p7Va0nVIjdNTLDmww9nKSZAeWsrJKicFotFt1KgPp+Puro6Tday05nr/7ZkyRJNx4pFTsbWJsPGjRu5/fbbNb1mXV0d1dXV0VzGuaxBMBiMRgIlih6b26kop2dyki/+8z+zbseOmO9/uGFDwteqr6/XLRUuGAzS0dGhqTUTkbhLDxGhoaFBs7HiseCVE+Cee+7h29/+NmvWrEnrOjabjaamJvr6+ujp6aGtrS3h9najo6MJJ3vrkZ2S7JqsKBTiCz/7GYVzlDXJTzC0sbCw0JAeqhGvuhb7sQ0NDXEtfW1tre4eYlgkyglTsaHbtm3jC1/4QsrX8Hg8tLa2Jm3ZItEwiYbftbe3a77JnozltIZCPPzCC7jn2S/e8sQT/OU//zMt86TIlZaW6rZFdC1erxe3251WN26HwzFnYIgR601YwGvOeKQSn1paWkpeXl7KAdc+ny/paaXWBYsTtpxK8fCbb9Jw+vS8hzr8fhz9/Wx84QUK77+f0t5eiru6GGhqonXpUk4VFhLIQB2giCMu1V4oNTU1cy5VjFhvwiJUzhtuuIE33ngj4Tt5JNMlnS2AkpKSpKveab01kOj1PnvyJMvDNYATxd3VxZ1PPRV9vvToUdYDn7Xb6bjpJp7eti2p62nB8PBwNGUvGYqLi+csVyMihpU2XTTT2ghut5u1a9fOeUxDQwOVlZU0Njbi9/vT7jGZSmSN1msa+9AQ//V//S/+0z/+I9/+2c8oiaGsd7W20vzii8ldeC3wbPhn7TVjTk6y5NgxijNU3jKV2kORfjbxqKmpMSwBfNEp5/j4eNxympEE6fb2dvr6+qLFhLWgr68vqYoCWgcjhCwW7IEABVevUt3Wxl/9z//JJ65exR4K0eD389DBgzOs37w8DBwA/pEppVwLfCf2oUsz0B/G7XYn3X8zkdrBRq03YRFOa8+dOxdzw9rlckVzKfVgbGyMoqKihOvnal0759pJvADb/umfyE912+avgXIgUsppCPiX2Ife9uyzNF+TgzleVsYv77svtbETIJX0tkRuiEatN2ERKueqVasYHBzk9ddfB6aUMhAIUFRUpHue4PDwcDR6SK9ws3jEihBKSTHXMmUhXwQeCj+uYEoxj8U+xd3VBdd4ficcDtBROb1eb7SkS6SUS4SIs81qtWK326Md6hIpjWpkK41Fp5wiwm233UZtbS15eXl4PB52797N4cMx00w1Z3BwMKGek4nunybK0nTD29YypZRF/GFtuTn1yzn9fgpDIUZ0rK43PDycUvuJeFRXV6e1RZMsi27NGWHZsmXU1tYiIqxbt87Qsb1eLyUlJXMe09/fr+ldui7dgPPvAJFAq7eJO4VNhsosquaeCEauN2ERK+d0qqurefTRR/nGN77BJz7xCUPGTCRaqK2tTbNKDO50WxBOn7bOMYVNhps//JDmgQGuGx+nNBgEjSpY6IWR601YhNPaeESySpqamjh79iy7du3SLVAbSLg6YFdXV1JpavEoSbdi3bHw41rgvwBfSe9yAOt27GD6nOVHf/d39BhY3SBZjG7daFrOGCxfvpzvfOc7uiRtRxgZGUkoRC8UCmliQa05NoXMNiorKzUprp0MpnLGwW638+Uvf5kbb7xRtzGSqYGTbkaHaLFv+hhT683H0r9UrmH0ehNM5ZwTm83GF7/4xbSzWbQgnfZzZYEARVo0Dz4GPIom681YZPOKMxPKmb0T/CzBYrHw2c9+ljNnzmjeis7n8yEiCVUpmJycnLXmiezXhUKhOYPyG65exZklbfRyFVM5sxSbzcaKFSv44IMPNL3u0NBQwhn882V2zOU0KtJ4z3SxUVxcrEtXu/kwp7UJctNNN+lyXa16sbS1tcWtBjdXwnRWoXGanFZUVc1ZL103TOVMkGXLlqXdjiAWkc5jWtDR0RGzvWB+mk2djMJqcEhjorjd7oyMaypnglitVt08t1pVlVdKRZW9oaGB+vp66uvrc0Y579+xA0uG0svmQq9awvNhKmcS3HbbbbpFELW2tmpSNCoQCNDW1kZ7e3u06HKejsEUWtJ48iR36ZQVlA6mcuYA5eXl/NEf/ZFurdD1qB0E5IxyAuRlofPKnNbmCB988IFuVdkBuru7NS+7mK9TSUo9sGXZlk9hYaGhmSjTMZUzSVauXKlrWcTInmVVVZUm5f5LAgEKcsVbC9g07CSuBZma0oKpnEmTl5dnSAB0b29vTM9rMtiCQb707LNYs9DJEg9rlilnpqa0YCpnStxwww2GjNPX15f6+lYpvrRnDx6N2yHqTbZNa03LmWMsW7bMkHHS6ay97ehRrjeouoOWZJvlNJUzx6ioqDAsnCuVIth3t7ay5uWXdZBGf6ym5YxiKmcKiIiuuZ7TGRgYoKamJuHj1w8MsPnpp3WUSD9CIgQNqgmbCPn5+YbncE7HDHxPkWXLlnHs2DFDxkq0NcOysTHu/fGPc8YBNJ6Xx3sPP8xAWRm9hYV02mwEdSz4lSyZtJpgKmfKXH/99Qmne6VLV1cXJSUlc5Y2cQcCPPDEE+Tp1ERWD7pWrmSvQTOQVMi0cmbPbSrHyM/PN7SmzFzV+vKCQR5++mmKBwYMk0cL2gyuepgspnLmMEZV6oMp6xkr+EFCIb60YwdV6RbwMpiQCO/rEKqoJZnc4wRTOdNi+fLlho01OTkZ0zH0xYMHWaJxErgRdH3iEwxarZkWY05ywnKKSKmIbBeRMyJyWkQ2iUi5iLwmIufCj2Vxzv2ZiPSKSOzuQTlMRUVFWrV9kuXaNef9Z8+ycu9ew8bXkvbm5kyLMCdOpzMj1Q+mk6jlfBzYpZRaAawBTgPfA15XSt0AvB5+HoufA/enKWdWopTSvG3CXAwODkat56beXm559lnDxtaa4wYFcqRKZWWl5g2Mk2Ve5RSRYuAO4KcASim/UsoHPABENtSeBrbFOl8ptR9Ir8FlljI0NGRYO/UIFouFG0dGuPvJJ8nOoh7z09vQQFcWF4+GzK83ITHLuQzoA54SkfdF5Cci4gKqlVJdAOHHtAqtiMijInJERI709fWlcynDOHPmjOFjdnZ2csf27dhzuEi0t6Ul0yLMS6bXm5CYctqA9cAPlVLrgFHiT2FTRin1hFKqRSnVkg1/mPkIBAIcOHAgI2O3btqUkXG14pSBXu5UyYbvYCLK6QW8SqlIFPV2ppS1R0Q8AOHHuWs3LjCOHTumaXu5ZNh//fWM61BszAh8FRWczwHZc0I5lVLdQLuIRPYNtgCngB3AI+HXHgF+o4uEWUgwGMyY1QS4arVy7q67MjZ+OnhbWrK2BGYEi8Uyb4tGQ+RI8LjvAs+IyHGm+kz9PfAPwL0icg64N/wcEakVkZ2RE0XkWeAQsFxEvCLyJxrKnxEsFktSfU704EBzc1bFoSbKOZ3q/2qJ0+nMuKcWEoytVUodA2Kt4rfEOLYT2DrtuQbN4rILEWHTpk28nMG0rB67nYsbNuRUzuaoy8WJBPqSZho96hOnQu7derOE1atXJ9QAV0/eveOOjI6fLB3r1qFywNqbypnj2Gw2Nm7cmFEZPnK56MgBz2eEi6tXZ1qEhDCVcwHQ0tKiayW+RDh+330ZHT9RJm02jmXBxn4iOLMk4dtUzjTIy8tj/fr1GZXh9+XlXMlQo51k6LzxRsZzYEoLpuVcMNx66626VYBPCIuFU/dnf+iyd+3aTIuQMKblXCCUlJSwatWqjMrwdlMToxmsdZMIJzLQfDZVTMu5gNiU4XA6v9XKR1tm7WplDb0NDfRkeaD7dEzlXEDU1NQYVo0vHgfWrCGQpQrQmeF1ebKY09oFxic/+cmMjn/ZauV8hmWIx9kc2u4B03IuOJYuXZpUfVk9OLx5c0bHj8VwcTFnCgoyLUZSmMq5wBCRjFvPi04nbTffnFEZrqVj3bqsD3S/FnNauwBZuXJlxrMZjmWZY6g1y24WiWBazgWI1Wrl1ltvzagM75eU0Ftfn1EZIkza7bxfXp5pMZLGVM4Fyvr16zP7zxXJmqCEzhtvZCJHooKmY05rFygOh4OWDNfIedvjYTgLkoXbcygqKILdbseaJfV0TeXUgY0bN2b0HxyyWjmTwYD4kAgXNmzgvRyKCoqQLVYTzEZGulBYWMjGjRs5ePBgxmR4e8UKVjscOA1sRjtpt3Pujjt459ZbabfbDRtXS7JlvQmmcurGPffcw8WLF+nq6srI+MNWKx/fcQc3GVAR/mpBAWfvvZe3V69mIEumhKmSTcppTmt1wmq18tBDD2HPoAU5dMsthHTcY/SVl3Poa1/jn/7mb9ixbl3OKyZkl3KallNHKioq2Lp1K7/5TWYKE3Y4HFxav55l772n6XV7Gxo4uXUrB6qrCeWgN3YuzDXnImLNmjVcuHCBEydOZGT8o3ffrZlytt10E8fuvZf3S0pyLuonUUzLuYgQET7zmc/Q3t6Oz+czfPwPCwvZfN11eM6fT+n8oMXCxVtu4d077+Sj/HyNpcs+sslyLqw5SZbidDp56KGHMlYL9WQK2yp+h4MP77uPn/7t3/LM/fcvCsUE03IuSurr67n77rt54403DB/7UGUlG9xuSvv75z3WV1HB+bvvZv/KlQwtsPVkIpjKuUjZvHkzFy5c4NKlS4aOqywWzvy7f8etzzwT8/2A1Urb+vV8+MlPcrSkBBahUkbIpmmtqZwGYrFYePDBB/nRj37E2NiYoWPvW7aMNfn55E8b93JNDefvvJN3PvGJBbENogWm5VzEFBcX8/nPf57nnnvO0HHHrVY+uucelu/dy6VNmzixYQOnCgoWrNc1VUzlXOSsWLGClpYWjhw5Yui4b27cyJ4NG7hq6Ki5RTZNaxfv4iLD3HfffVQZXAx6MBSioqHB0DFzjWyynKZyZgi73c5DDz1keCvBkZERQ8fLNUzLaQJAVVUV9xmc2jUwMJDxQmTZjKmcJlFaWlpYvnz5/Aea6E62NM2NYCpnhhERPv/5z1NUVGTYmN3d3ZSVlRk2Xq6QTVYTTOXMCgoKCnjwwQcNHTPTjX+zkWxyBoGpnFnD0qVLue222wwbr6Ojg/xFEi+bKKblNInLXXfdRV1dnSFjhUIhw7dysh3TcprExWq1snXrVsPG6+3tzZpKc9lATlpOESkVke0ickZETovIJhEpF5HXRORc+DGmh0FE7heRsyLysYh8T1vxFx61tbV4PB5DxhobG6O2ttaQsXKBnFRO4HFgl1JqBbAGOA18D3hdKXUD8Hr4+QxExAr8APg0sBL4iois1ELwhUxzc7NhYw0PDxs2VraTc8opIsXAHcBPAZRSfqWUD3gAeDp82NPAthin3wJ8rJS6oJTyA78Kn2cyB6tWrcLhcBgyls/nM61nmFxccy4D+oCnROR9EfmJiLiAaqVUF0D4MZZ3oQ5on/bcG35tFiLyqIgcEZEjfX19SX2IhYbT6TS0lX0wGDRsrGwm5ywnU5kr64EfKqXWAaPEmMLGIVa4hYp1oFLqCaVUi1KqpbKyMsHLL1yMbOnQ09OD2+02bLxsJRctpxfwKqUOh59vZ0pZe0TEAxB+7I1z7vQ0iHqgM3VxFw8ej8cwxxBk3xczE+Sc5VRKdQPtIhIJAN0CnAJ2AI+EX3sEiFWc9V3gBhFZKiIO4Mvh80wSwEjHkNfrXfRRQ9l2g0rUW/td4BkROQ6sBf4e+AfgXhE5B9wbfo6I1IrITgClVAD4D8Bupjy8zyulPtT0Eyxgbr75ZsMcQwDlOdhLU0uyzXImlEyolDoGxFoEzWqjrJTqBLZOe74T2JmifIsah8PBzTffzHsaV2yPR1dXFw6HA7+BzY+yiVy1nCYZwsip7eTk5KLO9cw2y2kqZ5bj8XgM3Ye8cuVKVuU0GompnCZJY6T1HBkZMSz4Ppuw2+1Ysqxeb3ZJYxITIyOGAMNr6mYD2bbeBFM5cwKHw8Hq1asNG+/y5ctUV1cbNl42kG1TWjCVM2cwcmoLZN0UT29My2mSMjU1NYauBbu6uigtLTVsvExjWk6TtDDaehpZdCzTmJbTJC1uuukmQ+/wi6nOkGk5TdLCaMfQYqozZCqnSdoYPbVdLHWGTOU0SZvq6mrq6+sNG2+x1Bky15wmmmC09VwMzY9My2miCUY7hgYGBgxN/M4EpuU00QS73W6oYwhAqZjVZRYMpuU00QwjawzBwm9+ZFpOE82oqqqiweAu1Qu5jIlpOU00xWjHUEdHBwUFBYaOaRSm5TTRlJUrVxr6pQqFQizUsqWm5TTRlEw4hnp6ehZcUEJpaWlWfqaECnyZZC/Nzc38/ve/N2y88fFxGhoaaG9vn//gLENEqKyspKamhurq6uijy+XKtGgxMZUzx4k4hoxUlqGhIcPGShWn0zlDCWtqaqisrMRmy52vfO5IahKX5uZmQ5VzcHCQ2tpaOjuzo3h/aWnpLEUsKSnJ+UJlpnIuAFauXMmuXbsYHx83bMxQKGTYWBGsVitVVVUzlLC6ujorPa1aYCrnAsBut7NmzRoOHz48/8Ea0d3dTXl5OVeuXNHl+gUFBbOsYUVFRVY6bvTCVM4FQnNzs6HKCeByuTRRzoqKilmKWFhYmPPT0nQxlXOBUFlZSWNjI21tbYaN6fV6cblcjI6OJnS83W6f4SWtqamhqqrK0LKfuYSpnAuI5uZmQ5VTKYXb7Y6pnMXFxbMUsby8fNFbw2QwlXMBEXEMGVkUuq+vL2oBpyviQg3zMxJTORcQNpuNNWvW8M477+g6ztKlS1m1ahUejyfn9g5zCfOvusBobm7WRTnz8vJYu3Ytzc3NZot6gzCVc4HhdrtpamqitbVVk+s1NDTQ3NzMypUrsdvtmlzTJDFM5VyANDc3p6WckRKcLS0ti65nSjZhKucC5MYbbyQ/Pz9px1BNTQ0tLS2sWrUqK1OoFhumci5AbDYba9eu5dChQwkdu2rVKlpaWqitrTW3OrIIUzkXKM3NzXMqp9vtpqWlhdWrVy+algu5RkLKKSKXgGEgCASUUi0isgb4EVAIXAL+WCk1K5dIRP4K+DNAgCeVUv9bE8lN5qSiooIlS5Zw6dKl6GsWi4WVK1fS0tJCY2OjaSWznGQs591Kqf5pz38C/I1Sap+IfAv4T8D/Of0EEVnFlGLeAviBXSLyqlLqXJpymyRAc3Mzly5dorS0lObmZtatW5e1icUms0lnWrsc2B/+/TVgN9coJ3Aj8I5S6iqAiOwDHgT+MY1xTRJkxYoV/PEf/zHXXXedaSVzkERrCClgj4i8JyKPhl87CXw+/PvDQKw6jSeBO0SkQkQKgK1xjkNEHhWRIyJypK+vL/FPYBIXm83G9ddfbypmjpKocm5WSq0HPg38hYjcAXwr/Pt7QBFT09YZKKVOA48xZVl3AR8AgVgDKKWeUEq1KKVaFmqFNxOTZEhIOZVSneHHXuAl4Bal1Bml1H1KqWbgWeB8nHN/qpRar5S6A7gCmOtNE5MEmFc5RcQlIkWR34H7gJMiUhV+zQL8HVOe21jnR45rBL7AlCKbmJjMQyKWsxr4nYh8APweeFUptQv4ioh8BJwBOoGnAESkVkR2Tjv/RRE5BbwM/IVSakDTT2BiskCRbOwe1dLSoo4cOZJpMUxMdEdE3lNKxexKZVZ8NzHJUkzlNDHJUkzlNDHJUkzlNDHJUkzlNDHJUrLSWysifYA2dTaSxw30z3tUdmHKbAx6yNyklIoZEpeVyplJRORIPNd2tmLKbAxGy2xOa01MshRTOU1MshRTOWfzRKYFSAFTZmMwVGZzzWlikqWYltPEJEsxldPEJEtZFMopIg0i8qaInBaRD8MVAae//zciokQkZhMQESkVke0iciZ8jU05Ivd/DJ93UkSeFRHd+7PHk1lE/ruIdIjIsfDP1jjn3y8iZ0XkYxH5nt7ypivzfP+jtFBKLfgfwAOsD/9eBHwErAw/b2CqOFkr4I5z/tPAn4Z/dwCl2S43UAdcBPLDz58HvpEpmYH/zlS1xrnOtTJVUWNZ+O/8QeTzZrHMcf9H6f4sCsuplOpSSh0N/z4MnGbqywvw/wL/makiZrMQkWLgDuCn4fP9Simf3jKHx0pZ7jA2IF9EbEABU0nxujKPzPNxC/CxUuqCUsoP/Ap4QB9J/0A6Mqf5eedkUSjndERkCbAOOCwinwc6lFIfzHHKMqAPeEpE3heRn4TLtRhKsnIrpTqA/xtoA7qAQaXUHiNkjTBd5vBL/0FEjovIz0SkLMYpdUD7tOdeNPqiJ0oKMs91blosKuUUkULgReCvmaoC+N+A/2ue02zAeuCHSql1wChgyFooQipyh79IDwBLgVrAJSJf01fSGeNHZVZTnQB+CFwHrGXqZvH/xDotxmuG7fWlKHO8c9Nm0SiniNiZ+uM9o5T6NVN/9KXAB+F2E/XAURGpueZUL+BVSkXuhtuZUlZDSEPuTwEXlVJ9SqlJ4NfAJzMkM0qpHqVUUCkVAp5kagp7LV5m1jWux4CpOKQlc8xztWBRKKdMVVX+KXBaKfV9AKXUCaVUlVJqiVJqCVNfjPVKqe7p54aft4vI8vBLW4BT2S43U9PZW0WkIHydLUythwyXOfy6Z9phDzJVcPxa3gVuEJGlIuIAvgzs0FPesGwpyxzvXE3Q2xOWDT/AbUxNj44Dx8I/W6855hJhrydT08Cd095bCxwJn/9vQFmOyP0/mKqOeBL4/wBnpmQOj38i/PoOwBNH5q1MeTzPA/8tk3/nRGRO5H+U6o8ZvmdikqUsimmtiUkuYiqniUmWYiqniUmWYiqniUmWYiqniUmWYiqniUmWYiqniUmW8v8De1HNfXVjd7UAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Create a figure with one subplot\n", "fig, ax = plt.subplots()\n", "\n", "# Plot polygons\n", "polys.plot(ax=ax, facecolor='gray')\n", "southern.plot(ax=ax, facecolor='red')\n", "\n", "# Plot points\n", "pip_data.plot(ax=ax, color='gold', markersize=2)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Perfect! Now we only have the (golden) points that, indeed, are inside the red Polygon which is exactly what we wanted!" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.12" } }, "nbformat": 4, "nbformat_minor": 4 }