{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Working with Geospatial Data\n",
"\n",
"This section of the tutorial discusses how to use `geopandas` and `shapely` to manipulate geospatial data in Python. If you've never used these libraries before, or are looking for a refresher on how they work, this page is for you!\n",
"\n",
"I recommend following along with this tutorial interactively using [Binder](https://mybinder.org/v2/gh/ResidentMario/geoplot/master?filepath=notebooks/tutorials/Working_with_Geospatial_Data.ipynb)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Coordinate reference systems\n",
"\n",
"The ``GeoDataFrame`` is an augmented version of a ``pandas`` ``DataFrame`` with an attached geometry:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"Most operations that will work on a `pandas` `DataFrame` will work on a `GeoDataFrame`, but the latter adds a few additional methods and features for dealing with geometry not present in the former. The most obvious of these is the addition of a column for storing geometries, accessible using the `geometry` attribute:\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 (POLYGON ((-74.05050806403247 40.5664220341608...\n",
"1 (POLYGON ((-73.83668274106708 40.5949466970158...\n",
"2 (POLYGON ((-73.8670614947212 40.58208797679338...\n",
"3 (POLYGON ((-74.01092841268033 40.6844914725429...\n",
"4 (POLYGON ((-73.89680883223775 40.7958084451597...\n",
"Name: geometry, dtype: object"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nyc_boroughs.geometry"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Whenever you work with novel geospatial data in a `GeoDataFrame`, the first thing you should do is check its **coordinate reference system**.\n",
"\n",
"A [coordinate reference system](https://en.wikipedia.org/wiki/Spatial_reference_system), or CRS, is a system for defining where points in space are. You can extract what CRS your polygons are stored in using the `crs` attribute:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'init': 'epsg:4326'}"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nyc_boroughs.crs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case `epsg:4326` is the official identifier for what the rest of us more commonly refer to as \"longitude and latitude\". Most coordinate reference systems have a well-defined EPSG number, which you can look up using the handy [spatialreference.org](https://spatialreference.org/ref/epsg/wgs-84/) website.\n",
"\n",
"Why do coordinate reference systems besides latitude-longitude even exist? As an example, the United States Geolocial Service, which maintains extremely high-accuracy maps of the United States, maintains 110 coordinate reference systems, refered to as \"state plane coordinate systems\", for various portions of the United States. Latitude-longitude uses [spherical coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system); state plane coordinate systems use \"flat-Earth\" [Cartesian coordinate](https://en.wikipedia.org/wiki/Cartesian_coordinate_system). State plane coordinates are therefore much simpler to work with computationally, while remaining accurate enough (within their \"zone\") for most applications.\n",
"\n",
"For this reason, state plane coordinate systems remain in use throughout government. For example, here's a sample of data taken from the MapPLUTO dataset released by the City of New York:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Borough
\n",
"
Block
\n",
"
Lot
\n",
"
...
\n",
"
Shape_Leng
\n",
"
Shape_Area
\n",
"
geometry
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
MN
\n",
"
1
\n",
"
10
\n",
"
...
\n",
"
12277.824113
\n",
"
7.550340e+06
\n",
"
POLYGON ((979561.8712409735 191884.2491553128,...
\n",
"
\n",
"
\n",
"
1
\n",
"
MN
\n",
"
1
\n",
"
101
\n",
"
...
\n",
"
3940.840373
\n",
"
5.018974e+05
\n",
"
POLYGON ((972382.8255597204 190647.2667211443,...
\n",
"
\n",
"
\n",
"
2
\n",
"
MN
\n",
"
1
\n",
"
101
\n",
"
...
\n",
"
3940.840373
\n",
"
5.018974e+05
\n",
"
POLYGON ((972428.8290766329 190679.1751218885,...
\n",
"
\n",
"
\n",
"
3
\n",
"
MN
\n",
"
1
\n",
"
101
\n",
"
...
\n",
"
3940.840373
\n",
"
5.018974e+05
\n",
"
POLYGON ((972058.3399882168 190689.2800885588,...
\n",
"
\n",
"
\n",
"
4
\n",
"
MN
\n",
"
1
\n",
"
201
\n",
"
...
\n",
"
6306.268341
\n",
"
1.148539e+06
\n",
"
POLYGON ((973154.7118112147 194614.3312935531,...
\n",
"
\n",
"
\n",
"
5
\n",
"
MN
\n",
"
2
\n",
"
1
\n",
"
...
\n",
"
2721.060649
\n",
"
1.008250e+05
\n",
"
POLYGON ((980915.0020648837 194319.1402828991,...
\n",
"
\n",
"
\n",
"
6
\n",
"
MN
\n",
"
2
\n",
"
2
\n",
"
...
\n",
"
2411.869687
\n",
"
8.724423e+04
\n",
"
POLYGON ((981169.004181549 194678.8213220537, ...
\n",
"
\n",
" \n",
"
\n",
"
7 rows × 90 columns
\n",
"
"
],
"text/plain": [
" Borough Block Lot ... \\\n",
"0 MN 1 10 ... \n",
"1 MN 1 101 ... \n",
"2 MN 1 101 ... \n",
"3 MN 1 101 ... \n",
"4 MN 1 201 ... \n",
"5 MN 2 1 ... \n",
"6 MN 2 2 ... \n",
"\n",
" Shape_Leng Shape_Area \\\n",
"0 12277.824113 7.550340e+06 \n",
"1 3940.840373 5.018974e+05 \n",
"2 3940.840373 5.018974e+05 \n",
"3 3940.840373 5.018974e+05 \n",
"4 6306.268341 1.148539e+06 \n",
"5 2721.060649 1.008250e+05 \n",
"6 2411.869687 8.724423e+04 \n",
"\n",
" geometry \n",
"0 POLYGON ((979561.8712409735 191884.2491553128,... \n",
"1 POLYGON ((972382.8255597204 190647.2667211443,... \n",
"2 POLYGON ((972428.8290766329 190679.1751218885,... \n",
"3 POLYGON ((972058.3399882168 190689.2800885588,... \n",
"4 POLYGON ((973154.7118112147 194614.3312935531,... \n",
"5 POLYGON ((980915.0020648837 194319.1402828991,... \n",
"6 POLYGON ((981169.004181549 194678.8213220537, ... \n",
"\n",
"[7 rows x 90 columns]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nyc_map_pluto_sample = gpd.read_file(gplt.datasets.get_path('nyc_map_pluto_sample'))\n",
"nyc_map_pluto_sample"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This data is stored in the Long Island State Plane coordinate reference system ([EPSG 2263](https://www.spatialreference.org/ref/epsg/2263/)). Unfortunately the CRS on read is set incorrectly to `epsg:4326` and we have to set it to the correct coordinate reference system ourselves."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'init': 'epsg:2263'}"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nyc_map_pluto_sample.crs = {'init': 'epsg:2263'}\n",
"nyc_map_pluto_sample.crs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Depending on the dataset, `crs` may be set to either `epsg:` or to a raw [proj4](https://github.com/OSGeo/PROJ) projection dictionary. The bottom line is, after reading in a dataset, always verify that the dataset coordinate reference system is set to what its documentation it should be set to.\n",
"\n",
"If you determine that your coordinates are not latitude-longitude, usually the first thing you want to do is covert to it. `to_crs` does this:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Borough
\n",
"
Block
\n",
"
Lot
\n",
"
...
\n",
"
Shape_Leng
\n",
"
Shape_Area
\n",
"
geometry
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
MN
\n",
"
1
\n",
"
10
\n",
"
...
\n",
"
12277.824113
\n",
"
7.550340e+06
\n",
"
POLYGON ((-74.0169058260488 40.69335342975063,...
\n",
"
\n",
"
\n",
"
1
\n",
"
MN
\n",
"
1
\n",
"
101
\n",
"
...
\n",
"
3940.840373
\n",
"
5.018974e+05
\n",
"
POLYGON ((-74.04279194703045 40.68995148413111...
\n",
"
\n",
"
\n",
"
2
\n",
"
MN
\n",
"
1
\n",
"
101
\n",
"
...
\n",
"
3940.840373
\n",
"
5.018974e+05
\n",
"
POLYGON ((-74.04262611856618 40.69003912689961...
\n",
"
\n",
"
\n",
"
3
\n",
"
MN
\n",
"
1
\n",
"
101
\n",
"
...
\n",
"
3940.840373
\n",
"
5.018974e+05
\n",
"
POLYGON ((-74.04396208819837 40.69006636010664...
\n",
"
\n",
"
\n",
"
4
\n",
"
MN
\n",
"
1
\n",
"
201
\n",
"
...
\n",
"
6306.268341
\n",
"
1.148539e+06
\n",
"
POLYGON ((-74.04001513069795 40.7008411559464,...
\n",
"
\n",
"
\n",
"
5
\n",
"
MN
\n",
"
2
\n",
"
1
\n",
"
...
\n",
"
2721.060649
\n",
"
1.008250e+05
\n",
"
POLYGON ((-74.01202751677701 40.70003725302833...
\n",
"
\n",
"
\n",
"
6
\n",
"
MN
\n",
"
2
\n",
"
2
\n",
"
...
\n",
"
2411.869687
\n",
"
8.724423e+04
\n",
"
POLYGON ((-74.01111163437271 40.70102458543801...
\n",
"
\n",
" \n",
"
\n",
"
7 rows × 90 columns
\n",
"
"
],
"text/plain": [
" Borough Block Lot ... \\\n",
"0 MN 1 10 ... \n",
"1 MN 1 101 ... \n",
"2 MN 1 101 ... \n",
"3 MN 1 101 ... \n",
"4 MN 1 201 ... \n",
"5 MN 2 1 ... \n",
"6 MN 2 2 ... \n",
"\n",
" Shape_Leng Shape_Area \\\n",
"0 12277.824113 7.550340e+06 \n",
"1 3940.840373 5.018974e+05 \n",
"2 3940.840373 5.018974e+05 \n",
"3 3940.840373 5.018974e+05 \n",
"4 6306.268341 1.148539e+06 \n",
"5 2721.060649 1.008250e+05 \n",
"6 2411.869687 8.724423e+04 \n",
"\n",
" geometry \n",
"0 POLYGON ((-74.0169058260488 40.69335342975063,... \n",
"1 POLYGON ((-74.04279194703045 40.68995148413111... \n",
"2 POLYGON ((-74.04262611856618 40.69003912689961... \n",
"3 POLYGON ((-74.04396208819837 40.69006636010664... \n",
"4 POLYGON ((-74.04001513069795 40.7008411559464,... \n",
"5 POLYGON ((-74.01202751677701 40.70003725302833... \n",
"6 POLYGON ((-74.01111163437271 40.70102458543801... \n",
"\n",
"[7 rows x 90 columns]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nyc_map_pluto_sample = nyc_map_pluto_sample.to_crs(epsg=4326)\n",
"nyc_map_pluto_sample"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Coordinate order\n",
"\n",
"`shapely`, the library `geopandas` uses to store its geometries, uses \"modern\" longitude-latitude `(x, y)` coordinate order. This differs from the \"historical\" latitude-longitude `(y, x)` coordinate order. Datasets \"in the wild\" may be in either format.\n",
"\n",
"There is no way for `geopandas` to know whether a dataset is in one format or the other at load time. Once you have converted your dataset to the right coordinate system, always always always make sure to next check that the geometries are also in the right coordinate order.\n",
"\n",
"This is an easy mistake to make and people are making it constantly!\n",
"\n",
"The fastest way to ensure that coordinates are in the right order is to know what the right x coordinates and y coordinates for your data should be and eyeball it."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Types of geometries\n",
"\n",
"Every element of the `geometry` column in a `GeoDataFrame` is a `shapely` object. [Shapely](https://github.com/Toblerity/Shapely) is a geometric operations library which is used for manipulating geometries in space, and it's the Python API of choice for working with shape data.\n",
"\n",
"`shapely` defines just a handful of types of geometries:\n",
"\n",
"* `Point`—a point.\n",
"* `MultiPoint`—a set of points.\n",
"* `LineString`—a line segment.\n",
"* `MultiLineString`—a collection of lines (e.g. a sequence of connected line segments).\n",
"* `LinearRing`—a closed collection of lines. Basically a polygon with zero-area.\n",
"* `Polygon`—an closed shape along a sequence of points.\n",
"* `MultiPolygon`—a collection of polygons.\n",
"\n",
"You can check the `type` of a geometry using the `type` operator:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"shapely.geometry.multipolygon.MultiPolygon"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"type(nyc_boroughs.geometry.iloc[0])"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"shapely.geometry.polygon.Polygon"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"type(nyc_map_pluto_sample.geometry.iloc[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Performing geometric operations\n",
"\n",
"The [shapely user manual](https://shapely.readthedocs.io/en/latest/manual.html) provides an extensive list of geometric operations that you can perform using the library: from simple things like translations and transformations to more complex operations like polygon buffering.\n",
"\n",
"You can apply transformations to your geometries in an object-by-object way by using the native `pandas` `map` function on the `geometry` column. For example, here is one way of deconstructing a set of `Polygon` or `MultiPolygon` objects into simplified [convex hulls](https://en.wikipedia.org/wiki/Convex_hull):"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 62.7 ms, sys: 2.64 ms, total: 65.3 ms\n",
"Wall time: 71.6 ms\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAFlCAYAAAAUB7oWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XeY3FXZxvHvQxOQ3pEmiBTpJSAdEhBBFJQuTQERRYSXKtIR6SAISi8iqChWQOkovfcivUUhVAGBUO/3j3PWXZbdze60M7+Z+3NdcwHJ7MyTkOw9pz0nJGFmZmbVMlHpAszMzGzkHOBmZmYV5AA3MzOrIAe4mZlZBTnAzczMKsgBbmZmVkGTlC7AzKyRImJG4D7gemB6YBZgVmBG4A3ghfwY1+ffB3r8Rz5na20s/OfTzDpJRBwPTCLpe/1+fCJgBlKgD+cxBfAiQ4f8/x6Sxjf712bWlwPczDpGREwOjAWWkfR0A15rZoYf+OMZXtiPA16R9GE99Zk5wM2sY0TEZsB2ktZq8fsGMA1pqn44YT8t8DIfD/fngYskPdDK+q2avAZuZp1kW+DMVr9pXit/LT8emdDzI2JS0pp8z/p8T7DPDVwZEXcCxwB/9zq8DcYjcDPrCBExD3AHMGeV16Pz1P2WwO7AW6Qg/52k94sWZm3HAW5mHSEiDgBm6b95raryprt1gT2BeYDjgTMlvVG0MGsbDnAzq7wcdo8DG0q6s3Q9jRYRy5FG5GOAM4ATJf2rbFVWmhu5mFknWIN0brvjwhtA0q2SNgVGkY633RcR50TEYoVLs4Ic4GbWCbYFzipdRLNJelLSLsD8pM1yl0fEpRExJu+Ety7iKXQzq7SImB54EviMpJdL19NKEfEJYAtgD+Ad0oa330p6r2hh1hIOcDOrtIj4DrB6nmLuSnkPwBdJG94+A5wAnC7p9aKFWVN5Ct3Mqm47umD6fCiSPpT0V0lrAF8DlgWejIijImLOwuVZkzjAzayyImIJUgOUK0vX0i4k3S5pc2AZYFLg3og4N/9eWQdxgJtZlX0TOEfSB6ULaTeSnpL0f6Qp9QeAv0bEZRGxlje8dQavgZtZJeUNXGOB5SU9Ubqedpd/vzYnbXj7gLTh7QJJ7xYtzGrmEbiZVdWXgfsc3sMj6R1J5wCLAXsD2wCPR8QeETFt0eKsJg5wM6uq7ShwcUnVKblU0prA+sBSwBMRcUxEzFW4PBsBB7iZVU4OmuWBP5Supcok3SlpC1KIB3BPRJwXEUsWLs2GwQFuZlW0NWn99u3ShXQCSc9I2h2YD7gXuDgiroiItb3hrX15E5uZVUpuWvIosJmk20rX04kiYjJgM9KGN4BjgV97w1t78QjczKpmVeBN4PbShXQqSe9KOhdYghTiW5Iaw+wdEdOVrc56OMDNrGq2Bc6Spw+bLm94u1zSWsCXgEVJG96Oi4h5CpfX9RzgZlYZ+bjTV4DzS9fSbSTdLWkr0qj8A+DOiPhVRCxduLSu5QA3syrZFLhS0oulC+lWkp6VtCdpw9sdwJ8j4qqIWMcb3lrLm9jMrDIi4mbgEEl/LV2LJRExKemD1R7AJKQNb7+S9E7RwrqAA9zMKiEiFgEuA+Zx7/P2k0ffY0hBvjjwU+DnvtK0eRzgZlYJEXEs8B5wOKnpyET5n33/ve+PjZP0fplqu1tELE5q17oWcDxwkoO88RzgZtb2IuIzwD35Pz/MD+XHh/3+2fNNbWLg+5J+09pqrUdELATsB6xNCvITHeSN4wA3s7YVEZMDBwHfAiYHPj3cDWwRsQzwG+BiYDcfOyvHQd4c3oVuZm0pIqYGrgQWBHYGbhzJ7nNJdwCjgI2ARZpSpA2LpH9K2hJYBViYdAvavhExTeHSKs0Bbmbt6qfA48CGwKeBu0b6ApL+A4wjjd6tsAGC/NGIWK5wWZXlADezthMRGwMrATtJ+pB0W9bdNb7crMALjarN6tcnyLcHLnIzmNo4wM2srUTEnMBJwBaS/pt/eElqCPB8tGkWHOBtSdJFwI7A3yJiidL1VM0kpQswM+uRbxo7Fzih56axvBb+KeDhGl5yGmC8pPGNq9IaSdIfI2ISUoivLumR0jVVhQPczNrJ7qTvS0f2+bElgPtrbN7i6fMKkPS73Of+8ohYWdLY0jVVgQPczNpCRCwF7AmM6hfWS1LDBrbM0+cVIemMiJieFOKrSHq5dE3tzmvgZlZcREwJ/ArYVdLT/X66ng1ss5B2oVsFSDoauAg4u3QtVeAAN7N2cDRwh6RfDfBz9YzAPYVePT8DvCt9GDyFbmZFRcR6wJdIa939f25S0nnh+2p8eU+hV8+/gVkiYhL3sh+aR+BmVkxEzAqcDmwl6bUBnrIw8LSkt2p8Cwd4xeTQfoF08sCG4AA3syLyGe2zgDMlXTfI05ai9ulz8Bp4VT0LzFW6iHbnADezUr4LzAwcPMRzamrg0ofXwKtpLA7wCfIauJm1XER8jnTL2IqS3hviqUsCf63jrTyFXk3PAnOWLqLdeQRuZi0VEZ8gHRn7gaRHh3heUP8I3FPo1eQp9GFwgJtZq/0YeIK0/j2UeYA3R3KFaF8RMRkwFfCfWr7einKAD4On0M2sZSJiTWAzYAlJmsDT62ngAml9/aV8m5lViwN8GDwCN7OWiIgZSR22vjnMNpmNmD73+nc1zYpnTibIAW5mTZfXs08FfivpimF+Wb1HyGbF699VNQa4qnQR7c4Bbmat8E3gs8APR/A1HoF3Lwf4MHgN3MyaKiLmJV0Purqkd4b5NTMC0wJP1vHWDvAKiojZSV3Y7ixdS7vzCNzMmiZPnZ8BHC3pgRF86ZLAPXVuQPMRsmoaDfy9xvvfu4oD3Mya6VvA1MBxI/y6ete/wV3YqsrT58PkADezpoiIuUhnvret4Vapete/wVPolZNnbNbEAT4sDnAza7j8jfg04KeS7q/hJeq5A7yHp9Crp+ce8IeLVlERDnAza4atgdmAI0b6hRExBTAf8GCdNXgKvULyh74jgSOG0eTH8C50M2uwiPgUcDSw9gQuKhnMosAjkt6to4YgjcBrasNqRXyRdIHJ6aULqQqPwM2sYXJwngycJqnWKfBGbGCbFnhb0vg6X8daICImIX3o26vGD31dySNwM2ukTYH5gU3qeI1GbGBzF7Zq+QbwMnBR4ToqxQFuZg0REbMAxwNfGW7DlkEsBfymznK8A70iImIq4BBgfa99j4yn0M2sUU4EzpV0a60vEBETA4sB99RZiwO8OvYgNW65rXQhVeMRuJnVLSK+Rpr6/kadL/VZYJyk1+p8HR8hq4DcNnVnYJnStVSRA9zM6pL7lp8EbCzp7TpfrhHr3+AjZFVxCHCmpKdKF1JFDnAzq9fxwO8k3dCA11qS+qfPIY3AR9J73VosIhYF1gcWKF1LVXkN3MxqFhHrASsysmtChzIH8FQDXsdr4O3vKODHkv5TupCq8gjczGoSEdMBpwBbS3qzQS/7DjBZA17Hx8jaWESsSRp5b1C6lirzCNzManUMcLGkqxv4muOByRvwOh6Bt6l80uAY4Af1dNszj8DNrAYR8QVgLdKRr0YaD0zRgNdxgLevLYG3gN+XLqTqHOBmNiIRMTWpX/UOkl5v8Mu/TZ0j8Ij4BPBJ4NWGVGQNExFTAocCm7hpS/08hW5mI3UkcJWky5rw2o2YQp8ZeNEB0ZZ2BW6WdFPpQjqBR+BmNmwRsTrwFRo/dd5jPDBjna/h6fM2lFvt7gZ8vnQtncIjcDMbloj4JHAG8F1JzZqernsKHQd4uzoQOE/SY6UL6RQegZvZcB0K3CLpL018j0ZMofsIWZuJiIVIN9QtVLqWTuIAN7MJiogVgc2ARZv8Vo3Yhe4RePs5AjhS0sulC+kkDnAzG1JETAGcBezcgm/AnkLvMBGxGrAE6QOgNZDXwM1sQg4E7pN0YQveqxEjcE+ht4mImIjUtOWHksaXrqfTeARuZoOKiFHAN4HFW/SWjVgD9wi8fWwKCLigdCGdyAFuZgPKDVHOBnaT1KoRrafQO0RETA4cTuqV/2HpejqRp9DNbDD7Ak8Av2rhe3oTW+fYGbhb0rWlC+lUHoGb2cdExJLAd4AlW9zRrK4ReEQEDvDiImJGYC9g5dK1dDKPwM3sIyJiUtLU+V6S/tXit693BD4d8JakdxpUj9Vmf+C3kh4uXUgn8wjczPrbi7SL+5wC7/02MGUdX+/Rd2ERMT/pxrHPla6l0znAzex/ImIR4P+ApQtdBvIS8GFEzCfpiRq+3kfIyjscOE6SP0g1mafQzQyAiJiE1LBlP0nPlKgh71b+GXBmnsofKY/AC8od+5YHji9dSzdwgJtZj12BN4HTCtdxaK7jhBq+1gFeSN5AeAzpA+BbpevpBg5wMyMiFgB+AGxf+syupA+ArwOrR8R3Rvjls+IAb7kc3tuSNiCeV7icruE1cLMul9tdngUcUuO6c8NJej0ivgJcHxEPSfr7ML90FuC+5lVmABExFTAKWBFYIT9eBL5e+gNgN3GAm9n38j9PKlpFP5Iei4ivAxdGxMnAyZL+PYEv8xR6g+XR9bykkO4J7AWAe4AbSR/+viXpuWJFdikHuFkXi4j5gAOAldpx5CTp6ohYidTV6/6IuBQ4QdItg3yJAzzLt8jNCczV5zE76ajev/s/etat89ctQ29Yrwi8TwrrG4FfAnf5rH15UeakiJmVlqfOrwT+Juno0vVMSERMR1pn3Zl0VOwE4PeS3u3znEeAL3d6A5Hcp34OeoO5f1DPBUwF/At4ts/jOdI69ezAp/q8xuzAh/nxCeAV4CnScsQtwL2koH+u7++3leUAN+tSEfFtYDtgRUnvl65nuCJiYuDLwC7AgsDpwOmSxkbEf4B5Jb1assZ65ONzs/PxQO4b1NOTwvhZYCwfDemex0v9Z1UiYjJgKT46Hf4J0sj6blLv+1eBGUkB3/OYI/9zVuA10geDj43i+zxeqNKfqapygJt1oYiYG7gTWE3SA6XrqVVELErq2b45cC3wJWDyvJO97eQPH7Mx9Mh5ZtIyQP9A7hvU44bza4yI2ejdZLYisCTwKHBTftwIPDHcpj151mZmPhruAz1mJDXlGSjoxwK3SXplOO9pg3OAm3WZvCnpb8D1kg4tXU8jRMTUwE7AIaSR6TnA2ZKeamENPeE21Mh5duBlBh81jyWtR4949Job8SzOR0fX0/HRsL5N0hs1/yKHX8ukpNH6QOE+D7A08Avg8BZeVdtxHOBmXSYivgl8H1hO0nul62mUiFiGNJ2+bX58nTQtfCbwR0nj63jtII0qhxo5zwG8zuCj5meBfzVqDTnf+NV3dL0s8AwpqHsC+5F23JyYZwb2IfVMPxU4usrLHqU4wM26SER8inT8Zy1Jd5eup5EiYldgVUlfy/89ObABKcyXBn4DnCnprn5fF8C0DD1ynhN4h8FHzc8CYyW93aRf20Sky0H67gyfDbiV3sC+pWohmJdyDgDWJ7VfPUHSf8tWVR0OcLMukYPqz8Ddkg4oXU8jRcRnSUG2cv8d6Hl6fTlgG2Bd4D3gSdLVpT2bxT5k8FFzTzi3LFgiYlpST/GewF6etC7eM7K+CXigXdf6Ryp3AjwYWAM4AjilnhmTbuEAN+sSuSnKD4FlOukMb177vRb4taQTI+KTpPXwLYG5gcn4aBhPQtq9/rn8dacCF5Waas4frBbgo13N5gVupzewb5b0Yon6WikiFgd+RJox+RFpH0PHLPM0mgPcrAtExKyks7zrSbqtdD2NFBE/BEaTjpZ9i7S2eh3wE+AR4JWBdllHxAykdfJtgRmAs4FzJD3d5HoHakP6Bh9du763m4MrIj5PutTm08CBwG86ZbahkRzgZl0gIn5HOi60d+laGikilgIuA44ltYS9CzhgpOv7+XW2JR1Hu5PUHvRP9U7jDtKGdEF625DeBNw0jBaxXSki1gB+DEwD7E/6f+LQyhzgZh0uIjYifRNcslmbrEqIiCmBh0lT5PcA+w/RYnW4r9mz8W07UsOTS0ndyB7KjyeGGgn2aUPaN7A/4KNr13d20hJGs+UPQeuS/gy/B+wLXOEgd4CbdbSImIkUQBtJuqF0PY2Qd2RvDPyctPlsQ0nXNuF95gHGAIsAC+fHbMBjpDB/ihTOM5HOf89BWld/kD6ja+AZh0398v/3jUhn/ccB+0q6vmxVZTnAzTpYRJxHaqm5a+la6pVHYl8mbW7qaRSyUCs3d+VR/4KkMJ8HmIh0jeaLwPPAPT2Xglhz5E2LW5HWxh8E9pN0Z9mqynCAm3WoiPgy6Wzt4pLeLF1PrXJwf4EU3JMBh5GOGu0i6aKStVk5+UKX7UlT6jcCe7XLffatMlHpAsys8fLNXScD21U8vFcjHfU6HjiadLxoLeBqh3d3k/SOpJ8B85M2Ht4aEQfmfQhdwSNwsw4UEWcB4yV9t3QttcjHiH4EzAccBPxK0gcRsR5wIrCEpNcLlmhtJnd1O460+fD7ki4pXFLTOcDNOkxErE1qTrJYKy6uaKR8nOsQYAlSgJ/Tcx46b8i7F9isGZvWrDPkP/8nktbHd23lhTat5il0sw4SEdMApwE7VCm8I2K6iDgfuAS4HFhA0ul9wjuAU4DzHd42FEmXAYuR+sTfHhH75eOBHccBbtZZjgSulHR56UKGK+/s/gfwKjC/pBMHaKCyBbAQqZmH2ZDy+vhhpBvalgHuj4h1CpfVcJ5CN+sQuWvVL4FFJf2ndD3DFRGnA1MAWw3S8nQu4A5g7f43iZkNR0SsC/yUtASzq6RnCpfUEB6Bm3WAfIHHGcCOFQvvzYHVgO8MEt4TkXqUn+DwtlpJ+iuwKOl++Dsj4gcRMVnhsurmEbhZB4iI44GZJG1Zupbhioj5SZ3KvjBYOEfE90i3iq0s6f1W1medKSLmI43GPwPsJOnqwiXVzAFuVnERsRJwIWnq/OXS9QxHbsJxI+m6yJMGec6CwPXASpIeaWV91tnypsivACeQPkTuXsULZTyFblZhuWnFWcD3qhLe2VHA08DPBvrJ3C7zXOBAh7c1mpI/k3rXPwHcGxG7RcSkhUsbEY/AzSosIo4CPi1pk9K1DFdEbEDqrLaUpFcHec7+wCqkjWv+JmVNFRELACeRLqvZSdJ1hUsaFge4WUVFxHLARaSGLS+Urmc4cres24D1Jd08yHOWAf4GLC1pbCvrs+6Vp9U3InVzu5rUW31c2aqG5il0swrKa8hnk47EVCW8JwV+DRw7RHhPQToKt6vD21opT6v/jjSt/gLp7Ph380mItuQRuFkFRcSPgMWBDaoyxRwRh5P6VK8r6cNBnnMsMBewaVV+XdaZImIRUkviiYBvSXqgcEkf4wA3q5jcL/wyYMmq7JzN/anPJE2LDzhjEBGrA+eTrj+t0oY861B59L0DqS//KcCPB+gSWEzbTg2Y2cflaeizSetzVQnv2YFzSJ3WBgvvafJzvuXwtnYh6UNJp5Au1/kccE++4rYteARuViERsR+wEmkauu3/8kbExKTLSa6TdNAQzzsLeE/St1tVm9lI5RMUJ5JmwPYc7BRFq3gEblYREbEosAvw7SqEd7YPMDFpCnJAEbE+qZ3q7q0qyqwWkv4ELAKMBx6IiE3z7vUiPAI3q4C8FncdcJ6kk0vXMxwRsSpwAbCspH8N8pxZgHuAjSVd38r6zOoRESsApwNPAd8tcUGKR+Bm1bA1MCnpru+2FxEzkTakbTtEePfc8f0Lh7dVjaSbgKVJrVjvjIhd8pJRy3gEbtbmImI64CHgK5JuK13PhORgvgh4UNJeQzxvG9K0+ShJ77SqPrNGy53cTgWmAraXdE9L3tcBbtbeIuKnwCeqssErInYDNgFWkfTeIM+ZB7gdWLNV3+zMmil/cP0mcATpfoJDJL3V1Pd0gJu1r4hYkrSLe+EqHK/K7V0vBpaT9NQgz5kIuBK4XNIRLSzPrOkiYlZSr//lSMcim3ZdqQPcrE3loLsWOFdS269956n+O4E9JP1hiOftQhqhryrpg1bVZ9ZKEbEuqS3w2pJub8Z7eBObWfvaCvgEqYNZW8vTh6cDf51AeC8M7Ads4/C2Tibpr8DvSSPxppikWS9sZrXLo9kjSLd2VSHovg3MT/rQMaDcRe6XwH6SHmtVYWYFPQbM16wXd4CbtaeDgYsk3Vq6kAmJiCVIjVpWmkCf6H2BF6nIUTizBngaGNWsF3eAm7WZHIibk3ovt7WImIrUrOX/JD0yxPNGAd8BlqpQFzmzej0NfLpZL+41cLM2kteSfwbsL+ml0vUMw8+AmySdN9gT+tzx/f2qXMBi1iDPkK7HbQqPwM3ay1bA5MAZpQuZkIjYkrRBZ9kJPPVw4C5JFzS/KrO2Mg6Yvlkv7gA3axMRMS1wJLBBu29cyy0jDyZdEfrmEM8bA2wELN6q2szahaQPIqJps06eQjdrHwcDF0u6pXQhw7A28IqkGwd7Qt5JfzawnaRXWlaZWXt5tlkv7BG4WRuIiMWBr1OBjWvZdqRz30M5gfSB5LIW1GPWrhzgZp2qz8a1A6qwcS0ipgHWArbN/z0xsAGwFDAT8AqpAc0qwGKFyjRrF55CN+tgWwJTMuERbbuYGnhL0msR8TngZmAP4D3S3d4TAzsBMwN/i4jvRcRsxao1K2vhZr2we6GbFZQ3rj0EfLUia989I+63gfWBc0kNWk6XpDyb8GfgftKa/hdIfc/XA+4GfgP8RtJrJWo3a6XcffAlSdM25fUd4GblRMRPgKklbV+6lpGIiJ6p/q9Kuq7Pj28LfJ90G9m7fX58CtLGty2ANUnXLf7Ym9usk0XECsDJkpZsxut7Ct2skIhYjBRo+5SuZSQiYhtgGmCvfuE9L+kY3FZ9wxtA0tuS/iRpY2ARYArg4YjYMyImb2H5Zq00GriqWS/uADcroM/GtQMlvVi6nuGKiJ1Jfc+vAt7q8+MTAecAR0m6b6jXkPRvSd8lbXJbkRTkW+bXMOskY4Cm3QfuvzBmZWwBfJKKXOwRyX7ALsCqwIN8tEXkrqTvJ8cN9zUl/VPSV0m/FzsBt0fEmo2r2qycvGy0HHBts97Dx8jMWiwfwzoS2LDdO67B/2YLjiatYa8i6bmIGEu+pCEiFiEtAyxfy69H0vURsSKwIXBKRDwK7C3p3kb9GswKWBG4V9IbzXoDj8DNWu8g4FJJN5cuZELyjvPTgJWB1SQ9l3/qWWCuiJiMdFHJPpKeqPV9lFxIamTzV+CKiDg7Iuas71dgVkxT17/BAW7WUnnj2lbAD0rXMiE5nH8FzAes1W/H+FhgTmB/UqOKMxvxnpLelXQisEB+3Xsi4rB83M6sSpq6/g0+RmbWMnkq+u/ABZJ+XricIUXElMCFpOYsm0oa3+/n5ySd634fWFLS802qY07gEOBLwKHAqf13uJu1m/yBcywwc/+/O43kEbhZ63yd1MXs1NKFDCWv0f+N1BJ1o0G+Af0HmBHYtVnhDSBprKRtSa1b1wUejIiN84chs3a1KnBLM8MbHOBmLZFD8Shgp3beuBYRM5Gm/R4Atpb03iBPPZx0jOymVtQl6V5J6wA7kjbM3RQRq7Tivc1q0PT1b3CAm7XKgcBlkloSeLWIiDlIR14uJ33Q+HCQ561FaqN6Px89StZ0kq4ElgVOAs6LiD9FxEKtrMFsGJq+/g0OcLOmi4hFafONaxHxGeA64BeSfqhBNsdExPSkNqjbAU+TNrK1lKQPJZ0HLAhcD1wXESf7whRrBxExCzA3cEez38sBbtZEea32JOBgSS+Urmcg+QPGP0hd1I4c4nmTAecBf5J0BfkoWWuq/DhJ4yUdAyxEms5/ICIOjIipStVkBqwBXCvp/Wa/kQPcrLk2B6YFTildyEAiYjngSlJf80FrzOfBfwl8AOyWf7jnKFlRkl6WtDtpan1B4JGI2CEi3KjKSmjJ+jc4wM2aJm9cO5o23bgWEWsAFwPfkvSrIZ43NfAnYDpgkz4b28ZScATen6QnJX0d+Arpg9O9EfFl71i3FmvJ+jc4wM2a6QDgckk3li6kv4j4MvBb0hnvi4Z43tzADcBzwHr9jsU8SxuMwPuTdDtpFLQncATw9zzTYNZUETEP6ajo/a14Pwe4WRPk/uDbAHuXrqW/iPg6cDrwJUnXDPG85UnHxM4Bvj3AkbK2mEIfSG7NegmwBHAu8MeI+E3erGfWLKOBawbbBNpoDnCzBmvnjWsRsSPpPPqakm4d4nmbAhcBO0o6bpBvSM8BM+XNbW1J0vuSziS1Zr0fuDUijo+IGQuXZp2pZevf4AA3a4ZNgelps41rEbE3sBfpUpIBp/jytaEHkEJ+raGm1/O6/jhg9mbU20iS3pR0KLAwMCnpDvK985WPZnXLH9xbtv4NDnCzhsobvo4hbVxr+jGS4cihfDhpSn8VSY8P8rzJScfE1gM+L+meYbx80aNkIyXpBUk7ASuR7mp+OCK2zrvszeqxIPAuUPOtfCPlADdrrAOAKyXdULoQgIiYCPgZqZf4qpL+NcjzZiGNHCbho9eGTkjbroMPRdLDkjYk7VbfEbgjIr5QuCyrtjHA1a1a/wYHuFnDRMTngG/QJhvXImJS0gauRYDRkl4a5HmLAreQzoNvLuntEbxNJQO8R/6gtRLpxrOTIuKyiFiicFlWTS1d/wYHuFlD9Nm4doikcW1Qz+Sk60BnANaR9Pogz1uHNPLeT9IBg/U/H0KlptAHknes/4H0QecvwGUR8YuIqPSvy1onL8GsDgx6qqMZHOBmjbEJKSxPLl1IXoe/BHgb2EDSWwM8JyJiN1Jf869KOr/Gt6v0CLwvSe9J+hlpx/qzwN0RcURETFe4NGt/SwDjJP27lW/qADerUw7MY2mDjWs5bK4gbaTZQtK7AzxnGlITl82BFepcr6/8CLw/Sa9L2g9YHJiJtNFtl3Y+LmfFtXT3eQ8HuFn99geuKr1xrU943wzsMFD71rzefRvwErCypKfqfNuOGYH3J+lfkrYH1gS+ADwUEZu6NasNoOXr3wDRwg1zZh0nb1z7B7BoybXvHN6XAzcC/zfQTtiI2AI4Hthd0rkNet+JSVP1Uw002u8kETGa1Nv+A2APSdcWLsnaQJ6ZeQn4tKRXWvnevq3HrEZ5JHYi8KMNNfT7AAAgAElEQVTC4T0tcBmp7enHwjsiPgEcRxpFjpF0b6PeW9IHEdHTzOXpRr1uO5J0dUSMAjYDzo2Ie4G9JT1UuDQrazngkVaHN3gK3aweG5PWSH9eqoAc3peTjoHtOkB4zw1cSwrYZRsZ3n107DR6f5I+zDe3LUSaebk2Ik6NiLbvRmdNU2T9GxzgZjXJLTiLblzrM/K+FdhlgPBeO//c74ANJb3WpFI6biPbhEgaL+lYUvet14H7I+LgvKHRukuR9W9wgJvVagfgdknXl3jzvJP8UuB24Pt9wzsiJsr9zM8iXRd6TJO7Q3XNCLw/Sa9I2hNYBvgM8EhE7Jib6FiHi4hPkv7fF/k+4AA3G6E8+t6b1L2rxPtPQxp53wns3C+8ZwQuJu2cXlbSP1pQUteNwPuT9JSkLUl95DcG7ouI1ctWZS2wEnCnpDdLvLkD3GzkdgBukXRXq9+4z8j7LuB7/cJ7FHAH8ABps9pw+5nXq2tH4P1JuoP04WlP4PdeG+94xda/wQFuNiIlR995ffVvwD30Ce/cVe3bpO5ru0vaU9J7LSxtLF0+Au8rt2a9CDib1CPAOlex9W/wOXCzEYmIXYDVJX21xe87NWnkfR/w3Z6e5RExJal969KkjWqPtLKuXMOcwK2SPtXq925neTnjYdLVrI+VrscaKyKmJx2dnKlUDwSPwM2GKZ+nPgj4UYvft2fkfT8fDe/PkrquTUQKiZaHd/Y8MJM3bn2UpJdJjXOK7JWwplsNuLFkAyMHuNnwCXgK+GNE/Cgilo+IpjZDyuH9V+BB4Dt9wvurwA2k0ffWpTbRAORjdOMAj8A/7nhgjYhYqnQh1nBF17/BAW42bJLelbQUsD4wJXAG8FJE/Dkivh8RizSyT3ZETEUK738CO0r6MCImiYijSMGwnqSTm3xEbLi8kW0Akv4LHAocVroWa7jROMDNqkXS3ZJ2l7QYqZHHb0g3V10M/Dsizo+I7SLi07W+R7/w/jYwaURsC9yb32sZSbfW9ytpqK4/SjaE04EFfaysc+TTBbOTToMU4wA3q4OkcZJ+LWl7SfMCKwLXkI4S3RIRj0fEafkWq5mH85o5vC8BHiEdR9qDdD3opsDOwDqSXmrGr6cOHoEPIq+R7g8c4ZvMOsYawD8GuvGvlXyZiVkDSXqSNLV+Rv5mvShprWxL4LSIeIp07OQq0ga0KYBZgVnyYw7gO8B40i1fj5M2sH1J0t0t/cWMzFhgntJFtLFfA3uRll/+VLgWq98YCh4f6+FjZGYtkndpL0v6yz8GGEXqo/0CaRPYy8AqpF3dP8//vE/S2CIFj0BEbAxsJmnD0rW0q4hYFzgGWKz0yM3qExFPkj5UP1i0Dge4WXm5p/LFpF3u2/XsNq+KiFgBOF7S8qVraVd5RuYfwFmSzilcjtUoIuYjnQD5VOkNpF4DNyssj8wvIjWF2L5q4Z15E9sE5G/2PwAOjojJS9djNRsNXF06vMEBbtYOtiL9XdyuwlOrbuYyDJJuBO4m7XOwamqL9W/wFLpZUbkRzEOkkXcrbg5rmoh4FlhZ0tOla2lnEbEoKQA+K+n10vXY8OVlkOeB5SU9Vbgcj8DNCtsUeK7q4Z35KNkwSLqf1Nd+j9K12IgtAvy3HcIbHOBmxUTERMC+pE5dncABPnwHAjtFxKylC7ERKd59rS8HuFk5XwPeAK4oXUiDeCPbMOUR3C9JH+CsOtpm/Rsc4GZF5LW0/YBD22E3a4N4BD4yhwFb5GNJ1ubyfpVVSZ0W24ID3KyM9fI/Ly5aRWN5BD4Ckl4ATgQOLl2LDcvSwFhJ40oX0sMBbtZiHTr6Bo/Aa3EssFZELF66EJugtlr/Bge4WQlrAVMDfyhdSIN5BD5Ckt4ADgd+XLoWm6C2Wv8GnwM3a7mIuBY4VdL5pWtppLxG+BbwSUnvla6nKiLiE8DDwJaSri9dj31c/n/0EjCXpP+UrqeHR+BmLRQRq5HuEb6gdC2NJul90sUss5eupUokvQMcgK8bbWcrAA+2U3iDA9ys1fYDDs9h14k8jV6b84HpgC+VLsQG1Hbr3+AAN2uZiPg88FnS+d9O5Y1sNcg98H8IHB4RE5euxz6m7da/wQFu1kr7AUd2+PqwR+C1u4h0P/zXSxdivSJiamAJ4MbStfTnADdrgYhYGlgKOLt0LU3mEXiN+lw3ekjeNGXtYRXgNklvlS6kPwe4WWvsCxwtaXzpQprMI/A6SLoOeBDYoXQt9j9tuf4NDnCzpsvXR64EnFa6lhbwCLx+PwT2zVO3Vl5brn+DA9ysFX4I/KQdp+CawAFeJ0n3AFcC/1e6lm4XETMC8wG3la5lIG7kYtZEEbEAcAMwX+661dHczKUx8gUntwELSXqxdD3dKiI2Ar4pqS2P93kEbtZc+wAndkN4g5u5NIqkJ4Bfk2ZvrJy2Xf8GB7hZ00TEvMBXgJ+WrqXFvJGtMQ4Fto6IeUoX0sXadv0bHOBmzbQ3cEq7tV9sAa+DN4Ck54GfAwcVLqUrRcScwIzAvaVrGcwkpQsw60T5L/8mwAKlayngWRzgjXIM8GhELCLpgdLFdJnRwDWSPixdyGAc4FZEbhc5ZZ/HFP3+e6BH/+cE8Dzw7/x4BrijTTZP7QmcJeml0oUUMBaYu3QRnUDSaxFxJOm60Q1K19Nl2nr9G7wL3fqJiIkYXpjWGro9j57dym8Bb/f598EeAz0HYDbgU/kxf/7v04HTJI1t6G/OMEXEbKRmHJ/L06BdJSI2ATaRtFHpWjpBREwOPAJsKumm0vV0g3wr3DPAGEmPlK5nMB6BV0QO1smpL0yH87zJqC1Qnx/Gc/o+3lUTPj1GxOeA7wL3RsQ1wM9I02Ct/KS6G3B+N4Z35k1sDSRpfEQcRLpudPUW/1nuVvOTZvgeLV3IUDwCbxMRsRiwDrAwsBAwAx8N3MmB8Yx8lDrS57zTCd8gcherrYDvAe+RRi//bMH7zkQaLS0h6dlmv187ioi5gJslzVG6lk6Rz9ffB+wm6W+l6+l0EbEjsIKkbUrXMhQHeGERsQ7prPB8wO9Jf0kfJp2l7Rus4zshWFstT4XtR/rLuG4L3u9QYBZJXdvL2s1cmiMivgYcACzdzhurOkFE/Ba4RNIvStcyFAd4IRExBXAc8AVSgP8hN8GwBouIyUhr0ttJ+kcT32c64HFgVG7E0bUiYiywoqRnStfSKfKH0ZuB4yX9unQ9nSovV44jfVBq61k0nwMvICIWIv1FnIH0h+S3Du/mkfQucCRpZ3gz7Qxc3O3hnfkseIP1uW700Pyh1JpjMeDVdg9vcIC3XERsDVxH2ly1maTXCpfULX4JLBsRCzfjxSNiKlKAH9aM168gb2RrAknXAI8B25eupYO1dfe1vhzgLRIRU0XEL0jT5aMlneY17dbJ93D/nLRDvBnWAu6S9HCTXr9qPAJvnn2A/SLik6UL6VBtf/67hwO8BSJibtLNQh8Ay0q6r3BJ3epkYKOImLUJr70acE0TXreqPAJvEkl3AtcCu5SupdNExKTAKlTk77IDvMkiYg7Sp7nTJW0r6c3SNXWrfC3jBaRz4o22GtC0DXIV5BF4c+0P7Jbvq7bGWRZ4siodFB3gTZQ7cl0FnCHpuNL1GAA/Ab4TEVM26gXz7vP5gTsa9ZodwCPwJpL0KHAhaVObNU5l1r/BAd40ETEzcCXwa0lHlK7HkrxGfROwdQNfdmXglrzb3RKPwJvvEGDbfHGONUZl1r/BAd4UETEDcAXwF9JfMmsvxwL/l897NoKnzz/uOWDmvKZoTSDp38Bp+LrRhsi9OZYj7S+oBAd4c5xG+kOwr3eat6XrgNeB9Rr0eqviAP+I3NfgBWD20rV0uKOA9XNvCavPisC9kt4oXchwOcAbLCI+DywP7O3wbk/5/8sxwB71vlb+1L4YcGu9r9WBPI3eZJJeBY4GDi1dSweo1Po3OMAbKrc6PAo4UNLbpeuxIf0emDsiRtX5OgsBj+dz5vZR3sjWGicBn4+IJUsXUnGVWv8GB3ijrQvMCJxbuhAbWp7iPQHYvc6Xmp82v3KwII/AW0DSW8A9+Pe6ZhExLbAIaYNrZTjAGyQiJgaOAPZxX/PKOANYKyI+XcdrzAs82ZBqOs+zOFRaZXrg1dJFVNiqpJMklZpJc4A3zpbAa8BFpQux4cmbVc6kvo5WnwaeakQ9HWgsnkJvlemB/5QuosIqt/4NDvCGiIjJgR/hjWtV9FNgm9yMpRaTAL7zemCeQm+d6fAIvB6VW/8GB3ij7ATcKemG0oXYyEgaC1wCfKvGl3iJtO/BPs6b2FrHU+g1iohZgLmpYCdFB3id8shtb+CHpWuxmh0L7FLjHcsvArM0uJ5O4WYuLZCPMgZQqfXbNrIGcG0V9y45wOu3N/AXSQ+WLsRqI+lu4GFgkxq+/EHAx3cG4GYuLTMd8KqX72pWyfVvcIDXJd80tgNuZdgJjgF2z2f5R+JmYJkaR+/dwOvgzefp8/pUcv0bHOD1Ooh0TejY0oVY3S4FJiNNpw2bpNeAx4GlmlFUB/BRsubzDvQaRcQ8wNTA/aVrqYUDvEYRsTCwAXBk6Vqsfnn68Thqa696PbB6QwvqHD5K1nzegV670cA1VV1+cIDX7jDgyNyL2DrD+cBSEfG5EX7d2aRNcDM0oaaq8xR683kEXrvKrn+DA7wmEbEisAypB7F1iNyF6WfAbiP8utuBC4BTalhD73Q+StZ8XgOvQf67Wtn1b3CAj1j+n34kcEDV2u7ZsJwMfC0iZhvh1+0DLAxs1fiSKs0j8ObzFHptFgLeBZ4oXUitHOAjtx7pL8wvSxdijSfpZeA3pOY8I/m68cAWwLERMW8zaqsoj8Cbz1PotRkNXF3V9W9wgI9InwtLfiDpg9L1WNP8BPh2REw5ki+SdC9wOHBeREzSlMqqx81cms9T6LUZQ4Wnz8EBPlJbk1pn/rV0IdY8kh4FbgC2qeHLjwfeBn7Q0KIqqk8zl5EuSdjwzUL6PbZhyoOx1XGAd4fcrvAQfGFJtzgW2C3/RR82SR8C3wB2jojlmlFYBfkoWZNExKrAKNJ94DZ8SwLjJP27dCH1cIAP3/eAWyXdXLoQa4kbgDeBlUb6hbmxz06kqfSpGl1YBXkjWxNExEbA74HNJT1Tup6KGU2Fj4/1cIAPQ0RMD+yFLyzpGnmW5VLSX/Ravv5C0oeA4xpZV0V5I1sDRcTEEXEEqf3v2pIqH0QFVH79Gxzgw/UD4A+SHi5diLXUNdQY4NkuwJiIWL9B9VSVR+ANEhEzkvbgLAssK+nOwiVVTr63YCXg74VLqZsDfAIiYi5ge+Dg0rVYy90ALD3S3eg9JL1OOhd+ag3nyjuJ+6E3QEQsCdwG3At8UdJLhUuqquWBRyS9UrqQejnAJ+wg4JSqb3awkZP0X9LmoBXreI0bgdOAs7u4S5s3sdUpIrYArgD2kbRnFe+ubiMdsf4NDvAhRcQipMYtR5WuxYq5H5i/ztf4ETADI2wO00E8hV6jiJg0In5CmgEcLemC0jV1gEq3T+3LzSaGdhhwRL4y0rrTOGDWel5A0nsRsSVwY0RcLenBxpRWGc8Bs0TEJB45Dl9EzAL8FngLGOWLk+oXEZ8k3WNxXelaGsEj8EFExMrAEsDPS9diRb1AapRRl9wcZh/g/Ij4RN1VVYik94AXgdlL11IVETEKuJ10Ve2XHd4NsxJwl6Q3SxfSCA7wAfS5sGR/Se+UrseKakiAZ2cCT5Om1LuNj5INU0RsC1wCfF/Sfm7b3Bj5ut/N6ZD1b/AU+mC+AkwF/Kp0IVZcwwJckiJie+CeiPibpGsa8boV4XXwCcjHm44nrdGuJumhwiVVWm6itDLpzPdo4LPAjXTQB2gHeD/5EorDgd39yddIAV7XGnhfkl6KiO2AX0TEEl00NeqjZEOIiE8BvyMtNSyXjyDaCOSlqc+TwnoMqV3qHaQR9y6kTprvlquw8RzgH7cNaePSpaULsbYwjsZNoQMg6dKI+BNwckRs3iW99X2UbBARsRJwAXAKcFjup28TkO8pWIYU2KOBFYCHSDvMDwZukPRWuQqbL7rje8fw5IYdjwBfk3Rr6XqsvIiYCHgH+GQjP73ny3FuBw6XdF6jXrddRcSmwEaSNi5dS7vIe212JPWa+KYk33I4hPz7tQi9U+Krkj4YXk0aZV8rqavuRfcI/KN2Bm5yeFsPSR9GxEvATEDDmvlIejs357g8Iq6X9FSjXrtNeRNbHxExOemEyyhgJUmPFS6p7eTAno/eKfE1gDdIgf1rYAdJ48pVWJ4DPMs7FPeghtunrOP1bGRraDc+SXdHxFHAuRGxRofvufAmtiy3Z/4D8ASwQu74Z/xvL0DPlPgYYFLS6Poy0lXOTxcsr+34GFmvfYALJT1SuhBrO3U3cxnCccD7pNvuOtn/mrmULqSkiFgduJXUoGWzbg/viJgxIjaMiJ9FxEPAfcBXSctLawNzSNpK0tkO74/r6r9MPSJibmBbYNHStVhbauRZ8I/IU/TbAHdExBWSbm/G+5SWu9H1NHN5tnQ9rZb3UuwC7A1sKenKwiUVkY92rULvOvb8pGY1VwFnAPd4E9/wOcCTg4GfS3qudCHWlpoW4ACSno2InUld2pbulC5RA+g5StY1AZ53Sm8K7Au8Bny+C/Y7/E9e6+97tGsJ0o1qVwPfA27LnfqsBl0f4BGxGLAu6ZC/2UCaGuAAki6IiC8BxwDfaeZ7FdRzlOym0oU0S+61vRTpeNMo0k7pZ4DdgMs7/chgXiLpOdo1hnR15wOkwD6QdLTr7XIVdpauD3DShSWHuXGCDWEcsFAL3mdn4O6IWE/SxS14v1brqI1s+djpkqTAWjY/5iXdYHcHcA3pe0vHXl6TlwYWpXfjWc8HlquAE0hHu3wZVJN0dYBHxKqkP3wbla7F2lrTR+AAkl6LiK2B30bEkh14RKayR8nyuf0lSCHdE9ifAR6k99KR44EHOq3bV1/5aNdn6F3DXoO0NHA1cB6wvaQXylXYXbo2wPtcWLKfLyyxCWhJgANIui4izgLOyiPxTppyHUtaD21beQp4flLDkEVIH/AXBT4N/JMU1jcBJwH3d8P3joiYg94p8dHAxKQR9t+APSU9U7C8rta1AQ5sAExBaghgNpSWBXh2EOnShR2Bk1v4vs1WrB96vihkJmDmPo/+//3Z/HiONA3+APAn4MfAQ508su4rImYCVqc3tGckLQdcTbon4pEO+2BZWV3ZSjV/yr4f2EXSZaXrsfaWp09fBaZo1TeuiFgAuAFYtVNupcrHNW+UVHeI5/8nM4/g8UngZdJlIYM9niQFdaeeAhhQRExNWrvuWceeD7iOFNhXA/f6aFd76tYA/xawGbCmP0nacETE68BcrdyQk8+H/5jUIfB3Ve/UFhGTAm8CU0p6v8+PBzA1Ew7hviPmSfho+L7E0OH8H4dQko92rUhvYC9Oai5zFSmwb/fRrmrougDPO0cfBTaQdFvpeqwaIuIxYB1Jj7b4fdcgnZSYhtSv4MJ2D6K8M3l6Bp6m3ge4gjQi7vvj7zF0APd//NcfvocnzzguS+8a9nKkGcieS0Bu8tGuaurGNfBdgOsd3jZCPevgLQ1wSddExIrAF0kBvn9EHAT8sVVBngNgRoY/XT0D8F8GDt6XgXtIG8H+N2p2gDRO/gC1GL1r2KuQlgeuBn5COtrlY7MdoKtG4BExI/Aw6QKBln4jtmrL93f/QtIfC9YQwJdIQT5J/udd+afV5zGh/xYwnE1dPY9pSHsAJjRN/b/p7MGmYCPiQuC3kn5bz++F9cp/Lj5L75T4GqT/Xz1T4n+X9GK5Cq1Zum0E/kPgAoe31aCZF5oMS54yvjgiLgG+DPyA1Fs8+jzo998D/RikC1QGCuS7BvixVxu4/t5RzVxKiYg56Z0SH51/+CrgYmB3SV3TrrabdU2AR8Q8wDdIZzvNRqrVR8kGlYP8L/lRNcWOklVZRMxMOtrVE9rTk452XQUcCjzmPQHdp2sCHDgEOEnS86ULsUp6AViwdBEdoO2bubSDiJiG3qNdY0iNZK4lTYn/nNREpq03M1rzdUWAR8TipLtlFyhdi1VPRIwC1iE1+LD6eAQ+gHyuvedo1xjSTOEtpMD+Nulo1/uDv4J1o64IcFL3oB9756UNVz5uuBnwXdIO7FOB04sW1Rl6biTravlM/Ch617BHAfeSAnsf0tGu8eUqtCro+F3oEbE6cBawULe0QrTaRcSCpOs8tyK1Mz0ZuKzqTVTaxWDNXDpdPtq1OL1r2CsDT9C7U/w6SW+Uq9CqqKNH4Pl4xVHAvg5vG0wOlfVJwb0I6QPfMpKeKllXJ5L0XkS8SNo937E7pfP3ngXonRJfnbTr/2rSn69tJL1UrEDrCB0d4MCGpF/jBaULsfaTj+J8C9geeJw02v5DN9wwVVjPUbKOCvDc6310n8eHpBH2n0n3LvyrYHnWgTo2wPOo6jBgJ+/WtB55KnMMabS9Ouk2urUl3V+yri7TERvZImIWUtOUnsCelt4LQA4BHvfRLmumjg1wYDvgaUlXlC7EyouIGUh9AHYE3iaNtrfxumMRldzIFhHTAqvRG9hzA/8gBXbP/eAObGuZjgzwiJgKOIDUrcq6VF6HHEUabX+V1KXqG6Qdvv5GW04lRuD5JMKK9G48+xxwM2lafHvgzm7aiGftpyMDHNgV+IekO0oXYq2Xv/FuTjoCNj1wCrCX+0G3jbZs5pKX3Zajd+PZssDdpBH2XsDN3h9h7aTjAjy3HNwVWL50LdZaEbEQaYp8K+AGYD/SETDvgWgvbdEPPe+HWJLeKfGVgMdIgX0k6WjXf8tVaDa0jgtwYF/g15IeL12INV+fI2DfJU1xngksLenpooXZUIpMoecllQXpnRJfnXRJzdXAGcCWkl5pdV1mteqoRi4RMS9wO7CwpBdK12PNk4+A7UBai3yU3iNgPu/f5lrZzCVfYtQzJT4aeJfe5inXSPp3M9/frJk6bQT+I+CnDu/OlKc81yRtSlsN+BXwBR8Bq5Y+zVxmI02nN0xEzEo62tUT2FORwvoq0sbWJ72B0TpFxwR4RCxF7/le6yARMSO9R8DeJI22t/L6ZKX1HCWrK8AjYjp6j3aNAeag92jX8cCDDmzrVB0T4KQLSw71ud7OkNcrlyN9IFufdARsa9JOYH9Drr6aNrLlEwYr07vxbCHgJlJgf5N0tMt9660rdESAR8QYYH58W1TlRcQnSUfAvgNMRzoCtof7RnecZxlGM5eImIz0Qa5nSnwZ4C7SlPgewC0+2mXdqvIBntdFj8QXllRaRCxMmiLfEriedJrgch8B61gDjsAjYmJ6j3aNITVSeYQ0wj4MuF7Smy2s06xtVT7ASR22AH5XtAobsbwbeQPSEbCFSEfAlpL0TNHCrBWeBZbLSyUL0zslvjrwHCmwTwU2l/RqqSLN2lnlj5FFxDnAjZJOK12LDU9EzEXvLWCPkDal/dEzKN0jIlYg3bf+PKk3fc8lINdIeq5kbWZV0QkB/hCwmaR7Stdig+tzBOy7wCqkI2CnSHqgaGFWTER8Hhgn6cnStZhVUaUDPB8heRaY3pcKtKd8BOybpPXtN4Cfkzrl+QiYmVkdqr4Gviy+Eajt5HXN5Uk7yb8CXETanHaLj4CZmTVG1QN8eeDW0kVYko+AfZ0U3NOS1rZ39xEwM7PG64QAP7d0Ed0uHwH7DrAFcB2wD3CFj4CZmTXPRKULqFWfadpbStfSjSJisojYJCKuIe0efo10BGwDSb7C08ysyao8Ap8b+JAGX4ZgQ8tHwHpuAfsnaZr8Tz4CZmbWWlUO8OXxpqiWyEfA1iIdAVsZOB8YI+nBooWZmXWxygd46SI6Wb8jYK+TRttfdytLM7PyKrsGTrrgwDvQGywipoqITSPiQuAxYFHS5rRlJJ3u8DYzaw+VbOSSe2i/CnxK0uul66m6iJgaWA/YmHSBxA3AhaS17VdK1mZmZgOr6hT6osDTDu/aRcQ0wJdJob0G6Qaw3wHb+fIIM7P2V9UA9/R5DSJiWlJntI1Itz5dSwrtb0j6T8HSzMxshKoa4N7ANky5X/z6pNBeFfg7aXp8a0mvFSzNzMzqUOUAP7F0Ee0qIqYnhfbGpJu/rgZ+A2zhZQczs85QuU1see3236QbyN4rXU+7iIgZgA1Iob0icBVpevwSh7aZWeep4gh8WeBuhzdExEz0hvbngSuAc4BNJL1RsDQzM2uyKgZ4V69/R8TMwFdJa9rLA5cBZwIb+o5tM7PuUcUAX460nts1ImIWUmhvDIwCLgVOA77qxipmZt2pUmvg+QayfwErSnqqcDlNFRGzAl8jhfbSwN9Ia9qXSnqrZG1mZlZe1UbgcwITA0+XLqQZImI2YENSaC8JXELabX+ppLdL1mZmZu2lagG+HHBrJ91AFhGfIoX2RsDipND+CXCZpPElazMzs/ZVtQDviA1sETEHvSPtRYGLgGOAKxzaZmY2HFUM8ENLF1GLiJiL3tBemBTaR5JC+52StZmZWfVUbRPbG8BcVejbnTfcLUm65Ws9YH7gL6SNaFdKerdgeWZmVnFVC/B/SlqodB2DiYgpSddx9oT2W6SR9sXA9Q5tMzNrlKpNoT9UuoD+ImJu4EukwF4FuJ0U2MdKeqRkbWZm1rmqFuBTlS4gIiYm7YbvGWV/Cvgr8AvSZSFtP71vZmbVV7UAn6HEm+Z7tNcmBfY6pMtULgZ2JB1r+6BEXWZm1r2qtgb+MrCcpMeb/D4TA0sAawFfBJYBriOF9iWSnmnm+5uZmU1I1QL8QOCzkrZs8OsG8DlgdH6sBjxPupLzMuAa9xw3M7N2UrUAnxp4FFhH0l11vE4An6E3sNcA/gtcA1xNCuzn6q/YzMysOSoV4AARsRXwY+CLkh4cwdfNRW9gjwYmIoV1T2A/1fhqzczMmqNyAQ7/C/GjSddp3jTAzwcwN7ACvePqujIAAAD0SURBVCPs6egdYV8NPNpJPdXNzKy7VDLAASJiHeBc4GDgH6SbypYj3Ze9HPABqW96T2g/IOnDMtWamZk1VmUDHCAiFgf2J21AGwfc2ufxL4+wzcysU1U6wM3MzLrVRKULMDMzs5FzgJuZmVWQA9zMzKyCHOBmZmYV5AA3MzOrIAe4mZlZBTnAzczMKsgBbmZmVkEOcDMzswpygJuZmVWQA9zMzKyCHOBmZmYV5AA3MzOrIAe4mZlZBTnAzczMKsgBbmZmVkEOcDMzswpygJuZmVWQA9zMzKyCHOBmZmYV5AA3MzOrIAe4mZlZBTnAzczMKuj/AVfT71aaBaHpAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%time gplt.polyplot(nyc_boroughs.geometry.map(lambda shp: shp.convex_hull))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can perform arbitrarily complex geometric transformations on your shapes this way. However, [most common operations](https://geopandas.org/geometric_manipulations.html) are provided in optimized form as part of the `geopandas` API. Here's a faster way to create convex hulls, for example:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 55.6 ms, sys: 6.45 ms, total: 62.1 ms\n",
"Wall time: 39.9 ms\n"
]
},
{
"data": {
"text/plain": [
"0 POLYGON ((-74.24712436215984 40.49611539517034...\n",
"1 POLYGON ((-73.94073681665428 40.54182008715522...\n",
"2 POLYGON ((-73.98336058039274 40.56952999448672...\n",
"3 POLYGON ((-74.02305574749596 40.68291694544512...\n",
"4 POLYGON ((-73.87830680057651 40.78535662050845...\n",
"dtype: object"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%time nyc_boroughs.convex_hull"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is beyond the scope of this short guide to dive too deeply into geospatial data transformations. Suffice to say that there are many of them, and that you can learn some more about them by consulting the [geopandas](https://geopandas.org/) and [shapely](https://toblerity.org/shapely/manual.html) documentation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Defining your own geometries\n",
"\n",
"In this section of the tutorial, we will focus on one particular aspect of `shapely` which is likely to come up: defining your own geometries.\n",
"\n",
"In the cases above we read a GeoDataFrame straight out of geospatial files: our borough information was stored in the [GeoJSON](https://geojson.org/) format, while our building footprints were a [Shapefile](https://en.wikipedia.org/wiki/Shapefile). What if we have geospatial data embedded in an ordinary `CSV` or `JSON` file, which read into an ordinary `pandas` `DataFrame`?"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"It is extremely common for datasets containing light geospatial data (e.g. points, maybe line segments, but usually not whole polygons) to be saved in a non-geospatial formats.\n",
"\n",
"In this case can import `shapely` directly, use it to define our own geometries, then initialize a `GeoDataFrame`. The `pandas` `apply` function is the best to do this:\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 POINT (-73.950057 40.767373)\n",
"1 POINT (-73.90903900000001 40.86267)\n",
"2 POINT (-73.961275 40.716507)\n",
"3 POINT (-73.987768 40.749788)\n",
"4 POINT (73.96049599999999 40.702401)\n",
"dtype: object"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from shapely.geometry import Point\n",
"\n",
"collision_points = nyc_collisions_sample.apply(\n",
" lambda srs: Point(float(srs['LONGITUDE']), float(srs['LATITUDE'])),\n",
" axis='columns'\n",
")\n",
"collision_points"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From there we pass this iterable of geometries to the `geometry` property of the `GeoDataFrame` initializer:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
LATITUDE
\n",
"
LONGITUDE
\n",
"
DATE
\n",
"
TIME
\n",
"
geometry
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
40.767373
\n",
"
-73.950057
\n",
"
04/16/2016
\n",
"
4:13
\n",
"
POINT (-73.950057 40.767373)
\n",
"
\n",
"
\n",
"
1
\n",
"
40.862670
\n",
"
-73.909039
\n",
"
04/16/2016
\n",
"
4:30
\n",
"
POINT (-73.90903900000001 40.86267)
\n",
"
\n",
"
\n",
"
2
\n",
"
40.716507
\n",
"
-73.961275
\n",
"
04/16/2016
\n",
"
4:30
\n",
"
POINT (-73.961275 40.716507)
\n",
"
\n",
"
\n",
"
3
\n",
"
40.749788
\n",
"
-73.987768
\n",
"
04/16/2016
\n",
"
4:30
\n",
"
POINT (-73.987768 40.749788)
\n",
"
\n",
"
\n",
"
4
\n",
"
40.702401
\n",
"
73.960496
\n",
"
04/16/2016
\n",
"
4:50
\n",
"
POINT (73.96049599999999 40.702401)
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" LATITUDE LONGITUDE DATE TIME geometry\n",
"0 40.767373 -73.950057 04/16/2016 4:13 POINT (-73.950057 40.767373)\n",
"1 40.862670 -73.909039 04/16/2016 4:30 POINT (-73.90903900000001 40.86267)\n",
"2 40.716507 -73.961275 04/16/2016 4:30 POINT (-73.961275 40.716507)\n",
"3 40.749788 -73.987768 04/16/2016 4:30 POINT (-73.987768 40.749788)\n",
"4 40.702401 73.960496 04/16/2016 4:50 POINT (73.96049599999999 40.702401)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import geopandas as gpd\n",
"nyc_collisions_sample_geocoded = gpd.GeoDataFrame(nyc_collisions_sample, geometry=collision_points)\n",
"nyc_collisions_sample_geocoded"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"In most cases, data with geospatial information provided in a CSV will be point data corresponding with individual coordinates. Sometimes, however, one may wish to define more complex geometry: square areas, for example, and *maybe* even complex polygons. While we won't cover these cases, they're quite similar to the extremely simple point case we've shown here. For further reference on such a task, refer to the `shapely` documentation.\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Joining on existing geometries\n",
"\n",
"Sometimes the necessary geospatial data is elsewhere entirely.\n",
"\n",
"Suppose now that we have information on obesity by state."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"We'd like to put this information on a map. But we don't have any geometry!\n",
"\n",
"We will once again have to define a geometry. Except that this time, instead of writing our own, we will need to find data with state shapes, and join that data against this data. In other cases there may be other shapes: police precincts, survey zones, and so on. Here is just such a dataset:\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAEXCAYAAAAgInRqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzsnXd0G+eV9p+ZQa8Ewd6rRFG9y5ZkW7Ic994dpzvF6zjJJpvNZjdf2maTOJvenGpvYseJ7dhxYjvuTVaX1Tt7r2BDbzPz/QESIgkQGAwKSfH+zuEBMAAGQwLEM/e+9z6XEUURBEEQBEHEhp3tAyAIgiCI+QAJJkEQBEFIgASTIAiCICRAgkkQBEEQEiDBJAiCIAgJkGASBEEQhAQUce6nnhOCIAhiocFE20gRJkEQBEFIgASTIAiCICRAgkkQBEEQEiDBJAiCIAgJkGASBEEQhARIMAmCIAhCAiSYBEEQBCEBEkyCIAiCkAAJJkEQBEFIgASTIAiCICRAgkkQBEEQEiDBJAiCIAgJkGASBEEQhARIMAmCIAhCAiSYBEEQBCEBEkyCIAiCkAAJJkEQBEFIgASTIAiCICRAgkkQBEEQEiDBJAiCIAgJkGASBEEQhARIMAmCIAhCAiSYBEEQBCEBEkyCIAiCkAAJJkEQBEFIgASTIAiCICRAgkkQBEEQEiDBJAiCIAgJkGASBEEQhARIMAmCIAhCAiSYBEEQBCEBEkyCIAiCkAAJJkEQBEFIgASTIAiCICRAgkkQBEEQEiDBJAiCIAgJkGASBEEQhARIMAmCIAhCAiSYBEEQBCEBEkyCIAiCkAAJJkEQBEFIgARzDvLykQ609tshiOJsHwpBEAQxjmK2D2Ah4/EHoeRYKLip5y2PvnUOoy4/jFolVpRlY0WFFSvKrajIM4JlmFk6WoIgiIUNCeYs8oe3G3CyYxhfvHElynON4e0cGxJFhyeA3ef6sftcPwDAqFViWWk2lpVlY1mZBTUF5gixJQiCINIDCeYsolFyaOwdwwO/3YWPbF+MmzdWgmUYcGx0EXR4Atjb0I+9DSEBVStY1JVYsLTUgmWl2agryYJerczkr0AQBLFgIMGcRXJNGgBAgBfwm9fOwKRV4YqVJXB5A5Ke7wsKONY2hGNtQwAAlgG+edd6rK/JS9sxEwRBLFQonzeHGLR74PIG4PIFZT1fEDEltUsQBEGkDhLMWSQoTK2CtTm86B/zyN6f1agOR60EQRBEaiHBnEV4Xphye8juxcGmAdn7W1JsAUNVtARBEGmB1jBnEYdn6lrlkbYh7GtMQjBLLMkeEkEQBDEDFGHOIu0255TbvgCf1P6WlGQl/Jwxtx/PHWjFyY5huGWunRIEQSwEKMKcRaxGNbQqDh5/ckIJACzDoKbAnPDzznaP4OFXTgMAGAC/vf9SlOYYkj4egiDSjz/Ig2PZcO82kV4owpxFHrhqGR59YBtWlGcnva+qfCPUSi7h5zX12sPXRQB5Zm3Sx0IQRPp493QvOm1OdAw68NlH9uDTv9uF4+1Ds31YCwKKMGcZi0GNb929Abf/4LWkUrKLihJPxwJAU99Y+HqOUSNLdOc7jb1jeO5AK/RqJfQaRfjSoFZCN37bMGm7SrHw/kbE3OHxnY3oHXEBCPViA8D3/34Mf3hwGxX9pRkSzDmAWslBwTLwJbGPumK5gnk+wizK1iVxBPOXTpsTrx/vlvx4JcdCr1FAo+SgUnBQKViox6+rFSyUCg5qJRu6reSg4liolByKLDpsX16cxt+EuNBp7bejbdARsb1/zIPuYRdKrLSckk5IMOcI/3LVUpzqHIHVoEbboBO7zvQikVkli2VEmGNuPwYm9X0WWfQJ72MCcXyyynw8w3VIdFaaIMALGHX5E34do1aJLUsKKEIlZPPWyZ4Z7zvUYiPBTDMkmHOEHStKsGNFSfh2+6ADKgWHnmEX/EEBi4rMaBtwwKBVornPjtNdIzjTOYKu4VBqRquSsX45KR0LACadSvbxj7r8+Piv3kGJVY9SqwGlOYbw9QKLDso5bBIv1YowWRyeAByeAKxGEkxCHhM+0tE43DyIG9dXZO5gFiAkmHOUCYu7Qsv5NKnVGHLxWVyUhWvWlAEIRYlnukbC9yVC86R0LAC8e6YXH9m+WNYIMZvDC4cngDNdozjTNTrlPpZhkJ+lRVG2HkUWHYqz9SjO1qPAokNBlnbWI65EI0y5KDkW2QZ1Rl6LuDBx+2du/TrWPoQgL9AEozRCgjnPMetU2LQoX9Zzpwtm/6gb/qAAjYzCn0H7zJZ+giiid8SN3hE3DkW5P8eoCYtnYZYOeVla5Jm1yDfrkGPSpD06dXoyI5i8IM7LlDUxh4ixTuPx8zjTPYrlZclX3RPRIcFcwLT0TxVMQQS6bE7UFCbezyknwp3A5vDC5vDiZEfkfcz4vh+8ZpnsE4N4ODMUYQqiCJcvQCPYCNkIYuzKhqf3NKOmwAStir7a0wHF7gsUf5BH15ArYvs3nz6EYac34f1V5ZvSEgmKCAmqUpG+j2qmBBMI9dARhBxEUYQ/GLv1bH/jAD79211onlafQKQGEswFSvugM+rZav+YB8/tb0t4f0qOxaevXpqCI4tOril9hgrTPX3TyavHulKynyGHFy8cak/JvhLB7vbjFy+fxMd+8Tb+ebgj7hc4kTp6ht1weuPbV3YNu/DZR/Zg1xk6OUs1JJgLlFhnoK8d7wIvCDPePxNXrCxFfpqcgnKSSPnGI5MR5rnuUfBCIg1DIRp7x3C01YbDLTa81zyIfQ39+Nv+1jQcYXQCvIDnDrTiI794G/842I6uYRd+8uIJfOTnb+Nv+1vhTdIHmYjPyc5hyY8N8AJ+8PzxmLUFROJQonsWONg0AJc3iK31hbPmAdnSH9n8PMGw04d/Hu7A9esqEtonxzL4z1tX45tPH8KQIxkbhqno1Qro1On7qLoknLWniqAgYsTpQ06Cc0ufO9Aa1Vzh+ffasaEmFzkmDTg2Pee/Z7pG8Pn/24NoOm9zePGrV0/jL7ubcOumKly3tjyt79VCpblvDI++eS6h59QUmJClp6rsVEKf7FngjRPdeOtkD14/0YUv37J6VopAXL7YUdXDr5xGc58dn7tuRUL7zdKpMZxCsQSSKyiKhyCKMUv104FPRhpzplafn790EkAoJV5TYMKSUguuXl2GsnEDfV4Qkz4pe+tkT1SxnMyoy4/fv3EWT+1pxk0bKnHj+goYtVTclCreOtmDEZf0/6tSqx5fvX3dnO5/no+QYM4C7YOhsV4HmwZxpMWGLUsKM34MQpxvQF4QZfnTNvaNJeRQJIV09i76ZyGVKMcoQa+JLT4BXsCZ7lGc6R7Fs/tasb4mF/6ggNOdIyix6rGuOhdb6wsTdoTiBQHvnJ7ZXWY6Dk8Aj73TgGf2tuDGDRW4ZWNlUoYYRIhNi/Lx9N4WSY/N0qvwrbs30AlLGqDTjwzj8gXQOWkOpj2DBSeTqY3TOmLSKrFjReK+py19dmxfViT3sKKSzrTSbKy99Y0mvq50SX1iJ1UHmwZxrG0IAV5A64ADT+9tiei7lcKR1iFZNoBufxB/3tWED/7sTTzyxlmMuRPfB3GeJSUWmCWceFiNanzn/RtRYFmYvtDphiLMDLPrTB8C/PmCmp+8eAJPvNsIuyeAgiwt7ri4GqIIlOYYwAsCeobdoQ+/KKK2KEuWqUA0rl9fgZ2ne3GmezTq/WadSpYDzwcvW4T+UQ/ejOF5mShZ+vRFKN4UzCJNFIcncfFYUpwFg0YhqUoyGmolh0uXJp7JePOEdFP6aHj8PJ7c04y/H2zDdevKceumSmQb0pdiv1DhWAYba/NiVllX5Brx33evpxF9aYQEM8O8EeULaNAe6ntsH3Tif/9+bMbnciyDqnwT6kssWFKShTVVuZLOOqOhYBkYYqRstiYY0UzAjNvgleUY0DEpkk6GdH7BzkaEKSdiYxgG2QYNnN7E/qabavNw15YaaFWKhNfKvQEee871JfScWPv6694W/P1AG35+3xZU5BlTst+FxEWL8mcUzDVVOfjKrWvipu6J5CDBzCADYx4cb5M/6JUXRDT2jqGxdwx/Pwisrc7F/9y9Xpbd2s7TvTjYNDjj/SvKrTjdNYKBMQ8uW5pYipVhGFyzpgy/evV0eBvLMPjaHWvRP+pG76gHvcMu9Izb5U2OuKNRapU/RSUe3kBmC34A4Gz3iKzn5Zo0kk9ClByLD21bhNs2VUn6fAyMeTAw5oEgiqgvsUDBsdh7rg+eFEfgAV7AC4fa8cBVS8kmMAbRirXWVOVApWDhD57/f2EA3HtJLe7eWjtrFfcLCRLMDPLWyZ6UFsQcah5E15ALpTmJjfTxB3n87o2zMR/zkxdPwGb3YlWlNWHBBIDr1pVjyOGFQaNEh82JlRXWqNZ2gijCZveid8SNgTEP+sc8GBhzo2/Ug/5RNwbGvChJ8PdLhNlIyR5pHYIvwCc8rFuqecOysmx8/voVKM6WdqLx9qke/PD54/AFeJRa9fjt/ZcCSD4dOxPPv9cOpYLFjesq0DvqxurKnLS8znxl77l+vNc8gAevWT5lu0alwJrKHOxrHAAAWPRqfOnmVfT3yyAkmBlCFEW8cSI1Li8TXLWqNGGxBEJfWJPnYEajd8QNALC75RUlKTkW9+1YEvdxLMMgz6ydcd2FF4S0RiJLSiz4zacugcsXhMsbgMsbhMsXgHP80u0Lwuvn4fEH4Q2MX/p5eALnt/sCQlyPz8kEeAFjbn/Ca02Lisx4+Whn3MdV5hkliWWQF/D7N87i2UkGCDtWlIBhGIy6fDiWRDYkHs/ua8Wz+1qxdUkhVlVYKdqcxKDdg70N/fj01csi/i4XLc4PC+YDVy0lscwwJJgZornPHm4nSQUGjQIfu7wu4ee9faoHv48TXU6msXcULf12VOWbEn6tVJCuZvwJ1EouPEotGYK8AF+Ahy/IwxcIXWcZQMGxUHAsOJYJXWcZcBwLtQxv3KWl8adQsAwTtwIaCH0ef/j8MTRNqpxlAGxfHqqMNutU+PWnLsWHf/5WwseZCO+e6cUP/sEh16TBHZuryTQcIaOOIYcPTX32iPcylKU5ASXHYm117uwc4AKGPp0Z4oldTSnbl0rB4t9vWpVwf1v3sAs//MexhKzZBBHoG3HPmmDOFyaEUY/0FV2ckbD2+eA1y3DlqtKYj/nec0fx9qmeiM9BtlGNvlE38sxaMAyDQosOP/rIxTjRPoxTncM41TmSFhvB146HMi+rKnOwssKa8v3PdQbGPDjVOYzVlTnI0quRM5563984ECGYWXo1Ll9eDLvHT45KswD9xTPAoZZB7D6bmmpDAPjI9jpsrE181NWZrhH4gol7xO4624eVlVYaSzXLHG6xxX1Mk4QpFfsbB6KeNA05fPjBP47h89evRF1xFtRKDvUlFtSXWABUgxdEHG2z4dWjXTjZMQy1ksPAmCdu0ZZUDjYNoDLPuGCMDl460oHHdzbCNl4l/9AHNmKRgsM3n34PALC/oR/3XlIb8bwv3LAiprUlkT5IMNNMgBfw8MunUrtPGaIHAJcvL8aTu5sTbvd440Q3Fhdn4cb1FRH3iSINRc4UH9m2GDvjjAeTYoWmVytmjBT7Rj3498f2QcEyWFSUhU2L8rBjRQmsRg04lsHaqlysrTqfCvzpP0/gxUNRBpnK4Om9LXjhUDtuWF+BWzdVyW6Zmi+4fcGwWGbpVcg2aKBUsOFe24beMQw5vBHWkBzLSkq7E6mHnH7SzLG2IXRGmTuZDHLTYg5vAP2jblnPPdg0gBcPteO1Y1147J0GtA040D/qxn/9+aCs/RGJI8UikI3TWtA34pbUexoURJzuGsEjb57DvT95Ex9/+B187BdvR4xC+8w1y/FvN6yMuz+pePw8ntzdjA/97E3883BqhHiusqEmL3y9vsSCshxDxP/ngaaBTB8WEQOKMNNMOmKvK2RY1gHA7rN9slKyQMhqbXLf5uM7G8PX95zrw8WLC2Ttdy4Q5AWwLDOjwflcQaNSQKdSxDSLf+dUD1ZX5GB9TW5E5L+voR8PPXcUbl9ivaeCKIazEvf/Zie+cec6VBecj3AmR0BKjkVxth5jbn9CZuHT8fh5PL6zAVevLr1gMxi5Zi3USg6+AA/L+MnQ9OzPG8e7cWl9Ea1XzhHoXUgzihRPCzBplbJaSYDQLMZ0MN//mf+8qwmP72yETqWATqOAXh1yxdFrQpcaFQetSgGNkoNWxUGj5KBRKaBVKaBWslArQ9tUCg5qJQe1ggPDhJrPg7wAXhDD1xUcm5TLjVGnjCmYQw4f/t9fDuKmDRW4/8rzA71FUcTDr5xKWCynM2j3YtDuRXWBGbwgoKXfAVEUUF9iwemuEdx+cRU+dNlieAM8/n6gFU/uboZL5msOOXxo6bdPEecLCY2Sw6X1hXj1WBfyxgt9jNOcek50DOOZfS34wKWLZuMQiWnM72+6eUBZjgG3X1SFY+1DaOodizsmKR71pdmyz7gbe+MXhMih1Jo+Y4FMYHOE1pHc/iDc/iDil9bIZ9uyIvzHzatlP19qFPzcgTZ4AzweuGopVAoODMPgtouq8POXkl9Pf/tUDzYtysc7p3rx0HNHoeRYlOUaUJ1vxF2bawCExOCS+iLJEzZm4r3mwQtWMAHgylWlePVYV9icfro5xbrqXNwQpXaAmB1IMNOMxaAON/C7vAGc6BjGsfYhHG8bQnOfPWHnn6WlFlnH4Q/yaB1IfWWdXq1I6/itTDBReJFuSq16fPHGVRl5LQB4+UgnrltbHi4Q0aWox/FQ8yDOdI3gJy+eABAqbFtRno1r15SH3YtEUcRDzx2JWPNMlINNg7hzXITnA7wgomPQgVyzFgYJvq6V+Uaw49kIIDTG7dL6QtQVZ2FJiQXVBSZZQxCI9ECCmUH0GiU2LcoPW8TZPX6caB/GsbYhHGsbQttgfEGTK5gt/faE+i+lUppjmPdrTEOOzAjmiMsHpzeQVPXn9evK4fIGwbIMlBwLlWLihwtfKrjQeizLMijIOj/myZGiHkq7J4DPPbonfDvPrMVHt9dN+WI/1TmCM13JLwGc6hyByxuY86birx7rxKtHu9DYOwZvgMe379kgyVhAr1bir198X7hly6hV4j9vXZPuwyVkQoI5i5i0KmyuK8DmulDBzKjLh1OdIzjVOYzTnSNo7B1DcJLIKTn55eQNPWlKx6bR5zVTDGYownR6g/j+P47hv+9aL3sft26qkv3cVK+nT7BlSUFEFHSmS57B/HQEUcThVhu2zsKQdam8crQTP3z++JRtzf12yU481N88fyDBnENk6dVTBNQX4NHQO4bT4y4rAGSnZ9ImmPN8/ZIXxLS418zEmIzRXqlCSo+mHBZFOYlrkjGseibeaxqcs4LZN+LGw69ErgsfbRvCHRdXz8IREemEBHMOo1ZyWF6WjeVl8T1E49HQm54K2TIZEebv3ziLd071IM+sRX5WyHi9IEsXvsw1adIWDU0n2arRRHF4Z08w37eyBGNuvyTz/URYVJgVsU2K45BUDjYPzEmDDEEU8YPnj0UdgdY+6AAvCGn3QiYyCwnmAsDjD6IzRcOcp1Oak/isyu4hJ/rHR3mdiNKbzjIMck0a5GdpkW/WwWpUI8ekgdUY+skxapClV6dk/p/Ll7noEgCcSRbBJAPDMLjj4mrcsL4Cv3v9DJ5/rz3pferVChRl66Zsc/uC6E6hWceQw4fWAcec8zP+5+EOHG8fjnqfze7F4RYb1k8yJyDmPySYCwCHJ4D1NXlo7rOHWyhSgZJjUWjRxX/gNGyO2A3tgiiGBRWI/oXEMkxYSHOMGiwpscha33NlMB0LhFyaBFGcVZMEjZLDp95Xj54RNw41zzxEXArRir66h10pnfsKhJym5pJg+gI8XjkSe9TarjN9JJgXGCSYC4A8sxbfHC80GXX50NxnR1OfHU19Y2jus6N7WF40UJStk5VySkVVqiCK4SZ6APAGeHmCmeGUrCAC/gAPzSyPsVJwLL5x5zo8s7cFT+xqgk+CXV40LPrIlqLRJBx+ZiJThVkz0WFz4kDjAJr7xtDUZ0fXkDNuT3W7hKr3WIiiiIExD1473o0lJVlTPHyJ2YEEc4GRpVdjbXXulAo+ty+Iln57+Mugw+ZEx6AzpqMMIK/gxx/kUxrlTpClk9cLmsmCnwkytT4bDyXH4q4tNegZceGVo/KGmxu1kRWeo2kobKorjlwnzRQOTwCfe2R3widXgzI/525fEN977igaekcxNJ6NuWZNGZaWWMCyDPVlziIkmAR0agWWlWVj2aTiIlEUYXN40THoRLvNiY5BBzpsTrQPOsLTFOQU/PSNpq7YZDJZenm9jS5vZiNMBkjJ2msqkXuyAQD+KN7E6fj1FhfNnmAatUpcsbIEzx1oS+h5ctPu/iCPvQ390Ko4LCvLRqFFh38e7kCWToUCiy7uvFMifZBgElFhGAa5Ji1yTdop0agoihhx+dAx6IwYOySFvhF501LiYZYpmO4MF/0oFeycq/asl2mGAQCj7sj0a6pNBnQqBYqtiReXpQpeECVNeJmOKMpbyTVolPjd/Zei2KoHyzAQRRGlVj30GiV2rCiRtU8iNZBgEgnBMAyyDRpkGxIXSwDoHUntqLMJ5KdkMxthpqsXMhlWVeZAq+KitkfEo30gsvpaTiFYLDbXFcxakVSQF/DVJ9+LWxxVlK2DWavCmUkDDoYcPgR4IeH3XMGxUwxBXL4geobduPWiqjmXnVhokGASGaV3rqVkMxxhblky98ag+QK8LLEEQnZ/Y27/FLu/EqseKgUbNV0rh5s3VqRkP3LgWAZnJ7kWKVgG5blGVBWYUFNgQnW+CVUFJujVSoiiiB+/eAIvj1fPCqKI3mEXynLlT6cBAK2Kwzune9Buc+B/7t4w520CL2RIMImM0pumlGxWlGpNKYgiUvrlHo+5OHnDpFWi1KqXPei8dcCOVRU54dscy6I815iS6TgryrNn9W/GMAw+vG0x9GoFqvJNKMkxzBgxMgyDj22vQ5FFj38cbIPNEariVik5vHS4A019dnzltjXQxqiQHnJ40dxnR65Jg8rxNhqOZVFfmg2PLwhvgCfBnEWYOHn21Lt1EwuaT/5qpyST+UR5/LPbI0YjJYI/yMPlDcLpDcDlC8A5ft3jD8LjC8Lt5+HxB+H2BeH1B+ELCvAF+PCPN8DDHwxFahOjmqKxtNSCH374YtnHmS4aekbx511N2HOuP+bjOJaBTh2aGWrUqmDUKnHbpqoI39QfPn9MduXtZP7fbWuwZY7a4sUiwAv4w1vncKTVhg6bE9euLcc9W2pgmsF4f19DP378wonw0O3r1pbhwWuWh+9v7htDZb5pzg85v4CI+oemCJPIGKIoonc0TUU/SUwAAUIevSoDB0uSo8pEUcTz77XjqT3N8Ph5KDgGgiDCPu7wk2+WL+rpZFFRFr5082q8eKh9XBCV0KkVMGgU0KmVMIwP01ZJLFpKhcmAWsFi3Txt/FdyLHasKEGeWYsNNXkoiLOua9Aow2IJACPOqcVUczEzsRAhwSQyxojLJ7tBPhYGjWLO9KYxDIMb1lfgunXl4AURSo6FIIp4ek8z/rK7GeskTrAYdfnC49gmkkCCKCIQFBDgQz/+IB++7Q8KCAQF+Hke/mDodpAXwAvi+I8AQRBRYjXgipXRKy01Si6paSiTScX7sboyBxrl3Hhf5VCRZ0RFnrT1y4o8I5Qci421ebhxQ8WUFi9i7kCCSWSMubZ+mU5YhgHLMeHrd26uwZa6Qjy9txmXS2gNuP8372LYmXrHHAXLYFWlNan0tRRS4fazOANmBb4AjxGXb8rc0ESYcM5yeAO4pL5QdspUEEV88n1LcN3a8jnXdkSchwSTyBjpEszsJNOomaLYqsfnrlsh6bH5Zm1aBFOl4JAjo382UWKt40oh36zFLRsrU3Q00emwOfHdZ49gy5IC3LO1NuZjeUFA15ALrf0OtPTb0TJgR3Offcp79OKhdnz33o2y7CJNWhWuX1eR8POIzEKCSWSMhRRhJstly4qm9PSlCm+Ahy8opD3VOZM9Hscy4YkzVqMGLl8AzX125Ju1yDVrkWvSoMiiw6rKnLT67R5tteFLj+8HAJTlnu95FMc9itsHHWgfdKJt0IHWfjvaB50I8LErqY+3D6Ol3yF7yDsx9yHBJDLGXI0wD7fY4A/y0KsnFbhoQkUvs1WV2JzCAcyTEUQR+xv6cenSorTsf4LNdQVYVWmFQa1Ell4Fs16NLL0KRo1ySspxtuZcZk+Ksk90DONHzx8PiaTNmdSM1C/83x589PI63LQhvdExMTuQYBIZYyQNUyyA5CPMR948O2PPoE6lgFbNQatUQKtWQKvioFUpoFFyUCm50KWChVrBQT2+rThbh421+bKPRxRFvJfk2K1Y/O/fj2HPuX78240r0+Y8dEm9tFaQdImlwxOAQaOYcf9FFh04lgEviLDZvXj5aOxRXVLxBYVZNYon0gsJJpExvn3PBtg9AXQPu9A95ApdDrvQM+xC/5gHDpnDlS0yXX4miDWxxO0Pjk9tSUzsN9bmYeuSQoy6fGAYBpX5RsnjmRiGwU0bKvHIm2cTek2pBHgBb5/qQV1xFm5O8zphJhFFEfsbB/D68W7sa+jHHz+zbUYLR8X4LNeuFA66BoBVFVYsmkWjeCK9kGASGYNhGJh1Kph1KtSXRBp+e/xBDIx5MGj3jl96MOTwwubwYcjuhc3hjSpuyfZOpmOI9P7GAexvHJiybfuyIty3Y4kk0/rVldaUH9N0fvXqaZRY9RfEkGN/kMcPnz+Ot072hLf1DLtjeh6XWg0pF0wFx0ZN408MLZDrwUzMDUgwiTmDVqVAea4R5TG8N30BHkMOL4acPgw7vBhx+VCZJ79JXhTFjBmwv3myB93DbvzLVfWoK449IaSm0Ay9WpH2AdffeOoQvnbH2nktmna3H19/6j2c6hyZsr1nxBWzn7E0x4C9DbGdjRJlJuP5pj47HvzdLly8OB93b62lwqB5ytwbnUAQMVArORRl67G8LBuXLi3CTRsqkZeEe47Hz0OQOYZJDud6RvHS4fjrZSzDZETEAryADlvkxJG5ypuqkevSAAAgAElEQVQnuqccb9eQE599dHeEWAJAd5zosTQn9SPDXjzUgYNNAxHb3z7VAxHA7nP9eO5Aa8pfl8gMFGHOcb7/j2NweAJYXGRGXbEFi4rMMJD5csqItX6ZLk53jYAXhLj9ep96Xz3eGf+iTScVSU7TyBT9o2788pVT8AcFXL26FGadCs/sa53xPeyJU5Vdak18AHo8BFHEN546hM9eu3yKo9Lkytvl5OIzbyHBnMOI4y0Adk8A+yaljkqteiwuzsLioiwsKjKjKt80Z6zh5hvpWL+MR4fNidePd+PKVaUxH2cxqFGcrUfXcHpmiALA4qIsrKnKifu4IC/A5QsZ0rt9Qbh8oUuHJwCHJ4A8s1ZyZSwQ+mzbPQHY7B7oNcoIp52jbTa81zSI3hE3BsY86Bt1h/14AeC5A21xX6Mnzt+tJA2CCSBqv+ZHti3Gma4RtA448Ie3G3DV6rK0vDaRXkgw5zDdw64pXxITdA650DnkwuvHuwGEmsErco2oLTJjcVEWagvNYW9KIjazEWECQGPvWFzBBIBsozotglmUrcPaqlzcsrEyovXiTNcIfvv6Gdjdfrh8Qbh8wbgewNevKw8LZs949fOoyw+Hxw+HJ4BRtx+jLh9GXD6Muvyw2b1hYbn9oirct2PJlP0dbxvG03tbkvode4bdMfs8jVolLHp1ytudrltbFuHXa9Kp8JHti/HVv7xHE0fmMSSYc5gzXdKcXnhBRHO/Hc399vDwWiXHoirfhC/fsnrGQgQCGSv4mY7UdcO7ttTgePuBlL724qIs/PRjm2e8f0mJBSsrrHji3SbJ+5xsHvHXfS148VCH5OdGswDMMSVfTer2BzHm9sfs0y3N0adMMDfV5mHLkkJsWxZpCjEw5sEP/nEcAKBU0InsfIUEcw5zuiuykEEqAV5AQ89o0mOvgJCPJsCAYy+8M+PZijCPtQ2hw+ZEWU7stODysmyoFSx8KRxw3TfqRv+oG/kxDMc/dNlinO4awdHWIUn7nCyYZm1in7khpzdiW24KBBMIZWmy9Go4PAEcaOyPML4vzTHgePtwwvtlGaAoW4/qfBOqC8yoKTShvsQy43Don790Muyvu3UezvckQpBgppFkbb/OJCGYQOgfWqdO/i3e3zCAbz97BMXZepRYQz+lOQbUFJgljy+aq6yvycX3PrAJLl8ALm9w6qUviDdPdMOfQrGazMmO4biCqVJwqC3KwsmOxL/UZ2LM7UeHzRlTMAHAJTH6vmdLDdZMMmUwJ2gkMeyIEmGmwCBeq+LCU1l0agVOdo5ECGa8vz8QmrVamW9Edb4JlXkmVOUbUZpjSKhu4LaLqtDcb8eQ3Ytr19D65XyFBDONHGgawMOvnMbysmwsL8/GijIr8rO0kkTU5QugbcCR1OvXFCQ/xBcIpQ8DvIC2QQfaBs8f09YlhfjKbWtS8hqzRZZeHTNlt+tMb9oEs1/iMG2TNrVV0bkmTdxCn3M9ozPaBU5nZaV1SmtPolmNIUdkhJmTgvFjHj8P9bjJPMcycHgiDeEnV8qatEqU5hjGe4HPX1r0alknvqIoonvYhaJsPVaUW/H4Z7bj1WNdcYdJE3MXEsw0crZ7FL0jbvSOuPHqsS4AobWZ5WXZWFFuxbKybJRa9VH/Gc91jyXdTpCqKe0zrbelo49tLiGKYlJG3PGQGv17/Kkdun39uooZW1qGnV4carbhyd3S1y8d7qlpbbMuMeclly8Irz84ZTqJQaOAWsklPXBcnNRj22lzgRfEKUsLdcVZ+P4HN6E0x5DyqTe7zvbhuQNt+P4HNwEIOV1JKfQi5i4kmGnkcIstYpvN7sVbJ3vCFl4Pf2IrqvIjI8Fk07EAUFOYuggzGgNjHhxttaE814gsveqCG3zrC/AQ0tgE2WmTVv165aoS5GdpQ0OpmdAXL8cyUHAsOIYBxzFQK0LG7yoFGzaD12uU0KsV4cvHdzYiwAu4c3N11NcRRREP/m43bFEivljYp0VuctbNh50+FGWf/zpiGAY5Rg26k6gQVitYmMaPRRiP9gbHPFMiPL1GieXl6bEhzNKpMJKGmabE7EGCmSZOdgzjbJx5hnq1YkYbuDPdKRDMFESYgiiicwbBfP14d7i1xahV4iu3rsGqyvg9ffOFdNvS7W3oQ/9obdy1xG3LirFtWXHSr/fpq5ehe2jm6lyGYXDFyhL8eZf06BKIHBadJcMMf8jpQ1H21IyF1ahOSjCLsvXhFo4hR6iNpXPImVBKdGJyTK5Ji9IcfULDoRcXZ2FgzIP9jQPYtEj+9Bpi7kCCmSbeONEd9zH1pZaolaeCKEpuKZmJHJMmJRWyg2MeeCWkxRyeAIwpXmubbTz+9Aqm0xvEN546hJ/dtzmhL2K5cCyDsjiuPv5g4inQ6VNmTHIizGjrmEkW/hRNEsaOwdCJQvewC+sRMmLoG3Wj0+ZC15ATS0osEb6zoijie88dxZvj2aDNi/PxX7etDf/P2t3+8IDppr7Qz6feV4+VFaGIVaXgsKjIjMfeacDG2rwLLgOzECHBTBPrqnPxz8Oxe9GWlka3yBJFIM+sTarlIRXRJSC9XxAACi0X1ppmidWAl75yDXwBHh5/EB4/D+/4pccfhD8owBfg4QvyocuAgAAvIDj+E5j4CU7cFhHkBQiiCF4Qw5fNffY5MRIqwAt480RP/AdOY3pKVsmxCRvHR+/FTK7wZ8uk9o3+MQ8A4Nn9rXjhvXb0jLjBT8q333lxdVgwO2xO/OujeyCMr2GrOAZ+XsTuc/3476cPIcALaOm3Rz3mE+1DYcEEQm1Bf9ndjLdP9aQkS0DMLiSYaWK1hNTkstLoEys4lsEXrl+BB3+/W7YxeCorZKWQbVCnpIVlrsEyDLQqxYz9dRcS3/rrYVlN/NHcqMx6VfKCaZRXhMMyDK5eUzrFQGCiLad/1BP1OROp32f3t+KxtxvGZ6CO749lUZKlwZDDi/JcA0QRMw74PjnNBH55uRV/2d2Mh/52FN3Dbtx7Sa2s34mYG1z43wKzRLzFfgXLxIwqagrNuGtzNZ5IcD1pguoMC+b09SdifiGKIo61RRapSSFau4ZZp0LPsLS2GSBUnTsdKXNDp1Ni1eMX922ZUnH7XvNg3CWS7mEX3j3di1+/ejrivq/ctgbra/Lg8oZ6c61GNeqKs/CNpw+FH1Oea0D7oBNnphnr15dYwDKAIAJ/3duMuzZXQ0GWlfMWeufSxFiUL5HJ1Baawz1iM3H31hqU58oziE5VSnamgp/pFFFv2bxmYMwju31lelsJkHhryVA08wIZKdnLlxeHxXJgzIOdp3vxveeOxn1ep82Jh6I8bl11bnjMml6jRJ5ZC45lsWlxPn758a340GWLAAArK6yoyDXCG+DR0n++V1mnVqCmwIxckwarK3Mi0tfE/IIizDRRW2iGWaeKqCCcYKmEET8qBYcv3LASn3tkd0LtDSatMiXWYqIoon1QmmAWU4SZduweP1r67OBFEYIQWv/kBXHKmqkvKMAf4MOXbn8QLm8ATl/ocl11Lu7ZGpkWNGiUYBlG1hJANBHISrDwJ1qEmWjRz32X1yF33EChsXcMn3tkN4IS/3FCj5v62FKrHv91a3RjDpZhUF1gQud41bHHz2NjbR7aBh042DQArz+IUZcfW+sL8f0PXRT35JiYH5Bgpgklx2LbsqIZxxBJ9WVdXJSF2y6qxlN7miW/dk2hOSUVeaMuv+TCI0rJph9BEPHQc0ejrvdJ5VTnCPQaJW5cXzFlu16jxNrqHBxsir42FwunNwBBFKdM4Ui0Qjva75RtVOOxz2zHsNOHn7x4Ai399hmfzzIMbtpYGZ7Q8+ibZyWLZTQULIMv3rQq7rr8n3Y2AgDePtmDL9+yGgDwh7cbAISmsGytLySxvICglGwauW5t+YzC+Oy+VrQPSrO++8CltSi1Shek6ihGCHJIpEK2ODv9KdnHdzbi288cxq9ePY0ndzfjtWNdONQ8iOY+O4Yc3ilVjxciWXo1vnbHOqiSnXYxQxT5L1culTWcXBAR4YiUaIuRwxOIaGlhGQZ5Zm3YjeeiRflTRFmnVoRN3wVRBDN+KYpi0lNoPnPtclTmGfHmie6YbkP37ajD5cuLYTGo0dJvx1dvXwuzTgWLXh0x4ouY/1CEmUZKcwxYX52LfY0DEfcFeAE/ev44fvDhi+NGmyoFh8/fsBKff3SPJLu8VK1f9o9JL9rIREvJ0VYbTsQwIf/+BzelxbXlmX0tON4+DI2Sg1bFQaNSQKPkQj8qDioFF3bXUSpYqJUcFCwDjmOh5NiQK8/4bcW4Q4+CY6c497Bs6PrE9pmoK87CNWvKJA1QnomZJp8UZevx7fdvwLf+ehgDY9GrSaNh1IachCbTMyL9swOEMjJCjBMevUaJr9+5Dn0jbvSPecAyoUxKz7ALP3z+OJaXW6HgWDzxbiPaB50JfXans6zUgmf2teCP7zRgxOnD8fYh3L2lJqrBxMbafGysDZkSTETZ9SUWaNWhzwhxYUGCmWb8UaavT3CmexT/ONiGmzdWxt1PfYkFN2+sxLP7W+M+NlWWeO9bWYpL6ovQOz4QuGso9NM55ESnzRluG8hUS0m0MVCTSbUX6AQNPWPY19Cfln1Pp644C/dfWY+64ugtRwDwwcsWwe0L4t0zvbIKdQ42DeCOi6Pb4y0uysKvPrEVH/zZW5LT8QxCEeJkw4JEe4gDvIAOmzNuP2qBRTfFqae6wIxffHxr+Pbehn409EgzjZ9MfpYWZp0KDT1jyDaop7SHvHSkEx4/H065zsTEiY7FkJ7PITH7kGCmmep8U1RP2QkefescNi3KlzTk+UPbFmNfY3/Mcn2tikvpeqJGyaEy34TKaWleURQx6vKja8gJRwZmSoqiGHUM1GTSJZjeNDv+TOZs9yge39mIb929YcbH6NVKfOGGlSEv3zZp8yonc7x9GHvO9uHiuoLo+9cosboyB3vO9UlKc9s9ARxuseGySX2Pt26qxM7TvZKP6bq1ZTDpVHh2XwtOdo6gLMeAuzZXT2kPkUJA5mSZZaXZ+OKNK/HCe+347RtnI+6/ZVP8k1riwofWMNNMNGP1yfgCPH7y4okpUxVmQqPk8PnrVsR9vVgpvVTBMAwsBjWWl1tx8eLoX7ypxO0PxrTo41gGBk16zv88SU7MSBQpA4aHHF4cb09cLCf4338cizm+6/4r6/HE5y5HicSTr9JpcyXLcqTNSd1Qk4snP78DD16zHAXjKc/dZ/vw511N+OJj+2asMo+GlP+hmagrzsKIy4e9jQNR1ywffuUU9p7rT+o1iPkPCWaakWIgcKTVhteOd0na3/JyK65fV57U681H4kWXZl36pqV4UzxeKx5vnozvQ5xtUMdM28bD7QviXx/dM6N9o9WoQZZeDaWEAqOqfFPE5+6NE9I+z5vrCsKZgYExD4YmVcs29Izhs4/sxl/3tsDlDWDY6YXXH8SJjmH89vUzU4qE3L4g/nm4A60yZsjmm7XQqTm8/8dv4NAMDj5nukbx9afewyd+tRPvnpEeORMXFpSSTTMlVj0senVcy7Ffv3oG66vzJK1/fOzyOhxoGohq85WqCtm5hkbF4faLqjDs9GHY6cOQw4sRly9s/J2udCyQfhP26RxtHcIvXz6F+3bUQaWIXjjCMAw+dWU9fvriCTT1zdxuEYsAL+DXr57G9uXFUQtUeoZdkgTovh11U26LoigpHbu01IKrVpehb9SN371+Fqc6hyPaS3pH3Pjt62fwx3ca4Avw4FgmnCbuHnLhgauXItekxW9eO42XjnTGfc1ohKpheyT1OnfYnHhqd7OkLABx4UGCmWY4lsVVq0vjjkxyegP4xcun8JXbojdKT0arUuBz167Al/+0P7xNp1LA7Q+mbGj0XCPXpMV9O5ZEbA/wAkacPllTNqSS7BBjOfz9YBvu2Vozo2ACoQKdj26vw38+cUD263gDPA42DUCj5MKONpPvi2dmYNGrI3yT9zcO4Hj7zNXMQEikvnvvRgChyDJe1DbxHkxeU93b0I+yXAM+ur1O9rq9kmNxSX0hfvTCCcnPaem3wx/kY743xIUJpWQzwDVryiDFp+DdM73Yc65P0j7XVOXg6tXnp7c/8sBl+H+3rUmZ6fp8QcmxyDNrUWKVZyEohUxHmBNIWYtORRr6W389jK/8+WDEmmZVvgn/cfOqmI47FXnGiON8OUakd/eWGlTkGnHD+gqoFBx8AR7/88xh2cdus4cqp+uK5U17uWF9OVoHHAjEqGafTlAQp9jfEQsHEswMkGfWYoXE/sCfv3QSLolVpx/fsQQ5Rg1YBtBrFNiypJBm7qUBKfNA04GU9zJVYzRZholqp3jp0iL85v5LsHSGyTrT24lEUcSR1lBV+ITrjkbJ4c7N1agvseDOzdV4+JNb8b6VJfAGePzq1dMYdcn3V911tg8dgw6sKLfiurVlCT9/UZEZzTEchGbihffa0T0kf7g1MT8hwcwUEnVsyOHDI29GlrVHQ69R4jPXLoNZp6b0UJrgBQF+ma0KySLl3IdL0QnSygrrjOvAerUS3//QRWFXnclMd6DiBRGb6wpwz9Ya/O7+S/Htezbgq7evxUe31+GHH74IWpUCYy4/FByLN453xZ0ZGw9fgA8vd9y0oTJhF6SH/nY0ZtvXTLx2vAsPv3oq4ecR8xtaw8wQFy8uwNmuUUnRyguHOrBtWXHEBPhobKzNx8cuT38f5ELF7Zud6BKQJpipyiiICFnKzbQ/lmFQlW/CsPN8Fen168ojjNwVHIt/v2lV+PZkkwGGYeDyBfBvf9yLD1+2GH8/2JaSY3/zZA/qSiy4cX0Fvv3+jfjPP+2XfJKTiJui1ahGRZ4JVXlGVOWb4raMERceJJgZ4sb1FbhiRQme2tMctwAIAH7y4gn89v5LJe2bPCvTx2ytXwIAIyEtwUo08Y/H0dYh/PSfJ/HR7XUz+sAWZ+txvH0Iy8uycaJjGDUFppjG4naPHwwY6DUK/HVvCwbGPGAYoGvIhW8lsW4ZjYdfPoXeETc+cGktvv3+jTjcMogn3pU3SzbHpEGp1YDSHD1KrQaU5RhQmW9K2FCeuPAgwcwgOrUCH7xsEfQaBZ7a3Rx1Uv0EHTYnfAGeJh3MMrMqmBK0MJUmFf883IGdp3vw0L2bUFMYWW39/ktqUV9qwWVLizDm9kdtRfEHefz9YBv6Rz14r3kQTm8AG2vz8Prx+L2lySACMKgV0KuVWF6WjfqSLLx1sge9UTxtdWoFirP1yDVpkGfWIsekQZ5JiwKLDqVWQ0ZsHon5CX0yMgzLMLj9omrY7N64BtrH24ciSv2JzDK7EWZ8pI6Jk4rTG8Tn/7AXP/rwxRFmBGadCpctLQpfj0bHoBO/e33qGny6xRIANtTm4Z5LzqeHOZbF9+7diB+9eAJGjRJrq3ORY1Sj0KJHoUVHxXGELEgwZ4n7r1wKjz+IV45Gd0SpzDOiMo/WSGab2VzDlBJiplgvAYQKaY602mS5RiUyEi5VfOzyOty6qTIi2s7L0uE779+Y8eMhLlyoSnYW+cCli8Kl99O5eWMlcqKU+ROZZa5HmOnyDZZjBNHYO4a/SZimk0rKcgy44+JqcKnqryGIGNCnbBbJNWlRmR9pUl2ZZ8T7qJBnTsCxDIqz9bAaQyPMMmFsP4GkNcx0hJgAXj3WhaYY5uzT4QUBX/i/PWhI4DnJUmjRxR25RRCphFKys8zly4uRZ9Kib9SNIYcPo+P+qCIkt24SaWTTonxsWpQfvi2KIoKCCF+Ahy/AwxvgEQgK8Ad5+IMCfEEePC8iyAsICqHLAC+AF8TwpSiKEAQRghgaOiyK49en9ThIWZ9Ml4D3jrjxmUd24zefumRGF6X+UTf2nOvHygorKvKMsJo0MUfPpRKTVomffmwzTFqqXCUyBwnmLHPThkrctOH8rL0gL8DpDWQ0kiGkwzAMlBwDJcfCoInefpFJEm3UTwReEDHm9qNkBpOq1gEHfvXqaQChAczRhgGkGpNWiR0rSnDbRVWzKpa8IKJ72IWWfjta+u1o7bejc8iFX3/yEqpsv4AhwZxjKDg2rZM3iAuLdFZ73rKpEouLZvZozTVpw9dTIZY6lQK5Zg20KgXUSg4aJQe1kkOhRYeqPBMq840oseqjrlcGeAGjrtAkm45BJ1oG7FhZbp2SHZCDKIoYtHvROeREp82J1n4HWvrtaBt0RDVHONY2hA21VNl+oUKCSRDzGJWCxea6gvFqWQYsExJRjg1FwUoFCwXHQsmxUCtYqMZFSK1gwbEsgoKAQDCUNg4EBbx8tBO9I25sXVKIT15RH/O1o/VhAqFZnQwDePw8irP1KLLoYDVqYNQqYdKpYNKqYNIqYdQqoVEpoFKwUCk4mLTKGU8ARl0+tPQ7cKh5EIN2L4advrBATh7zNhm9WilJMEVRhNMbxKDdg+5hFzptIXHsHApdT8RLeH9jPwnmBQwJJrGgaeodg1atgEGjhF6tgGKGquW5SpZeja/evjYl+xLGjdN7R9wzmq1P4PIG8PsZPI/v2lyNa9eWg2OZKQIoiCLcviDsbj/sngCGnT54A24Extd+B8Y86Btxw+MPrQ17A0H4AgKGnd6EDdrrirOwoSY36n1P721Gl82FQbsHg3YvBsY8KTPYP9A0GNNikJjfkGASCxZRFPHg73dPmfeoVXHQqhTQqRTQqhXQqjjoJlKEKg4apQKa8esqBQeVgoVayYUiOCUHBcdAwbFQsCwU42udHMtCObGdY8EyDBgG0y7PR4cMA4hiyN8VYmhbJtxnXjvWhVOdIzDrVFgcZ1xWY+8Yti0rwvXrysMnGwatEjoVhx+/cBLP7GvFqNuHnmE3ekdcGHX54fAEYs7WTBU6tQJfvHHljMVK/zzckbbipIExD9oGHKgkn9kLEhLMBc4jb5xFy4AdJVYDirP1KLXqUWI1wGpUX/Bnyd4AH/EF7vHz8Ph5DMM3S0cVSaFFh//79La0v84VK0tw5arSuI8TRRGLirLg9Abg8gbg9AUx5PCGb+883QPfLEx4UXIsvnjjSlw67kYUjV++fCpjlbzEhQcJ5gLnaNsQzvWM4mDT4JTtWhU3bkBtQIlVj9IcA0qtBhRl6y6YUWJOiXNHZ5vpI7x4QYTd7Q9FoAhFoxOIYqgnMsALCPKhVhZ/kEeQF+EPhlph/EEB3vG2GF+AR1AQoWAnIuCJCJmBSafCxYsLIo7n60++h32NA2n9nRMh36zFf4z3Y9aXRE8ldw+58ItXTuFIy2DU+1PF8rJsii4vYEgwFzCiKKJzBiszj59HQ+9YRCM6y4QGYpdYDbhsadG8npSSbdDgic9dDpcvCLcvAJc3CJcvCI8/CLcv9OPxB+H2B+ENr6vx8PqDcHmDaBt0ZOQ4p5sTBHgBd/3o9Yy89mOf2Y488/lqWFEUcaZ7NCOvLZWP71gyo1AOO7347t+Oom3AgTG3/EHVUrlxfUXaX4OYPUgwFyDnekax52wflpZmw52g9ZsgAn2jHvSNelAjw2t0LsGxDKxGDayRZktx6R1x48M/fyv1BxWF6RWgGiWHfLMW/WPp7Xu8anXpFLEEQmt0mRAeqXzmmmXYWl8Y9T5eEPHIG+dwrG0oI8eSY9Lg4rrk2liIuQ0J5gJkzOXHX3Y3A2hOaj9F2frUHNA8xJeiqkopjLh8cHkD0E8ySijJMaRVMA0aBT6xYwkcngDcvgDUSg5ZevWsmKtHw2pU477Ll2DbsujrlXa3H7957QxeOx59uEGqMetUeOCqpeRpe4FDgrkAGbSn5os2FYLZ3DeGX792BjlGDXJMGuSaNLAaNcgxhi6z9OqUj7BKBalqQ5BK97ALiyaZCJRa9TjUnL71OKc3iFv+99Xw7QeuWoob1lfg3CynYzfU5uHG9RVYXWmNKk6+AI/v/+MY9p7rR4CPXnikHo/Q9WoF3P4geobdMz42HhzL4Mb1FXj/JbVzwvmJSC8kmAuQVH3ZF2frkt5H15ArZsqMZRhkG9TINqqRrVcj26iBRR+6bdGrkaVXjd/WzNhInw4yGWECob/TFMHMid4ykQ4KLTpcvqIYh1tseGZfZqeRTCbXpMF/3bpmxvdZFEU8s68FO0/3Rty3Y0UxrltbjkKLDmadKlwB7g/yONkxgi//aT+UHIuCLC2KsvXwBwUcabXFPJ6qfBP+4+ZVKM+VkdMn5iUkmAuQkx3DSe9Dq+JgSYGF36DdG/N+QRRhc3hhc8R+3O0XVeG+HUuSPh6pZFowu4ddU26XztBjmGqUHIuffnQz9GolHn3rbMJr3qki36zFA1cvjXlSdKxtCH94uwFAyHM2z6xFoUWPWzdVYskMRUEqBYeVFVb88cFtyDFpw9kMURTx130teOSNsxCitI6yDIN/v3ElieUCgwRzAZKK/soiiz4l+0lVetisy6wRt0/GvMhk6BqaJpg5mVk/Ls81wKRTobF3DK39makKjsZVq0uxsTZ2QU1NoRm/+dQlyDNroVWd/2oTxVAbTv+YB32jbrh9QfCCiNpCM2oKTOBYBvlZU7MlDMPg9ouqUZ1vxneePQz7tMKrmzdWUPvIAoQEcwGypDgLu8/2JbWPohSkY4H4EaZUMm1YP9sRpkWvhl6tgMuX3ojvizeuAgD8bX+r7HW+RGEZwKRTIUunRkGWFgUWHVaUzzAyZRIGjXLKOuKhlkE89nYD2gYd8Pgj3y8GwMd21OH2i6oBhPpye4Zd6Blxo644CwVZOqypysHP7tuCbz51CM39dgChgqN7L1mUml+WmFeQYC5AdqwowWPvNCTlxlJkSU2EM5sRptsXxKjLB70MH9lMC6bdM7WVg2EYGDTKtAtmflaorSSZCF7JsTDrVDDrVDCNX56/roRZp4ZJp4RZq4LFoIZRq0qq0Ktn2IVfv3o6rrmCCOCZva1493QfekZcU9p3jFolvnzLaqytykVBlg4//MjFePTNs8gxaXBpfT3FH/4AABdhSURBVFFGrAqJuQe96wsMURTx+zfPJm1dVmxNjWDaUhZhJv6FfrTVhm88fSh8W61goVMroZvwkFUroFWFfjRhL1kOV6wsybhgRpuPmn5X1lALkkbJweULoDLPGPKN1Shh0CjCfx+NkoNawYHjGOjVypD46VQw69Qw61TQqriM2Sz6Ajy+9Ph+DEhsuRlxhaadTMfhCeArTxzAh7fV4Y6Lq6BRcrj/yqWpPlxinkGCucBoH3TitWPJ96YVWZJPyQZ4ASPO1Hi2yomAnL6p61K+oABfMPoX6GRWlFsz3lYSVTDTbGS+rjoXrxzthNsfhCCG3nOXLwib3YMu2/mJIl5/yF4PAG5YX44HrlqW1uOKxYGmAcliGQ9BBB558ywaekbxhRtWUlRJkGAuJERRxHMHUtMWkIoezCG7N2VRklnGGqYzygxFKWhUXMYjzEz74JflGBDgBTyxqymh5x1oHMADV8V/HC8IcHgCGHX5Yff44fYFkx72DABb6grwtdvX4q/7WnCqcyTp/QHArrN96LA58bU71s44AYVYGJBgLiBaBxx46Uhn0vtRKzlkG1LRUpKaSEA7nipNFKdX3vqfWsHBn+FpHJlOyS4ttcj6rDAMA1EU8dbJHoy4fBhz+THmjvyJZnz/xOcuh9WomXHfvgCPvlE3eobdGLB7YLN7YbN7YHN4YdGr8aWbV4NjGVxcV4BNi/Pxr4/uwdkUGS102Jx48Pe78cuPb0VhCrIrxPyEBHMBYdAowTJM0jMJiyy6lKxJpSp1JrcgxeWTH2HqM5yeiyaY6VRMuUU3KgULhmHw+vEuHGqJ3fg/naOtNly+YqqZ//eeO4rG3jEMO30xp8t8+ZbVU46ZZRh89trl+PTvdoGP1kgpg421eSjI0sZ/IHHBQsaHCwgFx8BiSL5fMVUesrPdUiJ3vJdayeH9l9TiSzetQn6WFmurcmTtJxEymZLVqRSyPWMnRr/J+YwcaY10fFJyLDpszrjv1Y9fOI4/7WycYgxflW/Cxy6vS/g4olGRa8Tnrl1+wc+IJWJDgrlAEEURX3/yEIYcyRfZFKdMMGc3wpSdklWGqj63Ly/GHx/cjv+5ZwOWlWXL2pdUorW8JJspiIZOrcCWJQU40S7PDUqlCB2nX4axw+HWwYhCpjUST0Y8fh5/fKcBP3r++JTtt26qwnfv3Rhuj5HL/VfVQ6OihNxChwRzgcAwDL73wU24c3N10mbmc820QHZKVmaEOX29lGEYfOf9G/DNu9bh2rVlKLXqkWfWoizHAGUCvZ2xUER5z9IhmNuWFeHVY12ys70TEaacNd4hhy9iPuvqyhwk8mnd29CPXWemesmurszBd+7ZmJQ5ehW5+hCgNcwFhUbJ4aPb67B9WTF+9tJJ2Z6yqTMtmH8pWZZBVBFUKThsrM2PsG/70QvH8fKRThi1yoi5lonARYswU7Q2N0FpjgF7zyXnAHU+wpRXFHWk1YaySf6sJp0KtUVmNPSMxXjWVH7y4gkYtEqsqjgfnRZb9fjt/ZfgG08dSrgQaElxFkzazFovEnMTijAXIBV5RvzvBzfhCzeskNXwP9dSsnJ+BwCyXHIm0rFSuXtLDf7v09vwh09vw2VLi8AAqCkwRQxmthrVuGxpEbbUFWB9TS4+dNkivH9rLQya0DltJiLMpSUWDDuTGw49IZhybfQOR1nHXFuVm9A+7J4AvvTYfvzhrXNTtmcbNHjw6mUJRaxV+SZ88+71Cb0+ceFCEeYChWUYvG9lKTYvLsCf3m3EcwfaJFUTqhUsso3Jt5R4A3xSEddksmSvYSb++hplYv8yBZNMvb98y2p84ool4FgGg3YvBsY8sBrV6Bl2Y9Oi/KiN8TdtrMCfdjaiZ8QdcZ+Qws4WtYJNyQnM+ZSsvD7V421DCPLClDXbtdW5+HOC/aAA8MSuJgiiiHsvXRTOCtQUmvHjj27GL146iYbe2FFrRa4R3713I0WXRBiKMBc4eo0Sn7iiHr/65CVYWx3/TL7Qoo/e4pAggylqKQHkmRbwggi3jAhTpUzuX2ZiKHZtoRmb6wpQV2zB9uXFM7rImLQq3H/lUnzpplUR9/EpjDA3LsrH4QTbQKIRjjBlpmTd/iBOd001HFhSnCW7jecvu5vxzN6WKdvqirPw0Ac2QauauXe3LMeAhz6wMeNTcIi5DQkmASD0BfE/d6/HN+5cF7MxOxVDo4HUrV8C8iJMOWIJAP2jHnzn2SMYc/sxFGdGZyqJVrCSSms8lklNW6cyyTVMIOQWNBkFx2KdhJO5mXhiVxP+frAN/KSQXKdW4Oo1ZVEfX5Ktx3fv3ZjxCTjE3IdSskQYhmGwaVE+1lTl4GDTILqGXOgZcYVGHg27YXN4U9iDmboIU84Xm9wKWQB4+1QP9pzrgz8o4KF7N2JVZfr7MKPx7zetgigCI04vfvnK6aT2larmfgWbXIQJAPsbByKGgW9alI93TvfO8IzY+AI8fvnyKfSPuvGJK+rD23csL8bbJ3swPMnP2KhV4rsf2BjTcYhYuJBgEhGoFBw21xVEbPcGeARTNBNxIsLkWAYbavJgc3gx5PBixOkDyzL4+h3rYPeEorgRlx8jTh+GnV6MuvwRri8mXeLtAnJNCyaYiKC++uR7WFeVg/W1eagryoIghqz6UnViEYutSwoBAA09ydu/pcp1aaJliU9igbXD5kTfqHvK+u+6mtykXapeOdqFu7fUwqgNfV6qC8z4yUc3476H3wl7A3/22uXINZGbDxEdEkxCMholB8jwbI3GRDoz16TB1+9cF94e5AUMO30w61RQx3gtf5CH3R3AmNsXLjRJhOmTSuTiC/DYfa4fu8/1h7flmbX4/b9cipZ+BxYXmdPuDiPXgGEy53rGUJajR4fNFf/BMWDHBTPZgPVA4wBuWF8Rvm3SqrC01IITMluhgNBJ0nMHWnHP1hpw45FwnlmLjbV52Hm6FztWFIdPQggiGiSYxKzwmWuW4SPbFsMxLdJTcGxEy0U0VAoOOSYOOSZ5qTNXCkRmJgbGPLj7R6/D6Q1i65JC3L2lGme7R7G6Mid8vC5vEC8e7kAgyEOjUqC13w61ksM9W2sTNvdOJr08GatRm7RgcmHBTE4xDzRNFUwA2FxXkJRgAsDjOxth1qmm7PuWjZWwGjX4wCW1Se2buPAhwSRmBYZhYNKpYJqlKsRkU7Lx9x8S5HfP9OLdcecZg0YBj5+HVsXB4+ejrhu+dbIHN2+sxN1baiTPX8wyqHHlqhJwLAsFx0DBslBwLJQcC6NWOelHBaMmdF3JseBFEaIYMkAQRBE/nGYrJwduPJpO1lThcIsNA2OeKSdPV64qxZ/ebUy6HWlvQz+uX1cejvyXlFiwpMSS1D6JhQEJJrEgSVVUlggTIhorhRrgBTy1pxn7Gvrx/9u78+AoyzsO4N/33ftgs1nIHXKQCwkhEEjAUIyiKBWleGBRBqntHx4d2ukltWo7Aw52nM7U2v5hD51a6k11amfQ8aQKCiKHB6CQICQcuYRkk+wme/aPJJSQTfZ5d9/d7G6+n5mMx7775hlN+O77HL+fw2rA0qo8XDd3+rj3rSpwoEqFWrZqNEiWL6xhRheY/kAQL3/UNKIZtdmgxS2LZuDvlxQkUGr/8U7cv3U3Hlu3iMXUSREGJk1Kaqz7xVJzZy+aO3uxoFT8OMXBE53w+gLw+YPw+QPwBQLw+AJwuj3ocXuHvgb/3uXxQZYkyJIESRqcSjUbIq+1Omx4bVCNKkSv72/BbfUlIzbhlGarU9N18cxshiUpxsCkSSnSXpjxpuTg/JZ/HRjR3koJs16LJ+9egj1H26I6j6lR6Qlz8B4BOF3eEYH52icno76v1agL+9ROFAoLF9CkFOs1TLUo6bBhMUb++dfl8eGTpo6I3z9MrU0/AFBXloWSS54oIy04cbEb5hfAxFZdFAEGJk1KiT4lO8znFw+eaNpXAUBTmxM31hZGdY8LganCE2ZaiPO1WQI7qMMpyU6L+h40OTEwaVKaiE0/kXB7xIN9SpSBua+pI+KWb8NkFZ8wQ03rVuRGH3aGKOsB0+TFnxyalJJlSlZJYFqiDEy9VoPjbT1R3WP4CTOSYhKX2tc0uhi8GhX8xiuIQTQeTuRTUvEHAlj/x/dgNmhhtxhgN+uRZtHDZho802kz6Yb+qkfa0D+H6mEZSS/MieD2iLfJinZK1hrFGuiw4ZKB6RZDxBuQhvW4R78/0vOSeq2MqgIHvm7vgUWF3cA0OTEwKaloZBn1Fdn4994TONnRK/QevVaG1aiDdejQvsWow7k4dhqJhltBsEczJWsz6XD4VPQ1abv6BguZ2616IMo9RBf3xBw2M8+OurLMUR1NwsmwmbBl7cLoBkSTHgOTks73rqrAnmNtaO0SKxju8Q3Wp724K0WycCmYkl1WnY/3j5zF2RDNpsNJMxvgVKGhd1ff4FNhugqtsfq9fhxqOYfK6SOLMtyzbBZ63V6YDVpYDFpYjDqY9BpoZRmyPHSuVJIgSRJkWYIsAdv3NyMYDPLsJUWFgUlJx2zQ4hffmYv7t+5WrS1Vonrz01NoqMwVquQzfZoVa5eU4Xevfaroe1iNOrR8I/a0Hs75oQ8l6VblgSlLEtKtejisRjisBjisoad186Za8Pu76hXde2lVHsOSosbApKQ0u8CBzWtqsenlfej3iq/zJZsBrx/b950ULn23sDwTGllS9EFi+lQLjpyOfjoWALpcQ4EZ4glTwmCol+emIc9hwdQp/w9Gh9UIm1l/YdOQ2tiyi9TAwKSkNb8kA4/duQgPPfexKtOJierjxg74/IGQa3qXspn0mFc8TVERgq4oN+eMuNcYU7ImvQb/2LB0wortE6mBx0ooqVXk2rFl7UJVCocnqt5+L3YfbQt/4ZCVCooPZNiMEa15jmV408+lU7Jujx9PvnkYQRXOZxJNFAYmJb2ynDRsWlMLvTZ1f5y3fXRcOGxqSzNRliN2wL/T2Y8cuxnpVv2IL7tFjwybEXkOM4oyrCjJtmFxRRbmFI4/Nez2+NHv9SPdMvpJ8p3PTwvvbCZKRKn7JwxNKlUFDjx4S81EDyNmlGzKkSUJRRlThK4NAjjb5cL5Xs+Ir64+Dzqc/Th9zoUTHb1oanVi7RXluLwiO+w9u3oHYA+xhqmVJWTbuZZIySt157Fo0llUnoXL8uyqbWBJJGU5dkW7PBtbu1X9/kadBoUZVryy53jYa7tcA6jItePFn14Dp8sDp9sLp9sDWZJgZNFzSmL86aWUsrK2CEdOH5zoYahuYVmm8LV9/V40d6o79Tng9aO5s1do9+35Xg8kSRqsxKTCeUyiRMEpWUopS2blqHJoPtE0VOYIX7vj0BnVz6cGATz43MdYs7gUuQ7zuNcmS51eIqX4hEkpRaeRcXlFFrbvb57ooahGK0vweAPC1+9Toa9lKOd6B9DhdGN9QwUeffUAACDXYcacgqkozpqC4kwbijOn8OgIpSwGJqWcutLMlApMXyAIp9uD7PTxn+yGzS5wYNdX4sdQlNj5ZSvubChHZpoJ6xrKcHVVfsyKDRAlGgYmpZy5xVOh08jw+sd/Ksu2m1CWY4fXH4DXH4Cr35uQG4ZWzC8QPiYCIKYba9440IJrq/Px9A+vhE6gkAJRKmFgUsox6bW4b3kl/vT6F2Ou5Rl0GmxeU4uCS45fvP3ZKTyx/QsMJFC5vZULioR3yLZ3u/HMjq9iOp7/fHISv7xJrFRfIvD6Awx3UgUDk1LS9TUFmD7Vgs3b9ocs4L3h27NHhSUAXDMnH6XZadi8bR9OfdMXj6GGpNPIqCvNwBWVucJTsQDw0odNF8rTGXQaaCRJUccTEZ4J+DAx4PWjs6cfnc5+dDrd6Ozph8NqxLLq/LDvXfv4O8i2m/GH79ezADtFRQpTPYR1rCipnevtxws7m7B9f/OFKdpl1fn4+crqcd/nGvDhvr9+oGrZOBEGrYx7l1fiilk5ETU6/rrNic6efhRMsyIjzQSfP4C9jR1474vT2H20/cJ/g1yHGRW5dlTk2tHt8uCFnY3Cv+y1pRl45PY6xWNT6rW9J/D+4bM42dETslawRpZw1excuAZ86Bvwwe3xIRgEAoEgAsEggkGgNMeGPcfa0eP24okfLEZFrj3m46aUEPKTFQOTJoW2Lhee/eAYjp7pxuN31Qut8+04dAaPvnIgDqMbZNRpsGlNLaqLpsbk/n39XjS2OkPuZN3b2I7fvnpQ6EiIVpbw8Or5WFSeFZNxDjt6pgsbntql2v1W1RXh3usqVbsfpTQGJpE/EBTe1RkIBvHjp3bh6Fl1q+aEYtZr8cgdtaOaJcdT63kXNm/bh8ZW56jXZAm4edEMuD0+vPv5adSWZsa8FOH53gFs/OduVBU4MLvAgXSrATuPtGLHoTPoueSJc06hA26PH8fG+X9l0muw9UdXY4pJ+ZM7TToMTCKlDp7oxMate2L6PSwGLbasrcPMvPSYfh8RA14/nn73S7RctH6rkSXc/q1SzMofHF/fgBdGnQYaObYbaYLBYMg1x/ZuN3717B60fNMHg06Du5ddhutrCiBJEl7ZfRx/fuvImPdcf2U57lhSFsthU2pgYBJF4oFn92D/8c6Y3NthNWDTmlpFx0ZosI3YSx82YUVNIfKmWka89vJHTfjb21+GfJ/NpMMzG5amdDs4UgUDkygSexvb8dDze1W9Z5pZj+8uLsGK+YUw6jSq3puAF3c14el3Q4fm8nnT8ZMb5sR5RJRkQgYmP2YRhbGgJAMF06yqFDS3mXRYXV+ClQsK2bkjhm6rn4GDJzpDzgy8caAFdaWZWDwzfKsyoovxNC9RGJIkYVVdUdT3WViWiWc2LMVt9SUMyxiTJAn3XDsr9GMCgCffPBzX8VBqYGASCVg+bzrmFkd+3KOuLBMP3VrDtbM4yhmn4EN7txuuAXULOlDqY2ASCdDIMh68uQZZdpPi9y4oycDDt9ZAr+VaZTy1drnH3YTx38Nn4jYWSg0MTCJBNrMev1m9AAat+K+NQafBxlVzGZYTwB2mJODzHzSGLdBPdDEGJpECJdk2rL2iXPj66+bmsz/kBCnPSUNptm3M19u63Xjr01NxHBElOwYmkULzFKxl3ji/MIYjofFIkoT7lldivMJOz+/kUyaJY2ASKTQjyyY8LTvNpnzNk9RTOd2BWy8vGfP19m43DrWci+OIKJkxMIkU0mpkVOSJdb3QK1jvpNhY11CG4szRrdyGne8diONoKJnxt5koAtVF08JeI0sStGxcPOH0Wg3uXzV3zCbSDEwSxd9mogisqCkY8w/gYToNmxUnihlZNqy/KvRmLZcn/g2xKTkxMIkikG41YFl1/rjXDPgC8AdYjjlR3LxwBuorRvfwVHJMiCY3/qQQRWhdQxmWzs4ds/waMNi0mRKDRpbw69Xz8cBN80b8+7JcdoohMQxMogg5rEZsvGke/nJvA66sDB2crV2uuI+LxiZJEhoqc3DZ0KYth9XA1mokjO29iFTS3NmL7fub8fZnp9DjHnyyvLY6Hz9bWT3BI6NL+QMBtHW5McWkxxSTbqKHQ4mH/TCJ4sHj8+PDr9pw7Gw30i0G3LKoGJLEDUBESYSBSUREJCBkYHINk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISAADk4iISIA2zOtSXEZBRESU4PiESUREJICBSUREJICBSUREJICBSUREJICBSUREJICBSUREJOB/8rJG4CHZJx0AAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import geoplot.crs as gcrs\n",
"gplt.cartogram(result, scale='Percent', projection=gcrs.AlbersEqualArea())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Save formats\n",
"\n",
"You can read data out of a geospatial file format using `GeoDataFrame.from_file`. You can write data to a geospatial file format using `GeoDataFrame.to_file`. By default, these methods will infer the file format and save to a `Shapefile`, respectively. To specify an explicit file format, pass the name of that format to the `driver` argument. For example:\n",
"\n",
"```python\n",
"nyc_boroughs.to_file('boroughs.geojson', driver='GeoJSON')\n",
"```\n",
"\n",
"The simplest and increasingly most common save format for geospatial data is [GeoJSON](https://geojson.org/). A geojson file may have a `.geojson` or `.json` extension, and stores data in a human-readable format:\n",
"\n",
"```\n",
"{\n",
" \"type\": \"Feature\",\n",
" \"geometry\": {\n",
" \"type\": \"Point\",\n",
" \"coordinates\": [125.6, 10.1]\n",
" },\n",
" \"properties\": {\n",
" \"name\": \"Dinagat Islands\"\n",
" }\n",
"}\n",
"```\n",
"\n",
"Historically speaking, the most common geospatial data format is the [Shapefile](https://en.wikipedia.org/wiki/Shapefile). Shapefiles are not actually really files, but instead groups of files in a folder or `zip` archive that together can encode very complex information about your data. Shapefiles are a binary file format, so they are not human-readable like GeoJSON files are, but can efficiently encode data too complex for easy storage in a GeoJSON.\n",
"\n",
"These are the two best-known file formats, but there are [many many others](https://en.wikipedia.org/wiki/GIS_file_formats). For a list of geospatial file formats supported by `geopandas` refer to the [fiona user manual](https://fiona.readthedocs.io/en/latest/manual.html)."
]
}
],
"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.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 1
}