{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spatial join\n",
    "\n",
    "[Spatial join](http://wiki.gis.com/wiki/index.php/Spatial_Join) is\n",
    "yet another classic GIS problem. Getting attributes from one layer and\n",
    "transferring them into another layer based on their spatial relationship\n",
    "is something you most likely need to do on a regular basis.\n",
    "\n",
    "In the previous section we learned how to perform **a Point in Polygon query**.\n",
    "We could now apply those techniques and create our own function to perform **a spatial join** between two layers based on their\n",
    "spatial relationship. We could, for example, join the attributes of a polygon layer into a point layer where each point would get the\n",
    "attributes of a polygon that ``contains`` the point.\n",
    "\n",
    "Luckily, [spatial join is already implemented in Geopandas](http://geopandas.org/mergingdata.html#spatial-joins), thus we do not need to create our own function for doing it. There are three possible types of\n",
    "join that can be applied in spatial join that are determined with ``op`` -parameter in the ``gpd.sjoin()`` -function:\n",
    "\n",
    "-  ``\"intersects\"``\n",
    "-  ``\"within\"``\n",
    "-  ``\"contains\"``\n",
    "\n",
    "Sounds familiar? Yep, all of those spatial relationships were discussed\n",
    "in the [Point in Polygon lesson](point-in-polygon.ipynb), thus you should know how they work. \n",
    "\n",
    "Furthermore, pay attention to the different options for the type of join via the `how` parameter; \"left\", \"right\" and \"inner\". You can read more about these options in the [geopandas sjoin documentation](http://geopandas.org/mergingdata.html#sjoin-arguments) and pandas guide for [merge, join and concatenate](https://pandas.pydata.org/pandas-docs/stable/user_guide/merging.html)\n",
    "\n",
    "Let's perform a spatial join between these two layers:\n",
    "- **Addresses:** the geocoded address-point (we created this Shapefile in the geocoding tutorial)\n",
    "- **Population grid:** 250m x 250m grid polygon layer that contains population information from the Helsinki Region.\n",
    "    - The population grid a dataset is produced by the **Helsinki Region Environmental\n",
    "Services Authority (HSY)** (see [this page](https://www.hsy.fi/fi/asiantuntijalle/avoindata/Sivut/AvoinData.aspx?dataID=7) to access data from different years).\n",
    "    - You can download the data from [from this link](https://www.hsy.fi/sites/AvoinData/AvoinData/SYT/Tietoyhteistyoyksikko/Shape%20(Esri)/V%C3%A4est%C3%B6tietoruudukko/Vaestotietoruudukko_2018_SHP.zip) in the  [Helsinki Region Infroshare\n",
    "(HRI) open data portal](https://hri.fi/en_gb/).\n",
    "\n",
    "\n",
    "\n",
    "- Here, we will access the data directly from the HSY wfs:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "from pyproj import CRS\n",
    "import requests\n",
    "import geojson\n",
    "\n",
    "# Specify the url for web feature service\n",
    "url = 'https://kartta.hsy.fi/geoserver/wfs'\n",
    "\n",
    "# Specify parameters (read data in json format). \n",
    "# Available feature types in this particular data source: http://geo.stat.fi/geoserver/vaestoruutu/wfs?service=wfs&version=2.0.0&request=describeFeatureType\n",
    "params = dict(service='WFS', \n",
    "              version='2.0.0', \n",
    "              request='GetFeature', \n",
    "              typeName='asuminen_ja_maankaytto:Vaestotietoruudukko_2018', \n",
    "              outputFormat='json')\n",
    "\n",
    "# Fetch data from WFS using requests\n",
    "r = requests.get(url, params=params)\n",
    "\n",
    "# Create GeoDataFrame from geojson\n",
    "pop = gpd.GeoDataFrame.from_features(geojson.loads(r.content))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Check the result: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geometry</th>\n",
       "      <th>index</th>\n",
       "      <th>asukkaita</th>\n",
       "      <th>asvaljyys</th>\n",
       "      <th>ika0_9</th>\n",
       "      <th>ika10_19</th>\n",
       "      <th>ika20_29</th>\n",
       "      <th>ika30_39</th>\n",
       "      <th>ika40_49</th>\n",
       "      <th>ika50_59</th>\n",
       "      <th>ika60_69</th>\n",
       "      <th>ika70_79</th>\n",
       "      <th>ika_yli80</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>MULTIPOLYGON Z (((25476499.999 6674248.999 0.0...</td>\n",
       "      <td>3342</td>\n",
       "      <td>108</td>\n",
       "      <td>45</td>\n",
       "      <td>11</td>\n",
       "      <td>23</td>\n",
       "      <td>6</td>\n",
       "      <td>7</td>\n",
       "      <td>26</td>\n",
       "      <td>17</td>\n",
       "      <td>8</td>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>MULTIPOLYGON Z (((25476749.997 6674498.998 0.0...</td>\n",
       "      <td>3503</td>\n",
       "      <td>273</td>\n",
       "      <td>35</td>\n",
       "      <td>35</td>\n",
       "      <td>24</td>\n",
       "      <td>52</td>\n",
       "      <td>62</td>\n",
       "      <td>40</td>\n",
       "      <td>26</td>\n",
       "      <td>25</td>\n",
       "      <td>9</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>MULTIPOLYGON Z (((25476999.994 6675749.004 0.0...</td>\n",
       "      <td>3660</td>\n",
       "      <td>239</td>\n",
       "      <td>34</td>\n",
       "      <td>46</td>\n",
       "      <td>24</td>\n",
       "      <td>24</td>\n",
       "      <td>45</td>\n",
       "      <td>33</td>\n",
       "      <td>30</td>\n",
       "      <td>25</td>\n",
       "      <td>10</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>MULTIPOLYGON Z (((25476999.994 6675499.004 0.0...</td>\n",
       "      <td>3661</td>\n",
       "      <td>202</td>\n",
       "      <td>30</td>\n",
       "      <td>52</td>\n",
       "      <td>37</td>\n",
       "      <td>13</td>\n",
       "      <td>36</td>\n",
       "      <td>43</td>\n",
       "      <td>11</td>\n",
       "      <td>4</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>MULTIPOLYGON Z (((25476999.994 6675249.005 0.0...</td>\n",
       "      <td>3662</td>\n",
       "      <td>261</td>\n",
       "      <td>30</td>\n",
       "      <td>64</td>\n",
       "      <td>32</td>\n",
       "      <td>36</td>\n",
       "      <td>64</td>\n",
       "      <td>34</td>\n",
       "      <td>20</td>\n",
       "      <td>6</td>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                            geometry  index  asukkaita  \\\n",
       "0  MULTIPOLYGON Z (((25476499.999 6674248.999 0.0...   3342        108   \n",
       "1  MULTIPOLYGON Z (((25476749.997 6674498.998 0.0...   3503        273   \n",
       "2  MULTIPOLYGON Z (((25476999.994 6675749.004 0.0...   3660        239   \n",
       "3  MULTIPOLYGON Z (((25476999.994 6675499.004 0.0...   3661        202   \n",
       "4  MULTIPOLYGON Z (((25476999.994 6675249.005 0.0...   3662        261   \n",
       "\n",
       "   asvaljyys  ika0_9  ika10_19  ika20_29  ika30_39  ika40_49  ika50_59  \\\n",
       "0         45      11        23         6         7        26        17   \n",
       "1         35      35        24        52        62        40        26   \n",
       "2         34      46        24        24        45        33        30   \n",
       "3         30      52        37        13        36        43        11   \n",
       "4         30      64        32        36        64        34        20   \n",
       "\n",
       "   ika60_69  ika70_79  ika_yli80  \n",
       "0         8         6          4  \n",
       "1        25         9          0  \n",
       "2        25        10          2  \n",
       "3         4         3          3  \n",
       "4         6         3          2  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pop.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okey so we have multiple columns in the dataset but the most important\n",
    "one here is the column `asukkaita` (\"population\" in Finnish) that\n",
    "tells the amount of inhabitants living under that polygon.\n",
    "\n",
    "-  Let's change the name of that column into `pop18` so that it is\n",
    "   more intuitive. As you might remember, we can easily rename (Geo)DataFrame column names using the ``rename()`` function where we pass a dictionary of new column names like this: ``columns={'oldname': 'newname'}``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['geometry', 'index', 'pop18', 'asvaljyys', 'ika0_9', 'ika10_19',\n",
       "       'ika20_29', 'ika30_39', 'ika40_49', 'ika50_59', 'ika60_69', 'ika70_79',\n",
       "       'ika_yli80'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Change the name of a column\n",
    "pop = pop.rename(columns={'asukkaita': 'pop18'})\n",
    "\n",
    "# See the column names and confirm that we now have a column called 'pop17'\n",
    "pop.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-  Let's also get rid of all unnecessary columns by selecting only\n",
    "   columns that we need i.e. ``pop18`` and ``geometry``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subset columns\n",
    "pop = pop[[\"pop18\", \"geometry\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pop18</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>108</td>\n",
       "      <td>MULTIPOLYGON Z (((25476499.999 6674248.999 0.0...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>273</td>\n",
       "      <td>MULTIPOLYGON Z (((25476749.997 6674498.998 0.0...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>239</td>\n",
       "      <td>MULTIPOLYGON Z (((25476999.994 6675749.004 0.0...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>202</td>\n",
       "      <td>MULTIPOLYGON Z (((25476999.994 6675499.004 0.0...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>261</td>\n",
       "      <td>MULTIPOLYGON Z (((25476999.994 6675249.005 0.0...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   pop18                                           geometry\n",
       "0    108  MULTIPOLYGON Z (((25476499.999 6674248.999 0.0...\n",
       "1    273  MULTIPOLYGON Z (((25476749.997 6674498.998 0.0...\n",
       "2    239  MULTIPOLYGON Z (((25476999.994 6675749.004 0.0...\n",
       "3    202  MULTIPOLYGON Z (((25476999.994 6675499.004 0.0...\n",
       "4    261  MULTIPOLYGON Z (((25476999.994 6675249.005 0.0..."
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pop.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have cleaned the data and have only those columns that we need\n",
    "for our analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Join the layers\n",
    "\n",
    "Now we are ready to perform the spatial join between the two layers that\n",
    "we have. The aim here is to get information about **how many people live\n",
    "in a polygon that contains an individual address-point** . Thus, we want\n",
    "to join attributes from the population layer we just modified into the\n",
    "addresses point layer ``addresses.shp`` that we created trough gecoding in the previous section.\n",
    "\n",
    "-  Read the addresses layer into memory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Addresses filpath\n",
    "addr_fp = r\"data/addresses.shp\"\n",
    "\n",
    "# Read data\n",
    "addresses = gpd.read_file(addr_fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>address</th>\n",
       "      <th>id</th>\n",
       "      <th>addr</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Ruoholahti, 14, Itämerenkatu, Ruoholahti, Läns...</td>\n",
       "      <td>1000</td>\n",
       "      <td>Itämerenkatu 14, 00101 Helsinki, Finland</td>\n",
       "      <td>POINT (24.91556 60.16320)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Kamppi, 1, Kampinkuja, Kamppi, Eteläinen suurp...</td>\n",
       "      <td>1001</td>\n",
       "      <td>Kampinkuja 1, 00100 Helsinki, Finland</td>\n",
       "      <td>POINT (24.93169 60.16902)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Bangkok9, 8, Kaivokatu, Keskusta, Kluuvi, Etel...</td>\n",
       "      <td>1002</td>\n",
       "      <td>Kaivokatu 8, 00101 Helsinki, Finland</td>\n",
       "      <td>POINT (24.94168 60.16996)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Hermannin rantatie, Kyläsaari, Hermanni, Helsi...</td>\n",
       "      <td>1003</td>\n",
       "      <td>Hermannin rantatie 1, 00580 Helsinki, Finland</td>\n",
       "      <td>POINT (24.97193 60.19700)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Hesburger, 9, Tyynenmerenkatu, Jätkäsaari, Län...</td>\n",
       "      <td>1005</td>\n",
       "      <td>Tyynenmerenkatu 9, 00220 Helsinki, Finland</td>\n",
       "      <td>POINT (24.92160 60.15665)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "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  Bangkok9, 8, Kaivokatu, Keskusta, Kluuvi, Etel...  1002   \n",
       "3  Hermannin rantatie, Kyläsaari, Hermanni, Helsi...  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.94168 60.16996)  \n",
       "3  Hermannin rantatie 1, 00580 Helsinki, Finland  POINT (24.97193 60.19700)  \n",
       "4     Tyynenmerenkatu 9, 00220 Helsinki, Finland  POINT (24.92160 60.15665)  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check the head of the file\n",
    "addresses.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to do a spatial join, the layers need to be in the same projection\n",
    "\n",
    "- Check the crs of input layers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'init': 'epsg:4326'}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "addresses.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the crs information is missing from the population grid, we can **define the coordinate reference system** as **ETRS GK-25 (EPSG:3879)** because we know what it is based on the [population grid metadata](https://hri.fi/data/dataset/vaestotietoruudukko). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define crs\n",
    "pop.crs = CRS.from_epsg(3879).to_wkt()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'PROJCRS[\"ETRS89 / GK25FIN\",BASEGEOGCRS[\"ETRS89\",DATUM[\"European Terrestrial Reference System 1989\",ELLIPSOID[\"GRS 1980\",6378137,298.257222101,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4258]],CONVERSION[\"Finland Gauss-Kruger zone 25\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",25,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",1,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",25500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"northing (N)\",north,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"easting (E)\",east,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"unknown\"],AREA[\"Finland - 24.5°E to 25.5°E onshore nominal\"],BBOX[59.94,24.5,68.9,25.5]],ID[\"EPSG\",3879]]'"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pop.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Are the layers in the same projection?\n",
    "addresses.crs == pop.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's re-project addresses to the projection of the population layer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "addresses = addresses.to_crs(pop.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-  Let's make sure that the coordinate reference system of the layers\n",
    "are identical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PROJCRS[\"ETRS89 / GK25FIN\",BASEGEOGCRS[\"ETRS89\",DATUM[\"European Terrestrial Reference System 1989\",ELLIPSOID[\"GRS 1980\",6378137,298.257222101,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4258]],CONVERSION[\"Finland Gauss-Kruger zone 25\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",25,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",1,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",25500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"northing (N)\",north,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"easting (E)\",east,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"unknown\"],AREA[\"Finland - 24.5°E to 25.5°E onshore nominal\"],BBOX[59.94,24.5,68.9,25.5]],ID[\"EPSG\",3879]]\n",
      "PROJCRS[\"ETRS89 / GK25FIN\",BASEGEOGCRS[\"ETRS89\",DATUM[\"European Terrestrial Reference System 1989\",ELLIPSOID[\"GRS 1980\",6378137,298.257222101,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4258]],CONVERSION[\"Finland Gauss-Kruger zone 25\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",25,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",1,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",25500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"northing (N)\",north,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"easting (E)\",east,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"unknown\"],AREA[\"Finland - 24.5°E to 25.5°E onshore nominal\"],BBOX[59.94,24.5,68.9,25.5]],ID[\"EPSG\",3879]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check the crs of address points\n",
    "print(addresses.crs)\n",
    "\n",
    "# Check the crs of population layer\n",
    "print(pop.crs)\n",
    "\n",
    "# Do they match now?\n",
    "addresses.crs == pop.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now they should be identical. Thus, we can be sure that when doing spatial\n",
    "queries between layers the locations match and we get the right results\n",
    "e.g. from the spatial join that we are conducting here.\n",
    "\n",
    "-  Let's now join the attributes from ``pop`` GeoDataFrame into\n",
    "   ``addresses`` GeoDataFrame by using ``gpd.sjoin()`` -function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a spatial join\n",
    "join = gpd.sjoin(addresses, pop, how=\"inner\", op=\"within\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>address</th>\n",
       "      <th>id</th>\n",
       "      <th>addr</th>\n",
       "      <th>geometry</th>\n",
       "      <th>index_right</th>\n",
       "      <th>pop18</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Ruoholahti, 14, Itämerenkatu, Ruoholahti, Läns...</td>\n",
       "      <td>1000</td>\n",
       "      <td>Itämerenkatu 14, 00101 Helsinki, Finland</td>\n",
       "      <td>POINT (25495311.608 6672258.695)</td>\n",
       "      <td>1514</td>\n",
       "      <td>515</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Kamppi, 1, Kampinkuja, Kamppi, Eteläinen suurp...</td>\n",
       "      <td>1001</td>\n",
       "      <td>Kampinkuja 1, 00100 Helsinki, Finland</td>\n",
       "      <td>POINT (25496207.840 6672906.173)</td>\n",
       "      <td>1600</td>\n",
       "      <td>182</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Hermannin rantatie, Kyläsaari, Hermanni, Helsi...</td>\n",
       "      <td>1003</td>\n",
       "      <td>Hermannin rantatie 1, 00580 Helsinki, Finland</td>\n",
       "      <td>POINT (25498443.209 6676021.310)</td>\n",
       "      <td>1904</td>\n",
       "      <td>275</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Hesburger, 9, Tyynenmerenkatu, Jätkäsaari, Län...</td>\n",
       "      <td>1005</td>\n",
       "      <td>Tyynenmerenkatu 9, 00220 Helsinki, Finland</td>\n",
       "      <td>POINT (25495645.995 6671528.068)</td>\n",
       "      <td>1550</td>\n",
       "      <td>1435</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Itäväylä, Vartioharju, Vartiokylä, Helsinki, H...</td>\n",
       "      <td>1007</td>\n",
       "      <td>Itäväylä 3, 00950 Helsinki, Finland</td>\n",
       "      <td>POINT (25506149.985 6678773.518)</td>\n",
       "      <td>3007</td>\n",
       "      <td>155</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "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",
       "3  Hermannin rantatie, Kyläsaari, Hermanni, Helsi...  1003   \n",
       "4  Hesburger, 9, Tyynenmerenkatu, Jätkäsaari, Län...  1005   \n",
       "6  Itäväylä, Vartioharju, Vartiokylä, Helsinki, H...  1007   \n",
       "\n",
       "                                            addr  \\\n",
       "0       Itämerenkatu 14, 00101 Helsinki, Finland   \n",
       "1          Kampinkuja 1, 00100 Helsinki, Finland   \n",
       "3  Hermannin rantatie 1, 00580 Helsinki, Finland   \n",
       "4     Tyynenmerenkatu 9, 00220 Helsinki, Finland   \n",
       "6            Itäväylä 3, 00950 Helsinki, Finland   \n",
       "\n",
       "                           geometry  index_right  pop18  \n",
       "0  POINT (25495311.608 6672258.695)         1514    515  \n",
       "1  POINT (25496207.840 6672906.173)         1600    182  \n",
       "3  POINT (25498443.209 6676021.310)         1904    275  \n",
       "4  POINT (25495645.995 6671528.068)         1550   1435  \n",
       "6  POINT (25506149.985 6678773.518)         3007    155  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "join.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Awesome! Now we have performed a successful spatial join where we got\n",
    "two new columns into our ``join`` GeoDataFrame, i.e. ``index_right``\n",
    "that tells the index of the matching polygon in the population grid and\n",
    "``pop18`` which is the population in the cell where the address-point is\n",
    "located.\n",
    "\n",
    "- Let's still check how many rows of data we have now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "28"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(join)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did we lose some data here? \n",
    "\n",
    "- Check how many addresses we had originally:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "34"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(addresses)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we plot the layers on top of each other, we can observe that some of the points are located outside the populated grid squares (increase figure size if you can't see this properly!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fa2f6d68d30>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmgAAAHfCAYAAADkyCQZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3df7RlZXng+e9TFDCBiVZRP4hdQAOjIGqQWNdf2NogK0h3ZwZnjY5UiMEfK8yYhGWnV6LiSppZSqOkTRPSq03GIAi2Fu0QR52YSJMW7HQg6C1BQ6wYbUJClShVqaI0kAjFfeaPs6+eunXr/Nzn7Hef8/2sVatO7V9n7/fsc+u57/vs543MRJIkSeVY0/QJSJIk6VAGaJIkSYUxQJMkSSqMAZokSVJhDNAkSZIKY4AmSZJUmLVNn8A0bdy4MU899dSmT0OSJIkdO3bszcxNq62bqwDt1FNPZXFxsenTkCRJIiL++kjrHOKUJEkqjAGaJElSYQzQJEmSCmOAJkmSVBgDNEmSpMIYoEmSJBXGAE2SJKkwBmiSJEmFMUCTJEkqjAGaJElSYQzQJEmSCmOAJkmSVBgDNEmSpMIYoEmSJBXGAE2SJKkwBmiSJEmFMUCTJEkqjAGaJElSYdY2fQKSJM2bU9/12UP+/dD7/0VDZ6JS2YMmSZJUGAM0SZKkwgwUoEXEuoi4LSL+IiJ2RsTLq+VXRMTXI+LPI+LXq2XHRMRNEfFnEfGViDiv6zh3VdvfX/3ZXC0/NiL+U0R8MyLujYhTu/a5LCK+Uf25rGv5adW236j2PaaWFpEkSWrYoD1o1wOfy8znAi8EdkbE+cDFwNmZ+XzgA9W2PweQmT8O/CTwGxHR/T6XZuY51Z9Hq2VvBfZn5rOB64BrASLiBOAq4KXAS4CrImJ9tc+1wHWZ+Rxgf3UMSZKk1usboEXEM4BXAR8GyMwnM/Mx4G3A+zPz+9Xy5WDrecB/6Vr2GLDQ520uBm6uXt8GXBARAbwGuCMz92XmfuAO4KJq3aurban2fW3/y5UkSSrfID1opwN7gJsi4r6IuCEijgfOAF5ZDTN+ISJeXG3/FeDiiFgbEacBW4GTu453UzW8+WtVoAWwBXgYIDMPAgeADd3LK7uqZRuAx6ptu5dLkiS13iBlNtYCLwKuyMx7I+J64F3V8vXAy4AXA5+IiNOBG4GzgEXgr4G7geVA6tLM3B0RPwr8HvBG4BYgOFyOsPwwEXE5cDnAKaec0vdiJUmatu6yG8OU3LBcx+wapAdtF7ArM++t/n0bnYBtF/DJ7PgisARszMyDmflLVY7ZxcA64BsAmbm7+vt7wMfp5JUtv8fJABGxFngmsK97eeUk4FvAXmBdtW338sNk5ocycyEzFzZt2jTA5UqSJDWrb4CWmd8GHo6IM6tFFwBfAz5FJw+MiDgDOAbYGxHHVUOgRMRPAgcz82vVkOfGavnRwE8BD1TH/Ayw/ITm64DPZ2YCtwMXRsT66uGAC4Hbq3V3VttS7fvpURtBkiSpJIPOJHAF8LGqlMWDwJuBx4EbI+IB4EngsszMqnTG7RGxBOymM4wJcGy1/GjgKOCPgN+t1n0Y+GhEfJNOz9klAJm5LyLeC3yp2u49mbmvev1O4NaIuBq4rzqGJElS60WnM2o+LCws5OLiYtOnIUnSyFbmnfViTlrZImJHZq5a6cKZBCRJkgpjgCZJklSYQXPQJEmaSW0uVbHyXFdey6jlO9Q8e9AkSZIKY4AmSZJUGAM0SZKkwpiDJkmaa03kZvUqlTHOVE+aHfagSZIkFcYATZIkqTAOcUqSNAF1DWP22rfNJULUmz1okiRJhTFAkyRJKowBmiRJUmHMQZMk6QimVcZimHy10qZvMg9uMuxBkyRJKowBmiRJUmEM0CRJkgoTmdn0OUzNwsJCLi4uNn0akqSW6pVvNUy+Wq+8snH0yv8a5tzNI5uOiNiRmQurrbMHTZIkqTAGaJIkSYWxzIYkSQOqa2hyUuU7plUWpNd7OjxaD3vQJEmSCmOAJkmSVBgDNEmSpMKYgyZJ0ohGLbMh9WMPmiRJUmEM0CRJkgpjgCZJklQYc9AkSWrYpHLZBq1JNk7tsl614ayJNjp70CRJkgpjgCZJklQYhzglSapBXdNA1Wkaw43DTPXk8Ofg7EGTJEkqjAGaJElSYQzQJEmSCmMOmiRJEzBMjtWo+Wp15b2NmkdW53F1KHvQJEmSCmOAJkmSVJjIzKbPYWoWFhZycXGx6dOQJKkW45TyqGu4sdc5OKTZW0TsyMyF1dbZgyZJklQYAzRJkqTCGKBJkiQVxjIbkqSZN408qWHywSaVmzWN6aYsnTEd9qBJkiQVxgBNkiSpMAZokiRJhTEHTZI0EyaRbzUp4+TEde/bb9tJ5IeZczYd9qBJkiQVxgBNkiSpME71JElqpV7lHqZVCmIaw6qTKgPiUGXznOpJkiSpRQzQJEmSCmOAJkmSVBhz0CS1Sr88mmlM6aNmjPPZ9zLOfTHM/TZoeYw6p4wa9Fh+N5phDpokSVKLGKBJkiQVxgBNkiSpME71JKk4w+TgWNtpfg1TB60Jkzo/7/H5YA+aJElSYQzQJEmSCuMQp6SZ5fDnbOlXVqOu8hR1lesY5vzqGv4sfdhXg7MHTZIkqTAGaJIkSYUxQJMkSSrMQFM9RcQ64AbgBUACb8nMeyLiCuAXgYPAZzPzHRFxDPB/AwvAEvD2zLyrOs5W4CPAjwB/UK3LiDgWuAXYCvwt8IbMfKja5zLgV6tTuTozb66WnwbcCpwAfBl4Y2Y+2es6nOpJaqc6p76RVio9N2uYfLlRc+v83jSjjqmergc+l5nPBV4I7IyI84GLgbMz8/nAB6ptfw4gM38c+EngNyJi+X1+G7gceE7156Jq+VuB/Zn5bOA64NrqxE8ArgJeCrwEuCoi1lf7XAtcl5nPAfZXx5AkSWq9vgFaRDwDeBXwYYDMfDIzHwPeBrw/M79fLX+02uV5wH/pWvYYsBARzwKekZn3ZKfb7hbgtdU+FwM3V69vAy6IiABeA9yRmfsycz9wB3BRte7V1bZU+y4fS5IkqdUGKbNxOrAHuCkiXgjsAN4OnAG8MiL+DfAPwC9n5peArwAXR8StwMl0hi1PpjPcuavruLuALdXrLcDDAJl5MCIOABu6l6/YZwPwWGYeXOVYh4iIy+n02nHKKacMcLmSmjbMUEyvUgbDDNtYkqN8JQx197v/2sR7vGyDDHGuBV4E/HZm/gTwOPCuavl64GXArwCfqHq2bqQTMC0CvwncTSdHLVY59nIC3JHWDbv88IWZH8rMhcxc2LRp06oXKEmSVJJBArRdwK7MvLf69210ArZdwCez44t0esg2ZubBzPylzDwnMy8G1gHfqLY/qeu4JwHf6nqPkwEiYi3wTGBf9/IV++wF1lXbrjyWJElSq/UN0DLz28DDEXFmtegC4GvAp+jkgRERZwDHAHsj4riIOL5a/pPAwcz8WmY+AnwvIl5W9bT9LPDp6pifAS6rXr8O+HyVp3Y7cGFErK8eDrgQuL1ad2e1LdW+y8eSJElqtUHLbJxDp8zGMcCDwJvpDHXeCJwDPEknB+3zEXEqncBqCdgNvDUz/7o6zgI/LLPxh8AVVZmN/wH4KPATdHrOLsnMB6t93gK8uzqVf5OZN1XLT+eHZTbuA35m+YGFI7HMhlSOYaakGWa6mknk1Zif1py6pl2a1HtOY8qoXpr4Pqg+vcpsDDQXZ2beT6eu2Uo/s8q2DwFnHr4pZOYinVpqK5f/A/D6I+xzI51AcOXyB+mU3pAkSZopziQgSZJUGAM0SZKkwgw0xClJdetVT6q0vJnSzkerq2tKpFG2G3ff7vMb5lz7vUfJ3yv1Zg+aJElSYQzQJEmSCjNQmY1ZYZkNafZYAmO2lTaVkveX6tSrzIY9aJIkSYUxQJMkSSqMAZokSVJhLLMhqdXMCdJqepVxGYc5j5oWe9AkSZIKY4AmSZJUGAM0SZKkwpiDJkmaCd35YOPknNV1HGkc9qBJkiQVxgBNkiSpMA5xSpKK1UQZC4c1VQJ70CRJkgpjgCZJklQYAzRJkqTCmIMmTYHTw0hqE39mNc8eNEmSpMIYoEmSJBXGIU5pAvo9pt+9vt/QwTDbzqNeQzEO00ij8bvSPHvQJEmSCmOAJkmSVBgDNEmSpMJEZjZ9DlOzsLCQi4uLTZ+GZlSvvLOV+RzmldVnmHZXffrlWdbV9nV9V6YxfVOv73m/bZtQ+vnNg4jYkZkLq62zB02SJKkwBmiSJEmFMUCTJEkqjHXQpJoMk7Nhfofablr38DDvM2pO1TD5ab3yzMbJc2s6L9WfSeWxB02SJKkwBmiSJEmFscyGJKlYo5avaZvShxibHoKdVZbZkCRJahEDNEmSpMIYoEmSJBXGHDRJmiMr87TalE/U5hyzfqbxOTi1U3nMQZMkSWoRAzRJkqTCGKBJkiQVxqmeJGnGDJOr1eactGnoV2utrqmopJXsQZMkSSqMAZokSVJhHOKUhlDXdCdOm6JJqmsKpBKGP0sbFhynDUr+rg/TziVfxyyxB02SJKkwBmiSJEmFMUCTJEkqjDloUk2GyStrYloX80bKm+qmhLyfYfLVmsidrKuMxTDXOS/flV7XOeh90O84Gp09aJIkSYUxQJMkSSqMAZokSVJhzEGTeuiVazFMHkYJORvzWHuthCmP6qrjNc759bpvmzifSel1ncN8l+vKT5vGPVVCu2sy7EGTJEkqjAGaJElSYRzilHroNXwwqSGwaZVTaJNeQzrjlKoofRqhSZxfXffBtO6nUT/fJs5vUsOfk7pPS7v/dSh70CRJkgpjgCZJklQYAzRJkqTCmIMmTUGb87+aMMw0M3VpIvdvnFIQs6q0fKs6p20bNJdy5boSSnv4M2z67EGTJEkqjAGaJElSYSIz+28UsQ64AXgBkMBbMvOeiLgC+EXgIPDZzHxHRBxdbfsiOkOot2Tm+6rj3AU8C/j76tAXZuajEXEscAuwFfhb4A2Z+VC1z2XAr1bbX52ZN1fLTwNuBU4Avgy8MTOf7HUdCwsLubi42Pd6pUmbRCXwfkM4bRqimNZw3qBDiOO03TjX0qbPbBh1fb7jDP1N4z1LMOh1t+26ZkVE7MjMhdXWDdqDdj3wucx8LvBCYGdEnA9cDJydmc8HPlBt+3rg2Mz8cToB1/8REad2HevSzDyn+vNoteytwP7MfDZwHXBtdeInAFcBLwVeAlwVEeurfa4FrsvM5wD7q2NIkiS1Xt8ALSKeAbwK+DBAZj6ZmY8BbwPen5nfr5YvB1sJHB8Ra4EfAZ4EvtvnbS4Gbq5e3wZcEBEBvAa4IzP3ZeZ+4A7gomrdq6ttqfZ97QDXK0nSISKX2Pj4fhhgREmalkF60E4H9gA3RcR9EXFDRBwPnAG8MiLujYgvRMSLq+1vAx4HHgH+BvhAZu7rOt5NEXF/RPxaFWgBbAEeBsjMg8ABYEP38squatkG4LFq2+7lkiQNLHKJ7dvfzT0ffBO3br+SyKWmT0kCBiuzsZZOPtkVmXlvRFwPvKtavh54GfBi4BMRcTqdocingX9Urf/jiPijzHyQzvDm7oj4UeD3gDfSyT2LlW9Kpydu2OWHiYjLgcsBTjnllAEuV5oN85JTMs70SJPKO+t1nDaXy5hE7mQ/o77PoO284YkDbN29k6OXnmbr7p1seOIAi//+Z454nGncM+MYZ+ozlWWQHrRdwK7MvLf69210ArZdwCez44vAErAR+Gk6+WpPVcOefwIsAGTm7urv7wEfpxPMLb/HyQDV0OgzgX3dyysnAd8C9gLrqm27lx8mMz+UmQuZubBp06YBLleSNC/2HreOHVvO4qk1R7Fjy1nsPW5d06ckAQMEaJn5beDhiDizWnQB8DXgU3TywIiIM4Bj6AROfwO8OjqOp9PD9hcRsTYiNlbbHw38FPBAdczPAJdVr18HfD47j5feDlwYEeurhwMuBG6v1t1ZbUu176dHbANJ0ryKYNu2a3j5z3+ES7a9D2K1ARpp+gadSeAK4GMRcQzwIPBmOnlmN0bEA3QeBLgsMzMi/gNwE53gK4CbMvOrVbB2exWcHQX8EfC71fE/DHw0Ir5Jp+fsEoDM3BcR7wW+VG33nq58tncCt0bE1cB91TEkSRpKxhr2Hr++/4bSFA1UB21WWAdNpSg9j2VQk5pWZpj2aXNbTqr9xjGJullNTM/V5ly/Ou8D66CVrY46aJIkSZoSAzRJkqTCDJqDJqlGg04x1DbTeMS/zvZqomxEtxKHldp0b/YaIm5zeZNpDX2XeP81remfCd3sQZMkSSqMAZokSVJhDNAkSZIKYw6a1ICm8xymVfZg1Pdpc/7QvGhziZCVet1vTdyLdZYwKeFzaZOS2sseNEmSpMIYoEmSJBXGAE2SJKkw5qBJDWhz3tmodbJKyu3QYMbJv2pT3mAJ+XTDTFs1ao6a38EyPutB2YMmSZJUGAM0SZKkwjjEKc2JYYZQehl16KqJcgBtGmZrSq82WvmZlDwcNI5hymzU9T1qQpuG9yalTddsD5okSVJhDNAkSZIKY4AmSZo7kUtsfHw/ZDZ9KtKqzEGTNHI5hTqnwZlEOYB+5zfq+7Q5l2dSeVKl5191i1xi+/Z3s3X3TnZsOYvTSDLG76+oa4qmfm056lRUbbpPZQ+aJGnObHjiAFt37+TopafZunsnG5440PQpSYcxQJMkzZW9x61jx5azeGrNUezYchZ7j1vX9ClJh3GIU2qRaQ2tjTpbQGnGGSoa9X1KHEaaxmdY53D3JKw8v9NINjxxoBOcRfTcdhomdQ+VeD9qMAZokqS5k7GGvcevb/o0pCNyiFOSJKkwBmiSJEmFcYhTKlwT+U2D5g8NM31TEzlJk8qLKj2vp7T8ryaUWH5iEt+rYdapXexBkyRJKowBmiRJUmEM0CRJkgpjDppUuEHzuNqee1Ly+Y+T0zWN6yr9/JrQ5um4VmrzuWt09qBJkiQVxgBNkiSpMA5xSi0yqaGOYYZORx1mHefcmyjXUVebTOoc5kVdn32vIc9ZGg7V7LAHTZIkqTAGaJIkSYUxQJMkSSpMZGbT5zA1CwsLubi42PRpqEWaKGsxTJ7NNPJoSivtUWdeVq/z7zXF1qTKWpSQczaJz3TS1xW5xIYnDrD3uHUQMdA+bcozm8YUZeNML6XRRcSOzFxYbZ0PCUiSWityie3b383W3TvZseUstm27hgwHh9R+3sWSpNba8MQBtu7eydFLT7N19042PHGg6VOSamGAJklqrb3HrWPHlrN4as1R7NhyVmeYU5oBDnFKPUyr3tWoep1PXTlp08o96ZXz1a3fZ1JXjtes5tz0a78mcg7HqnUWwbZt1wydg9attDpoJdTVa7oNZA+apFm2tATf+Q7M0cNQ8yhjDXuPXz9ScCaVygBN0mxaWoLzz4eTToLzzuv8W5JawiFOzb1Bh9YGWT8rhmmTaes3TLO8fuPj+7nnv/0JRy89zVP/7U94+ds/Dsevr+UcplUqYxJTXJU2TA/TKSMxif1mWWnDvvPIHjRJM8nkcUltZg+apNlUQ/K4JDXFAE3SzPpB8rgktYwBmtSlielOms6/mfSxjmTUEiHDlNkYNF+tVJM4v9JLxwxjnPt00OsuYUo3zSdz0CRJkgpjgCZJalzkEhsf319vzTrr4KnFHOLUzJvWUMKo71NXZfwShmd7Xcsw7z/OtvM4dDSNCv+T1D3h+dH/5BVw552wpn//Qc/7b7kO3t13w7nnDnzMw44zp2yD5tmDJklqVPeE59x9N+zZM/5B9+zpHOvgwfqOKU2RAZokqVHdNes491zYvHn8g27e3DnW2rX1HVOaIoc4JUnN6qpZt/hbl9ZTsy6iM6y5Z08nOLMOnlrGAE0zb5wcpWnk9pSeM9XE+Y0zzcy0pmGahkGve1LX1eu44+ROruZIE54Pcw6HWbMGTjxxtH3nnFM9Nc8hTknqMpGnCSVpSAZoklRZfprwng++iVu3X0nkUtOnJNXOX0LawQBNkirdTxNu3b2TDU8caPqUpFr5S0h7RM5RBL2wsJCLi4tNn4amYKy8lR7HKTm/aVo5IpNq27qM9Rllcuv2K9m6eyc7tpzFJdveN1ByeV1t38Q9VNdnNo3jjNM+TeRQNf0zYTUbH9/PPR98E0cvPc1Ta47i5T//kbHnqzU/bXQRsSMzF1Zb50MCkrSs62nCvcet88k/zZzlkibLv4TsPW5d06ekIzBAk6QuP3iaUJpF/hLSGgZomglNTOe0slt/UqUzeh13GkMLw1zLOENg0yg9Ute0UKUN6fQ7n173bRPpAMO8T9NlXsb5rCd1v43bJv4S0g4+JCBJUqEinfB9XhmgSZJUoOUnLjnpJDjvvM4E8JobBmiSJBVoueyLE77Pp4Fy0CJiHXAD8AIggbdk5j0RcQXwi8BB4LOZ+Y6IOLra9kXV8W/JzPdVx9kKfAT4EeAPgLdnZkbEscAtwFbgb4E3ZOZD1T6XAb9ancrVmXlztfw04FbgBODLwBsz88kx2kItU1fuzKDHrNOo+Vbj5Oo0kTc1zHWWVmKirjykYYyaazTMZ19avmav/SapDW17yBOXz3oul/y7L9Wa1F9CWRwd2aAPCVwPfC4zXxcRxwDHRcT5wMXA2Zn5/YjYXG37euDYzPzxiDgO+FpEbK8Crt8GLgf+lE6AdhHwh8Bbgf2Z+eyIuAS4FnhDRJwAXAUs0AkMd0TEZzJzf7XNdZl5a0T8TnWM3x6zPSRJKoNPXM61vkOcEfEM4FXAhwEy88nMfAx4G/D+zPx+tfzRapcEjo+ItXR6yp4EvhsRzwKekZn3ZKc67i3Aa6t9LgZurl7fBlwQEQG8BrgjM/dVQdkdwEXVuldX21Ltu3wsSZJmwpEmkdfsGyQH7XRgD3BTRNwXETdExPHAGcArI+LeiPhCRLy42v424HHgEeBvgA9k5j5gC7Cr67i7qmVUfz8MkJkHgQPAhu7lK/bZADxWbbvyWJIkSa02yBDnWjr5ZFdk5r0RcT3wrmr5euBlwIuBT0TE6cBLgKeBf1St/+OI+CNgtfB/+bnhI60bdvlhIuJyOsOqnHLKKattIg2srtyZYfLn2pz70aZz73Wuk6prV6JJnF8T9fvG+R41kbN3pPcf9hyG2beuPMsmPt95MEgP2i5gV2beW/37NjoB2y7gk9nxRWAJ2Aj8NJ18taeqYc8/oZNDtgs4qeu4JwHf6nqPkwGqodFnAvu6l6/YZy+wrtp25bEOkZkfysyFzFzYtGnTAJcrSYeLXGLj4/utRyVpKvoGaJn5beDhiDizWnQB8DXgU3TywIiIM4Bj6AROfwO8OjqOp9PD9heZ+QjwvYh4WZVD9rPAp6tjfga4rHr9OuDzVZ7a7cCFEbE+ItYDFwK3V+vurLal2nf5WJJUq+V6VPd88E3cuv3KTvHQOWSQKk3PoE9xXgF8rHqC80HgzXTyzG6MiAfoPAhwWVUy4z8ANwEP0BmKvCkzv1od5238sMzGH1Z/oPMAwkcj4pt0es4uAcjMfRHxXuBL1XbvqfLZAN4J3BoRVwP3VcfQGNrcTT2px8WnUQ5gloY0e+l3nU0PK/WyXI/q6KWn2bp7Z+epuhmdKudIbb8cpC5Psr1t2zVkDFZKc1pTbNV139RVwqeusimjDlv2O59h3rPk7+esGihAy8z76QxTrvQzq2z7d3RKbax2nEU6tdRWLv+HHvvcCNy4yvIH6eS7SdJEHVKPastZnZIHc2aeglSpBE6WLkn9WI/KIFWaMgM0SRrAD+pRzSuDVGmqDND0A3U+Zn2k49aZ5zbovqVPgTTOftPIBWlbjtyopQNKmH5opUlMCTaOSQWpg/6MmFaeVF05rXWV9qirzMY4uXWTmFpPvTlZuiRJUmHsQZMkqQRLS7BnT6eMiUPIcy9yjurZLCws5OLiYtOn0Rq9hopKGG7pNqmSF3UZZ2i3rmrfkzCpCvuTGlKc1uc7iXMofWh5pWE+z7o+l1F/Zk0qvWMY3WVMjv4nr4A774Q1a2p9j2HU1Sal/cwqTUTsyMzVqmQ4xClpdllYVW3RXcaEu+/u9KRprhmgSZpJVv9XmyyXMXlqzVFw7rmweXPTp6SGmYMmaSatVlhVKlZXGZPF37rUHDSZg6YjayKfYxjTyDsbJnes6by8ab1/07lQA597Jpx3Xme46Nxz4a67Dv9Pbzkpe/Pm2v5DbPo+GPYcVqqr5EWv82kiL2lS59D096EuTfxMMD/NHDRJ8yiik2i9a9eRg7Pzz4eTTuoEckvzMwRqbp6WeS+UywBN0uxaswZOPHH13rE9ezq9awcPzlVStrl5WtZ9L8zbLyltYIAmaT5t3twZ+ly7dq6Sss3N0zKfHC2bDwnMmWFyx3rlB9RVy6muHIRJ5TKMOlXQtPJS2pT/0i+PcdT8vpFzviI47WW/zIazq7klr/yDWs5nqHOYkJ7n/s7fP2TS85UJ6YO2ddvyh9qcczape3H5ydGtu3ey41nP5ZJ/96Uj5mJO4udbCbmJJTNAkzS35nIC9BWTnj/k04Lza8W94JOjZTFAk6Q5M5eBqVblvVAuAzRNVZ1d2KV1hw/T7d/0UFG/oYVpPDY/jbIMdV1XncNaTX/2dRnnHmrz0FYTpXZKG1bVdPiQgCRJUmEM0CSpBtaTklQnhzglaUzL9aSWn4zctu0aMvz9V9LoDNDm2Dj5OfOSnzBMrswwj6HXVaZkUup6pH7Ue2qc+2vUcx+nrMzCFf/xsNpiR0q8LiHncNB1K5VedmYSxmmfEnIemzArn33T/BVPksa0XE/qqTVHsWPLWZ2SBZI0BnvQJGlc1pOSVDMDNEmqgfWkJNXJAG3GTareUJvqGA2TD9Zv22Gmyhp03aTy1cY55jTyRqYxzdcw19nrnq7zfm9zHbRhzn1WrrOfuqbEG+ZnyyS+n3X+TDDvrB7moEmSJBXGAE2SJKkwDnHOuHGm3hl1OG8Y0xiC7XdddXXH1zU9Ul3nV8KwSOlGHXoepy2n0WZ1fda97qF5HcYq4Z6vq4xFCdeiI7MHTZIkqTAGaJIkSYUxQJMkqTDO7Spz0DSQuh4X73fcUfcdpkTCNM5nGHWdX7/3H56cgVQAABgcSURBVCcfcVSlTWs0jCbuoSaMUwJmHk2jxNDKuV1PIxuf23Ua0zdNI++5TexBkySpIBueOHDY3K6aPwZokiQVxLldBQ5xzrU6K9iP2uU9qXIAoxrnfEatnD7OuTdRsqG0oYbSh+HG+R6N2tbjDFuW3p6TUNyQXZ+5XUsv1RK55Ly0NTBAkySpMG2d23Vl/ty2bdc0nj/XVraaJEmqhflz9TFAkyRJtTB/rj4Occ6ZYfIKRs2pWmlSj2dP45HstpVMGFRduXYrNZGfNqkyEZPKqxw1J6203L9+2na+R1JnLu40cmWnUe6n33ueRpqDVgMDNEmSVJu25s+VxiFOSZKkwhigSZIkFcYhzpaYRD7YtDSdtzVMHlm/bSfRfk3Uguu37aDnVHqe0TDXOUwu2zi5ib3es9dxh7kXp7VtmzVxLW3OKexllu6LktiDJkmSVBgDNEmSpMI4xNlCTUy9U3oXdl1DOMMcd161qU2GKXtQ+j3eyzDTN03q82vTfdFP09fS9PsPa9Rh8bZd57TZgyZJMyByiY2P74fMpk9FUg3sQZOklnP+Q2n2+A2WpJZz/kNp9tiD1qBRc15KHLdv4tH8SUxpMq22LT3fqcR7bBT9ymGU/jkMavG3LoWv/i7cfTdHn3vuUPMfzkobNGWY9uv1vWpTCSTzeKfDHjRJarsIuPNO2LUL7rprqvMfmvsmTYYBmiTNgjVr4MQTpx6cbd/+bu754JvgvPNgaWlq7y3NOgM0SdJIunPfuPtu2LOn6VOSZoY5aFM0Tq5Hr3pN0xjjn9TUQNPIfxnmPafVtnXVDTJ/aDgltN8s5eTsPW4dO7ac1QnSzj0XNm9u+pSOqITc09JMqn6k6mGAJkkaTQTbtl3DhicOdB5UmOLw6qgil9jwxIHOgxQtOF/NLwM0SdLIMtaw9/j1rQh2rBenNjFAa4nSh7JG7f4eZsipri72Xu85r934vYYzZnmoY9Cpnvpd8zBTRg3zPoMOQY3zmYw67NtvvxLvk9Xqxe09fn1t6Sf9tGmYtbTzmUf+6iBJmgvLOXNPrTmKHVvOGqpenDRt9qBJkuZDV86cOWgqnQGaJGlu/CBnTipc5BxVf15YWMjFxcXG3r+uKUHGeY/S8wpmJR+sXz5OaeUdpnFv9tOm/Jxe6vxsJ/VzYBJK/4zqmtaotKnr1G4RsSMzF1ZbZw6aJElSYQzQJEmSCjNQDlpErANuAF4AJPCWzLwnIq4AfhE4CHw2M98REZcCv9K1+9nAizLz/oi4C3gW8PfVugsz89GIOBa4BdgK/C3whsx8qHrvy4Bfrba/OjNvrpafBtwKnAB8GXhjZj45QhsUwe7tjtLbodcQbNND2P3Op/RSLbOirrILK9cPc9xpDMtN67ta19BkXWa1rIbKM2gP2vXA5zLzucALgZ0RcT5wMXB2Zj4f+ABAZn4sM8/JzHOANwIPZeb9Xce6dHl9Zj5aLXsrsD8znw1cB1wLEBEnAFcBLwVeAlwVEcvZndcC12Xmc4D91TEkSTMscomNj++HOcqf1nzqG6BFxDOAVwEfBsjMJzPzMeBtwPsz8/vV8kdX2X0bsH2A87gYuLl6fRtwQUQE8Brgjszcl5n7gTuAi6p1r662pdr3tQO8jyRpCCUFRMszAdzzwTdx6/YriVxq+pSkiRmkB+10YA9wU0TcFxE3RMTxwBnAKyPi3oj4QkS8eJV938DhAdpNEXF/RPxaFWgBbAEeBsjMg8ABYEP38squatkG4LFq2+7lkqSalBYQrTYTgDSrBslBWwu8CLgiM++NiOuBd1XL1wMvA14MfCIiTs+qbkdEvBR4IjMf6DrWpZm5OyJ+FPg9OkOgtwCrVQvMEZYfJiIuBy4HOOWUU/pdq9RTr7yRYab76cVcMfNzSnGkqZEGMYn8q+WZAJbn0hxmJoBp5OW1vcyRyjJID9ouYFdm3lv9+zY6Adsu4JPZ8UVgCdjYtd8lrOg9y8zd1d/fAz5OJ69s+T1OBoiItcAzgX3dyysnAd8C9gLrqm27lx8mMz+UmQuZubBp06YBLleSBAVOjVTNBPDyn/8Il2x7nzMBaKb17UHLzG9HxMMRcWZmfh24APga8N/p5IHdFRFnAMfQCZyIiDXA6+nkrlEtWwusy8y9EXE08FPAH1WrPwNcBtwDvA74fGZmRNwOXNP1YMCFwJXVujurbW+t9v30OA0hSVqhwKmRnAlA82LQqZ6uAD4WEccADwJvBh4HboyIB4AngcuWhzfpBGa7MvPBrmMcC9xeBWdH0QnOfrda92HgoxHxTTo9Z5cAZOa+iHgv8KVqu/dk5r7q9TuBWyPiauC+6hiSpBoZEEnNcKqnCeo33Y9m1zi5KHXl7oyTYzNoPl2J93Sbpgur6zOq87ijvucw5zCpc6+rvmCv+7+J+oeaXU71JEmS1CIGaJIkSYUZNAdNA5rU8MqoUww1MR2LQwCHG+e6mmiTQe+pOofxB/3u9HvPNt1D/c51EtMcNfF9HGeaqmkYtHyONE32oEmSJBXGAE2SJKkwBmiSJEmFscxGoYbJK5tGWYESph+apVyQYT6zSXy+deUhTWtqm0nkYpVoGvlhTZSxmCV1/Sy0LQWW2ZAkrRC5xMbH98Mc/ZIutYlPcUrSnIlcYvv2d/9g0vFt264hw9/XpZL4jZSkObPhiQNs3b2To5eeZuvunWx44kDTpyRpBXvQClVa3aBxaptN6j2bzuGYVj2pSdUdm4a6avSNUxetZON8j0bNYwTYe9w6dmw56wc9aHuPW7fqsUrIPW3CvOQ8qmwGaJI0byLYtu0aNjxxoBOcRTR9RpJWMECTpDmUsYa9x69v+jQkHYEBmgbSxKP5kxpeGWZIrK5zmNa1dKtrWLpXew3z2U6qDdo8rDSt1IFh7vF5Gdacl+tUe/mQgCRJUmEM0CRJkgpjgCZJklQYc9BaYtRyBXXlKI2z7TAmkRdSeq7JtMpEDFM+ode2k/rsNTnTmA5uXtmemhR70CRJkgpjgCZJ82RpCb7zHefglAoXOUdf0oWFhVxcXGz6NEYy6DBdiY/t1zUE22vfaZxPnUooTzGoaQ1v63Cjfu8PsbQEe/bA5s2QyZ+e+sJD5uD8q2v/54Hfs02fbwnf+7raq82zZai3iNiRmQurrbMHTZJm1dISnH8+nHQSnHcefOc7zsEptYQBmiTNqj174O674eDBzt8R7NhyFk+tOeqwOTgllcWnOCWpxSKX4JFHOvNpnnjiofNqbtwIxx8PBw50/t60yTk4pZYwQJsx/XIVRs21mNS0QdMwTC5bE4YpedFr3zrzUtqchzRPJSUil9j+8StZ+vU/B+CLJz2fbT/9PjI6gyMbH9/PPd/9HkcD/N3fwd69h+WcrTRoiZV5zU3sdd2zdJ1qnkOcktRSG544wNbdO1lD54f5woq8sr3HrWPHSc+DtWvhFa/oPCggqRUM0CSppfYet44dW85iCVgCFrc879C8sgi2bbsGdu2Cu+5ySFNqEYc4JamtItj20+9j4+P7SYK9x68/LAjLWNPJTZPUKtZBa4lBcz/65YWMmiMxTJ7WpGqx1VXrrN++0zCNz2FSSs+zqSsnqLTconE++155ZZP6rpTQZqOqK1d3Uu/f5rbVoayDJkmS1CIGaJLUEpFLbHx8v9M0SXPAIc4ZMI3H24fpcp/lYbhpDNdO6j2Hef82l0yYxLlPashpmM82cont2999yDRNy+U06jTotczL9EN1DS2P856z2rZyiFOSWm+5pIbTNEnzwQBNkgrRawhzuaSG0zRJ88EyG5JUgL5DmFVNM6dpkuaDOWgzYBrlAEovY1HCdY9agqCu/L46p7Rqc85L0zmQo5aZ2fj4fu754Js4eulpnlpzFC//+Y906ppNWZs/+2mY1P1lu88nc9AkqXAOYUrq5hCnJJXAIUxJXQzQJKkQGWsaGdaUVB4DtBkwaj7TqO8xz+rKPxnmOINua+2k/iY1Ddk03nNa517aFFelqWv6MKkfc9AkqVDOHCDNL3vQJKlAK8tu8L5/Dmv8nVqaFwZoM6ZtQxLDlK0YtbTHOG3SxJBYr3Po1SYr1znkebhJDS/WVX6l28qZA9izB048caRjDWPQ+03Ds/00DH8dk6QCdZfd+OqPPQc2bWr6lCRNkQGaJJUogp++5Gq++mPP4exH/hLOPx+Wlpo+K0lTYoAmSYU64e+/x9nf/gZH5xLcfXdnmFPSXDAHbc7UlUM1TC5FE6U92pxHM6myDBrdOPdQr8+z32e9+FuXwld/txOcnXsubN488PuOen6DrpM0WQZoklSqCLjzzk7P2ebNzi4gzREDNEkq2Zo1U3l6U1JZzEGTJEkqTOQcVaheWFjIxcXFpk+jUdPIQRumVlebc1wmVRNtmJwl1Wecz9PPZX7M0s8wNS8idmTmwmrr7EGTJEkqjAGaJElSYRzinGNNTFvUz6jTOfU6Tp1DEpNoM4dIyjDMZ+tnJqkODnFKkiS1iAGaJElSYQzQJEmSCmOh2jlW15RCw7xPv3ywukpKTKsExrTfX2Ww1MJ4bD+pP3vQJEmSCmOAJkmSVBiHOFW7XsMVdQ5ljDrcOKnhFYc1281hNkklsQdNkiSpMAZokiRJhTFAkyRJKsxAUz1FxDrgBuAFQAJvycx7IuIK4BeBg8BnM/MdEXEp8Ctdu58NvCgz74+IrcBHgB8B/gB4e2ZmRBwL3AJsBf4WeENmPlS992XAr1bHujozb66WnwbcCpwAfBl4Y2Y+2es6nOppOKPmVI1TvmMaeUA+4i9JKkEdUz1dD3wuM58LvBDYGRHnAxcDZ2fm84EPAGTmxzLznMw8B3gj8FBm3l8d57eBy4HnVH8uqpa/Fdifmc8GrgOurU78BOAq4KXAS4CrImJ9tc+1wHWZ+Rxgf3UMSZKk1usboEXEM4BXAR8GyMwnM/Mx4G3A+zPz+9XyR1fZfRuwvTrOs4BnZOY92em2uwV4bbXdxcDN1evbgAsiIoDXAHdk5r7M3A/cAVxUrXt1tS3VvsvHkiRJarVBetBOB/YAN0XEfRFxQ0QcD5wBvDIi7o2IL0TEi1fZ9w1UARqwBdjVtW5XtWx53cMAmXkQOABs6F6+Yp8NwGPVtiuPJUmS1GqD1EFbC7wIuCIz742I64F3VcvXAy8DXgx8IiJOr3rHiIiXAk9k5gPVcWKVY2efdcMuP0xEXE5nWJVTTjlltU3UEnXljplzJkkq3SA9aLuAXZl5b/Xv2+gEbLuAT2bHF4ElYGPXfpfww96z5eOc1PXvk4Bvda07GSAi1gLPBPZ1L1+xz15gXbXtymMdIjM/lJkLmbmwadOmAS5XkiSpWX0DtMz8NvBwRJxZLboA+BrwKTp5YETEGcAxdAInImIN8Ho6T1kuH+cR4HsR8bIqh+xngU9Xqz8DXFa9fh3w+aon7nbgwohYXz0ccCFwe7Xuzmpbqn2XjyVJktRqg071dAXwsYg4BngQeDPwOHBjRDwAPAlctjy8Seehgl2Z+eCK47yNH5bZ+MPqD3QeQPhoRHyTTs/ZJQCZuS8i3gt8qdruPZm5r3r9TuDWiLgauK86hqak1zDhpKY8cmhSkjQvBgrQqjIZq9Xp+JkjbH8Xndy0lcsX6dRSW7n8H+j0uK12rBuBG1dZ/iCd0huSJEkzxZkEJEmSCmOAJkmSVJiBpnqaFU71NDmTyDsz50ySNMvqmOpJkiRJU2KAJkmSVBgDNEmSpMIMWgdN6mnUumjmmUmSdDh70CRJkgpjgCZJklQYhzg1ccMMY3YPhzr8KUmaV/agSZIkFcYATZIkqTAGaJIkSYUxB01FMe9MkiR70CRJkopjgCZJklQYhzjVqF6zDKzk8KckaV7YgyZJklQYAzRJkqTCGKBJkiQVxhw0Ncq8MkmSDmcPmiRJUmEM0CRJkgpjgCZJklQYAzRJkqTCGKBJkiQVxgBNkiSpMAZokiRJhTFAkyRJKowBmiRJUmEM0CRJkgpjgCZJklQYAzRJkqTCGKBJkiQVxgBNkiSpMAZokiRJhTFAkyRJKowBmiRJUmEM0CRJkgpjgCZJklSYyMymz2FqImIP8NfVPzcCexs8nTawjXqzffqzjXqzffqzjXqzfforuY3+cWZuWm3FXAVo3SJiMTMXmj6PktlGvdk+/dlGvdk+/dlGvdk+/bW1jRzilCRJKowBmiRJUmHmOUD7UNMn0AK2UW+2T3+2UW+2T3+2UW+2T3+tbKO5zUGTJEkq1Tz3oEmSJBWplQFaRJwcEXdGxM6I+POIePsq25wXEQci4v7qz79esf6oiLgvIn6/a9l/6tr+oYi4v1p+akT8fde635n8VY5ngm10TkT8abX9YkS8pGvdlRHxzYj4ekS8ZrJXOJ5pt4/30CFt9MKIuCci/iwi/r+IeEbXOu+hI7TPPN5D1c/hP1v+PnUtPyEi7oiIb1R/r+9a15p7CKbfRm27jybYPq+vjrcUEQsrjlfGPZSZrfsDPAt4UfX6R4G/BJ63YpvzgN/vcYx/BXz8SNsAvwH86+r1qcADTV93CW0E/Gfgn1Wv/zlwV/X6ecBXgGOB04D/DhzVdDsU1D7eQz9c9iXgn1av3wK813tooPaZu3sIeAjYuMryXwfeVb1+F3BtG++hhtqoVffRBNvnLOBM4C5goWt5MfdQK3vQMvORzPxy9fp7wE5gy6D7R8RJwL8AbjjC+gD+d2D7+GfbjAm2UQLLPR7PBL5Vvb4YuDUzv5+ZfwV8E3gJhWqgfVpngm10JvBfq9d3AP9b9dp7qONI7dM647ZRDxcDN1evbwZe27W8NfcQNNJGrTKp9snMnZn59VVWFXMPtTJA6xYRpwI/Ady7yuqXR8RXIuIPI+L5Xct/E3gHsHSEw74S+E5mfqNr2WnVUMQXIuKVNZz61NTcRv8S+LcR8TDwAeDKavkW4OGu7XZRzw+ZiZtS+4D30LIHgP+lev164OTqtfdQx5HaB+bvHkrgP0fEjoi4vGv5iZn5CHT+Awc2V8tbew/B1NoIWnof1dw+R1LMPdTqAC0i/kfg94B/mZnfXbH6y3SmUHgh8O+BT1X7/BTwaGbu6HHobRzae/YIcEpm/gTVkER05c2UbAJt9DbglzLzZOCXgA8vv9Uq2xb/iPAU28d76IfeAvxCROygM2Tx5PJbrbLtPN5DR2qfubqHKq/IzBcB/4xOm7yq31utsqz4ewim2katvI/m8h5qYly1jj/A0cDtwL8acPuH6MzH9T46EfFDwLeBJ4D/2LXdWuA7wEk9jnUXXWPWpf6ZRBsBB/hheZYAvlu9vhK4sutYtwMvb7oNSmkf76FDv2dd258BfNF7qH/7zNs9tMry/wv45er114FnVa+fBXy9rffQtNuojffRJNrnSNdf0j3UeMOP+GEFcAvwmz22+bGu/yhfAvzN8r+7tjmPFYmFwEXAF1Ys20SVJAicDuwGTmi6HZpoIzrj/+dVry8AdlSvn8+hiZUPUnBybgPt4z30w39vrv5eUx3/Ld5DA7XPXN1DwPHAj1bLjwfuBi6q/v1vOTQB/tfbeA811Eatuo8m1T5d+97FoQFaMffQWtrpFcAbgT+LqhQG8G7gFIDM/B3gdcDbIuIg8PfAJVm1fh+XcPjDAa8C3lMd62ng/8zMfeNfxkRNqo1+Drg+ItYC/wBcXh3vzyPiE8DXgIPAL2Tm0zVfU52m2j54D3XbFhG/UL3+JHBTdTzvoY5V24c5u4ci4kTg/+08s8Va4OOZ+bnqGO8HPhERb6Xzn/Hrq+O17R6CKbcR7buPJtI+EfG/0hkO3QR8NiLuz8zXlHQPOZOAJElSYVr9kIAkSdIsMkCTJEkqjAGaJElSYQzQJEmSCmOAJkmSNKSIuDEiHo2IBwbY9rr44WTufxkRj/Xdx6c4JUmShlPNSvB3wC2Z+YIh9rsC+InMfEuv7exBkyRJGlJm/lfgkBpyEfE/RcTnqrk//zginrvKriunk1xVWwvVSpIkleZDdIr/fiMiXgp8EHj18sqI+Md0Zij4fL8DGaBJkiSNqZrQ/Vzg/6lmL4DOlFHdLgFuG2R2AgM0SZKk8a0BHsvMc3pscwnwCz3WH3IwSZIkjSEzvwv8VUS8HiA6Xri8PiLOBNYD9wxyPAM0SZKkIUXEdjrB1pkRsauamP5S4K0R8RXgz4GLu3bZBtyaA5bPsMyGJElSYexBkyRJKowBmiRJUmEM0CRJkgpjgCZJklQYAzRJkqTCGKBJkiQVxgBNkiSpMAZokiRJhfn/AZCbY49xhqslAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Create a figure with one subplot\n",
    "fig, ax = plt.subplots(figsize=(15,8))\n",
    "\n",
    "# Plot population grid\n",
    "pop.plot(ax=ax)\n",
    "\n",
    "# Plot points\n",
    "addresses.plot(ax=ax, color='red', markersize=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also visualize the joined output:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-  Plot the points and use the ``pop18`` column to indicate the color.\n",
    "   ``cmap`` -parameter tells to use a sequential colormap for the\n",
    "   values, ``markersize`` adjusts the size of a point, ``scheme`` parameter can be used to adjust the classification method based on [pysal](http://pysal.readthedocs.io/en/latest/library/esda/mapclassify.html), and ``legend`` tells that we want to have a legend:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAegAAAGoCAYAAACNCHZgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeXzU1b3/8ddnshH2XcFAgYLIvkWQ1h0BRVu0WoVrW+py23pbe1t769bfvbYudamtWmu9bV1qvVZU6tpalRb3BQWkCKiIiBDZZQkQQpb5/P74nsQhZJmEhJkk7+fjMY/MnO855/v5zkzyyTnf78wxd0dERETSSyzVAYiIiMj+lKBFRETSkBK0iIhIGlKCFhERSUNK0CIiImlICVpERCQNKUFLi2ORe81sm5m9Wc32c83suST7+qaZvdLAOI43s4Jatv+vmf13Q/puSlWP2cx2mdmAJNpdaWZ3NW10++3zp2b2fwdznwc7BjNbbWYnNVX/ScbQN7wPMlIZR2ujBN0KmdkLIXnlpDqW6hxIUgyOBiYDee4+vupGd3/A3accQP+Nwt2/4+7XQN3JvL4a4Tms5O7t3X1VEvV+7u4XNsY+01Vjv07V9P9HM7u2qfpvKHdfE94H5XXVNbN+ZuZmlnkwYmvJlKBbGTPrBxwDOPDllAbTdD4HrHb33akORESkoZSgW59vAG8AfwRmJW4I/73/1sz+HqazXjWzQ83s1jDifs/MxiTUHxJG49vNbJmZfTlh2wtmdmHC46rTpm5m3zGzD0Lfd4Sp6SHA/wITQwzbqzsIM+ttZk+a2VYzW2lm/x7KLwDuSmj/s2raJhVLlTY3h20fmdkpCeXnmdm7ZrbTzFaZ2ber2d+VZrYlTFWeW+X5vtbM2gF/B3qHmHeF4xtvZq+H53e9mf3GzLIb+hya2TQzWx5i/cTM/qu657aa+N3MBprZUWa2IXGa08zOMLMl4X7lVG/CKGqWma0Jx/+ThHa5ZnZfiPtdM7u0tpGpmQ0zs7nh9d5oZlfWUO/L4b24PbwHhyRsuywc904ze9/MJoXymJldbmYfmtmnZvawmXWtpu9qX6ewOdvM/hT6XmZm+QntepvZX8xsc3j/fL+G2L8FnAtcGvp+KmHzaDNbYmY7zOwhM2uT0O40M1scjvk1MxtZy/PoZvb98F7dYma/MLNYwvPw/8zsYzPbFI6nU9i2z6g4PLfXWPQ3YqeZPWdm3cNuXgo/t4fjmFhTPFIHd9etFd2AlcB/AOOAUuCQhG1/BLaEbW2AecBHREk9A7gWeD7UzQp9XQlkAycCO4HBYfsLwIUJfX8TeCXhsQN/BToDfYHNwMnV1a3hOF4EfhviHB3aT0qmfQNiKQX+PTwHFwHrAAvbTwU+DxhwHFAEjA3bjgfKgF8BOWH77oTn6I/AtQl1C6rEOQ44CsgE+gHvAj9o6HMIrAeOCfe7VMSZ5PMzMNz/EJicsO0R4PJw/6fA/4X7/UK7PwC5wChgLzAkbL8hvIZdgDxgSdXjT9hHhxD7j8Lr3QGYUM0+Dw/P72Si9+elRO/RbGAwsBbonRDf58P9HxD905oXXqffAQ/WEEt1r9NPgWJgWniPXA+8EbbFgIXA/4Q4BgCrgKk19F/5nkgoWw28CfQGuob3wXfCtrHAJmBC2PesUD+nhv4deD700xdYQfg9Bc4Pz9cAoD3wKHB/ldczM+H3+8PwnOeGxzdUV1e3ht80gm5FzOxoounfh919IdEv2L9VqfaYuy9092LgMaDY3f/k0bmnh4CKEfRRRL/EN7h7ibvPI0oWM+sR0g3uvt3d1xD90Rid5HH0ITrPfJm7F7v7YqJR89frse/6xPKxu/8hPAf3Ab2AQwDc/W/u/qFHXgSeIzqFkOi/3X1v2P434OxkAgqvwxvuXubuq4kSx3H1iLuqUmComXV0923uviiZOKp4kPAam1kHoqT0YC31f+bue9z9X8C/iBI1RM/Bz0McBcCva+njNGCDu/8yvN473X1+NfXOAf7m7nPdvRS4mSh5fAEoJ0q+Q80sy91Xu/uHod23gZ+4e4G77yVKuGdZ/c6hvuLuT4f3yP0Jx3kk0MPdrw6/J6uI/mmZUY++AX7t7uvcfSvwFJ+9zv8O/M7d57t7ubvfR/SP0FG19HWju28N75lb+ex39lzgV+6+yt13AVcAM2p5Hu519xXuvgd4mCR/fyV5StCtyyzgOXffEh7/mSrT3MDGhPt7qnncPtzvDax193jC9o+Bw+oRz4aE+0UJfdelN7DV3XcewL7rE0vlNncvCnfbA5jZKWb2Rph63U6UsLontN3m+54L/zjEXyczO9zM/hqmlQuBn1fpu664qzozxPexmb3YwKnHPwNfsegCw68Ai9z941rq1xRfb6IRbYXE+1X1Ifpnsi69iZ5fAMJ7cy1wmLuvJBop/xTYZGazE6anPwc8FqaItxONUMsJ/4QlqepxtgmJ7XNEU+LbE/q/sp59V9d/xfP4OeBHVfrvQ+3vscTnOvH9uM/zF+5n1hJrQ39/JUlK0K2EmeUSjVqOC3/wNwA/BEaZ2ajaW1drHdCn4vxV0Bf4JNzfDbRN2HZoPfqua4m1dUDXMIKrbt8HRUhSfyEaqR3i7p2Bp4mmuyt0CecuK/Qlir+q6o75TuA9YJC7dyT6w27V1KvOfv25+1vuPh3oCTxONOqpF3dfTvSH+xSi2Zc/17ePYD3RlHKFPrXUXUt0GqEu64gSFhB93C70+wmAu//Z3StmkRy4MaH/U9y9c8KtjbtX936q7/J/a4GPqvTdwd2n1VC/If1fV6X/tu5e26xG4nOd+H7c5/kL28rY95/0ZGiJxEaiBN16nE40KhhKNBU1GhgCvEx0jrm+5hMl4UvNLMvMjge+BMwO2xcTjbTamtlA4IJ69L0RyLOEC6ISufta4DXgejNrEy6KuQB4oAHHcSCyiaZNNwNlFl08Vt3Ht35mZtlmdgzRdO0j1dTZCHSruCgn6AAUArvM7Aii89/J2uc5DPs/18w6henfQqL3Q0P8Gfg+cCzVH0syHgauMLMuZnYY8L1a6v4VONTMfmBmOWbWwcwm1NDnqWY2ycyyiM5Z7wVeM7PBZnZi+KeqmGg2qOL4/xe4zsw+B2BmPcxseg2xVPc61eZNoNCiC9RyzSzDzIab2ZG19F/nZ84T/AH4jplNsEg7Mzu1yj+vVf04PO99gP8kOnUF0amKH5pZfzNrTzRj85C7l9UjHoh+H+L1PA6phhJ06zGL6JzRGnffUHEDfgOcW8/zbbh7CdHHtE4hurDst8A33P29UOUWoIToD8591C95zgOWARvMbEsNdWYSXYyyjuhc+VXuPrc+x3CgwhT794kSwzaiEeWTVaptCNvWET0H30l4jhL7eo/oD+SqMFXZG/iv0OdOoj/ED1VtV4vqnsOvA6vDdPl3gK/Vo79EDxJdLDUv4XRJfV0NFBBdhPgPYA5RMt1PeJ4nE/0DuAH4ADihmnrvEx3T7UTvyS8BXwrv1RyiC9O2hD56Es1IANxG9Lo9Z2Y7iS4Yq+4fgJpepxqFc9JfIvqH+KOw/7uAmhL83UTnybeb2eO19R36X0B0Hvo3RO+zlUQX+dXmCaIL1xYTXRNxdyi/h+j8+Ush1mLg4rpiqCamIuA64NVwHLWdD5daVFyJKiKSMmZ2ETDD3ateBCeNyMyc6JTJylTHInXTCFpEDjoz62VmXwyfvR1MNB39WKrjEkkn+io2EUmFbKKPjfUHthNdu/DblEYkkmY0xS0iIpKGNMUtIiKShlrVFHf37t29X79+qQ5DRERagYULF25x9x4Nbd+qEnS/fv1YsGBBqsMQEZFWwMxq+5a9OmmKW0REJA0pQYuIiKQhJWgREZE01KrOQYuIHGylpaUUFBRQXFyc6lCkibRp04a8vDyysrIatV8laBGRJlRQUECHDh3o168f0QJb0pK4O59++ikFBQX079+/UfvWFLeISBMqLi6mW7duSs4tlJnRrVu3JpkhUYIWEWliSs4tW1O9vpriFhFJJ+5QXgoeB4tBRhYowbdKGkGLiKSL0mLYUwgle6B0b/RzT2FUfgDOP/98evbsyfDhw/cpf+SRRxg2bBixWGyfL3F64IEHGD16dOUtFouxePHi/frdunUrkydPZtCgQUyePJlt27ZVbrv++usZOHAggwcP5tlnn61XvGvXruWEE05gyJAhDBs2jNtuu63R93nzzTdjZmzZ8tmS5gcSc5Nw91ZzGzdunIuIHEzLly9PrmLJHvfd22u+lexpcAwvvviiL1y40IcNG7ZfbO+9954fd9xx/tZbb1XbdsmSJd6/f/9qt/34xz/266+/3t3dr7/+er/00kvd3X3ZsmU+cuRILy4u9lWrVvmAAQO8rKws6XjXrVvnCxcudHf3wsJCHzRokC9btqzR9rlmzRqfMmWK9+3b1zdv3twoMVf3OgML/ABylkbQIiKp5h6NmGtTujeq1wDHHnssXbt23a98yJAhDB48uNa2Dz74IDNnzqx22xNPPMGsWbMAmDVrFo8//nhl+YwZM8jJyaF///4MHDiQN998M+l4e/XqxdixYwHo0KEDQ4YM4ZNPPmm0ff7whz/kpptu2ufc8YHG3BSUoEVEUq28tHHrNaKHHnqoxgS9ceNGevXqBURJddOmTQB88skn9OnTp7JeXl5eZYKtr9WrV/P2228zYcKERtnnk08+yWGHHcaoUaP2KW/MmBuLLhITkRbHd+3AF74MsQws/zgst22qQ6qdxxu3XiOZP38+bdu23e/cdV28mpF+Q6503rVrF2eeeSa33norHTt2POB9FhUVcd111/Hcc881qP3BphG0iLQovmsH5df8B/FH7yY+5w+U//xivHhPqsOqnSX5pzjZeo1k9uzZNY6eAQ455BDWr18PwPr16+nZsycQjT7Xrl1bWa+goIDevXvv03b+/PmVF6E9+eST+/VdWlrKmWeeybnnnstXvvKVRtnnhx9+yEcffcSoUaPo168fBQUFjB07lg0bNiTV/mBTghaRFsUXvgzFRVBaEp233bUdX5bmy8xmJPkVkcnWawTxeJxHHnmEGTNm1Fjny1/+Mvfddx8A9913H9OnT68snz17Nnv37uWjjz7igw8+YPz48fu0nTBhAosXL2bx4sV8+ctf3mebu3PBBRcwZMgQLrnkkkbb54gRI9i0aROrV69m9erV5OXlsWjRIg499NCk2h9sStAi0rJUN8qMpfmfOjPIyqm9TlZOgz8PPXPmTCZOnMj7779PXl4ed999NwCPPfYYeXl5vP7665x66qlMnTq1ss1LL71EXl4eAwYM2KevCy+8sPIjWZdffjlz585l0KBBzJ07l8svvxyAYcOGcfbZZzN06FBOPvlk7rjjDjIyMpKO99VXX+X+++9n3rx5laPsp59+usH7TIy5Jgcac1Ow6ubdW6r8/Hyv60USkebN9+ym/OcXw64dUUG3Q8m47BYsKzsl8bz77rsMGTIkucqlxdVfzZ2VA1ltGjcwaVTVvc5mttDd8xvapy4SE5EWxXLbkfGT3+BLF0Asho0Yn7LkXG9ZbSAzR98kJoAStIi0QNamLZZ/bKrDaBgzyGwm/1BIk0rzEzMiIiKtkxK0iIhIGlKCFhERSUM6By0ikka8rATf+DHs3QM5udghn8N0TrpV0ghaRCRNxD9aQvzlOfi7r+OrFuPvvk785TnEP1rS4D6b49KNtcUMcPvttzN48GCGDRvGpZdeCsCnn37KCSecQPv27fne975Xa//VtT/QmJuCRtAiImkg/tES/MP911ymvAz/cDFxINZ/ZL37zczM5Je//CVjx45l586djBs3jsmTJzN06FBuuOEGJk2axOWXX84NN9zADTfcwI033sjy5cuZPXs2y5YtY926dZx00kmsWLGi2i/uWLt2LXPnzqVv376VZfVpX9+Yn3/+eZ544gmWLFlCTk5O5WIZbdq04ZprrmHp0qUsXbq0xr5ran+gMTcFjaBFRFLMy0rw1TUnFQBfvRQvK6l3381x6cbaYr7zzju5/PLLycmJvnmt4ru427Vrx9FHH02bNrV/oUtN7bXcpIiI7Mc3fgzlZbVXKi+L6h2A5rh0Y9WYV6xYwcsvv8yECRM47rjjeOutt+rVX03ttdykiIjsb2+Sq22VNHxVrua4dGN1MZeVlbFt2zbeeOMN3nrrLc4++2xWrVqVdP81tddykyIisr+c3OTqZSdZr4p0XLqxoctN5uXl8ZWvfAUzY/z48cRisX0uTqtLTe213KSIiOzHDvkcZNQxoZmRGdWrp3RdurGhy02efvrpzJs3D4imq0tKSujevXvSz0dN7bXcpIiI7Mcys7F+w2uv0294gz4P3RyXbqwt5vPPP59Vq1YxfPhwZsyYwX333Vc5Fd2vXz8uueQS/vjHP5KXl8fy5cv3i7mm9lpuMsW03KSIHGz1WW4y/tGS6GruxAvGMjKxfsMb9BErOXi03KSISAsW6z8S73NEdLV2yR7I1jeJtWZK0CIiacQys7HDBqU6DEkDOgctIiKShpSgRURE0pAStIiISBpKKkGbWWczm2Nm75nZu2Y2MZRfbGbvm9kyM7splGWZ2X1m9k6oe0U1/T1pZksTHueY2UNmttLM5ptZv4Rts8zsg3CblVDeP9T9ILTVVRQi0uz5nt3EX32W+N9nE3/1WXzP7lSHJCmS7Aj6NuAZdz8CGAW8a2YnANOBke4+DLg51P0qkOPuI4BxwLerJNyvALuq9H8BsM3dBwK3ADeGul2Bq4AJwHjgKjPrEtrcCNzi7oOAbaEPEZFmK/732ZRf+Q3iD/ya+FP3E3/g19Hjv89ucJ/NcenG4uJixo8fz6hRoxg2bBhXXXXVfnWqW+JyyZIlTJw4kWHDhjFixAiKi4sPWsxNwt1rvQEdgY8In5lOKH8YOKma+jOBp4iuEO8GrAC6hm3tgVeAocDShDbPAhPD/UxgC2Chr98l1PtdKLNQJzOUTwSeretYxo0b5yIiB9Py5cuTqlf+9INeetG0Gm/lTz/YoP2vW7fOFy5c6O7uhYWFPmjQIF+2bJm7u8+bN88nTZrkxcXF7u6+ceNGd3fftWuXv/zyy37nnXf6d7/73Rr7rqn9smXLfOTIkV5cXOyrVq3yAQMGeFlZWdIxx+Nx37lzp7u7l5SU+Pjx4/3111+v3L5mzRqfMmWK9+3b1zdv3uzu7qWlpT5ixAhfvHixu7tv2bKl2n02VczVvc7AAq8jL9V2S2YEPQDYDNxrZm+b2V1m1g44HDgmTDO/aGZHhvpzgN3AemANcLO7bw3brgF+CRRV2cdhwNrwD0MZsCMk98ryoCCUdQO2h7qJ5fsxs2+Z2QIzW7B58+YkDldE5ODyPbuJP/dIrXXizz2C76n6p7NuzXHpRjOjffv2QPSd3KWlpfssXFHdEpfPPfccI0eOrFxZq1u3btV+E1hLW24yExgL3OnuY4iS7+WhvAtwFPBj4GGLnq3xQDnQG+gP/MjMBpjZaGCguz9WzT6qWzLEG1C+f6H77909393ze/ToUcthioikhi96BfZWPx1baW8x/vYrB7Sf5rR0Y3l5OaNHj6Znz55Mnjy5MuaalrhcsWIFZsbUqVMZO3YsN91000GPubEl80UlBUCBu88Pj+cQJegC4NEwjH/TzOJAd+DfiM5XlwKbzOxVIJ9o1DvOzFaH/fY0sxfc/fjQVx+gwMwygU7A1lB+fEIsecALRNPbnc0sM4yi84B1DXoGRERSrXBbcvV2bK27Tg2a29KNGRkZLF68mO3bt3PGGWewdOlSBgwYUOMSl2VlZbzyyiu89dZbtG3blkmTJjFu3DgmTZp00GJubHWOoN19A7DWzAaHoknAcuBx4EQAMzscyCZKnGuAEy3SjmiE/Z673+nuvd29H3A0sCIkZ4AngYortM8C5oXE/ywwxcy6hIvDphCda3bg+VCX0PaJBj4HIiKp1bFL3XUAOnVtUPfpuHRjXctNVujcuTPHH388zzzzTJ1LXB533HF0796dtm3bMm3aNBYtWtSoMR9syV7FfTHwgJktAUYDPwfuAQaEj0vNBmaFxHkH0cVgS4G3gHvdfUkd/d8NdDOzlcAlRCN0wrnra0I/bwFXJ5zPvgy4JLTpFvoQEWl2bOzRkFP7+V5y2mBjjq53356mSzfWttzk5s2b2b59OwB79uzhH//4B0cccUStS1xOnTqVJUuWUFRURFlZGS+++CJDhw5t1JgPtqS+i9vdFxNNU1f1tWrq7iL6qFVt/a0Ghic8Lq6pjbvfQ/TPQNXyVUTnu0VEmjXLbUdsyleJP3V/jXViU76K5batd98VSzeOGDGC0aNHA/Dzn/+cadOmcf7553P++eczfPhwsrOz91u6sbCwkJKSEh5//HGee+45hg4dyoUXXsh3vvMd8vPza2yfuHRjZmZmvZduXL9+PbNmzaK8vJx4PM7ZZ5/NaaedVmubLl26cMkll3DkkUdiZkybNo1TTz0V4KDE3BS03KSISBOq13KTf58dXc2deMFYThtiU75K7JQZTRShNAYtNyki0oLFTpmBHf/l6GrtHVuhU1dszNENGjlL86cELSKSRiy3LfaFKakOQ9KAFssQERFJQ0rQIiIiaUhT3CLSbHm8nPjKxbBjE3TuSezzo7FYaq+8FWksStAi0mzFl70Km9dCvBx2biVesoeMYfX/rHA6KSssZNtfn6Z04yayDulJl9OmkRm++UtaF01xi0jztaUgSs4Q/dy8tvb6aW7drbfzrzHjWf2jS/nkpptZ/aNL+deY8ay79fYG99lcl27cvn07Z511FkcccQRDhgzh9ddfr3Wfc+fOZdy4cYwYMYJx48ZVfhlJTao75nRbblIjaBFpvrJzYc/Ozx7nNN+PI6279XY+ufEX+5XHi4oqy3v/4OJ695uTk8O8efNo3749paWlHH300ZxyyikcddRRQLRe9Ny5c+nbt29lm7KyMr72ta9x//33M2rUKD799FOysrL26/v555/niSeeYMmSJeTk5LBp0yYAli9fzuzZs1m2bBnr1q3jpJNOYsWKFfX64o///M//5OSTT2bOnDmUlJRQVFRU6z67d+/OU089Re/evVm6dClTp06tcbGL6o65MWJubBpBi0izFRtxLGRmQywDsnKIDT821SE1SFlhIetvv6PWOutvv4PynTtrrVOd5rh0Y2FhIS+99BIXXHABANnZ2XTu3LnWfY4ZM6byu7OHDRtGcXExe/furbb/6o65uS43KSKSlqxjN2LHnk3sC6cTO+arWIckF51IM9v++jTxotrXeo4XFbH1r083qP/mtnTjqlWr6NGjB+eddx5jxozhwgsvZPfu3bXuM9Ff/vIXxowZU5nEE9V0zOm43KQStIg0axaLYW3aYbHm++esdOOmRq1XVcXSjQUFBbz55pssXbqUoqIirrvuOq6++ur96lcs3fjAAw/wyiuv8Nhjj/HPf/6z2noVSzf+4he/4Oyzz8bdD3jpxrKyMhYtWsRFF13E22+/Tbt27bjhhhtq3WeFZcuWcdlll/G73/1uv35rO+ZmudykiIg0raxDejZqvZqk09KNtS03mZeXR15eXuVI/6yzzqrcf21LZBYUFHDGGWfwpz/9ic9//vP7xVvXMTfX5SZFRKSJdDltGrG2tV/gFmvblq6nTat33+m6dGNty00eeuih9OnTh/fffx+Af/7zn5X7r2mf27dv59RTT+X666/ni1/8YrXPRW3HnI7LTSpBi4ikWGbHjvS6+Lu11ul18XfJ6NCh3n2vX7+eE044gZEjR3LkkUcyefLkei3dOHr0aMaOHbvP0o0VqwKef/75rFq1iuHDhzNjxoxql248+eSTG7R04+233865557LyJEjWbx4MVdeeWWt+/zNb37DypUrueaaaypH5hVXeCfGXJPGiLmxablJEZEmVJ/lJtfdejvrb79jnwvGYm3b0uvi7zboI1Zy8Gi5SRGRFqz3Dy7mkAu+ydaEbxLretq0Bo2cpflTghYRSSMZHTrQY+Y5qQ5D0oDOQYuIiKQhJWgREZE0pAQtIiKShnQOWkQkjezZsYNFc55gx/oNdOp1KGPPmk5up06pDktSQCNoEZE08fR1v+Dyw47g/gu/x5P/fS33X/g9Lj/sCJ6+bv9VruqjOS7deP7559OzZ0+GDx++T/k555xT+Tnnfv36MXr06MptNS2R+ZOf/IQ+ffpULhpSk5piXrhwISNGjGDgwIF8//vfr/ZrQZtExfemtobbuHHjXETkYFq+fHlS9f527U3+bTrUePvbtTc1OIZvfOMb/oc//MHd3ffu3evbtm1zd/d58+b5pEmTvLi42N3dN27c6O7uixYt8k8++cTd3d955x3v3bt3jX2vWbPGp0yZ4n379vXNmze7u/uyZct85MiRXlxc7KtWrfIBAwZ4WVlZvWJ+8cUXfeHChT5s2LAa61xyySX+s5/9zN3dS0tLfcSIEb548WJ3d9+yZUvlPl9//XVft26dt2vXrsa+aov5yCOP9Ndee83j8biffPLJ/vTTT+/XvrrXGVjgB5CzNIIWEUmxPTt28Oz1v6q1zrPX/4o9hYX17ru5Lt147LHH0rVr1xq3uzsPP/wwM2fOBGpfIvOoo46iV69ete6vppjXr19PYWEhEydOxMz4xje+weOPP16vY2koJWgRkRRbNOcJ9oblFGuyd/duFs15ot59t9SlG19++WUOOeQQBg0aVHksySyRWZOaYv7kk0/Iy8vbr/xgUIIWEUmxHes3NGq9RC116cYHH3ywcvQMyS+RWZOaYj4Yx1ITJWgRkRTr1OvQRq2XKF2Xbqxtucm6lJWV8eijj3LOOZ9941qyS2TWpKaY8/LyKCgoqPVYmooStIhIio09azo57drVWienXTvGnjW93n2n69KNtS03WZeKJTMTp56TXSKzJjXF3KtXLzp06MAbb7yBu/OnP/2J6dPr/zo0hBK0iEiK5XbqxNQrLqm1ztQrLiG3Y8cG9d8cl26cOXMmEydO5P333ycvL4+77767ctvs2bP3md6G2pfIvPTSS8nLy6OoqIi8vDx++tOfAtE59P/5n/+pM+Y777yTCy+8kIEDB/L5z3+eU045pV7H0v8a+YQAACAASURBVFBablJEpAnVZ7nJp6/7Bc9e/6t9LhjLadeOqVdcwrSf/LipQpRGoOUmRURasGk/+TEnXPzt/b9JrIEjZ2nelKBFRNJIbseOfPH8r6c6DEkDOgctItLEWtOpxNaoqV5fJWgRkSbUpk0bPv30UyXpFsrd+fTTT2nTpk2j960pbhGRJlTxOdrNmzenOhRpIm3atNnnI1+NRQlaRKQJZWVl0b9//1SHIc2QprhFRETSkBK0iBwc7hCPRz9FpE6a4haRphePw95dnyXnnLaQkZXamETSXFIjaDPrbGZzzOw9M3vXzCaG8ovN7H0zW2ZmN4WyLDO7z8zeCXWvCOVtzexvoY9lZnZDQv85ZvaQma00s/lm1i9h2ywz+yDcZiWU9w91PwhtsxvnKRGRRleyZ9+R894ijaRF6pDsFPdtwDPufgQwCnjXzE4ApgMj3X0YcHOo+1Ugx91HAOOAbyck3JtDH2OAL5pZxReaXgBsc/eBwC3AjQBm1hW4CpgAjAeuMrMuoc2NwC3uPgjYFvoQkXTk8VRHINLs1JmgzawjcCxwN4C7l7j7duAi4AZ33xvKN4UmDrQzs0wgFygBCt29yN2fr+gDWARUXJc+Hbgv3J8DTLJowc2pwFx33+ru24C5wMlh24mhLqHt6Q18DkSkqVWdzrYYHKQ1dUWaq2RG0AOAzcC9Zva2md1lZu2Aw4FjwjTzi2Z2ZKg/B9gNrAfWEI2atyZ2aGadgS8BFatpHwasBXD3MmAH0C2xPCgIZd2A7aFuYvl+zOxbZrbAzBboc4giKZKVE91iGVGybtM+1RGJpL1kEnQmMBa4093HECXfy0N5F+Ao4MfAw2FkOx4oB3oD/YEfmdmAis7CyPpB4NfuvqqiuJr9egPK9y90/72757t7fo8ePeo6VhFpCmaQ1SZKzDltNXoWSUIyCboAKHD3+eHxHKKEXQA86pE3gTjQHfg3ovPVpWHa+1Ugcbmt3wMfuPutVfbRByoTeCdga2J5kAesA7YAnUPdxHIREZEWoc4E7e4bgLVmNjgUTQKWA48TnQfGzA4HsokS5xrgRIu0IxphvxfqXUuUfH9QZTdPAhVXaJ8FzPPoi2ufBaaYWZdwcdgU4Nmw7flQl9D2iXoeu4iISNpK9nPQFwMPhI8yrQLOI5rqvsfMlhJdCDbL3d3M7gDuBZYSTUXf6+5LzCwP+AlRsl4UzYbzG3e/i+gCtPvNbCXRyHkGgLtvNbNrgLdCHFcnnM++DJgdkv7boQ8REZEWwVrTCiv5+fm+YMGCVIchIiKtgJktdPf8umtWT1/1KSIikoaUoEVERNKQErSIiEgaUoIWERFJQ0rQIiIiaUgJWkREJA0pQYuIiKQhJWgREZE0pAQtIiKShpSgRURE0pAStEgr5+XleHFRqsMQkSqSXSxDRFqg+NuvEv/jzVBeDv0Gk/G9q7E2uakOS0TQCFqk1fKiXVFyLi2BeDms+YD4X/8v1WGJSKAELdJaFW6DWMKfgLJS2Lg2dfGIyD6UoEVaq+6HQpu2nyXp7BwYNTG1MYlIJSVokVbKMrPI+K+bsaFj4bD+2Je+TuyLJ6c6LBEJdJGYSCtm3Q4h4z9+luowRKQaGkGLiIikISVoERGRNKQELSIikoaUoEVERNKQLhJLE4WvvMqO518gp29fepw7E8vUSyMi0popC6SBrX/9Gx99/4fE9xQTy81lx4svMeieP6Q6LBERSSEl6DSw4be/I76nGID4nj3smPtPynfuJKNDhxRHJiJy4HzPLuJLX4E9hdC1F7EhE7EMpZ+66BlKA7H27fYrs6ysFEQiItL44m//A4p2Ag6bPiYeyyBj6BdSHVba00ViaaDvT/+HWIf2xDp0IJbbhsOuvIxYmzapDktE5IB5vByKCgGPCuJx2LYxpTE1FxpBp4G2Q4cw8rWXKVryDtmH9SZ38OGpDklEpFFYLANy2sLeoooC6NgttUE1ExpBp4ms7t3odOLxSs4i0uLExpwEuR2i5NypB7EjJqQ6pGZBI2gREWlS1r4zGV88I9VhNDsaQYuIiKQhJWgREZE0pAQtIiKShpSgRURE0pAStIiISBpSghYREUlDStAiIiJpSAlaREQkDSlBi4iIpCElaBERkTSkBC0iIpKGkkrQZtbZzOaY2Xtm9q6ZTQzlF5vZ+2a2zMxuCmVZZnafmb0T6l6R0M+4UL7SzH5tZhbKc8zsoVA+38z6JbSZZWYfhNushPL+oe4HoW124zwlIiIiqZfsCPo24Bl3PwIYBbxrZicA04GR7j4MuDnU/SqQ4+4jgHHAtxMS7p3At4BB4XZyKL8A2ObuA4FbgBsBzKwrcBUwARgPXGVmXUKbG4Fb3H0QsC30ISIi0iLUmaDNrCNwLHA3gLuXuPt24CLgBnffG8o3hSYOtDOzTCAXKAEKzawX0NHdX3d3B/4EnB7aTAfuC/fnAJPC6HoqMNfdt7r7NmAucHLYdmKoS2hb0ZeIiEizl8wIegCwGbjXzN42s7vMrB1wOHBMmGZ+0cyODPXnALuB9cAa4GZ33wocBhQk9FsQygg/1wK4exmwA+iWWF6lTTdge6hbta99mNm3zGyBmS3YvHlzEocrIiKSeskk6ExgLHCnu48hSr6Xh/IuwFHAj4GHw8h2PFAO9Ab6Az8yswGAVdO3h581batv+f6F7r9393x3z+/Ro0f1RygiIpJmkknQBUCBu88Pj+cQJewC4FGPvAnEge7AvxGdry4N096vAvmhfl5Cv3nAuoR99AEIU+OdgK2J5VXabAE6h7pV+xIREWn26kzQ7r4BWGtmg0PRJGA58DjReWDM7HAgmyhxrgFOtEg7ohH2e+6+HthpZkeFkfY3gCdCn08CFVdonwXMC+epnwWmmFmXcHHYFODZsO35UJfQtqKvFq1kwwZ2vvkWZdu2pzoUERFpQpl1VwHgYuCB8FGmVcB5RFPd95jZUqILwWa5u5vZHcC9wFKiqeh73X1J6Oci4I9EF4/9PdwgugDtfjNbSTRyngHg7lvN7BrgrVDv6nA+G+AyYLaZXQu8Hfpo0bb+7Wk+uvgHWGYW7s7gRx6k/ehRqQ5LRESagEWD0dYhPz/fFyxYkOowGmzh4UOJ79xV+Th32FCG/+OZFEYk9eEb1sLundBnAJbdJtXhiEgTM7OF7p7f0PbJjqAlxdydeNGefcrKCwtTFI3UV/lj9+AvPAUZGZDbjozLbsU6dqm7oYi0Wvqqz2bCzOh2xulYbjTyiuXm0nPW11MclSTDP90YJefSEijeA4XbiD/zUKrDEpE0pxF0M9L/1ptpN3oku99ZSsdjj6HbGdNTHZIkY08RxBL+Fy4vj6a6RURqoQTdjFhGBodccF6qw5D66tUHuvaETeuhvBSysol98eS624lIq6YELdLELCOTjP+6mfg/H4Od24kdeQI2cFiqwxKRNKcELXIQWG47Mk77WqrDEJFmRBeJiYiIpCElaBERkTSkBC0iIpKGlKBFRETSkBK0iIhIGlKCFhERSUNK0CIizUzxzp38+uQz+G52N6783DDWLFqc6pCkCShBi4g0M3MuuZIVL7xCeWkpW9es5bbJ04nH46kOSxqZErSISDPz8cLFlO3dW/m4eNduirZtS2FE0hSUoOvJN35M+fLXiK9Zjus/VhFJgcEnHktWbi4AFovRoWd32nXtmuKopLHpqz7rIb7mPXzlQoiX47EMfMdmMkYcl+qwRKSVOf3nV1G8o5AlT/2dbv36cv4Dd2FmqQ5LGpkSdD14wXsQL48exMth4xp8WByLaSJCRA6ezOxsvvaH21MdhjQxZZb6yMze93EsBvqvVUREmoASdD3EjpgAGVnRLZaBDZmoaSUREWkSmuKuB+vYjdjRX4Hd26FNe6xNu1SHJCIiLZQSdD1ZVg50PiTVYUhjcw/XFzjEMnXqQuRgcY9+6nduP0rQIu5QsgfKS0OBQW57MJ0BEmkyVX/vstpAVk5qY0oz+gsk4vGE5AzgULq3xuoi0gjKS/f9vSst/uxTMgIoQYt8NsW2T9nBD0OkVakuGbu+/CmRErRILGP/6eys7OrrikjjyMjav8wyDn4caUznoEXMoE17KCuJ/oPPzI6Stog0nYxMyMqFsr3R72BWm+i7JaSSErQIhD8QukBF5KDKytZsVS3074qIiEgaUoIWERFJQ0rQIiIiaUgJWkREJA0pQYuIiKQhJWgREZE0pATdCmy89z4Wj87nX0dOZNszz6U6HBERSYISdAu348WXKLjm55Ru3ERJwSes+o/vsWflh6kOS0RE6qAE3cLtXryEeEnJZwUZmRQtW5a6gEREJClK0C1c26FHEMtJ+Kae8jJyBw9OXUAiIpIUJegWrvPkkzj0Py4i1q4dGZ078bkbr6ftEUrQIiLpzry6pfZaqPz8fF+wYEGqwxARkVbAzBa6e35D2yc1gjazzmY2x8zeM7N3zWxiKL/YzN43s2VmdlMoO9fMFifc4mY2OmybaWbvmNkSM3vGzLqH8hwze8jMVprZfDPrl7DvWWb2QbjNSijvH+p+ENrqG9dFRKTFSHaK+zbgGXc/AhgFvGtmJwDTgZHuPgy4GcDdH3D30e4+Gvg6sNrdF5tZZujnBHcfCSwBvhf6vwDY5u4DgVuAGwHMrCtwFTABGA9cZWZdQpsbgVvcfRCwLfQhIiLSItSZoM2sI3AscDeAu5e4+3bgIuAGd98byjdV03wm8GBFV+HWzswM6AisC9umA/eF+3OASaHOVGCuu291923AXODksO3EUJfQ9vSkj1pERCTNJTOCHgBsBu41s7fN7C4zawccDhwTpplfNLMjq2l7DiFBu3spUVJ/hygxDyUkfeAwYG2oVwbsALollgcFoawbsD3UTSzfj5l9y8wWmNmCzZs3J3G4IiIiqZdMgs4ExgJ3uvsYYDdweSjvAhwF/Bh4OIxsATCzCUCRuy8Nj7OIEvQYoDfRFPcVFdWr2a83oHz/Qvffu3u+u+f36NGjjkMVERFJD8kk6AKgwN3nh8dziBJ2AfCoR94E4kD3hHYz+Gx6G2A0gLt/6NGl4w8DX0jYRx+AcK66E7A1sTzIIxp9bwE6h7qJ5SIiIi1CnQna3TcAa82s4sOzk4DlwONE54Exs8OBbKLEiZnFgK8CsxO6+gQYamYVw9jJwLvh/pNAxRXaZwHzQhJ/FphiZl3CxWFTgGfDtudDXULbJ+px3CIiImkts+4qAFwMPBA+yrQKOI9oqvseM1sKlACz/LMPVR9LNOpeVdGBu68zs58BL5lZKfAx8M2w+W7gfjNbSTRynhHabDWza4C3Qr2r3X1ruH8ZMNvMrgXe5rPz2SIiIs2evqhERESkCRyULyoRERGRg0sJWkREJA0pQYuIiKQhJWgREZE0pAQtIiKShpL9mJWIiKQxd2fLw3PY+eprtB0xnEPO/yaWkZHqsOQAKEGLiLQA62/5Nevv+C3xoj1s/evTFC1dzoDbfpnqsOQAaIpbRKQF2PR/DxAv2gOA79nD1sceT3FEcqCUoEVEWoDMTp32eRxr1y5FkUhjUYIWEWkB+v3qF8TatyejQwdiubkM+PWvUh2SHCCdgxYRaQHajxnNqDdfo/jDD8np24esnj1THZIcICVoEZEWIrNLZ9rnj0t1GNJINMUtIiIClOzZQ+GmzaTLIlJK0CIi0uq9es/9XNK5D1f2HcKNEydRvGtXqkNSghYRkdatcOMmHvzuJZSVlFC2t4SCxe/w92t/keqwlKBFRKR1K9ywkcysrMrHZXv3smnlqhRGFFGCFhGRVu2QwYPIad8ei0UpMbttW8Z99fQUR6WruEVEpJXLatOGy974J4/88Ap2bNjI0RfOIv+cM1MdlhK0iIhI1759+PZf/i/VYexDU9wiIiJpSAlaREQkDSlBi4iIpCElaBERkTSkBC0iIpKGlKBFRETSkBK0iIhIGlKCFhERSUNK0CIiImlICVpERCQNKUGLiIikISVoERGRNKQELSIikoaUoEVERNKQlpuU5svjUFYKZpCRFf0UEWkhlKCleYrHoXgX4NFj2wtt2itJi0iLoSluaZ7KSqhMzhCNpuNlKQtHRKSxKUFLC6LRs4i0HErQ0jxlZbNPQo5lRDcRkRZC56ClebIY5LaH8jLAICNT559FpEVRgpbmy2KQmZ3qKEREmkRSU9xm1tnM5pjZe2b2rplNDOUXm9n7ZrbMzG4KZeea2eKEW9zMRodt2Wb2ezNbEfo6M5TnmNlDZrbSzOabWb+Efc8ysw/CbVZCef9Q94PQVn+ppVnzglXEn3uE+Px5eLw81eGISIolO4K+DXjG3c8KibCtmZ0ATAdGuvteM+sJ4O4PAA8AmNkI4Al3Xxz6+Qmwyd0PN7MY0DWUXwBsc/eBZjYDuBE4x8y6AlcB+USX7C40syfdfVuoc4u7zzaz/w193HkgT4ZIqviKdyi/46poyj4zC/vXa2R86/+lOiwRSaE6R9Bm1hE4FrgbwN1L3H07cBFwg7vvDeWbqmk+E3gw4fH5wPWhftzdt4Ty6cB94f4cYJKZGTAVmOvuW0NSngucHLadGOoS2p6e3CGLpJ/yZx+C0r0QL4eSYnzpAnzH1lSHJSIplMwU9wBgM3Cvmb1tZneZWTvgcOCYMM38opkdWU3bcwgJ2sw6h7JrzGyRmT1iZoeEssOAtQDuXgbsALollgcFoawbsD3UTSzfj5l9y8wWmNmCzZs3J3G4IgefVT2X7h5d+CYirVYyCToTGAvc6e5jgN3A5aG8C3AU8GPg4TCyBcDMJgBF7r40oZ884FV3Hwu8DtxcUb2a/XoDyvcvdP+9u+e7e36PHj1qPVCRVIl96euQkxvdsnOwY0/F2ndMdVgikkLJ/IteABS4+/zweA5Rgi4AHnV3B940szjQnWi0DTCDfae3PwWKgMfC40eIzhtX7KMPUGBmmUAnYGsoPz6hjzzgBWAL0NnMMsMoOg9Yl8SxiKQly+tPxlX/i69chnXuhg0cnuqQRCTF6hxBu/sGYK2ZDQ5Fk4DlwONE54Exs8OBbKLESbgA7KvA7IR+HHiKzxJuRT8ATwIVV2ifBcwL9Z8FpphZFzPrAkwBng3bng91CW2fqM+Bi6Qb69ydWP5xSs4iAiR/FffFwAPhCu5VwHlEU933mNlSoASYFRInRBeVFbj7qir9XAbcb2a3Eo20zwvld4fylUQj5xkA7r7VzK4B3gr1rnb3rQl9zTaza4G3Qx8i0tK5Q2lxdFIrM0vn6qXFss9yasuXn5/vCxYsSHUYItJQ7tEqZh7/rCynnZK0pCUzW+ju+Q1tr+/iFpHmI16+b3KGsLKZSMujBC0izUd137eu72CXFkoJWkSaj1jGvt+/bgZZOamLR6QJ6cSNiDQv2blRknaPErZG0NJCKUGLSPOjtb+lFdAUt4iISBpSghYREUlDStAiIiJpSAlaREQkDSlBi4iIpCElaBERkTSkBC0iIpKGlKBFRETSkBK0iIhIGtI3iYnIAfO1H1L+59th107smFOITT4T01dwihwQJWgROSC+Zzflt14Be3ZHj5/+M/GOXcg4alKKIxNp3jTFLdKIVr0+n1smncYvT5jG+y+8nOpwDo6NBdHCFRVK9sK7i1IXj0gLoRG0SCPZtPJDbp08nZLdRQDc8eZZXP7mC/QeNiTFkTWxrj2hvOyzx1k5cFj/1MUj0kJoBC3SSD546bV9Hsfjcd6f92KKojl4rGMXYuddCjm5YDFsxJHEJp2e6rBEmj2NoEUaSZc+h+3zOCMjky59+6QomoMrNnoiNuoRcMdi+r9fpDHoN0mkkQw56QS+cN7XycjKJCMri3Fnn8GoL09LdVgHjZkpOYs0IvPEiztauPz8fF+wYEGqw5AWbk9hIbiT26lTqkMRkRQys4Xunt/Q9priFmlkuR07pjoEEWkBNB8lIiKShpSgRURE0pAStIiISBpSghYREUlDStAiIiJpSAlaREQkDSlBi4iIpCElaBERkTSkBC0iIpKGlKBFRETSkBK0iIhIGlKCFhERSUNK0CIiImlICVpERCQNKUGLiIikISVoERGRNJSZ6gBE0onHy/ENH0FpCdazD5bbIdUhiUgrldQI2sw6m9kcM3vPzN41s4mh/GIze9/MlpnZTaHsXDNbnHCLm9noKv09aWZLEx7nmNlDZrbSzOabWb+EbbPM7INwm5VQ3j/U/SC0zT6wp0JaO/c48YXP4u+9ia9cRPyNp/DdO1Idloi0UslOcd8GPOPuRwCjgHfN7ARgOjDS3YcBNwO4+wPuPtrdRwNfB1a7++KKjszsK8CuKv1fAGxz94HALcCNoW5X4CpgAjAeuMrMuoQ2NwK3uPsgYFvoQ6ThCrfCru0QLwOPQ3kZ8bXvpToqEWml6kzQZtYROBa4G8DdS9x9O3ARcIO77w3lm6ppPhN4MKGv9sAlwLVV6k0H7gv35wCTzMyAqcBcd9/q7tuAucDJYduJoS6h7el1H65ILWIx8GrKRERSIJm/PgOAzcC9Zva2md1lZu2Aw4FjwjTzi2Z2ZDVtzyEhQQPXAL8EiqrUOwxYC+DuZcAOoFtieVAQyroB20PdxHKRhmvfBboeAhmZ0S0rh1jfIamOSkRaqWQuEssExgIXu/t8M7sNuDyUdwGOAo4EHjazAe7uAGY2AShy96Xh8WhgoLv/MPEcc2DV7NcbUL4fM/sW8C2Avn371nKY0tqZGbFRJ8Knn+BlpVjXXlh2m1SHJSKtVDIj6AKgwN3nh8dziBJ2AfCoR94E4kD3hHYz2Hf0PBEYZ2argVeAw83shYR99AEws0ygE7A1sTzIA9YBW4DOoW5i+X7c/ffunu/u+T169EjicKU1MzOsex6xQ/srOYtIStWZoN19A7DWzAaHoknAcuBxovPAmNnhQDZR4sTMYsBXgdkJ/dzp7r3dvR9wNLDC3Y8Pm58EKq7QPguYF0bizwJTzKxLuDhsCvBs2PZ8qEto+0S9j15ERCRNJfs56IuBB8JHmVYB5wG7gXvCx6VKgFkV09tEF5UVuPuqJPu/G7jfzFYSjZxnALj7VjO7Bngr1Lva3beG+5cBs83sWuDt0IeIiEiLYJ/l1JYvPz/fFyxYkOowRESkFTCzhe6e39D2+gyJiIhIGlKCFhERSUNK0CIiImlICVpERCQNKUGLiIikISVoERGRNKQELSIikoaUoEVERNKQErSIiEgaUoIWERFJQ0rQIiIiaSjZxTJEWiR3x1f9C/9kBWRmExv6Baxzz1SHJSKiEbS0br7+Q/zjZVBSDEWFxN/+B15SnOqwRESUoKV1820bIV6eUGJQVJiyeEREKihBS6tmHbtBLOOzAo9DbofUBSQiEihBS6tmeYOh14AoSWflEBt5HJaTm+qwRER0kZi0bmZGxpCJMGRiqkMREdmHRtAiIiJpSAlamh93vLQEL/gI37Yl1dGIiDQJTXFL8+EOJXvwwq2U3/4/sHMHlJdj02aScfI5qY5ORKRRaQQtzUdZCZSXEp/7KGzbAnuLoawUf/pBfLtG0iLSsihBS/MRj0c/C7dDecJnlzMyYPeu1MQkItJElKCl+ciIzsjY+OMhK/uzso5d4JDDUheXiEgT0DloaT4ys8DbEBsyGr55Cb7gZejcndgpM7DMrFRHJyLSqJSgpXnJyom+UGTMMTDmmFRHIyLSZDTFLSIikoaUoEVERNKQErSIiEgaUoIWERFJQ0rQIiIiaUgJWtLGqtfnc/u0M/nNqWex+q2FqQ5HRCSl9DErSQsb3lvBrSdNp6SoCIAVL77Cfy95nR4D+qc4MhGR1NAIWtLCu/94Ho8nfH2nw4rnX05dQCIiKaYELWmh82G9iWV8NqFjZnQ+rFcKIxIRSS0laGkS7k7ppk2U70puEYvRp5/GqDNOIzM7m4zsbPJnnsnQqSc1cZQiIunL3D3VMRw0+fn5vmDBglSH0eLFi4t5f+bX2f322xB38q64jEMv+lZSbXdu3oKZ0b57tyaOUkSkaZnZQnfPb2h7jaCl0W286x52L/4XvrcELy2l4KabKf5odVJtO/ToruQsIoIStDSBvWvW4sXFlY9jWVmUbtyYwohERJofJWhpdF2mnUIsNzd6EIthWZm0HTY0tUGJiDQzStDS6DodfywDfnMbHY6aQOepUxjytyfI6NAh1WGJiDQrSSVoM+tsZnPM7D37/+3de7BVZRnH8e8Pj0CCF0RMFBMvXBJGwbwkJqGO9xJHMzU1TcuxUbNhUGHIS1l5ndGa0RrG+2RiYGnZKJI6aZM35BbIIIiopCYIeE8Rnv5Y79HF8XA47LPb+93w+8ysOWu977ve9T5n7Xc/Z629z97SXEn7p/LzJc2TNEfSNansFEkzSstqSUMkbSbpr6mPOZKuKvXfRdI9khZIelpS31Ld6ZLmp+X0UvnOqe38tG/nav1SrON6HHUEA/80kX63jqdr3771Ho6ZWcNp7xX0r4CHImIgsCcwV9JBwEhgj4gYBFwHEBF3RcSQiBgCnAYsiogZqZ/rUh9DgQMkHZnKzwKWR8RuwPXA1QCStgYuA/YD9gUuk9Qj7XM1cH1E9AOWpz7MzMw2COtM0JK2AIYDtwBExMcRsQL4IXBVRHyUyt9sZfeTgbtT/QcR8VhzH8A0oE9qNxK4I61PAg6RJOBwYEpELIuI5cAU4IhUd3BqS9r32PUJ3MzMLGftuYLeBVgC3CZpuqSbJXUD+gMHptvMf5e0Tyv7nkhK0GWStgK+CTySinYAXgWIiE+At4Ge5fJkcSrrCaxIbcvlZmZmG4T2JOgmYC/gNxExFHgfGJPKewBfBS4E/pCubAGQtB/wQUTMLncmqYkiaf86IhY2F7dy3Kig/HMknS1pqqSpS5YsWXuUZmZmGWlPgl4MLI6Ip9P2JIqEvRj4YxSeAVYD25T2O4lWrp6B8cD8iLihxTF2hE8T+JbAsnJ50gd4DVgKbJXalss/JyLGR8Tek5XtXQAACIdJREFUEbF3r1692hGumZlZ/a0zQUfEG8CrkgakokOA54H7KF4HRlJ/oDNF4kRSJ+AEYEK5L0k/p0i+P25xmD8Dze/Q/hbwaBSfQToZOExSj/TmsMOAyanusdSWtO/97YzZMrFq5cp6D8HMLFvtfRf3+cBdkmYBQ4BfArcCu0iaTZGIT4/PPth7OMVVd/MtbCT1AcYBuwPT0r9gfT9V3wL0lLQAGEVxC52IWAZcATyblp+lMoCLgVFpn56pD2sASxa+xCX9hnBul20Yt/Ng/jN/Qb2HZGaWHX9ZhtXclfuM4JVpM4jVq0Fi+8Ff5tJZT9V7WGZmVeUvy7CG8+aCF4vkDBDB0oWL6joeM7McOUFbzQ046ECaunQBoKlzZ3Y7cFidR2Rmlh8naKu5M+4cz1dOPI5eu+7M0OOP4Qf33F7vIZmZZadp3U3Mqqtr9+6ccftveevlV1j9ySd09RdpmJl9jhO01VxEcNt3z2b6vfcDot/wYZz7wEQ2afLD0cysmW9xW83NnfIoM+97gJUf/peVH37Ii/94iqkTJq17RzOzjYgTtNXcu0uWrrG9auVK3n1z6Vpam5ltnJygreYGHvx11GkTmj+6vdOmTQw++vA6j8rMLC9+0c9qbsve2zHmmcd48BfXsmrlSg4d/SO2G9Cv3sMyM8uKE7TVxXYD+vG9O8fXexhmZtnyLW4zM7MMOUGbmZllyAnazMwsQ07QZmZmGXKCNjMzy5ATtJmZWYacoM3MzDLkBG1mZpYhJ2gzM7MMOUGbmZllyAnazMwsQ07QZmZmGXKCNjMzy5Aiot5jqBlJS4CX12OXbYCl/6fh1IpjyINjyINjyMPGEsNOEdGr0gNsVAl6fUmaGhF713scHeEY8uAY8uAY8uAY2se3uM3MzDLkBG1mZpYhJ+i2ja/3AKrAMeTBMeTBMeTBMbSDX4M2MzPLkK+gzczMMuQEbWZmlqENNkFL2lHSY5LmSpoj6YJW2oyQ9LakGWm5tEX9JpKmS3qgVLanpCcl/UvSXyRtkcoPlfRcKn9O0sGNFkOp/kuS3pM0uhFjkLRHqpuT6rs2UgySNpV0RyqfK2lsR8ZfjRgkLUrjmSFpaql8a0lTJM1PP3uU6sZKWiBpnqTDGy2GHOd0Jech1Wczpyt8LGU1pyt4LFU2pyNig1yA3sBeaX1z4AVg9xZtRgAPtNHHKOD35TbAs8DX0/qZwBVpfSiwfVofDPy70WIo1d8LTARGN1oMQBMwC9gzbfcENmmwGL4DTEjrmwGLgL71jCGNYZtWyq8BxqT1McDVaX13YCbQBdgZeLHe56GCGLKb0+sbQ6k+mzldwXnIbk5XEENFc3qDvYKOiNcjYlpafxeYC+zQ3v0l9QGOBm5uUTUAeDytTwGOT8eYHhGvpfI5QFdJXSqPoPYxpH2OBRZSxNBhdYjhMGBWRMxMx3wrIlZVHkFdYgigm6Qm4AvAx8A7FQdAx2Now0jgjrR+B3BsqXxCRHwUES8BC4B9O3KgWseQ45xuw9rOQ3Zzug1riyG7Od2GtcVQ0ZzeYBN0maS+FH8NP91K9f6SZkp6UNKgUvkNwEXA6hbtZwPHpPUTgB1b6fN4YHpEfNSRcZfVIgZJ3YCLgZ9Wa9xlNToP/YGQNFnSNEkXVWv8ULMYJgHvA68DrwDXRcSyqgRAxTEE8HC61Xt2qfyLEfE6FE96wLapfAfg1VK7xVTnCRCoWQxluczp9Yoh0zm9vuchxzm9vjFUNqc7ersj9wXoDjwHHNdK3RZA97R+FDA/rX8DuClauc0BDAQeTn1eBrzVos9BFLfzdm20GIDrgG+n9cupwu2wOsQwGniJ4nNyNwOeBA5psBgOAO4CNqWY4POAXeoVQ9puvtW7LcWt6+Fpe0WLPpannzcCp5bKbwGOb6QYSttZzOkKz0NWc7rCGLKa0xXGUNGcrsqJynVJv4zJwKh2tl+UHgRXUvzFvwh4A/gA+F0r7fsDz5S2+1C8lnFAI8YAPJHaLwJWAMuA8xoshpOA20t1lwAXNlgMNwKnlepuJT3J1iOGVsovJz3Rpyea3mm9NzAvrY8Fxpb2mQzs30gxpO1s5nSF5yGrOV1hDFnN6QpjqGhOV+UBl+MCCLgTuKGNNtvx2Ye17Etx60Et2oxgzauebdPPTqn/M9P2VhR/SVXlKqEeMaztQddIMQA9gGkUf2k3AX8Djm6wGC4GbkvH7QY8D+xRrxjSGDZP5d2AfwJHpO1rWfNNMdek9UGs+SaxhXT8jT21jiGrOV1JDC36vZw6z+kKz0NWc7rCGCqa01V50OW4AF+jeJ1gFjAjLUcB5wDnpDbnUbxxYibwFDCslX5GsOaT6gUUf1G/AFxVOoE/oXiNYUZp2baRYmixT7Umc81jAE5N/c2mlSeq3GOguO02MfX3PNW5Wqg4BmCXVDYz1Y8r9dsTeASYn35uXaobR3FreB5wZKPFQGZzutLzkNOc7sBjKZs5XeFjqaI57Y/6NDMzy9BG8S5uMzOzRuMEbWZmliEnaDMzsww5QZuZmWXICdrMzKxE0q2S3pQ0ux1try99ocYLklZUbRx+F7eZmdlnJA0H3gPujIjB67Hf+cDQiDizGuPwFbSZmVlJRDxO8alrn5K0q6SH0udvPyFpYCu7ngzcXa1xNFWrIzMzsw3YeIoPMZkvaT/gJuDT7wiXtBPFp+Y9Wq0DOkGbmZm1QVJ3YBgwUVJzccuvHj0JmBQd/CrMMidoMzOztnWi+KaqIW20OQk4t9oHNTMzs7WIiHeAlySdAKDCns31kgZQfKnHk9U8rhO0mZlZiaS7KZLtAEmLJZ0FnAKcJan5SzJGlnY5GZgQVf63KP+blZmZWYZ8BW1mZpYhJ2gzM7MMOUGbmZllyAnazMwsQ07QZmZmGXKCNjMzy5ATtJmZWYb+B9cFNxb/uf/kAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create a figure with one subplot\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "\n",
    "# Plot the points with population info\n",
    "join.plot(ax=ax, column='pop18', cmap=\"Reds\", markersize=15, scheme='quantiles', legend=True);\n",
    "\n",
    "# Add title\n",
    "plt.title(\"Amount of inhabitants living close the the point\");\n",
    "\n",
    "# Remove white space around the figure\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a similar way, we can plot the original population grid and check the overall population distribution in Helsinki:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg8AAAGoCAYAAADB1AZpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeZwU1b3//9cHEBAiooKKguIOIrK7JNEbNRqzEjUq7lGMmkRvlm+M3J83MXGFRBO9EfHmBgwuCe5KXEOiudeYiOASFwzigmEAFQREQYRhPr8/qqbndNOnqZruYQZ4Px+PftBUV1VXVVfPnPl8zjkfc3dEREREsmrX2gcgIiIiGxc1HkRERCQXNR5EREQkFzUeREREJBc1HkRERCQXNR5EREQkFzUeZJNlZp8xs7oqtr/RzH5Uy2NqS8xsFzP70Mzat/axyPqZ2f9nZr+p8PpcM/vshjwm2Xyp8SAbRPqD7aP0l9U7ZnaTmX2itY+rkZl93cz+Gi5z9/Pc/bIWeK8zzOwZM1tuZnVm9jMz6xC8vq2Z3WtmK8zsLTM7OXito5ndlV5PN7PPlOy7U9roecfMlpjZH8xs53LH4e7/cvdPuPvaZpzDQWY2LX2PRWZ2p5n1Cl7/iZmtST/vxsfuwet9zexxM1tpZv+s1S+9DNf2L2a2Kjim2SXbH5Eez8r0+HatxXHVgrtf6e5nt/ZxiIAaD7JhfdndPwEMBUYA/9nKx9NaugDfBXoABwJHAD8IXh8PrAZ2AE4BJpjZgOD1vwKnAm+X2fd3gIOB/YGdgGXAr2p8/ADbAL8G+gK7Ah8AN5Wsc3vaOGl8vBG89nvgOWA74GLgLjPrWYPjWt+1BTg/OKZ9GheaWQ/gHuBHwLbATOD2GhxT1cIG0MZqUzgHaaLGg2xw7j4feBjYD8DMdjKzqelfsa+Z2Tca103/gr3LzG43sw/M7FkzGxS87ma2Z/D/35rZ5eXe18zGmNnr6X5mmdkx6fL+wI3Awelfo8vK7cvMvpEe35L0eHcqOY7zzGyOmS01s/FmZpHzn+DuT7j76vRa3AZ8Kt1PV+A44Efu/qG7/xWYCpyWbrva3a9Nl5eLGOwGPOru77j7KmAKMKDMeo1//XvjD/X0r/LLzOzJ9Br9Mf2FWu4cHnb3O919ubuvBK5vPIf1MbO9SRqQl7j7R+5+N/Biet7l1n/IzK4J/n+7mU2KHFf02mZwLPByel6rgJ8Ag8ysX+S4/mJml5vZ39L75g9mtp2Z3ZZGPmaYWd/Ym5nZ6Wlk6T0z+5EFaYfgvr/VzJYDX0+X3Rpsf1qw/cWVTszMvpDe8x+Y2Xwz+0Hw2oVmttDMFpjZWeF3Kj3Hs4N1iyJ0Znadmc1Lz/cZMzskeK3cObQLvofvmdkdZrZtun7ndN33zGxZev12qHRe0nrUeJANzsz6AF8g+csTkr9C60j+Uv4acKWZHRFsMhK4k+Svwd8B95nZFs1469eBQ4CtgZ8Ct5pZL3d/BTgP+Hv612j3Msd8OHAVcALQC3iL5Bdz6EskEZVB6Xqfy3hchwIvp8/3Bta6+6vB6/8g0gAoYyLwqbRB1oUkcvFwxm0BTgbOBLYHOrLuX+0x4Tk0+nLa0HrZzL4ZLB8AvOHuHwTLKp3jWcBpZna4mZ1Cco2/U8VxXWVmi9NG0mdKjusfjf9x9xUk90ylaz+KpGG3M7AH8HeSCMy2wCvAJeU2MrN9gRtIPp9eJPdkaXppJHAX0J2kEVS6/YT0vXciieD0rnCcE4Fz3X0rkkb7Y+l+jib5jI8E9gLypo9mAINp+m7eaWadK5zDvwNfBf4tPe6lJJE2gDNIrkOf9HzOAz7KeTyygajxIBvSfelf9X8F/pekkdAH+DRwkbuvcvfngd+Q/qWdesbd73L3NcAvgM7AQXnfPP2LcoG7N7j77cAc4ICMm58CTHL3Z939Y+A/SCIVfYN1xrr7Mnf/F/A4yQ/ViszsTGA4cHW66BPA+yWrvQ9slfE4XwX+BcwHlgP9gUszbgtwk7u/6u4fAXeQ7Rz2B34MXBgsviN9757AN4Afm9lJ6Wu5ztHd3yb5RTIZuA44vaThETuu0msLcBGwO8kv6l8DfzCzPZpzXKmb3P11d3+fpJH2urv/yd3rSRq8QyLbfQ34g7v/1d1Xk1y/0kJDf3f3+9L7tfSX6NeAB9z9/9L78UdAQ4XjXAPsa2bd3H2puz+bLj8hPYeX0sbSTyrsYx3ufqu7v+fu9e5+DdAJ2CdYpfQczgUudve69Lh/AnwtjX6tIWk07Onua939GXdfnud4ZMNR40E2pK+6e3d339Xdv5X+MNkJWFLyy+Ativ8Km9f4xN0baIpS5JKGiZ9PQ6LLSP4CKxuWL2On9Lgaj+ND4L2S4wz7IKwk+WVU6Xi+CowFPu/ui9PFHwLdSlbtRtKnIIsJJI2r7YCuJDn8PJGHvOewZ7r/77j7E43L3X1W2lBb6+5/I/ml/7X05eac4wNAe2B2mrKpKHJtcffp7v6Bu3/s7pOBJ0miYM09rneC5x+V+X/s+u1E8X29kuR+Cs0jrnT7FWW2Dx1Hcp5vmdn/mtnB5fZDcI9nYWb/z8xeMbP30+/U1hR/p0rPYVfg3uA7+ApJ+m0H4BbgUWBKmkL5WTMjjLIBqPEgrW0BsK2ZhX/d7ULyl3OjPo1PzKwdSXh2QbpoJUknuUY7lnsTS3rN/w9wPrBdmpp4CWjsl7C+8rILSH7wNe6vK8kv6PnRLSpIw8X/Q9KJ9MXgpVeBDma2V7BsEOuG3mMGAb919yXpX3a/Ag6I9V2oRnpN/wRc5u63rGd1p+lavwzsXvKZr+8cryD5RdMriGDEjit2bbMcV9ifpitJKiLrtc9jIUGawcy2JLmfSo+t0vbh96JLme2bduQ+w91HkqSj7iOJDK2zH5LvXmgFke9X2r/hIpLoxTbpd+p9mq5nuXOYR9Kg6x48Orv7fHdf4+4/dfd9gU+SpAFPj52TtC41HqRVufs84G8keejOaQh8NMU53mFmdmwa2vwu8DHwVPra88DJZtY+/aXxb5G36kryg2wRFELa+wWvvwP0NrOOke1/B5xpZoPNrBNwJTDd3efmO+NC/4nbgOPc/enwtfQvyHuAS82sq5l9iiRvfEuwfacgr9wxvW6NP7BnAKeb2dbpX23fAhaEf33XgiXDPx8Dxrv7jWVeH2lm21jiAJJc9/3pOb5K8rldkh77MSSjQ+6OvNehJP0wTk8fv7LI8NNK19bMupvZ59L37JD2nziU5K9dgHuB/czsuPT6/hh4wd3/mefaZHQXSZ+QT6b33E8p/qWbZfsvmdmn0+0vJfLz3JLhvaeY2dZp6m85TZ1t7yDpyLhv2gAp7aPxPHCsmXVJo0yjg9e2AupJvlMdzOzHrBu5KXUjcEXa8MTMeprZyPT5YWY20JJ5R5aTpDFyDyOWDUONB2kLTiIZ8reA5Af4Je4+LXj9fuBEks5VpwHHpj8EIek492WSIYmnkPxVtQ53nwVcQ9Kh7R1gIEnIutFjJH9hvm1m6/yidfc/k+SV7yb5a20Pks5yzfEjkvDuQ9Y030CYWvgWsCXwLkln0m+6e/jX72ySkPjOJL/4PqIpKvIDYBVJf45FJKHqY5p5nJWcTdJ34JLgHD4MXh8FvEYS8r8ZGJemCcLXh5N8pmOBr7n7otI3MbNu6fbnp3+d/pWk899NQYMpVOnabgFcTnJdFgMXkKTSZgOk738cSZRjKclQz+Z+xhWln+cFJJ1uF5Jcp3dJGsZZt/82SaN2YXq8lSZEOw2Ya8moh/NIhvri7g8D15Lc/6+l/4Z+STJs+B2SPidho/5RkpTVqyTpjlVUTrVAkr6aCvzRzD4g+SPgwPS1HUkaRctJokz/C9xabifS+sx9fdFakdZjZj8h6UB1amsfi0hLsWTCtGXAXu7+Zisfi6fH8VprHoe0bYo8iIi0AjP7cpoO6EoyIuRFYG7rHpVINmo8iIi0jpEkqboFJHMsjHKFgmUjobSFiIiI5KLIg4iIiOSyWRUq6dGjh/ft27e1D0NERKTNe+aZZxa7e9mCdZtV46Fv377MnDmztQ9DRESkzTOz6IyjSluIiIhILmo8iIiISC5qPIiIiEgum1WfBxERyWfNmjXU1dWxatWq1j4UaSGdO3emd+/ebLFF9iKmajyIiEhUXV0dW221FX379qV8ORHZmLk77733HnV1dey2226Zt1PaQkREolatWsV2222nhsMmyszYbrvtckeW1HgQEZGK1HDYtDXn81XaQkREascd1q4BbwBrB+23ADU+NjmKPIiISG2sWQUfLYfVH8Gaj5N/P1qeLK/CWWedxfbbb89+++1XtHzJkiUceeSR7LXXXhx55JEsXbq08NpVV13FnnvuyT777MOjjz5adr/Vbt+cY2509dVXY2YsXrw413teeOGF9OvXj/33359jjjmGZcuW1eyY81DjQUREqrdmVdJgKPvax1U1IL7+9a/zyCOPrLN87NixHHHEEcyZM4cjjjiCsWPHAjBr1iymTJnCyy+/zCOPPMK3vvUt1q5dW/Ptm3PMAPPmzWPatGnssssuhWVZ3/PII4/kpZde4oUXXmDvvffmqquuqtkx56HGg4iIVMc93nBotObjZL1mOPTQQ9l2223XWX7//fdzxhlnAHDGGWdw3333FZaPGjWKTp06sdtuu7Hnnnvy9NNP13z75hwzwPe+9z1+9rOfFfU1yPqeRx11FB06JD0ODjroIOrq6mp2zHmo8SAiItVZu6a262X0zjvv0KtXLwB69erFu+++C8D8+fPp06dPYb3evXszf/78mm/fHFOnTmXnnXdm0KBBRcub856TJk3i85//fIsfcznqMCkiItXxhtquVyUvE+HIM6Kg2u1jVq5cyRVXXMEf//jHqt/ziiuuoEOHDpxyyinN2r5aijyIiEh1LOOvkqzrZbTDDjuwcOFCABYuXMj2228PJH91z5s3r7BeXV0dO+20U023nz59OoMHD2bw4MFMnTo10/G+/vrrvPnmmwwaNIi+fftSV1fH0KFDefvttzMfM8DkyZN54IEHuO222woNhDzb14IaDyIiUp32Gac1zrpeRl/5yleYPHkykPxCHTlyZGH5lClT+Pjjj3nzzTeZM2cOBxxwQE23P/DAA3n++ed5/vnn+cpXvpLpeAcOHMi7777L3LlzmTt3Lr179+bZZ59lxx13zHzMjzzyCOPGjWPq1Kl06dKl6FyybF8z7r7ZPIYNG+YiIpLdrFmzsq24+iP3Fcuij4bl73nD+4u84f1FuY9h1KhRvuOOO3qHDh1855139t/85jfu7r548WI//PDDfc899/TDDz/c33vvvcI2l19+ue++++6+9957+0MPPVRYPnr0aJ8xY0azt6/2mEO77rqrL1rUdD2yHPMee+zhvXv39kGDBvmgQYP83HPPrckxl/ucgZke+X1q3szerxuj4cOH+8yZM1v7MERENhqvvPIK/fv3z7ZybLjmFp3wjz4s/Ne69ajR0UmtlPuczewZdx9ebn11mBQRkdrYojN06FR+hsmg8SAbv0x9Hsysu5ndZWb/NLNXzOzgdPkFZjbbzF42s5+lyzqa2U1m9qKZ/cPMPhPs5y/p+s+nj+3T5Z3M7HYze83MpptZ32CbM8xsTvo4I1i+W7runHTbjjW5IiIi0nxm0KFj2pDoqKmpN1FZO0xeBzzi7v2AQcArZnYYMBLY390HAFen634DwN0HAkcC15gVdbE9xd0Hp49302WjgaXuvifwS2AcgJltC1wCHAgcAFxiZtuk24wDfunuewFL032IiIhIC1tv48HMugGHAhMB3H21uy8DvgmMdfeP0+WNDYF9gT8Hy5YBZXMmgZHA5PT5XcARlow/+Rwwzd2XuPtSYBpwdPra4em6pNt+df2nKyIiItXKEnnYHVgE3GRmz5nZb8ysK7A3cEiaOvhfMxuRrv8PYKSZdTCz3YBhQJ9gfzelKYsfWdMMFjsD8wDcvR54H9guXJ6qS5dtByxL1w2Xr8PMzjGzmWY2c9GiRRlOV0REam3lmwsKD39rTuGRhS9fXHhI25Cl8dABGApMcPchwApgTLp8G+Ag4ELgjrQxMInkl/lM4Frgb0DjL/lT0nTGIenjtHR5uaSYN2P5ugvdf+3uw919eM+ePddzqiIiUg2vX03D/Dk0vPECDfPn4PWrW/uQpAVkaTzUAXXuPj39/10kjYk64J50OOjTQAPQw93r3f17aZ+GkUB3YA6Au89P//0A+B1JP4bG9+gDYGYdgK2BJeHyVG9gAbAY6J6uGy4XEZFW0vDmCzQ8cRf+yt/xN57HX/k7DU/cRcObL1S137O+/e/ssEd/Bh50SNnX22J561hJ7n/84x8cfPDBDBw4kC9/+cssX74cgPfee4/DDjuMT3ziE5x//vnR/bZkGfE81tt4cPe3gXlmtk+66AhgFnAfSb8DzGxvoCOw2My6pGkNzOxIoN7dZ6VpjB7p8i2ALwEvpfucCjSOpPga8Fg6QcWjwFFmtk3aUfIo4NH0tcfTdUm3vb+5F0FERKrT8OYL+OvPw9r64hfW1uOvP8+W3aDroP3pOmh/bNe9Co+YMFXx9ZNH8fDdU9ZZDm23vHWsJPfZZ5/N2LFjefHFFznmmGP4+c9/DkDnzp257LLLuPrqq9fZJtSSZcTzyDra4gLgNjN7ARgMXEmSntjdzF4CpgBnpL/UtweeNbNXgItoSk10Ah5N9/E8MB/4n/S1icB2ZvYa8H2StAjuvgS4DJiRPi5Nl5Hu+/vpNtul+xARkQ3M61fjc1+qvM7cl5qdwjj0U59k2222KftaWy1vHSvJPXv2bA499FAgabzcfffdAHTt2pVPf/rTdO7cueJ+W7KMeB6ZJoly9+cpP2Li1DLrzgX2KbN8BUnnyXL7XwUcH3ltEklDpXT5GzSlPUREpJX4O2+tG3EotbYef+ctbOd4tCGvSuWtDzrooML/s5a3PvHEE5u9fVb77bcfU6dOZeTIkdx5551FxayyqFRGvKWOuRwVxhIREVj5ftMjr48/yrbe6ozrAda1e9OjWw9sq22hXfvC6ytXruSKyy/n0kt+vE7DpVzZhbZS3nrSpEmMHz+eYcOG8cEHH9CxY23mN2zJYy5HjQcREalOpy2zrdcx43oZvP7mXN6cO5dBQ4fRd489N3h56+aU5Abo168ff/zjH3nmmWc46aST2GOPPfKcdtVlyGtFjQcREamK7bArtF9PFrx9h2S9Ghk4YF/eXbiAua+/xtzXX9vg5a2bU5IbKKQZGhoauPzyyznvvPNynXe1ZchrRY0HERGBLls3PXKyDh2xvvtVXmeXfTFrl6QYwkfEipdmFR4nfe1YDj7oQGbPnk2ffQcz6fZ7sK7d8RXLCo/QgAEDOOGEE9h33305+uijGT9+PO3bJymPs88+m8bqyueffz4ffPABRx55JIMHDy78Iq+0fVYnnXQSBx98MLNnz6Z3795MnJj06f/973/P3nvvTb9+/dhpp50488wzC9v07duX73//+/z2t7+ld+/ezJo1a51jHjNmDNOmTWOvvfZi2rRpjBkzpmbHnIdKcouISFSektwNb76QjLoIGwXtO2B996PdLvuW3ygSsVjxj6a5Ibrs1hR+t67dC8/DRoPKfFdHJblFRKRVtNttf7xPv2T0xeqPoOOW2A67Yh06rn80hmxU1HgQEZGaMWuH7bhb8cK19evvE1Gi66D9i7dvFOxH0YbWoz4PIiIikosaDyIiIpKL0hYiIhLX0FC2v4LXvVl4br2b0hRFIx+Cok1hHYuiUtzBtNNhGiJcp1INjBYRSZNIE10VERGpGV+1En9xBixfBu07Yv0GYVknkZKNhtIWIiJSEw2P3kHDld/F756ET7sHf2QKDeMvoeFvf6xqv61S3nrsOPbcp1+zylvPmzePww47jP79+zNgwACuu+66wmsnnnhiYWbKvn37MnjwYACmTZvGsGHDGDhwIMOGDeOxxx7Lf8xtqSS3iIhsxtq1S0L3JeH7sKy2r1jG2j/cTMMDt8Hqj4u3X7Maf+IhGh6eUli0ctlHTY83FxQe/tacwiP09dNO45EHH1jn0IrKW4/8Cj8fNw7W1heXt/aG6IRU0fLWs19lyh13NLu8dYcOHbjmmmt45ZVXeOqppxg/fnxhwqfbb7+9MDPlcccdx7HHHgtAjx49+MMf/sCLL77I5MmTOe2008rue2MryS0iIlKWr1qJ/2XdX+6hhj/eiX+0sln7P/TQQ9Zf3vqzn+Xue+8FWr+8da9evRg6dCgAW221Ff3791+nwqW7c8cdd3DSSScBMGTIkEItigEDBrBq1So+/rikIdaCx5yXGg8iIlIVf3HGuhGHUh+vwp/7a03ft7G8NcCdd91V0/LWffr0KaxXTXnruXPn8txzz3HggQcWLX/iiSfYYYcd2GuvdTuD3n333QwZMoROnTq1yjFnocaDiIhk05i+aN8BX7648GD5svVvC4WURJfuWxYeoX8t+qDwCFMbvmIZvnI5NKwtqosx8b+uYfx11zJs8CCWL1lMxy22WKfORXN4w9qm90o1p7z1hx9+yHHHHce1115Lt27dil77/e9/X4g6hF5++WUuuugi/vu//zvfMaskt4iIbFS6dV//OgCf6Lb+dXLot/dePHrfncz8vz9z0teOZY/d+ubaPlreeufezJtXV1ivOSW516xZw3HHHccpp5xS6NfQqL6+nnvuuYcTTzyxaHldXR3HHHMMN998c7RUt0pyi4jIJsEGjoCO64bYi2zREdtncE3f991Fi4CkvPUVP/8F5551Rq7to+Wtv/wlptxxe7NLcrs7o0ePpn///nz/+99f533/9Kc/0a9fP3r37l1YtmzZMr74xS9y1VVX8alPfSr/Maskt4iItHXWrUfh0W67nWh31Ncqrt/uc8fTbo/+yYRS22zT9Mjg5HPO55NHfZHZc16j9659mThpEgBTpj7MPsM/Sf8DPk2vHXfgzFNPLmxTKG89+WZ679o3V3nrfbfuwvFHfZZ9+/dvVnnrJ598kltuuYXHHnusEJ146KGHCq9PmTJlnZTF9ddfz2uvvcZll11W2KaxP4NKcrcyleQWEcknU0nutG9Aw6N30PDHu4o7T3bshH3mS7T/8umFRWG/hJVvLig8X7ymqY9Bjy2ahoZ23S9SzjsQLc/djNkiW3V2y1aiktwiItIq2n3uBBh+SDL64oNlsFV3bOAIrLNmmNzUqPEgIiJVWfHSrOIFPfaEHtBlt/Id9sLIQJegevcuwWyJK5d9VHgeG0ER7idanrsZtSmKog2qc1GW+jyIiIhILmo8iIiISC6KwYiISFW6Dtq/8DzsbPivV/9VeL5Lz/Llua1rMEdE8Lxr0yjGovLfYTojTHlE0xbNkSFV4csXt8x7byTUeBARkZqp//BDlj72f6xZvISPOnWh06c+SbuuXVv7sKTGlLYQEZGaWHDtr3jhyyfy1pXXsODXN/HBr8bz3tdHs+L2O5u9z3nz5nH4qNPY94ij2e/ILzD+ttsKr436+tkM+fRnGPLpz9S+vHUVJbkbrV27liFDhvClL32paPmvfvUr9tlnHwYMGMAPf/hDAG677bbC/A6DBw+mXbt2PP/88/mOeQOW5NY8DyIiEpVpngeShsP8cT+Pvr7zRRey03cvAIpD/jGN6YyFCxey8N1FDB06lA8++IBB++3Hr39+DXvtvnvR+v/1+9vYeuut+fGPf8xzzz3HDjvswE477cSLTz3B0ceeQN0/X1wnvfDDH/6QbbfdljFjxjD2yitZunQp48ZexcvPTOfk0ecy/bFHWfjhaj772c/y6quv0p7g92WGkRe/+MUvmDlzJsuXL+eBB5Kqo48//jhXXHEFDz74IJ06deLdd98tTDHd6MUXX2TkyJG88cYb6+yz6JjHjk2Oedw4Zs2axUknncTTTz/NggULmo4540RReed5UORBRESqUr98OQt/Nb7iOgt/NZ61H3yQe9+l5a336Lsbby96t2idiuWt+/dj1aqP11/e+vTTuC+tUXH/gw9z4rFfraq8dV1dHQ8++CBnn3120fIJEyYwZsyYQsXM0oYDxItmrXPMKsktIiIbq6UPPETDypUV12lYuZIlDzxUcZ31mTt3LrNm/5PBA/YrWv50GmkoW976/j8wZP+B+cpbL1xIn947F9ZrTnnr7373u/zsZz+jXbviX7OvvvoqTzzxBAceeCD/9m//xowZM9bZ9vbbb482HlSSW0RENglr3nl3/SvlWK+cxvLWP/7+D9jqE58oem3qHx+Jlrcec8ll3Hjt1bneq9ry1g888ADbb789w4YNW+e1+vp6li5dylNPPcXPf/5zTjjhhKL3mz59Ol26dGG//fZbZ9uWPOa81HgQEZGqbLHDuqH3SuuFRbViD9p3KDzWrPqI4449llNOGsWpnzuKHlt0KDy6G/zx8cc54VMH4G/NKQwVLZS3vvVW9hw8ouxwyqLy1u8uSlII7TvQe5ttmDfnNVi6FF++mLq5b9Br6y5Fx1SpJPeTTz7J1KlT6du3L6NGjeKxxx7j1FNPBZKIwLHHHouZccABB9CuXTsWL27qA1KuaFb0mFWSW0RENlbbfOkLtOvSpeI67bp0YdsvfSH3vt2d0d/4Bv379+P73/veOq8/Pn06e++2G7177VhYVnV56yOP4PY/PMjHH6/mzblvMef1Nzlg2NCibSuV5L7qqquoq6tj7ty5TJkyhcMPP5xbb70VgK9+9auF0R+vvvoqq1evpkePpGHT0NDAnXfeyahRo/Ifs0pyi4jIxqRDt270uuDbFdfpdcG3ab/VVrn3/eSTT3LLrbfx2OOPM3jYMA4+4QQefeKJwut3PfIIxx99dNE2VZe33nsvjv/SFxhw5Of5/HEncv01Y2tW3vqss87ijTfeYL/99mPUqFFMnjy5kF74v//7P3r37s3uJSNJVJK7lWmopohs7mb02qXwfOgtlxeet//s6eVWz1WSe8F/jWfh9ROKOk9a5870OPF4+o69otnHHJvNccU/Xig8D4twNWvGx828AJZKcouISKvY6d+/zQ5nns6SBx9mzbuLaFizhm6f/jTtu1ZOacjGR40HEUK1E5EAACAASURBVBGpmfZbbUXPUScAZUp1yyZDjQcRkU3RyvcLT+t/cHLh+ZBjBhae+z23N60fSVtkERauCnXpvmXTf8K0QCiWIgjWj6UhwoJcVdsMUxXVUIdJERERyUWNBxEREclFcRoRkU3E2j/dXHgejp6wY08su7wasZEOq+rb8ezd9/P+wnfo3NDA/p/9DFtu9Qm69o4ddGSUQ940wmY+WmJDU+RBRERq4uGfj2fMLgO45ZzvMPWSK7njp2O59Khj+NP/TK563+stb73/IH54UTLnQWuXt543bx6HHXYY/fv3Z8CAAVx33XWF1370ox+x//77M3jwYI466igWLFgAwNNPP1043kGDBnHvvfcWtlm9ejXnnHMOe++9N/369ePuu+8u+76xY37mmWcYOHAge+65J//+7/9edirr3Nx9s3kMGzbMRUQ2VfXTJhceWZZnMWvWrLLLP3z+H4VHw/uL/MH//Imfy1bRx4OXjnWvX1P5UcE111zjJ510kn/xi18sLHvsscf8iCOO8FWrVrnXr/F3FsxfZz8vvPCC77bbbmX3eeGFF/pVV13l7u5XXXWV//CHP3R395dfftn3339/X7Vqlb/xxhu+++67e319feZrtmDBAn/mmWfc3X358uW+1157+csvv+zu7u+//35hveuuu87PPfdcd3dfsWKFr1mzprB9z549C///8Y9/7BdffLG7u69du9YXLVq0zntWOuYRI0b43/72N29oaPCjjz7aH3rooXW2L/c5AzM98vtUsR0RWUcs/B324C/SZesWPiKJCj6T2GfV/pMjC8/rv/XFsrvpcMODud42TFV89P5yHv3lhIrrPzruWg77zjfZslu3ouXhBFAh69q98Lxu4ds8+OCDXHzxxfziF78oLJ9www2MufBCOnVIZlIslLcOUhjrK2/9l7/8BUjKW3/mM59h3Lhx0fLWBx98cMVzbNSrV69C5cutttqK/v37M3/+fPbdd1+6Bee/YsWKwuySXYLpvVetWlVU1GrSpEn885//BKBdu3aF6axLz6XcMfft25fly5cXjv3000/nvvvu4/Of/3ymc4lR2kJERKry7P0P8/GKyiW5P16xgmfvur9Z+4+Wt57zKk/89a8cePAn+bfDDm+T5a3nzp3Lc889x4EHHlhYdvHFF9OnTx9uu+02Lr300sLy6dOnM2DAAAYOHMiNN95Ihw4dWLZsGZCkO4YOHcrxxx/PO++8s877xI55/vz59O7de53l1VLjQUREqvL+O4uyrbfw7dz7fuCBByuUt17L0mVLeepvT/LzcWM54aST21R568Yy4tdee21RxOGKK65g3rx5nHLKKVx//fWF5QceeCAvv/wyM2bM4KqrrmLVqlXU19dTV1fHpz71KZ599lkOPvhgfvCDH2Q+5lqdSymlLUQ2Z0Foe+3fmv4qDMPc0RSGtD1hWilMJQXLO1z9u5q8VZhW6L7Lrpm22XqHnvHJokr22Thi4smnnmLq1Kk89NBDrFq1iuXLl3Pqqady66230nvnnTn2q8esU966Z8+eAEy5/Y5M5a179eqVu7z19OnTOffccwG49NJL16msuWbNGo477jhOOeUUjj322LLvf/LJJ/PFL36Rn/70p0XL+/fvT9euXXnppZcYNmwYXbp04ZhjjgHg+OOPZ+LEievsK3bMvXv3pq6uruK5NIciDyIiUpWhx42kU9euFdfp1LUrQ4/7SsV1yqlY3nrkV3js8ceBSHnru+9usfLWlUpyuzujR4+mf//+fP/73y96bc6cOYXnU6dOpV+/fgC8+eab1NcnDau33nqL2bNn07dvX8yML3/5y4W+GX/+85/Zd999y55LuWPu1asXW221FU899RTuzs0331w4z2oo8iAiIlXZcutufO6i7zL1x/HKmZ+76LvrdJas1llnnslZZ5/NfoMG07HjFkyeNCkob/0EvXfeuWx56/POO4/hw4czZswYTjjhBCZOnMguu+zCnXfeCRSXt+7QoUPu8tZPPvkkt9xyCwMHDmTw4MEAXHnllXzhC19gzJgxzJ49m3bt2rHrrrty4403AvDXv/6VsWPHssUWW9CuXTtuuOGGQkNo3LhxnHbaaXz3u9+lZ8+e3HTTTUDS+Jg5cyaXXnppxWOeMGECX//61/noo4/4/Oc/X3VnSVBJbpHNWrQ8c5i2CNMZsZEXGm3RNgSfyYw9mmpYxD7b2OcWjsiY8+2f0b/fPkBxoauwboXtuhcAD13xcx696hd8vGJF4bVOXbbkc986ky+MbTqGokmcgjSGr1jWtM9YWe3YZFB59yNFVJJbRERaxRcuvpDDvjWaZ++eyvsL32Hrzu0Z8vkj2HKrT7T2oUmNqfEgIiI1s2W3bnzqzFOBeLVN2fhlajyYWXfgN8B+gANnufvfzewC4HygHnjQ3X9oZh2B/waGAw3Ad9z9L+l+hgG/BbYEHkpfczPrBNwMDAPeA05097npNmcA/5keyuXuPjldvhswBdgWeBY4zd1XN/9SiGx+Riz8V+F5OKoiLOEc1kUoolRF2xN8JuFnm0WYwpr4dlPY/5QvzWPuR6swjF2HNJXAfuu5ptoWu/YuP4LiX4s+KDzvsawp5bF4TdP6PbZo+jUUK7EdTiSVKQ0RTDON0hbr1ZzuC1lHW1wHPOLu/YBBwCtmdhgwEtjf3QcAV6frfiM9mIHAkcA1Ztb4PhOAc4C90sfR6fLRwFJ33xP4JTAOwMy2BS4BDgQOAC4xs23SbcYBv3T3vYCl6T5ERKSGlrz2Bivr1+JsPv3jNifuznvvvUfnzp1zbbfeyIOZdQMOBb6evtFqYLWZfRMY6+4fp8vfTTfZF/hz4zIzWwYMN7N5QDd3/3u635uBrwIPkzRCfpJufxdwvSVdZj8HTHP3Jek204CjzWwKcDjQ+OfR5HT7yvOjiohILo//ZCz8ZAzb7rk7H/1zdmH5e4ubJoZaGSwPhessbd/0t+qHaxvKLu/0yitl9+OrPiw8t86RCakamvbpS4NIxcr4vBKS6Ny5c9EslFlkSVvsDiwCbjKzQcAzwHeAvYFDzOwKYBXwA3efAfwDGJn+gu9DkoroQ5LCqAv2WwfsnD7fGZgH4O71ZvY+sF24vGSb7YBl7l5fZl9FzOwckmgHu+yyS7lVRDZfkfoHBKMqwnB2UwYxY1hcIzJaXnCNw3RTKEvdihGvv1h4PqR0P/98EP5ZvJ8Zn20a7vdhsGqY8sjrRl9e/oUM95EmM9uwsqQtOgBDgQnuPgRYAYxJl28DHARcCNyRRgsmkfwynwlcC/yNpE9EufkwG+NgsdfyLl93ofuv3X24uw9vnHVMREREmi9L46EOqHP36en/7yJpTNQB96SVO58miSz0cPd6d/+euw9295FAd2BOun4YF+kNLAjeow+AmXUAtgaWhMtLtlkMdE/XLd2XiIiItKD1pi3c/W0zm2dm+7j7bOAIYBbwOkm/g7+Y2d5AR2CxmXUhmXxqhZkdCdS7+ywAM/vAzA4CpgOnA79K32YqcAbwd+BrwGPpKIxHgSuDTpJHAf+RvvZ4uu6UdNvmlWsT2QycZ00z+0VDw0E4OFw/FN02JpaqKC3trZRGZmF4HsDvub3serlrWASfQbhtUSok+NxiaashwQRTz93blAoJJ6p69rR86a8s6RilKjasrPM8XADclg7DfAM4kyR9McnMXgJWA2ekv9S3Bx41swZgPnBasJ9v0jRU8+H0ATARuMXMXiOJOIwCcPclZnYZ0Fhn9dLGzpPARcAUM7sceC7dh4iIiLSwTI0Hd3+eZN6GUqeWWXcusE9kPzNJ5oooXb4KOD6yzSSSfhSly98gGb4pIiIiG5BqW4hsZsLRE7GQcZi2uP6bhxSeZ+m1Ly0rrDtRqihVkSVlFFknyz1SjdLUS6O8qYfwWujerL1KtS1UkltERERyUeNBREREclHaQkQq00RPbUrW0RahMJ1xXtem0e9hSiocGRFO9DR6x+6F5y2RwpC2S2kLERERqRk1HkRERCSXrPM8iMjmSqmKNi02oVNs5EWYhigqtx6kLcJ1Qi09CkM2Hoo8iIiISC5qPIiIiEguSluIiGxESkdX1Af/j6UqwsmUwhoTYRn2iW+fX3ge1jAJR3dkGdkhmwdFHkRERCQXNR5EREQkF6UtREQ2IrWs4RCbMCqcGCxMbayt2TvLxk6RBxEREclFjQcRERHJRWkL2ayFPcnzlgMW2VicP+GJwvPRwWRQ46ddX3b9cLKpcCKpZ0/7z8LzEQs3wu9LWKcFNAFaFRR5EBERkVzUeBAREZFclLaQzU5pSeOCLKWnN8fy1LFz3hyvxUYqHEkRpiHCSZ+K6lwEwlTFkGMGtsDRbUC6T2tGkQcRERHJRY0HERERyUVpC9ksFI2qCCa9iYbhYzb3sKdSFbmEJaxD1ZSzLk27ZRklFE4sFda5iKUwngtGZERF7oVM37UWFp5jrN6HVEeRBxEREclFjQcRERHJxdy9tY9hgxk+fLjPnDmztQ9DZOMSS+coBNy2BZ9bOOlTrDZGmG4IR1iEwnRLmBoIhRNShWJlvlticrYwXVRNimhzZ2bPuPvwcq8p8iAiIiK5qPEgIiIiuWi0hYhUpvRE21JpVFDwWYWpilCYbghHVUx8e1nh+egdu5ddPiJ473Db2ORRN66YV3i+IevIZE5VaPRQsynyICIiIrmo8SAiIiK5KG0hIlJJWwttlx5DOKoiZ0oiTDeMjqQhwuVr/3Z/2X2ODyabCtePlfZuEc0ZFdQWPs+NlCIPIiIikosaDyIiIpKL0hYiIlAU9p6xR1PYPgzPh6WtY5MtbWjhsQ695fKmF4L0QXjcYTojNqFTmHoIRy6cZ90Kz8NJn0IjFjaNpAgnaxra0mkL2aAUeRAREZFc1HgQERGRXJS2kE1LFT3jN+QkNtIGBfdLUfg/qPMQGzGwQe+dklEFI15vSjGEKYxQmJ4IUxhE0hYxsVRFTJuoKxGp8QHFn6e+8/ko8iAiIiK5qPEgIiIiuagkt2xaImmL1kxJbBbpkOC6hxMJbYjzLbq+nxxZfqWcKaxwVEEoDNuHEzKFaYHx065vOp5W/LzD6xIKjykcDRGOKsmbnmgTsqQsS1I+sUmsNtnvaU4qyS0iIiI1o8aDiIiI5KLRFrLxi4Urg+XRMOSGrlvQ1uokVCNSS+DZYHTC0FuallcdCo70mo+NgIiVpM4yuVMYto+F/2OKUiet+HmH1zs8h/B5OBoiLLcdpjPCOheZJsaq1Tnn3E+WdNl5XfsU/T+s86FJrPJR5EFERERyUeNBREREclHaQjZ+sZBmNaHOlgo3t+FURTh6oCg8HRlJ4ffcXn79FhK+d5iqCI+D4DOsVQnoWAg8VhciDI23xKiF8HOC+HmGx503ZRROkhWmoUbc0LRObBRRNd+pTCNnIsL1w/0U3R9SM4o8iIiISC5qPIiIiEguajyIiIhILpphUqScWvZ5aMPDM2MzKcaEQ9tiMxJG+060lMiQ0ZaeYTTc57ePPL/sOrXq8xC+V9gHAYqHUoZi1z78fIpmVQz7DFQxS2iWmS2r+U6Ew0jDoabhvRzep0VDTa/+XfHO2tj3sa3RDJMiIiJSM5kaD2bW3czuMrN/mtkrZnZwuvwCM5ttZi+b2c/SZVuY2WQzezFd9z+C/fwlXf/59LF9uryTmd1uZq+Z2XQz6xtsc4aZzUkfZwTLd0vXnZNu27E2l0REREQqyZS2MLPJwBPu/pv0l3QXYAhwMfBFd//YzLZ393fN7GTgK+4+ysy6ALOAz7j7XDP7C/ADd59Zsv9vAfu7+3lmNgo4xt1PNLNtgZnAcMCBZ4Bh7r7UzO4A7nH3KWZ2I/APd59Q6TyUtpBKahluL9pXGCptY2HSWKg3TEmEhZ5i4flYMahoiLw5qYNwBsQ9mkLRI15/sfz6bexaZxKZRTO8h8JzLxWG8fO+x3P3Nl3HcKjmhh6Sm0eW+3ejLPLVRlSVtjCzbsChwEQAd1/t7suAbwJj3f3jdPm76SYOdDWzDsCWwGpgfZ/eSGBy+vwu4AgzM+BzwDR3X+LuS4FpwNHpa4en65Ju+9X1nYuIiIhUL0vaYndgEXCTmT1nZr8xs67A3sAhaergf81sRLr+XcAKYCHwL+Bqd18S7O+mNGXxo7QRALAzMA/A3euB94HtwuWpunTZdsCydN1wuYiIiLSwLDNMdgCGAhe4+3Qzuw4Yky7fBjgIGAHcYWa7AwcAa4Gd0tefMLM/ufsbwCnuPt/MtgLuBk4Dbgas9E1JIhh5l6/DzM4BzgHYZZddyq0iUnNtOVURE4aqCXr0hymG0TsW9/QvKzL6oepRDsF1DHvQx4pk1WpURW7VjCQIUhKxURSZUhOl+w1GKIQh/ZgR4YyRwfJMIylaWuT+Uqpiw8oSeagD6tx9evr/u0gaE3UkfQ7c3Z8GGoAewMnAI+6+Jk1lPEnSZwF3n5/++wHwO5KGRuN79AFI0x1bA0vC5anewAJgMdA9XTdcvg53/7W7D3f34T179sxwuiIiIlLJehsP7v42MM/M9kkXHUHSCfI+kn4HmNneQEeSX+r/Ag63RFeSyMQ/zayDmfVI198C+BLwUrrPqUDjSIqvAY950pPzUeAoM9vGzLYBjgIeTV97PF2XdNumgckiIiLSYrIWxroAuC0dafEGcCZJv4ZJZvYSSafIM9zdzWw8cBNJw8CAm9z9hbQh8WjacGgP/An4n3T/E4FbzOw1kojDKAB3X2JmlwEz0vUuDfpPXARMMbPLgefSfYg0W1GBoWondmrpVEWNJp4Ke6iHkwSNWFh+xESY2sgSqm6p3vmxwl0bNEUUvG9YDCs28iQmTCmEwvuxUjojHCURS2mE24wP9htOOBWG/UuLb633+IJziE5aVTpBU6O8n1mWgnctMCmYFMvUeHD350lTDyVOLbPuh8DxZZavAIZF9r+q3Dbpa5OASWWWv0FT2kNEREQ2EM0wKSIiIrmotoVIGRu8PkNORWHZSE2CUCx0m6UX/vXfPKTw/PwJTxSex3q0b44h4/B+yXKNYtsWpSCCya/CESWVZBnlE34+4QRQ4XvHUg8xRSm/QLj/mGpGJsUmiWrOiJSNUgun7FTbQkRERGpGjQcRERHJJetoC5FNXhjO3dClpGMTHYXC9ET4PAxXhmHi2H5CYXh3YqQ8dxjOziKsf3Gjbx5pi6gMYeXY9Q1TUGFovzQ1FX7msbRVLH0Q3iMjgns+lraL1WyJva9FRnYUTUiWIdweG/2x2U8G1YoT0CnyICIiIrmo8SAiIiK5KG0hktoQIwNitQEyCUKUWUaD1Op8stRCCIW93jcXRZ/BhKb0T1Hp8CBFFJsYKhzlUBzmb1qnNAURK9cdjtYI61OE2xelv4IUS+yeKirbfWy+FEl4brF7MzZSJ0x/FO2/tSYI29Da4Hkq8iAiIiK5qPEgIiIiuShtIbIBRVMJQfg4HHlRFCYOQtexkRRFPeCrKJ8cph7CcHM4AVIYSg1724dpjhHrfaeNTKQcdBhKzjICIEwphKmG2ERNRROBlewrvF/yTu4Ufm7hfRFOvhSeT2zypdhoiCwjMkJ5R+rE0jRtJbRfM23wfBR5EBERkVzUeBAREZFclLYQaQuCsGSYkpgYpAmGBL3Mw97wsToXURl6boeph9FBKL1oJEVk2+ZM3BOOPmhrdQnCkHx0Aq8MqaBYTYlQmHaIhfwr1YsoSntQPv0Vvke4rxnB6I68qplULTbyJLxPw+tbHxzzJp2qaOMUeRAREZFc1HgQERGRXFSSW6S5NsDELZlKDsdGAATCHvmxEsjR9wpCw+d17VN4Hpbqjk5UlDEFEYb0w8mRWiuFkSUVFBupkEUsFdI+MuomVJo6iaUxYjUzwmtaNKpixbymlVrgfg7TE0XpmZwpj82x5HtrUUluERERqRk1HkRERCQXpS1E8tgAqYrY5E7hBDqxlEF4TGGYOBw9EQtPx9YP3yucJCocVREec/FEPxtPyeQsoyGKrnUo570Qm1QpWsMhkLdEOmSrq9ESKaIwLRK7j2IpMo2eaH1KW4iIiEjNqPEgIiIiuWiSKJE8wlBqpVEOVYRcw7B/rLx1Ufrg6vL7CXu03xhJbYSG3nJ503+CEQ+xEHPeYw6VpmaKajeEoxjCkQ4tEcYOPsMsaYIRN9TmGGKTTVUrlm4K0yQbcgTL+GnXF57Hrm9RLRSlKjYaijyIiIhILmo8iIiISC5KW4g0Vy1DrEH4vGgkRTCBzpBYjYXIceSdfCeccGfIMU0h5uuDdWLh9lgJ71BsFAlQdA6xegt5w+3RSYlylokuGgHQBpReh9joiTDdUjTa5urIiKEWGEkUfpbPRWqktLXrK9ko8iAiIiK5qPEgIiIiuShtIRu9Fp/0JkPtiHXCvHlDwME6Yboh7CUfmwyqpRWFv4tGXjSlFLLUKqg0ciJMaYTnOSIYSVJ0LWIpmeC6h8cUpluKam8EIfP6CqWucx1DBmGKKJZeKRrlEdYXCSZeguIUQDgBWHiembTAPRVeoxE3ZNgg+PzC44+lwmLvFYp+h0ATUVVBkQcRERHJRY0HERERyUVpC9nohemJWM2AllZp0qNqerFXExqPyVLSuOh9720Kq8dqVYSh92h4usoQcTUTK4VlvsMQeJbRFkUppSDsHbvX8n5mRamToAx3lnoUpWIpmSHBfqOlxDdA3ZZyYuXYYymcaAojOP7wHLOMTJL8FHkQERGRXNR4EBERkVxUkls2GlnC7S39vjGl8/ZXVbq5hcPHWUYMZAmThxMPhWpZhrs+MjFW9POPXLss5bZjspTJznRsEbFy5jGltUPCzyGsJbEhvyPVCEePxOqiFNVdCRSlBwOxdFSYFoENW+djY6SS3CIiIlIzajyIiIhILhptIRuN2AiGaO/xKhSlSLKUiy4Nn1aTbmjh0tNFPdcpn8KI9fQPl08smayo3HuFvf/D8HppeLpW9Q3CzydMMcRGaoTppdikROE9tTbYNhYyzytLqiKsd7LOuQSh+Lz3f3SCtfDzqNX9GNwX53Xts97Vs0wMFY5OCUvHx4TXUaqjyIOIiIjkosaDiIiI5KLRFtLmxHrGZ5kkKBa2DXvtxyacifaYzzD6oeIkUaEahYBjI0+iIylqlOaJ9Ywvqr0QhKRjvedLFZWMDuozhLIcd3TipmpGv2TZfxWTecXuzVDWUQF5a2/ERhK1xEiN8N4pKjsffDaxeyeWwsiSqii6TzW6IheNthAREZGaUeNBREREclHaQtq2INwe9qzOW546S82LMMwbG20Rfa8sZbsrbR8TSTeEYuWdwxBtrH5AqJqQbvi+scmjQpXSGbEQdZaJmPJOJJYlnZVX7N6MHVvs2sXC7aUjJMIJlNryxFClpcQbhRNbhSNPYmmL8PqG34nwvo7dg+F7Qdu+Xm2B0hYiIiJSM2o8iIiISC5KW0ibUE0v9iyh6qJe6JGwZ15hOiNMqUDLlNKuasREhvRHOLIl7MWepbd6eH2z9IAvVase8bGROtWkJKpJkeXdNjppU4VrkiVlVDS6oSXuzVDOyaBiwnui0gRjjWIlvCul0WKjPlS6O6G0hYiIiNSMGg8iIiKSi2pbSOsJwptFPemrqBkQS2HERlKEwjB3LLWRaUKmGspyrPWR8HwY6g7Dx2HJ7Nj+Yz3gR0SOMzz/0feWD6OHIeLSyZBaevKeWG/9TJN+RUZMhKKjPyK1NmL3S5aRE6WpjSwh/aL74oayqxSrpix8sH4s9RAT3ptFqcwMqbDwcx16bNN1HxL5fpQKa5toMqn1U+RBREREclHjQURERHLJNNrCzLoDvwH2Axw4y93/bmYXAOcD9cCD7v5DM9siXXcoSVrkZne/Kt3PMOC3wJbAQ8B33N3NrBNwMzAMeA840d3nptucATTO/nG5u09Ol+8GTAG2BZ4FTnP31ZXOQ6MtWl/N5tKPlH3OUko4S2/4osmWmlGeONbrP3d6IzL5VOnojnL7j17rSFooVt64ViMhYhNYtZTY+ccmE8o6KVM5eSfkKqrlEYTqY2I1RUr3FaaGws8wTENlmaAqJjZqI7af8H2zjH6JjZKI1TvJ8j2oJDzuotRWS49I2UjUYrTFdcAj7t4PGAS8YmaHASOB/d19AHB1uu7xQCd3H0jSGDjXzPqmr00AzgH2Sh9Hp8tHA0vdfU/gl8C49MC3BS4BDgQOAC4xs23SbcYBv3T3vYCl6T5ERESkha238WBm3YBDgYkA7r7a3ZcB3wTGuvvH6fJ3000c6GpmHUgiDKuB5WbWC+jm7n/3JNxxM/DVdJuRwOT0+V3AEWZmwOeAae6+xN2XAtOAo9PXDk/XJd22cV8iIiLSgrKMttgdWATcZGaDgGeA7wB7A4eY2RXAKuAH7j6D5Bf6SGAh0AX4nrsvMbPhQF2w3zpg5/T5zsA8AHevN7P3ge3C5SXbbAcsc/f6MvsqYmbnkEQ72GWX9YcfZeMQhiuLJv3JkKqI7SdMYcT2GYY5w57x4YgEiPd6z9TTPRCmZGLh3aL9X92U5ohNYrU2WD8MJYfCMHE1E2nF6hlMLFmepVd+3vRJeP7NKRNe7niyTLZUei80CkP4E4N18tbjKL0+pfUaCoK0RZZURezcws8wTIXcGNzLsfRPeC3Ce2rihGD0T3gdg5E64fEUiXzHw2O7nvwjoYq+U/du2BTbxihL2qIDSf+FCe4+BFgBjEmXbwMcBFwI3JFGBA4g+fm0E7Ab8P/MbHfAyuy7scNF7LW8y9dd6P5rdx/u7sN79uxZ/gxFREQksyyNhzqgzt2np/+/i6QxUQfc44mngQagB3AySf+INWkq40mgMerQO9hvb2BB8B59ANJ0x9bAknB5yTaLge7puqX7EhERkRa03rSFu79tZvPMbB93nw0cAcwCXifpd/AXM9sb6EjyS/1fwOFmditJ2uIg4Fp3X2hmH5jZQcB04HTgV+nbTAXOAP4OgMulNgAAIABJREFUfA14LB2F8ShwZdBJ8ijgP9LXHk/XnZJu2/zY6qaqmoleWkim8taBWLohb3nmMCQZphuKet5PKD+RUngMsboIpaHwLCmGvDLVjIhMdBOdGCn4PMIUCfc2PS8eGbBhSxjH0hmxyapCYaqm0gRVjcL7IhQtGX1v+TRoLNwepohiox9io3RCpfda7P6/cUX5UTVZRlWEx5ElXRQ7nyHBdye8v8J9xkbh1Ee2DT+n8L6+8YbIz5MKPwfzjqqRJllnmLwAuM3MOgJvAGeSpC8mmdlLJJ0iz0h/qY8HbgJeIkkv3OTuL6T7+SZNQzUfTh+QdMa8xcxeI4k4jAJI+0pcBsxI17vU3Zekzy8CppjZ5cBz6T5ERESkhWVqPLj78ySph1Knlln3Q5LhmuX2M5NkrojS5asqbDMJmFRm+Rsk/StERERkA1JJbmkSmZAIKAr3tXhNh9gEUDl7PUcngAlD+KEMIzVitSMqXbuiER1VXK9YL/lYaL9o8qFI6DYWJi/quR4J+ddqwijIFkqvZrKq2KiPTPuM3I8x4aRi4WdfdO1iE49FPqdKE2xlmnwrQ5nsoomYcqYUY+cWXq/YvRnuJ0xHhvdjNE1Z6WdWowrnEhs9lXvSuk2USnKLiIhIzajxICIiIrmo8SAiIiK5ZB1tIRu72HClYHnpTIKxgk416+cQee9Y4aa8M/FlyotnyJmG5zuE8kNHK+VVa3W9wvMpnaGxrAznlnU4YDnV9H0p/WzCoZfRWSlzDtUMRWf8DPLzsf4fsWJN4fDEcMhvrJ9DbNbS8H2LZmQM+gVUupdjr8W+L7FiVdUM5w77C8T6VITCzzjWNyecCbVI5L6O/gzJOFQzpsX7eG2kFHkQERGRXNR4EBERkVw0VHNzEQndxYZJQQsNV6qQJim8bxByLApFxo4ny0yaWUKdtTzf2HDT2BC9DGJDNWPCoXexIXOxokptZahaluGAtZo9NQzzx0RnSM2QFswyFHiDhMVrNfNssJ9YqiYUu2dj6Z8wFRJLr5X+zCqnWfdy5PsbO4ei1MsmREM1RUREpGbUeBAREZFclLbYFMR61bdQMaxqwqx5R0y0iBa4XqUh71hv/aKZLnNeu0oppnLvG+4/y3WPrRMrGBSbMTC2fqWe7bERFqFY8aWWFr12kRRcGHrPPeKnjRSwi4ndg2E6J0xhxNYJR2RkmeUzllKKjswKPo+i70Tp7LKRWTxj6Ykso5CqmXm1rVHaQkRERGpGjQcRERHJRWmLTUBRWDXWG3xzF+kZHgurVjMSArIVroqF4WOTBoWqKRJVJEPIPEwpxN43lnYIi2qFaZTmXK9weTXpn7xihaeypHNCscmQ8obw26TI+cQ+/5jw3i9K/2QoNlY0CVdEpWOIfZ7hfmMjSUIb1ee2HkpbiIiISM2o8SAiIiK5KG2xKYiEnmPpjDDst0nM1Z6zt3qt0jyxECvEQ6jnT3ii8DxvL/NYCiMMh8fSB1WFwzNMgBROpBNqTgg3liaJje5oCeGoglAsbF3VpFUbeLRUXkWjq0pHKzSK3BdZftZU+h6tTyytFU4wlfU7nmWUUGy0SWyU08ZOaQsRERGpGTUeREREJBeV5N5YZQjVF80HH5lMpS1OUBMLGccmPQpFJ4QJzi1TGDNDKDk2336pMLw5PhzdEZnoKDz/MM0RitWkiKmqB3jsGgXLYyMvRmRIqUH8WkRreLTwfZulV314bBOrmPSordQRCUVD+FcH/8lwX+RNVcTuo1j6KiZrqiL2s6bosw2Oo2iEVPgesZ+vmzBFHkRERCQXNR5EREQkF6UtNlYZQrWxMH9RXYRqQ741Ch/H6mWExx0+j9WOKHJv0wiAcPSDZUgdxELsmd6XeMokmg6JiJXVzhLqzpvaqJXYSIjKI0fKjyQJFdUVaIFURXh81UzyFd7LsfB88ciRtpe2CBWNHKrRdQ8//9h9mncUUTjCIutxxiaZGh0sDz+r4pLhlF2+KU0SVYkiDyIiIpKLGg8iIiKSiyaJ2lw0J70Qma8+lKkUc2QCnSyljqNaeLKi8FyKwqGBWKgaNq2JYkJFdR6CzzW87mGYN1bboFIJ7thohViNiVrJMjFU7H1jI2RiIfkwbN/SE141R2ykQ61qisTSOeH+Y9cxNpFUtccW+/zzHkeoqsnD2gBNEiUiIiI1o8aDiIiI5KK0xSYmywRLebeF+GRFoTBE/Vykt3KWEG0s/RGmD8JRD0WjRyI21lLlsVEorXUM4X0QGwkSUxTyL0lNZUmLZdESky/F0iVZRlJkmdyoraQtYqm6WFohywRQofDzqJS2atQS1yVruiyWPomtn2XiqrbyOeehtIWIiIjUjBoPIiIikosmidoERCeACjVjtEWYesgyYUssvBemM2IlbcOQZlGaIwiNr/3b/YXnsUmYYiHTMKzeViZxiZUGj6UAhh5bfrKqlk5nhPu/8YbgheCeyhK2HVKhZkcsBJylTHKouAd8y07cE35mQ45p+jxGZ6iL0RZD2OF1Hxosz3LtYj8Tspxn+POh6OdXcH+F3/3c6aiME7NFUzIT8tXYyJvO2Fgp8iAiIiK5qPEgIiIiuSht0UZU00u8Vr3KS1MeE4PQ8sQMc8vHJgeKvUdRSevgHIpCnUHIMTpiIlwnci1GLMx3jWL1LGKK0iilxxdRzefWJiahCs4x/MwyjRCZUHw/haHx2L1WFFaPfF9iJZ1j64dioyrCcHZsnaLzvLdpndjIi7YiNhlU+ypKTIf7iV2vLHVXwlRFNep/cHLZYysVptKuZ/2lujONpAnSrpsaRR5EREQkFzUeREREJBelLdqIotB4S09oFOnFXDpXeyz1EIb3wslkhpZbuYIsk8xkCu23wDWKTUKVJYVRS2HIvFbpqVqJTSqUtYZFNYpLIJe/LmFYubiWRNP6sdE54XGHKZkRN5Rbu1hRPYM2OMonFAu9j6jiO5UlfRlq6fs6TClOnFA8GVl4/kUphjAlF3zm1Ywq2dQo8iAiIiK5qPEgIiIiuShtsQEU9T4v7ZWfKhrp0Eq1FyqFVcMQdRjqC8PHsdoToSyjBIpC0q83pXCqmigmg1htj9jc/i2lLYa3y4nVqqhUIyWvakYrhOmDiZHaGcWTk51fdp1QLGwdhufDtF5bHGERkyXFEIpN4BWKfd9bOlURqztS8fPI8HN3c0xPxCjyICIiIrmo8SAiIiK5qCR3CwlDt89F5rov6pUdmfRog6YwSuaAD9MEsVoKeUPUeSc3itV/KNLC1yiadqpQU6OtjYzYkGLh7EopiNjkTjGZwsfB/RyW/I5tm6UnfXhueesWtJWQd3ivFo9CyXd8se9FOClTm5jMLKDvaT4qyS0iIiI1o8aDiIiI5KK0RQuJhfPDFEas534sNN7WxSbcqVlosI2kc8odQ+nnnTtcGzm3amqetAXV3hNZUhhhnYQs+43V3sg7AVDu9Epr3b8lYp9JmMKIjR4JZUlzhNc6HFG2Md7LmyOlLURERKRm1HgQERGRXJS2aCFZRiHkDW0XlbeNjdRoTW0kLLs+sbB1ljLcRaHXYPRH2MO8dL3YSJVQbFTJjBaojRArw9wWJ6eKTfYTio3iyHs+4X0RC9WHKZLweLKMvGgroy3y1hsJjzuW8shS80Vpi42P0hYiIiJSM2o8iIiISC6Z0hZm1h34DbAf4MBZ7v53M7sAOB+oBx509x+a2SnAhcHm+wND3f15M/sL0Av4KH3tKHd/18w6ATcDw4D3gBPdfW763mcAjfHBy919crp8N2AKsC3wLHCau6+udB5tYZKoqkLDYf2ANjwRS6uKpE5ik/tUE0qulOYo+kxikxUFJYDD9EQoHJFTq895UysrHPtsq0lhxCZSyqtFrmnJ6J/YvRNLbWZJz4TynkOLj7qSDaYWaYvrgEfcvR8wCHjFzA4DRgL7u/sA4GoAd7/N3Qe7+2DgNGCuuz8f7OuUxtfd/d102WhgqbvvCfwSGJce+LbAJcCBwAHAJWa2TbrNOOCX7r4XsDTdh4iIiLSw9TYezKwbcCgwEcDdV7v7MuCbwFh3/zhd/m6ZzU8Cfp/hOEYCk9PndwFHmJkBnwOmufsSd18KTAOOTl87PF2XdNuvZngfERERqVKWkty7A4uAm8xsEPAM8B1gb+AQM7sCWAX8wN1nlGx7IknDIHSTma0F7iZJQziwMzAPwN3rzex9YLtweaouXbYdsMzd60uWr8PMzgHOAdhll/WXkK2VcKKnocfeX2HN7MLwZFj+WgKRUR5h2Hpizt7m0UmiAkUl1UuEKaZw8p1weRhiDj/nopTXDes/1Lw21lRFKPwuTAxSDGEKY0TOfVaTqsiiRdKalIz6iIzUCdNfo+9df8oni2gNmjY86kqqkyVt0QEYCkxw9yHACmBMunwb4CCSPg53pBEBAMzsQGClu78U7OsUdx8IHJI+Tmtcvcz7ejOWr7vQ/dfuPtzdh/fs2TN+liIiIpJJlsZDHVDn7tPT/99F0pioA+7xxNNAA9Aj2G4UJSkLd5+f/vsB8DuSfgyN79EHwMw6AFsDS8Llqd7AAmAx0D1dN1wuIiIiLWy9aQt3f9v+//buP1bSqjzg+PeBFdKoKCBYKihoBatWWETQNm6ppEKtKf5Cl4olYkpqkNhYqmz/0AatFaXBn0CsrAUrbCnR1GoVMbJqAlVEMa6SFVxpdyv+wAW1KJLF0z/mvdwzs3Pufc/O75nvJ9nw3nfeeWfOnLm7h+c55zwR2yPiqJTSVuAk4NvAd+nMO9gcEUcC+9D5R52I2As4jc5cCZpza4BHp5TujoiHAS8EPtc8/AngTOAm4GXA51NKKSKuA96eTZJ8PrCheeyG5tpNzXOHkxsYQNcmTlmYcFizjKdxE59Zl688KIXw8/TCqDa6KW0GlX+nSt+vWgOlKqZwI7A2KYY2n12plHitNp9vaZVDKb3SZgMz6F7B06Z/8s9ibYvN0/LvfOnzWvvi1TdF0+xrM+cB4FzgoxGxD7ANeDWd9MXGiNgCPACcmZbXfa6jE63Ylt1jX+C6ZuCwN52Bwz81j10OfCQi7qATcVgPkFLaGRFvBZbmUlyQUtrZHL8J2BQRbwO+3txDkiSNWKvBQ7PUst9azzMK12+mMxciP3cfnX0c+l1/P51IRb/HNgIb+5zfxnLaQ5IkjYm1LQY1pDBuKaw6zlnMvWHIrvBu1s4Hb+yfIZqVTWBKtS1KBi2LXarP0LUBVLY6J/98i3UxBvjetWl/7Wc0jYa1GVZ16e3Ke87qipdS2mJYm7Bp8qxtIUmShsbBgyRJqmLaYpJKmw9lYehBQ+arvW4pRN5WKdyem5V0RptUwDD7o1QDIFd6jVLIeJBVGKUZ/aUy5NOy2iKXp1ty+eqG2lD6qGtBzJL8e5dvzlUql16qr6HZYNpCkiQNjYMHSZJUpe0+DxqFFmG8kYT8s9cthadXeqxrxQCrb6A0ktRLC6VVDsVSxS3KnOfvf5grEkqrMCh8Xnkb8jLfxw5QDrnrmvx4CjeGyuX9kH8HS5tHtVnpkH938lRFXpsk7/NhbTA1LUoptXwlxdrs7wc3sFs8Rh4kSVIVBw+SJKmKaYtF1KLENPSWg14O0T7rrtVD4MW9+CvTFm1SHm028emS15HIwv/5ypFavbP8S2WWS6WYc3novViTIUsfFGty5O9pkHTRFKYqcqNIH+TftcvS8nHXd+3S5ePaEtZTo7DyqrR6Ik/baLEZeZAkSVUcPEiSpCqmLaZFIXxYnCXf5vrKe/ZuElU7gzoPGecrBnpXcaymNHu+Vh5Kzje0yZU2yWpTenvF91ZISeQh9rX0T2Hkr3f5D5ZXDKwtpDyc6b4sXw0wiHmo7dFG/v3PV6fkvzulGizTns7SaBl5kCRJVRw8SJKkKqYtplxptn3XhkZ5KLFF6ez8mhVDspW1HrrCmwOUem4VGm25YmRJafZ4m/ec69pg6L7tDx33ftb5ffPnfOD696/6PvLPNF9Jkd+nzez+eQ63l5Q+r1r5Z5d/Z0v379rYawqVfu+6U3LLaYs8/ZN/Y4v1a4ZU82RSG8qpnpEHSZJUxcGDJEmqYknuaZGF4W9+0uoh0K4Z9m1SFZlS3YndriuEIosbQGUGSVsUtUhV5HUecqXNbYqrKgr9UVrZ0BsiH6gs8xTXkihuWjWFSmmLPHXUpuZJG/nnkof8S2mOPOU1jj5uU8+jlBYrtWcUK1JMW0wXS3JLkqShcfAgSZKqOHiQJElVnPMwLQp57mKucpC8eGHuQO98idIcgJJ8+WirHRAr29CVD83nYxTkcxXypXT5e6tdqll6P6WloL2v0TXHImtDaZ5LvjPmNOSA2+TOiyY4l6PNss18LkTxu9/ifbd5rYHmxIxDizlYQ/tulv5umbL5PovIOQ+SJGloHDxIkqQqpi2mXL4ELDfQMrm2YcJSYa0WKYNcq+cOKUS5q1A8KpcvycuXdrYJJef3X0mb126ze2arImljNKzlhisuEa5sZ94nr7v0S1XPrZUvYWzzOzgPBbaKfwdl3+VBvqeztPx30Zi2kCRJQ+PgQZIkVTFtMeVGEtLbk1nvtbtYZumJYjGdIYXh24SGa2fAtwnPr5QiKfVV6TmDrPoYp9IuhPnM+3zlSb47YVv5bqClAlX5Z7cnr7GnSu9tUZR2l+0qZpfJV2AVd3PV1DJtIUmShsbBgyRJqmLaYoZ0hefzMOE4NlNpkbaoXYUxrPedp3ZKm0HlSmHurpB09vmWQq8r6QrpZuHtQYqKTYM2q3/aFFhaSZsVDaVCVCWlgk61qzPy95ZzlUDBFK8c0upMW0iSpKFx8CBJkqqsmfQbUHuTDGdXr5ioXNExyGY6ecg4Twusza5pE57uqsdxUf9r2tY8uDlbcXDsS7KaHPnzZzB0m6eFcqV0RsluG3K1+L7kfTvICovSd6FNiiU/n7ehq+7KDPbryGR96ecyX4w8SJKkKg4eJElSFVdbzIHSBkilMOy0b3QzSAqjzWZQudJnVCrPXKrJsFtIdoLlp0eqMHs+3xgqr3mQl3NeKdVQW6K6VBq8Tcnw0oqXNqtzpr6UtjRErraQJElD4+BBkiRVMW0xB9pselMKz68Uhi3OIJ+RkHy+AiCvvXDOH72u7/V5qqIrDJ+t5piHEsvjVFqF0ZvCqC11rcnw+79YTFtIkqShcfAgSZKqmLaYM6UURh6Sz8P2K22MUwrj5ymAUg2HklIqYVhpkbz9eR2KUh2J3NjrhSyAtjUoXMUgTR/TFpIkaWgcPEiSpCqmLeZYKWTctkxyaaOkYlqhRboh38Qn36yqtEHPSGbeZ++zVG7bffiHI1+xU1rlAq62eMiMrGTSYjBtIUmShsbBgyRJqmLaYo7lIeN8tURJb7nl2k1ghlWWuHX9iAr5KoxcaUWGG+BoEiztrWli2kKSJA2NgwdJklSlVdoiIh4NfAh4OpCAs1JKN0XEucDrgF3Ap1JKb4yIVwJ/kz39GcCxKaVbI+KZwD8DvwH8J/D6lFKKiH2BK4FnAj8BXpFSurN57TOBpZj721JKVzTnjwA2AQcAXwNelVJ6YKV2LFraIpeH7UsrG3qVajrkz89THSMJ9Tv7XJImYhhpi/cAn0kpPQU4GrgtIv4QOBV4RkrpacBFACmlj6aUjkkpHQO8CrgzpXRrc59LgbOBJzd/TmnOvwa4J6X028DFwIXNGz8AeAtwAnA88JaI2L95zoXAxSmlJwP3NPeQJEkjturgISL2A9YBlwOklB5IKd0LvBZ4R0rpV835H/V5+unA1c19DgH2SyndlDrhjiuBFzXXnQpc0RxfC5wUEQGcDFyfUtqZUroHuB44pXnsec21NM9dupckSRqhNS2ueSLwY+DDEXE0cAvweuBI4LkR8ffA/cB5KaWbe577CjoDA4DHATuyx3Y055Ye2w6QUtoVET8FDszP9zznQODelNKuPvfqEhFn04l28PjH9y8PvAi66jZ8/Hf7XpPXmuiVpyq6NpbK6mdcdknhyYOkHkxVSNLUaZO2WAMcC1yaUloL3Aec35zfH3g2nTkO1zQRAQAi4gTgFymlLUun+tw7rfJY7fndT6b0wZTScSml4w466KB+l0iSpAptBg87gB0ppS83P19LZzCxA/hY6vgK8GvgMdnz1tOkLLL7HJr9fCjw/eyxwwAiYg3wKGBnfr7nOXcDj26u7b2XJEkaoVXTFimlH0TE9og4KqW0FTgJ+DbwXTrzDjZHxJHAPnT+USci9gJOozNXYuk+d0XEzyPi2cCXgT8H3tc8/AngTOAm4GXA55tVGNcBb88mST4f2NA8dkNz7abmuf8+yAcx7x68cfnjKZbC7lGqjVEq1V1k6kGS5kqbOQ8A5wIfjYh9gG3Aq+mkLzZGxBbgAeDMtLzucx2daMW2nvu8luWlmp9u/kBnMuZHIuIOOhGH9QAppZ0R8VZgaS7FBSmlnc3xm4BNEfE24OvNPSRJ0oi1Gjw0Sy37rfU8o3D9ZjpzIXrPf5XOXhG95++nE6nod6+NwMY+57fRWb4pSZLGqG3kQTOuzT75vTUl8g2g3p+dL5VWfpYbOknSQnB7akmSVMXBgyRJqmLaQg/ZLbWR/ZzXtrjsvmzfLtMTkrRwjDxIkqQqDh4kSVIV0xZqpVW5bVdbSNJCMPIgSZKqOHiQJElVTFtoeExVSNJCMPIgSZKqOHiQJElVTFtoMNkKi13n/dlDx2suumr5GtMZkjRXjDxIkqQqDh4kSVIV0xYaTJaSaLWRlCRp5hl5kCRJVRw8SJKkKg4eJElSFQcPkiSpioMHSZJUxcGDJEmq4uBBkiRVcfAgSZKqOHiQJElVHDxIkqQqDh4kSVIVBw+SJKmKgwdJklTFwYMkSari4EGSJFVx8CBJkqpESmnS72FsIuLHwH83Pz4GuHuCb2fcFq29YJsXxaK1edHaC7Z5Up6QUjqo3wMLNXjIRcRXU0rHTfp9jMuitRds86JYtDYvWnvBNk8j0xaSJKmKgwdJklRlkQcPH5z0GxizRWsv2OZFsWhtXrT2gm2eOgs750GSJO2ZRY48SJKkPeDgQZIkVZnJwUNEHBYRN0TEbRHxrYh4fZ9rToyIn0bErc2fN/c8vndEfD0iPpmd+9fs+jsj4tbm/OER8cvssctG38rd2jOqNh8TEf/VXP/ViDg+e2xDRNwREVsj4uTRtnC3toy1vXPex0dHxE0R8c2I+I+I2C97bGJ93Lz+WNs8D/3c/N30zaXvcHb+gIi4PiJub/67f/bYzP4u17Z3zvv4tOZ+v46I43ruN94+TinN3B/gEODY5viRwHeAp/ZccyLwyRXu8QbgqtI1wD8Cb26ODwe2zGObgc8Cf9wcvwDY3Bw/FfgGsC9wBPBdYO85bu889/HNwB80x2cBb52GPp5Qm2e+n4E7gcf0Of9O4Pzm+Hzgwmno5wm0d577+HeAo4DNwHHZ+bH38UxGHlJKd6WUvtYc/xy4DXhc2+dHxKHAnwAfKjwewMuBqwd/t8MxwjYnYOn/RB8FfL85PhXYlFL6VUrpe8AdwPGMyQTaO3EjbPNRwBeb4+uBlzbHE+1jmEibJ27QNq/gVOCK5vgK4EXZ+Zn9XV5Bqb0TN6o2p5RuSylt7fPQ2Pt4JgcPuYg4HFgLfLnPw8+JiG9ExKcj4mnZ+XcDbwR+Xbjtc4EfppRuz84d0YRGvxARzx3CW99jQ27zXwHviojtwEXAhub844Dt2XU7GM4vfLUxtRfmt4+3AH/aHJ8GHNYcT00fw9jaDLPfzwn4bETcEhFnZ+cfm1K6Czr/eAEHN+enpp/H1F6Y3z4uGX8fTzK0M4TQ0COAW4CX9HlsP+ARzfELgNub4xcCl6wUNgIuBf46+3lf4MDm+JlNJ+03D20G3gu8tDl+OfC55vgDwBnZdZcvXTen7Z3nPn4KnXTNLcBbgJ9MUx+Puc0z3c/Nz7/V/PdgOqHqdc3P9/bc455p6ucxtndu+zh7fDPdaYux9/HYP8whdsrDgOuAN7S8/k46hUb+gc6o7E7gB8AvgH/JrlsD/BA4dIV7dXXcLLcZ+CnL+30E8LPmeAOwIbvXdcBz5rW989zHPdcfCXxlWvp43G2e9X7uc/7vgPOa463AIc3xIcDWaenncbZ3nvu41KZJ9PFYP8whdkoAVwLvXuGa38z+kTge+J+ln7NrTqQn8gCcAnyh59xBNJNPgCcC/wscMA9tppOLO7E5Pgm4pTl+Gt0TcLYx3klW427vPPfxwc1/92ruf9Y09PGE2jzT/Qw8HHhkc/7hwI3AKc3P76J7AuE7p6GfJ9Deue3j7Lmb6R48jL2P1zCbfh94FfDNaJZTAn8LPB4gpXQZ8DLgtRGxC/glsD41n/Iq1rP7RMl1wAXNvR4E/jKltHPwZlQZVZv/AnhPRKwB7gfObu73rYi4Bvg2sAs4J6X04JDbtJKxtpf57uPTI+Kc5vhjwIeb+026j2HMbWbG+zkiHgt8vDOnmzXAVSmlzzT3eAdwTUS8hs4/RKc195t0P4+1vcxxH0fEi4H30RkgfSoibk0pnTyJPnZ7akmSVGXmV1tIkqTxcvAgSZKqOHiQJElVHDxIkqQqDh4kSZojEbExIn4UEVtaXHtxVpzrOxFxb6vXcLWFJEnzIyLWAf8HXJlSenrF884F1qaUzlrtWiMPkiTNkZTSF4GuvS0i4kkR8ZmmXsaXIuIpfZ56Oi0LQs7qJlGSJKm9D9LZMOv2iDgBuAR43tKDEfEEOrtTfr7NzRw8SJI0xyLiEcDvAf/W7FwJna2sc+uBa9vuTOngQZKk+bYXnSqkx6xwzXrgnBUe3+2GkiRpTqWUfgZ8LyJOA4iOo5cej4ijgP2Bm9re08GDJElzJCKupjMQOCoidjTFw14JvCYivgF8Czg1e8rpwKaWxSM7r+FSTUmSVMPIgyRJquLgQZIkVXHwIEmSqjh4kCRJVRw8SJKkKg4eJElSFQcPkiSpyv+Gdk2JAAAABUlEQVQDysH8Jui1+kAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create a figure with one subplot\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "\n",
    "# Plot the grid with population info\n",
    "pop.plot(ax=ax, column='pop18', cmap=\"Reds\", scheme='quantiles', legend=True);\n",
    "\n",
    "# Add title\n",
    "plt.title(\"Population 2018 in 250 x 250 m grid squares\");\n",
    "\n",
    "# Remove white space around the figure\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-  Finally, let's save the result point layer into a file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Output path\n",
    "outfp = r\"data/addresses_population.shp\"\n",
    "\n",
    "# Save to disk\n",
    "join.to_file(outfp)"
   ]
  }
 ],
 "metadata": {
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}