{ "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": 3, "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": 4, "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": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 5, "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": 6, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 6, "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": 8, "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": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Does polygon contain p1?\n", "poly.contains(p1)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 10, "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": 12, "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": 13, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 13, "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": 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": [ "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": 15, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 15, "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": 18, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 18, "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": 19, "metadata": { "collapsed": false, "deletable": true, "editable": true, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 19, "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": 20, "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", "
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": 20, "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": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'AeronavFAA': 'r',\n", " 'ARCGEN': 'r',\n", " 'BNA': 'rw',\n", " 'DXF': 'rw',\n", " 'CSV': 'raw',\n", " 'OpenFileGDB': 'r',\n", " 'FlatGeobuf': 'r',\n", " 'ESRIJSON': 'r',\n", " 'ESRI Shapefile': 'raw',\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", " 'SEGY': 'r',\n", " 'SUA': 'r',\n", " 'TopoJSON': 'r'}" ] }, "execution_count": 21, "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": 22, "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": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'AeronavFAA': 'r',\n", " 'ARCGEN': 'r',\n", " 'BNA': 'rw',\n", " 'DXF': 'rw',\n", " 'CSV': 'raw',\n", " 'OpenFileGDB': 'r',\n", " 'FlatGeobuf': 'r',\n", " 'ESRIJSON': 'r',\n", " 'ESRI Shapefile': 'raw',\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", " 'SEGY': 'r',\n", " 'SUA': 'r',\n", " 'TopoJSON': 'r',\n", " 'KML': 'rw'}" ] }, "execution_count": 23, "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": 24, "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": 25, "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": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Check the data\n", "print(\"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": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM8AAAD4CAYAAABL/rJKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3HUlEQVR4nO2deXhjV3nwf0e7ZUuW5H23ZzJLZl88noRAWMKSBpoAAQoBmrI0LXu3r1/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": 27, "metadata": {}, "outputs": [], "source": [ "# Select data \n", "southern = polys.loc[polys['Name']=='Eteläinen']" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "# Reset index for the selection\n", "southern.reset_index(drop=True, inplace=True)" ] }, { "cell_type": "code", "execution_count": 29, "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": 29, "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": 30, "metadata": { "deletable": true, "editable": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOYAAAEYCAYAAABIhYpmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+wklEQVR4nO29eXBb55Xg+zsAsZAguBPgApKSLMmyZVmiRDuyZct2FDu2Omm7o9ixM+5WJp2Xqk5PKpN0v0le93uvZv6Yrp433emkupJ0HHsSx1kcW1HseIklx21Z3mRL8iJrtXbu+yIS4Abge3+QwIAkAGK5uADI+6tikQDuvd8BgXPP+c53vnNEKYWBgUFuYcq2AAYGBgsxFNPAIAcxFNPAIAcxFNPAIAcxFNPAIAcpyLYA0aiqqlIrVqzIthgGBhnn6NGj/Uqp6vnP56RirlixgiNHjmRbDAODjCMil6M9b7iyBgY5iKGYBgY5iKGYBgY5iKGYBgY5iKGYBgY5iKGYBgY5iKGYBgY5SEKKKSJlIrJHRE6LyCkRuUlEKkTkZRE5O/u7PM75ZhF5X0Se1050A4OlS6IW8/vAS0qpdcBG4BTwHeAVpdQa4JXZx7H4xuw5BgYGCbCoYopICbAdeAxAKTWllBoG7gUenz3sceC+GOd7gD8BHk1fXAMt8fv9+Hy+bIthEIVELOYqoA/46aw7+qiIOAC3UqoLYPa3K8b53wP+CxCMN4iIfFVEjojIkb6+voTfgEFyBINBLl68yLPPPsu//uu/4vV6sy2SQRQSyZUtADYDX1dKvSMi3ye+2xpGRD4D9CqljorI7fGOVUo9AjwC0NLSYtQ70RClFN3d3Xz00UccP36c0dFRAGw2G8XFxVmWziAaiShmO9CulHpn9vEeZhSzR0RqlVJdIlIL9EY5dxvwpyKyE7ADJSLyC6XUw1oIbxCfsbEx3nvvPT766CP6+/sXvD45OYnP56OwsDAL0hnEY1FXVinVDbSJyNWzT+0ATgK/B3bPPrcbeDbKuf+XUsqjlFoBPAj8u6GU+uD3+/nVr37Fq6++GlUpQzz11FNMTk7qKJlBIiQalf068EsROQZsAv4B+EfgThE5C9w5+xgRqRORFzMgq0ESvPTSS3R1dS16XG9vL7/4xS8IBAI6SGWQKAntx1RKfQC0RHlpR5RjO4GdUZ4/ABxISjqDlPjoo484evRowse3t7fT09NDXV1dBqUySAYj82eJ0dfXx3PPPZf0eYbFzC0MxVxCTE1N8dRTTzE9PZ3UeTabDZcr1mqXQTYwFHOJoJTi+eefjxvoicXk5CQnTpzIgFQGqWIo5hIhtCySKi+88AJtbW0aSmSQDoZiLgG6urr4wx/+kNY1gsEgTz/9NGNjYxpJZZAOhmLmORMTEzz99NOaBG9GR0c1u5YeBINxszzzGkMx8xilFM8++yxDQ0OaXbO1tZV9+/Zpdj0tOHPmDK+//jqRnekmJib4xS9+weuvv55FyTJHTtaVNUiMQ4cOcfr0ac2ve/jwYVwuFy0t0ZauM0cwGOT111/nqquuIhgM0tjYiM/nY9++fQwNDdHb28vdd9/Ne++9x4kTJ+jp6WFiYoKbb74Zs9msq6yZRnKxP2ZLS4syCj7Hp62tjZ/97GcZc+dEhC984QtcffXVix+sAV6vlz179nDp0iVMJhPBYJA777yTI0eOzPEIzGbzAld769atfPrTn9ZFTq0RkaNKqQV3QMNi5iE+n489e/ZkdI6llOK3v/0tu3fvpr6+PmPjhLBareEtaKH39fLLLy84Ltr8t7+/n3fffRe73Y7VasVms2G1WnE4HJSVlWVU7kxhWEyN8Pl82O12TKbMT9t/+ctfcu7cuYyPA+BwOPja175GUVFRxsfq6+vjJz/5SdIJErFwOp1861vf0uRamSKWxTSCP2kQDAa5fPky+/bt44c//CG//vWvMz6mz+fjwoULGR8nhNfr1S0YVF1dzT333KPZ9UZHR/H7/ZpdT08MVzYBlFKMj49TVFSE3+/nwoULnD59mjNnzswpzXHu3DnefvttbrzxxowEIy5evMj7778fnoPpxbFjx7juuutYs2ZNxsfatGkTFy5c4Pjx42lfS0R48sknMZvNlJSU4Ha7aWxsJBgMEggEUEpRX1+PiGggubYYrmwCnDhxgj179lBWVobP52Nqairu8cXFxVitVurq6rjuuus0CaCMjY3x4x//OGsJACUlJXzta1/DZrNlfKyJiQl+/OMfMzw8nNFxbrvtNm6//faMjrEYhiubIoFAgFdeeQWA4eHhRZUSZpRocHCQ48eP8+STT/Lmm2+mJUMwGGTv3r1Zzcq5cuVK1GBMJrDb7ezatSuj8/W1a9dy2223Zez66WIo5iLMD9enwquvvsqhQ4dSSjAHOHjwIBcvXkxLBi04evQoly5d0mUsj8fDHXfckZFrV1ZW8md/9mc56cKGMBQzDpOTkxw8eDDt6wQCAfbt28cPfvADfvvb35JMFcCzZ8/y2muvpS2DVjz33HOaRU0XY9u2baxcuVLTa1qtVr7whS9gt9s1va7WGIoZh1OnTmled/X48eP88Ic/5Omnn6a3N1r9sv/N4OAge/fu1XT8dBkcHOTAgQO6jCUi3H333Zpe85577qG6ekFn9ZzDiMrGQYvIYCxOnjzJyZMnWbduHbfeeuuCsh5TU1P85je/YWJiImMypMrbb7/NunXraGhoyOg4SinN12vffvttTp06xfbt25mensbtdudklUBDMWOg13rh6dOnOX36NKtXr2b79u00NDSglOL3v//9ohY1Wyil2Lt3L1/96lcz9qVWSrF//34OHTpEeXk5TqeTiYkJRkZG0qrq19vbS29vLx9//DEADz30EGvXrtVKbM0wFDMGp06dQs+lpHPnznHu3DlWrFhBY2NjzlcUGB4e5plnnuHBBx/UPIhy5coVXnvtNbq7u6mrq2N8fJzW1tbw606nk4KCgoSCciLCqlWrOH/+fNTXq6qqNJNbSwzFjMHJkyezMu6lS5d0C66ky8cff8wbb7zBrbfeqsn1JiYmePPNNzl06FDcjJ3R0VHKy2M2l5uD2WzmoYce4oknnuDy5ctzXisoKMBqtaYlc6YwFDMKXq83q8sT+ZRGdubMGUpLS7n++utTvobf7+fw4cO8/vrrjI+PL3p8YWFhwktY9fX1mM1mmpqa5iim0+nks5/9bM62iDAUMwonT57U1Y2dT75YTLfbTUdHB8888wwWi4VrrrkmqfOVUhw7doxXX32VkZGRhM8TETweD93d3VFvYtdeey1FRUVYLBY2btxIb28vp0+fxuPxYLPZsNvt3HHHHVRWViYlr54YijmPqamprO+Kz5eWBQMDA9jtdiYmJtizZw8PPPAAtbW1TE1NMTk5yfT0NMFgEKVU+EdE8Pv9jI2NceTIEXp6epIe1+fz4fP5sFgs4cjw4OBgWNGsVit33XUXFoslfM5f/dVfafOmdcJQzHm8/fbb4W5Y2SIXlkhcLhd+v5/p6WnMZjNOpxO/3z+n7UJtbW24sl4wGOTJJ59c9LqNjY1zAjnpMD09HR6/sLAwfN3W1lb6+vr44he/qMt2tUxgKGYEXq+Xt956K6syWCyWnHBl/X4/g4OD4cfDw8MUFBTgcrnCAZNcKnc5f27a0dHBY489xsMPP5xwoCiXyMvMn0wlcx88eDChJPVMkiupYtGCIn6/n97eXtrb22lvb0/6mm63m87OTi3ES4jBwUGeeuqplGTNNnmpmAcOHOAnP/kJhw8f1szt6+rqIhe2moW2VVVVVeHxeFi5cmVWkq0z0QLebreHgzWR64dmszkjVs1isdDT08Pjjz+e8+vC88lLV9blcnH06FE6OzvZv38/69evZ8uWLXg8nqS/xJOTkxw4cIB33303J+qUTk5OUlVVFd6JUltbS1VVFRMTE7rOffv7+6MWvkqHyEh3aLnDarUSDAaZmJjQfAN4VVUVXV1d+P1+9uzZw9DQENu2bcvpXSUh8lYxQ/j9fj788EM+/PBDqquraW5uZuPGjYtO+pVSHD9+nP379+dU9fHR0dE5Cmi1Wunq6sJqtVJfX09HR4cuclRWVjIwMJCRa1ut1vA8OjR18Pv9OBwOvF4vpaWlTE5Opu0Nzd/P+corr9DT08OOHTtyvkhX3itmJH19fezfv58//vGPXHPNNTQ3N7Nq1aoFd8iBgQFeeOGFnNjjuBiRX+COjg4aGhro7u7OeIAoEzmwIkJNTQ0WiyVq4KiiogKLxcLw8DBVVVW4XC4GBgbC1fOSweVyRb2JHT9+nFOnTnHDDTdw66235mzUNi8Vs6ioKHx3jUYwGOTEiROcOHGCsrIympubaW5uprCwkDfeeIM33ngjb9oAzM9waWtro7S0lIKCgoxYNLPZTG1tbcYCJt3d3TFfi1TWkCtvMplwu91Jr3fGc4kDgQCHDh3i/fff5+abb2br1q05l5qXtzV/fv7znydl8ULJzBcvXsyJuWQyFBcXL3C3RSS8JqjVZygiVFZWplxpYTHMZjM2my3pwFJhYWHcVL26ujp6enrCN1uPx5PUjcXhcLB9+3a2bNmie0X3JVfzJ9lGq0opzp8/T3V1tS4FpbSktLR0wXNKKS5fvkxlZaVmEU2Px5MxpYQZS5XKbo5YbrXdbsftdtPd3Y3b7QZmLGwy6X0ws379hz/8gR/84AdJVZfIJMtGMUP09PRgt9spKSnRWKLMEe8u3t/fz/DwME1NTWlHG/WI+qaSbjg4OEhTU1P4cV1dXXh/Zk9PD8FgkM7OTmpqavB4PCm/j6GhIV544YWs5kmHSEgxRaRMRPaIyGkROSUiN4lIhYi8LCJnZ38vuG2LiF1E3hWRD0XkhIj8N60ET6c8xMjICNPT03lRYgJYNOkhZD3TqShQV1eX8XKRQMpzuba2trBVLCgoiKp8gUAg7QSGy5cvp9UAWCsStZjfB15SSq0DNgKngO8Aryil1gCvzD6ezyTwSaXURmATcLeIbE1batJTTJhJ4RoYGMDj8WghTkaJTI2LRyJbpqJhMpl0y3hKdQkkGAzS09MT9/Mym82abJnbv39/1vOVF43KikgJsB34EoBSagqYEpF7gdtnD3scOAB8O/JcNeMThKIWltkfTfyEkDt65cqVlK8RDAZpb2+nvr4epVQ4adtut2OxWMLBhJBro5Sak8StF1NTU5SWli46d+rr60t6/dFsNlNdXR03WqoVWgSW4gV1Cgq0WWTwer28+uqrmrZrSJZE3skqoA/4qYhsBI4C3wDcSqkuAKVUl4hEnfSJiHn2nNXAD5RS72giOTPzzHQUM0Sii/bZDKk7nc6EghpFRUUJK6bFYqG8vFwXpYQZxcnk/E3LiOrhw4fZtGkTtbW1ml0zGRJxZQuAzcCPlFLNgJfobmtUlFIBpdQmwAPcKCLXRTtORL4qIkdE5EiikTG954iR+/v0JtGq5J2dnQnfQKqrq3Ur+FVZWZnS3stk0NIdV0olHAjq7e1dsC6ebq5xIhazHWiPsHR7mFHMHhGpnbWWtUDcT1gpNSwiB4C7gQV1IZVSjwCPwMw6ZiLCpxqZTZVsKmaimT6BQID6+vq4ex5NJpPuy0Z6lIjUeoN5R0cH77//Pps3b455jFKK5557jqmpKW655Ra6u7s5d+4cfX19fOtb30q5dMmit2GlVDfQJiKhzjg7gJPA74Hds8/tBp6df66IVItI2ezfhcCnAM16k+ttMbWaw6RCMsGIaEocKXtJSQk9PT1MTU3R2NhIfX09jY2NGXt/hYWFumz3SjX4FY8//vGPca3foUOHaG9vp7e3l7179/LWW2/hdrv5xje+kVY9oUQ/ia8DvxQRK3AB+I/MKPVTIvKXQCtwP4CI1AGPKqV2ArXA47PzTBPwlFLq+ZSlnYce4f1I9M4KiSSZRPvBwUEKCgrCEUqXy8Xo6CilpaUMDQ2F/2/z59ZWq1XTCgMhKisrM74nUkQyopjj4+Ncvnx5QT0jn8/HsWPHwg2nIrl48WLac+mEFFMp9QGwIG2IGes5/9hOYOfs38eA5jTki4rX6+Wll17KaKX0aGRLMW02W1JuWuTWscLCQsbGxhgfH1/0izs1NaV5cnwo6T4ebrcbm81GR0dHyjnMRUVFKSW7J0Ko1f309DQff/wxx44d49y5czFTO8fGxti3bx8PPPBAykkfeZnEfuHCBd2tJSQegNEah8OR9PypsLAwXKsnmQCP1vNopdQCZQ/l5Nrtdqanp8NBofLyckZHRykuLqa0tHRBHdh42O32jCim0+kMZ4lduHCBPXv2LHqO2+3mT/7kT9LKxMpLxdywYQMbNmxgaGiIffv2cebMGV3GzZZiplJuZHR0lKqqqqQjoVpbzPb2dqxW65yIaVVVVdSc1KGhIUpKSrDZbLS2tiZV/yhTgbny8nKGh4ex2WwcPnx40eNNJhMPPfRQ2vVq8zZXFmb+acnWMs1HUgnKDA8Pp7Q8kW4v0GjMT7KPNyUYHr7C+fNjBIMqqTXETK0xh9rFP/HEEzHbLEQSDAZ55530l+rz0mJGomfjnXzbLpYs6WZSxcLn81FSUhJWnlhzzmAQHn98N21tDTQ0tPF3f/fHhMfIROKCyWTC6/Umna2kRVpgXivmmTNnsl5uUg/02tRtt9uTVkyTyRSutTQ4OBg1udxutye0ncrnc9DW1kAwaKatrYGzZ4dZv74Jn8+36PmJ5hMvRuj9hK6ZrFKKiCbtCfNWMYeGhnjmmWd0HVNEcLvdjIyM6JrkrFdl9mTmaVarlZqaGkZGRsLLK6H54Xx5E00ucDi8NDS0hS2mw+Hl8mUvTqeTxsZGYMZFn5iYwOFw4PP5mJyc1GyJp7i4mPLy8rSupZTid7/7HV6vl61bU9+vkbeK+dJLL+m+AyCy9EVhYWHYPRMRzdf+IsnUMsB8FrPMBQUF4errU1NTC97z2NjYApeytLQ04TVMEdi9+3F8PgcOh5dQUHN+gTKYWdoJRW/Tnc5UVVVht9tpb2/XpDCb1WpNe86bt4qZ7W05keuCmUw3M5lMGVk4j0Ysi1lVVYXZbKa3tzdu9fVgMLhgd0uiyfchTCYoLk7sRjQ2NsbY2FjKFrOuri5cxDpdXC4XN954IzabjWuuuSbtNe+8Vcxc2GUewm63Z0x5gsFgRrJxotHV1YXZbEZEsNlsVFRU4PP50tqqpccSU7KlRAoKCnA4HJqmCV5//fVs2bJFs+vl9XJJrpDpZHC9LKbf76e0tJTa2lrGx8dpb29PuhLfwMAAjY2NYYuRzLpoMAhjYw6SveemsqtE6+jzqlWrNL2eoZgakOnK3nrWJ7Lb7XR1dREMBlOu7tDa2orNZqO+vj7hJabQUsl3v/tNfvaz3SSzMhWtWFk8MnEj1dozMBQzx6mrq8vI2mI0Ghsb6ezsDK/DpZMF5PP58Hq9CSc5zF8q8fkcCY+VrKL5fD5NmzeF3H8tyVvFzKU5ZqZksVgsXLlyhb6+voxvcWtqalowj+3u7k4rtWx4eDjh0pqhpRKTKRBeKkmUZP//SqmUymjG4s4779R8b3DeBn9yiUwpZm1tbVhZUo3yuVyucHOgWNFHh8MRM7hUUVGR1hKC0+lMKM0v1lJJpvB6vUkXho7FG2+8wcaNGzW1wnlrMZc6brd7jrJ0d3en1AgnNGfs7e3F4/FEtYBVVVUxby7ppiEm4w6HlkqSVcpUa9W2t7eHExfSYWxsjF/96lea9sLJW8X8zGc+kzNFm/v7+6mtraWxsZHa2lrKy8sXDQbEW4CuqKiIamWSDXLU1tbOWUpob29nYmKCxsZGKioqqKpyYbU20N8fO/Kabn+UxdYIU43EhjCbzWlVT29vb9fEDW1ra+PQoUNpXydE3rqybrebr3zlK/z617/OSknJSAKBwAIZLBZL2MKJCCaTKfwjIlgslqj7DUtKSpiYmIi6BJBoECiUENDV1bUgMOL3+2ltbV2QML579+NEu5eMj49TXV2d8pff4XDElDtRGeLhcrnS+vyDwaBmUdqLFy/i9/s1KdGSt4oJM/OXL33pS+zdu1e3PZmJMj09HTex2mQyLWgW5HA4UErFrDEzNDREeXl5zDlbaWkpxcXFc0qGxEpZixYFjZVxk05mU1FRUUzFTEYGgJqaGvr7+3E4HNjtdqxWqyZKpVVriOnpaUZGRqisrEz7WnnryoawWq088MAD3HjjjdkWJSmCweAc1zRUZHqxL0k0972wsJDGxkauXLmyoI5PrC9uMlHQZDNrQlit1jnLJfPd1kgZPJ42lCKuSxuqYzQyMkJPTw9tbW2cO3curXmizWbTdA1Sq62BeduGbz4jIyN873vfy4xAGaSuro7BwUGKiooS2rpUVlYWLqtisVioq6ujo6Mj5h7AeOl8wSAJR0GLioqSrpUamTcby20NBsHrdbBnz+fjurSLpSU2NjYyODiYcAS5sbERpVS4WmCq2Gw2brnlFlwuFy6Xi9LS0qTWNGO14ctrVzaSTGffZIrJycmElRJm1gYrKytxOBz09fUtWhcnXlQ0MmF8MSUtLi5OWjEjq8LHcltNppmlkngubSK5wqHX3W532FLHUri6ujrNco937tzJ9ddfr8m1Isl7VzZEturxpENRURGBQCDpTb5DQ0O0trYmlEObyJaxRNLhUqmpMzIygtPpBOK7ztFeq66upqqqipqamqSUKOTiTk35cTpXU1NTu+CmbTKZNFkmgZkCXZnAsJhZoqmpiYGBgZSq/SUzj7ly5Uo4wSAWiwVhioqKUpo7hcaur69naGgoZgJBtOSC6enpRf83saz8zI3mz8Ou8V/91dO4XFUMDAxQXFxMe3s7RUVFSb+faGh1nfksGcXMJ4tZW1ubVGnGdDGZTHEVM1rlACCsVN3d3SkvSQQCATo6OnA6nTgchZhM0S14pFttMplwOp1YLJaYyzTxllrm32j6+4WJiRmrG3JvfT4ftbW1XLlyJaWN6CLCPffcww033JD0uYmwZBQzHyxmRUUF4+Pjuq+7LlYGUgT+4i8ep7+/murqvrD1WawHSiyiWbLR0dGE95UGg8HwhuxY58Sz8rFuNDB3zt3V1UVTUxMTExNJ1VWyWq18/vOfZ82aNQmfkyxLSjFFJKeS2yNxuVx4vV7d9lZGsljgJhiEn/98ofXp6+ujsLAwKZnjWbJUPpvW1tao1fviKV8yebft7e0EAgGcTifl5eV0dXXFvYmVlJTw0EMPUVNTk/R7SYYlo5gHDx5EKYXFYqGgoIDq6mpddv0nis1m07XUZiSLJVfHsj7j4+NJV0+IZ8lSqfYXq7TKYsqXaImSkEyhukIWi4XGxkZGR0fnJHKUlZVRWlrKrl27wgGtTJI/E7M4vPXWW+EyltPT04yPj9Pa2kpdXV2WJZtJjwsVsMoW0W4IkYv98SKmvb29Se1siXWtpqamlEp5uN3uqBYsmTXYZJienqa1tZWhoSFKS0vxeDw0NDQwPDzM2NiYLkoJS8Bi9vb28vLLL0d9rbu7G5fLlTVLBTN3/NCcsrCwMLyvUk9rPjExgcViQSmF1WplbMwXdjc9njbuv39PTOszMTGx6Pao+UoSq9JdKkRbpkkkx1YLxR0ZGZmT9ZTo3lItyHvFfO2112K+FgwGGRkZiZtfmmmsVms4L7asrCyskA0NDeGE5+7ubioqKjLacbm0tBSr1Up/f/8cd7O1tYnvfvebNDbOfMGjfYnjbauKVJKrrurh29/ex+TkOHa7nclJB729vTidzpSbQEVb411seUeL5Pho6KmYee3K9vX1cfLkybjHTE5OhmuQZoOOjo7wnCXSSra1tdHV1UVbWxvT09PhbWButzsjnbInJiYwm81MTU3NcTcBlJpbzmN+TmtfXx+VlZWYzeYFubeRSnL+vJvTpwfo6+ujra2N3t5eysrKKC8vT2kdtLKyMmqK3WJ5vumUKYmHYTETJJYLOx+v10tZWRkOh0O34skw80H6/X4uXbq06LHt7e3U1NTQ29sb7seoJWNjY+EvVsjd9HodPP3052lv/9+RzVjWpri4mIGBAerr6+nr68PhcGC1WjGbC2hs7KC1tT6qkoQ6ZaVCeXl51P2giwV+4kVs00HP/b95qZhKKY4fP87Zs2cTPifkSjkcjnAFda/Xi81mIxgMYjKZFuzMSJWSkhJKSkqSKluhlKK7u5v6+vqMurQhTCZwOr186Utzv+Beb3Q3MRS9DBXrioyU/sVf/DTufG5kZAS3253U+2psbIy7BzRe1DVTZUpOnTrF+vXrtbnYIuSlYl65coW9e/emdK7X641qNc1mc1L9GOMRCARSqiVTXl6u2c0hGtGWHeZ/waNZG4vFEu7QFW0Xy2JLExMTE7jd7riyRXozkYkNJpMpJTc4mYruiXLixAluvfXWRd+LFuTlHNPpdGo+ZwwEApr8w0tKSlJ2l7M1D44kZG2+9a1/4UtfmgkGud3utFvLxVvDbGhoYGJigqamJtxu95w2ffGS59MtS5IKBw4c0GWcvFRMk8nEX//1X3PDDTdomoo3OTlJbW1tWm3UQqlnsWr6iEjMhqxavJfCwkIqKyujKnmi65Hzi2Jp0QYw1Fl6PsXFxXR3dxMIBLh8+TI9PT1zxotVgCyVAtFaKPLp06d1SanMS1cWZrJZdu7cSXNzMy+++KImZQgj5zTJzolCKKVobW2lsLCQhoYGRkdHcTqd+Hw+BgYGwhubm5qa5iSyFxYWzlnrM5lMWCyWBUsV5eXl4UXuYDDI1NRUuEZQKLki5LKGumGJSLgvSSqk2yAHZvKEoy19zC+vMp+BgYE5N0qv18vg4GDSZUm0XEI5cOAADz30UGonJ0hCiikiZcCjwHWAAr4MnAF+A6wALgEPKKWG5p3XAPwcqAGCwCNKqe9rI/oMtbW1fPnLX+bYsWPhUh179+5Nu45LaHkhVWsxPj4ezvYJBZ4qKyvDc8iQuxvaF9jR0cH4+DjFxcVUVlYyNDQUzg8N1bcJBAIMDQ0lvCYb6oYFMymBqSpYe3s7DQ0NaWUvRXNJE7mm3++fc0yoUkGykddkFTkeH3/8Me3t7Sm3kEiERC3m94GXlFKfFxErUAT8HfCKUuofReQ7wHeAb887zw/8jVLqPRFxAkdF5GWlVPzFxyQRETZu3Bh+vH79+rRLCY6MjGjeZSsy9O9wOHA4HGGrWVlZSUlJCRcvXlxgQSYmJtJuOzg5OZlWPZq2tjasVmvKZTh6enrmeAkOhyMlj0REcDqdjI6OJhV51XoJ5cCBAzz88MNpXSMeiyqmiJQA24EvASilpoApEbkXuH32sMeBA8xTTKVUF9A1+/eoiJwC6gFNFTOKzJpcJ5PZQtPT0+Gd9KFWd6FOWZlK10v3/ZSXl4eVKZWUt8uXL4fT+0pKSlKaq12+fJmamhpGR0eTirxqsYQSuiEAnD9/ntbWVs0qIcwnES97FdAH/FRE3heRR0XEAbhnFS+kgHHTVURkBdAMvBPj9a+KyBEROZJOAd9ZedI6P0QokJMJOjs76evrw+v1zuk/qUWgJRplZWVJ1+yZT6heajqduTo7O7nqqqvSCqD09PSkVDkg1UrvMOPRjI6O4nA4aGhooL6+XpO4RiwSUcwCYDPwI6VUM+Blxm1NGBEpBn4L/GelVNQio0qpR5RSLUqplmQb6AQCAa5cuUJvby9nzpzRdIE+k3fFYDA4pziwyWTStJlqJFrsigjJmk7KWzAYjFuVIJGoqVJKt10eIULv3ev10tbWRkdHR1rR+8VIRDHbgXalVMjS7WFGUXtEpBZg9nfULRwiYmFGKX+plEotKyAOSin8fj9msxm73c5rr72maQ8JIGNbyFwu15ybiNlszthGby1Kr7S2tlJTU5NWZ65YbQWTtcI9PT1UVVVlNAAToq6ubsHNPlQ6NFMs+mkppbqBNhG5evapHczMEX8P7J59bjfw7PxzZWay9xhwSin1XU0kXjgGNpstnGq3e/duVq5cqfk43d3dae1a93g8FBQUhEtPVldXLyilPz09rUkV70yhlMJsMtEyPMSeTd+ak4QQi/lWMOSqz38+FSvc398fjhhnCpPJFDVhxOPxaLKMFItEo7JfB345G5G9APxHZpT6KRH5S6AVuB9AROqAR5VSO4FtwJ8DH4nIB7PX+jul1IvavYW52Gw2vvjFL/LMM89w4sQJza4bDAbp7u5OKjhTV1dHQUEBSqlwVHN0dJRAIBAzOyiyFquWpGuJncEg6/v7Wbd3L03HjhEwm/nj329gUube2yODQkrNXTv85jefo6enJ+qaYjpR03SzkuLh8Xiift5NTU0ZGxMSVEyl1AfAgmrRzFjP+cd2Ajtn/34D0L1KVkFBAbt27aK4uJh33okaa0qZWLms9fX1+P1+fD4fZWVliAj9/f3Y7fbwwnoiSw2Dg4M5V7uoNBBg97/9G+URc0NzIMC6sTE+jNhxMV/hPv/5PXOsYDBYCcRODkg1appKzdtEsNlsMeMVK1asyMiYIfIyJS8RRIRPf/rTfOpTn9L8upFYrVY8Hg8dHR309PQwOjpKW1sbra2t+Hw+xsbGktpf6fV6M1KrNFVFX+vz8SevvDJHKUPUzouszlc4IDwXXbGik9HR80Ds/ZSpRk1DzZa0pqamJuom8VBZz0yStyl5iSAibNu2Da/Xy9tvv63JNf1+Px6PB7/fj8ViCc9zYjE1NUVfX1/CLrDZbM5IJb1UFHP92Bj3fu97WKK4igGzmXdXrw4/DgZnGgJ5PG3h/Z2RVnDt2jI6O2ciOlpvywoV0aqtrdUsj7W0tDTm5xWKF2SSJa2YIe666y42bNjAW2+9xfHjx1O+jtVqpba2lp6enqQycZRSTE5OYrfbFz1PKRVzXpMOySpm/dQU9/zbv0VVSphxZW2z14x0YT2eNr75zX8JWz4RuOqqYjo7504BtN6WNT09TVdXl2YJGg6HI2aXs0zPL2EJu7Lzqa2tZdeuXXz2s59N+RrV1dVcvnw56fS4xsbGhJU5GAzS2tqqeXQ2GcUsCQb53KOP4likc9af//M/8/Uf/pCyEWvYhW1vbwgrZIh0umklS2jdOZ3sL5fLFXc92VDMDJBom7ZIKioqcLlcKW9ijtxfmCgOhzZ1akIkmidbEAzywJNPUpFAZcHC8XEqenv5/Es/Z737PGbx01L0IbvajrLW54NgkJKSEt0LobW2toa7fqVCvOwrk8mU0eWZEMvClY3k6quv5uDBgwmnvpWWluL3+5PuyBVJsh2rQLu0whAJKaZSPPDyy9R//HFi10Too5q6M2fYx3YEhXusF/kNbAXGCws598lPsjdD/T3i0d3dTXl5edLWerFSnfX19RmLAkey7Cym2+3muuuui3tMY2MjLpeLhoYGfD5fzFblidLa2pq0+6N1L5bqCxf4u//+3/nb//k/2f2731E4X1GV4r4jR1iTYJBsCjM38g51tFPOEI208QV+g4pYHSscH6dRo6BbKiT7uRUUFCxq3fVwY2EZWkyv18upU6eivma323G73RnpxDU4OJjWtql0CZpMWKansUxP4/jwQ/7m+HG+9/d/jx+onZ7mlv37WXX0aELX8mOiigFGKQEk/PstttFHNe6I7MzSwUEcwSBenbuxpRIEqq+vX/SzNxQzQ5w5cyaqcpSWljIyMpKx9nijo6OUl5dTUFCQ9i6PVAjMs8DmQICv/dM/UZikLEGEN7glrIwz++Znft/AYVxRUqb/wy9+gX9eCcvedet4PmIPrZaYTKakraXD4Vg0hiAiuswvYRkqZnNzMyMjIxw8eBCYKW0xPT2dcu3TZBgaGqK6ujorihltzpqMUgYRenDxIE/yFtsw4ydAAWYCBFF8gnd5g1uipnnVRum6bPb7IUOKGQwGGR4eDtc+8vl8c6ofhqYJBQUFWCwW/H4/Nptt0WoKtbW1unxPYBkqpohwxx130NDQQHFxMW63m+eee473339fl/ETTTbQel9mbYotCmBGKe/g33mTbQQxozBRwDSvcju38DoDVOOiN6ncy0IdIrUDAwOa5h3r5cbCMgz+hFi9ejU1NTWICM3NzbqO3draumgov6OjQ9O0r+o0Ft37qOYtthHAggLMTHMzb3IbBylA4U5SKQEcIyPoWndSAzKdHxvJslXMSDweD1/5ylfYvXs3q1at0mXMWGUZI+no6NCsj0l5FHcyUVz0chNvMlNPDW7gMP/OJ9PanWCdnub29nY2Dg/TNDlJcRr1iPQiUxvmo7HsXNloiEjYOjU1NXHy5En279+f9jJJPBJd8B8aGtIkzaw8jQrvAvyah/DQAZh4h6304qKW9CpF3PbYY+G/x5xO/vlv/iat62WSmpqaRRsAa4lhMechIqxfv56vfe1rGb1D9vf3U1VVtehxoUaqiRwbD3OarR/M/O8byUwcNr/c0HTRc34JhmLGxGaz8fDDD3PVVVdlbIxE0+4sFkva1lvSnM+56eVWDmJmmls5OGetcjlgKGYOYbFYePDBB1m3bl1Grp9o2t309HS4mHUqrB4fx5KmxRTgAJ+kAw+vcYf+u9+zjN6KacwxF6GgoIB7772Xs2fPar6Ekcy+y4KCgrBrHapwEFqP8/v9cRfHPb3JR02jYZqNwC43qqurM7J5PR6GYiaA3W5n9erVnDlzRtPr9vX14XK56E1gJ8diG4Dr6+tjKqcjzXYRy51MVyuIhuHKJkimGpZqtRO+q6srZhvBwhgbfnMJpXHSvpZotWSVDIZiJsjatWszUk6is7NTk+hvMBhkcHAw6gZrex4opikQIKmS7jqSbAFyLTAUM0FsNhtr1qzJyLVbW1s1KVw8PT2Nz+ejoaEh/OPxeLCnkY6nFw6vl/veey/bYkQl3aWqVDDmmElw++23MzExoXmld5hpdVdXV5d2i4TI9n8h8sGVBdj4/PMcX7+ec4WF2RYljMViSSsiniqGxUwCl8vF/fffn7Hr9/T0ZGQ+U5QnigngzNJ+1VhUV1drvmk9EQzFTJJM7kIJBAIMDw9rGgW0BALY0+ytqSfprrdqTTbcWDAUM2k2bNigSYOeWExNTYV3lmhxp16ZR0oJxCyXmS2yEfgBQzGTxul0Ultbm/FxtNj2Ver3c/fPfqaNQDphzTGLaShmHpGp6Ox8YjUeSoSCQIAHnnwyamuDXKYgx+aYhiubR+i1Z3NoaCi11n9K8YX9+6k7d057oTJMurtgtMRsNmekJ0oiGIqZAvX19brVfkllPnvve++xWuMuZ3phySGLWVVVldF4QjwMxUwBk8mkW5mJzs7OhKodhLjj8mU2Pfdc5gTKIH6zmUAGm8EmS7bcWDASDFJm1apVmie1x8LpdDKcQPZO89AQ2x5/PPMCacRATQ0n7rmHofJyuouK6DaZZroN5QjZCvyAYTFTZu3atbqN1dnZuajrvGpigrt+/GPMOZpvGo1L27bxalMTH5SU0F1QkFNKCYZi5iVlZWW67ToIBAJxg0CVfj9/+sgjeZVIAHBKx5tbKmTTlTUUMw30tJqhFvDzsQUC3P/EE5Sm0fQoGwxVV3M+xW5ceiAimrdCTAZDMdNAT8UcHR2lrq5uznMSDPLgc8/hzlBbh0zS3tJC2q2kM0hlZSXmLAaiElJMESkTkT0iclpETonITSJSISIvi8jZ2d9RF3xE5H+JSK+IpN7KOUfRqyVbCP+8dLXPHTrEig8+0G18LTl77bXZFiEu2XRjIXGL+X3gJaXUOmAjcAr4DvCKUmoN8Mrs42j8DLg7TTlzEr/fP6cnRqbp6ekJL53cdfYs1+3fr9vYWjJWXMxxjRvzak02Az+QgGKKSAmwHXgMQCk1pZQaBu4FQrH5x4H7op2vlDoI5NcEKEHSaWabKiUlJXyir49P/OpXuo+tFR2bN6NyLAI7n5xXTGAV0Af8VETeF5FHRcQBuJVSXQCzv9MKUYrIV0XkiIgc6cuT/M7Tp0/rPmZnRwe3PPEEpjzr+xHJheuvz7YIi5IPrmwBsBn4kVKqGfAS221NGaXUI0qpFqVUS7bvVokwOTnJO1lIe/MHAlzaulX3cbViymrlg4qKbIuxKPmgmO1Au1Iq9C3cw4yi9ohILcDs72VVcPTdd99lIkvrhm+1tBDMcVcwFh0bNjCV47KXl5frGtSLxqL/IaVUN9AmIlfPPrUDOAn8Htg9+9xu4NmMSJiDTE1NcejQoayN32W1cmnLlqyNnw6XdW55mArZ2lESSaK3rq8DvxSRY8Am4B+AfwTuFJGzwJ2zjxGROhF5MXSiiPwaeBu4WkTaReQvNZQ/K5jN5qyucQEcve22rI6fCgGTiffnrcXmIoU5UAwsoSR2pdQHQEuUl3ZEObYT2Bnx+KFUhctVzGYzW7du5eWXX86aDCeLi+m66ipqz5/PmgzJ0rV2LVdy3I0FdNvSF4/c/y/lKFu2bMn6B/jRXXdldfxkac8T91vPPpixMBQzRWw2G1uy/EU7VF3NcJajh8lwTOeOWamS7RsuGIqZFlu3bs3aDncAZTJx6tOfztr4ydBfV0dXDietR2JYzDzH6XRyfZYXyw+uWsV4DgQrFqNj8+Zsi5AwhmIuAW666aasjj9hNvPxjgUxuJzjbIaa/2YCQzGXAC6Xi9WrV2dVhjc3bsypWjnzGSsu5oTOjV/TwZhjLhFuvvnmrI7fZ7Fw4ROfyKoM8ejctCnnyobEw7CYS4QVK1boUp09HodvvTWr48fj0oYN2RYhKQzFXCKISNat5tnCQtpzcB7nLyjggzzYlBCJ4couIa699tqk6r9mgmN33pnV8aPRuW4d43nkxoKhmEsKk8nE1ixvxzpcUcFAKi0VMkj7pk3ZFiEpbDZbVvphzsdQTA1pbm7ObgK0CCfvzq0qLsd1qlivFbkwvwRDMTXFarXS0hIt118/Xm9owFtcnFUZQvR5PHQV5Fexf0Mxlyg33nhjVreETZvNnMmRuWY+ZfuEyIX5JRiKqTnFxcVs3LgxqzK8sX490zlgqc5effXiB+UYhsVcwmR76WSooIDz27Zlbfxpi4WTO3ZwMg9yeOeTK4qZ/dvqEqSyspJNmzbxQRaLMb9z882se+01Xcf0ORycufNOXt+wgaEcThGMR664soZiZoh77rmH1tbWrNSeBbhks3H5+utpOnYs42MNV1ZycudOXl+5kok8W7Ocj2ExlzhWq5Vdu3bx2GOPEcxSa7wPd+zIqGL2NjTw0c6dvOl253wB50QxLOYyoK6ujh07dmStNtD7JSVsbWjA1dam6XVb16/ngzvv5P3S0pxuDJQKhsVcJtx0001cuHCB89komiXCyU9/Gtejj6Z9qaAIFz7xCQ7fdhsf52FQJ1EMxVwmiAj33XcfP/rRj/D5fLqP/3ptLVtKS3GOjKR0/pTVytnbbuPNG27Im9Ig6ZArruzSmBjkOMXFxdx3331ZGTtoNnM6hWp6XoeD9z73OX747W+zZ9u2ZaGUYFjMZceaNWv4xCc+kZV+J6+vW8f1Viu2qalFj+31eDj7yU/y+ooVTC6RgE4yGIq5DPnUpz7FpUuX6Onp0XXcUbOZc9u3s/6Pf4z6+qTVysWbb+bDG27gdI73rcw0ueLKGoqpIwUFBezatYtHHnlkQXfoTPP2jTdyzSuvzGnf171yJee2b+ftpiZ8y9A6RsOwmMuU6upq7r77bp5//nldx+2wWrm0eTOuM2e4eOutvL9pExdzxDrkCiaTiYIcyDEGQzGzwubNmzl//jynTp3Sddx/v/9+un0+AoZ1jIrdbs+JTdJgRGWzgojw2c9+lpKSEl3H7ZiYoDLHKhzkErkyvwRDMbNGYWEhn/vc53QfN1dctVwkV+aXYChmVmlqauJWnctOdnZ26m6p8wXDYhqEuf322/F4PLqOme1qfrmKYTENwphMJnbt2qXr3bqzszOnrEOukEv/E0Mxc4CysjI+85nP6Dae3++nxggCLcBQTIMFXHfddbrWCurv789qb89cxHBlDaJyzz33UFFRoctYXq+Xuro6XcbKFwyLaRAVm83GnTqWnhwfH9dtrHwg7xRTRMpEZI+InBaRUyJyk4hUiMjLInJ29nd5jHPvFpEzInJORL6jrfhLj7Vr1+J0OnUZa2BgALfbrctY+UA+urLfB15SSq0DNgKngO8Aryil1gCvzD6eg4iYgR8A9wDXAg+JyLVaCL5UMZlMNDc36zqewQx5ZTFFpATYDjwGoJSaUkoNA/cCj88e9jhwX5TTbwTOKaUuKKWmgCdnzzOIw2YdK5h3dXUZ65qz5JvFXAX0AT8VkfdF5FERcQBupVQXwOxvV5Rz64HISlDts88tQES+KiJHRORIX19fUm9iqVFaWsqaNWt0G08v1znXySuLycwOlM3Aj5RSzYCXKG5rDKKl6qsoz6GUekQp1aKUaqnOs0anmWDLli26jdXR0ZHdLmU5Qr5ZzHagXSkVqomxhxlF7RGRWoDZ370xzm2IeOwBOlMXd/mwZs0a3SxZMBjE5Yrm8Cwv8spiKqW6gTYRCXWI2QGcBH4P7J59bjfwbJTTDwNrRGSliFiBB2fPM1gEvYNAvb29We1Slm1EBIvFkm0xwiQakvs68EsROQZsAv4B+EfgThE5C9w5+xgRqRORFwGUUn7gPwH7mInkPqWUOqHpO1jCbN68WbeNu+Pj48s64SCXNklDghUMlFIfANE6su6IcmwnsDPi8YvAiynKt6wpLS1l9erVnD17Vpfxrly5oss4uUguubFgZP7kPHoGgUZGRpat1TQU0yAp9AwCAQQCAd3GyiVyKSILhmLmPCaTSdeEg56eHiorK3UbL1cwLKZB0jQ3N+samFiOa5qGxTRIGr0zgdrb2ykuLtZtvFzAsJgGKaGnOwvoti80VzAU0yAl1qxZo2t1u+7ubqzLpMMXGK6sQYronQk0NTW1rOoCGRbTIGX0zAQCGBoayqlsmExiWEyDlCkpKdE1CDQ6Okp9fdRdeksOw2IapIWemUAAk5OTuo6XLQzFNEiL1atX6xoE6uvrWxZbwgxX1iAt9A4CATm1HSpTGBbTIG30DgJ1dHQs+UZEhsU0SBu9g0Cw9BsR5dqaraGYeYreQaCurq6c+/Jqhc1my7llIUMx8xS9g0DT09NLNuEg19xYMBQzb9F7OxjA4OBgzlkWLci1wA8YipnX6L0dbGxsbEkmHBgW00BTSkpKWLt2ra5jTkxM6DqeHuTiTpqEinEZ5C5btmzhzJkzuo3X39+Py+WitzdaGeHcxmKx4Ha7cbvd1NTUUFNTg8vlysmglqGYec5VV11FaWkpIyMjuo2ZDwkHTqeTmpqaOUpYXl6eN02UDMXMc0KZQAcOHNBtzFDCQS6UuxQRqqur5yig2+3G4XBkW7S0MBRzCdDc3Mxrr72GUlHbwmSEsrIy3RXTZrMtUECXy0VBwdL7Gi+9d7QMCQWB9JxrdnZ2YrVamZqaysj1S0tLF7iiZWVlS3K5JhqGYi4R9A4C+f1+mpqauHz5clrXMZvNVFdXz1FCt9u9LCv1RWIo5hIhG0GggYEBRCRhF7qwsHCBFayqqlrWzYxiYSjmEiGUCfTqq6/qNubY2Bgej4f29vYFr1VUVCxQQqfTuWxc0XQxFHMJEYrO6hkE8vv91NfXL1gbzMU0t3zCUMwlhNPp5Oqrr+b06dMZG8NkMrFu3TrWrVtHbW0tFRUVebM2mE8YirnE2Lx5c0YUs7S0lM2bN9Pc3Kxrk6PliqGYSwytg0Br1qyhpaWF1atXG5ZRRwzFXGJoEQRyOBw0NzezZcuWJV+5IFcxFHMJkmoQaMWKFbS0tLBu3TpjCSPLGIq5BEkmCGS329m0aRNbtmyhqqpKB+kMEsFQzCXKli1b4iqmx+Nhy5YtrF+/Pi92iyw3ElJMEbkEjAIBwK+UahGRjcC/AcXAJeA/KKUWZDWLyDeA/wMQ4CdKqe9pIrlBXKIFgSwWC9dffz0tLS1Ltn7PUiEZi3mHUqo/4vGjwN8qpV4TkS8D/yfw/0SeICLXMaOUNwJTwEsi8oJS6myachssgoiEg0But5uWlhY2bNhgLPznCem4slcDB2f/fhnYxzzFBK4BDimlfAAi8hrwZ8D/l8a4BgmyZcsWVq5cicfjMVLh8oxEF6YUsF9EjorIV2efOw786ezf9wMNUc47DmwXkUoRKQJ2xjgOEfmqiBwRkSN9fX2JvwODmDgcDhoaGgylzEMSVcxtSqnNwD3AX4vIduDLs38fBZzMuKpzUEqdAv4HMxb1JeBDwB9tAKXUI0qpFqVUS3V1dfLvxMBgCZGQYiqlOmd/9wK/A25USp1WSt2llNoC/Bo4H+Pcx5RSm5VS24FBwJhfGhgswqKKKSIOEXGG/gbuAo6LiGv2ORPwfzMToY12fui4RuBzzCixgYFBHBKxmG7gDRH5EHgXeEEp9RLwkIh8DJwGOoGfAohInYi8GHH+b0XkJPAc8NdKqSFN34GBwRJE9Ny7lygtLS3qyJEj2RbDwCDjiMhRpVTL/OeN7QIGBjmIoZgGBjmIoZgGBjmIoZgGBjmIoZgGBjlITkZlRaQPSK+ScOpUAf2LHpVbGDLrQyZkblJKLUh1y0nFzCYiciRa+DqXMWTWBz1lNlxZA4McxFBMA4McxFDMhTySbQFSwJBZH3ST2ZhjGhjkIIbFNDDIQQzFNDDIQZaFYopIg4i8KiKnROTEbOW+yNf/VkSUiEQtrCoiZSKyR0ROz17jpjyQ+Zuz5x0XkV+LiD1bMovIfxWRDhH5YPZnZ4zz7xaRMyJyTkS+k2l5tZB7sc8pZZRSS/4HqAU2z/7tBD4Grp193MBMIbHLQFWM8x8HvjL7txUoy2WZgXrgIlA4+/gp4EvZkhn4r8xUVIx3rpmZKhirZv/HH4beb47LHfNzSudnWVhMpVSXUuq92b9HgVPMfHkB/gX4L8wUHFuAiJQA24HHZs+fUkoN57LMsxQAhSJSABQxs5k9oywi82LcCJxTSl1QSk0BTwL3ZkbSuaQjd5rvOSbLQjEjEZEVQDPwjoj8KdChlPowzimrgD7gpyLyvog8OltiRTeSlVkp1QH8E9AKdAEjSqn9esgaIlLm2af+k4gcE5H/JSLlUU6pB9oiHrejwRc8WVKQO965KbOsFFNEioHfAv+ZmWp9fw/8v4ucVgBsBn6klGoGvICe85+kZZ79At0LrATqAIeIPJxZSeeMH5ZZzVTn/xFwFbCJmRvFP0c7Lcpzuq7lpSh3rHPTYtkopohYmPnH/VIptZeZf/hK4MPZFhAe4D0Rmd87oB1oV0qF7oJ7mFHUXJb5U8BFpVSfUmoa2AvcnCWZUUr1KKUCSqkg8BNm3Nb5tDO35rAHHdzvEGnIHfXcdFkWiikzFY8fA04ppb4LoJT6SCnlUkqtUEqtYOaLsVkp1R157uzjNhG5evapHcDJXJaZGRd2q4gUzV5nBzNzH91lnn2+NuKwP2OmEPh8DgNrRGSliFiBB4HfZ1LeCPlSljvWuWmjR9Qr2z/ALcy4RceAD2Z/ds475hKzEU5m3L8XI17bBByZPf8ZoDwPZP5vzFQwPA48AdiyJfPs+B/NPv97oDaGzDuZiWqeB/4+29+PRORO5HNK5cdIyTMwyEGWhStrYJBvGIppYJCDGIppYJCDGIppYJCDGIppYJCDGIppYJCDGIppYJCD/P8xGGlANWTQZAAAAABJRU5ErkJggg==\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": 31, "metadata": { "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 31, "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": 32, "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": 33, "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, 2, Kaivokatu, Keskusta, Kluuvi,...1011Rautatientori 1, 00100 Helsinki, FinlandPOINT (24.93985 60.17038)
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.92477 60.16488)
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, 2, Kaivokatu, Keskusta, Kluuvi,... 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.93985 60.17038) \n", "30 Urho Kekkosen katu 1, 00100 Helsinki, Finland POINT (24.93312 60.16909) \n", "31 Ruoholahdenkatu 17, 00101 Helsinki, Finland POINT (24.92477 60.16488) \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": 33, "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": 34, "metadata": { "deletable": true, "editable": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAAEYCAYAAACnR+FYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7pUlEQVR4nO29eXBb55Xg+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", "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.6" } }, "nbformat": 4, "nbformat_minor": 4 }