3_Vector-based geospatial analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## This is a notebook for the Vector processing in Python session\n",
    "You will be introduced to how to interact with vector objects in Python.\n",
    "\n",
    "### Learning Objectives:\n",
    "- to Create geospatial data \n",
    "- to Read / write spatial data from/to different file formats\n",
    "- to Deal with different projections\n",
    "- to Reclassify your data based on different criteria\n",
    "- to Do queries and subsetting \n",
    "- to do Spatial overlays \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## Creating vector objects \n",
    "\n",
    "We assume that you are familiar with the two major types of data models that are used in GIS, which are vector and raster. Don't worry if you are NOT. Watch [this](https://gisgeography.com/spatial-data-types-vector-raster/) and/or read [this](https://www.youtube.com/watch?v=hbOGp51nYGI).\n",
    "\n",
    "Basic geometric objects in vector are ``Points``, ``Lines`` and ``Polygons``, which consist of coordinate tuples where:\n",
    "\n",
    "``Point`` -object represents a single point in space. ~ either 2D (x, y) or 3D (x, y, z).\n",
    "\n",
    "``LineString`` -object represents a set of points linked to form a line. So, a list of at least 2 coordinate tuples.\n",
    "\n",
    "``Polygon`` -object represents a filled area that consists of a list of > 3 coordinate tuples that forms the outerior ring and possibly a list of hole polygons.\n",
    "\n",
    "\n",
    "### Useful functionalities of ``Shapely``:\n",
    "\n",
    "- Create ``Points``, \n",
    "- Create ``Lines`` or ``Polygon`` from a Collection of ``Point`` -geometries,\n",
    "- Calculate ``areas``, ``length``, ``bounds``, etc. of input geometries,\n",
    "- Do geometric operations e.g., ``Union``, ``Difference``, ``Distance``, etc.,\n",
    "- Conduct spatial queries between geometries such as ``Intersects``, ``Touches``, ``Crosses``, ``Within``, etc.\n",
    "\n",
    "Here is a detailed [manual](https://shapely.readthedocs.io/en/latest/manual.html) to it, but you may read it right before [spatial analysis methods](https://shapely.readthedocs.io/en/latest/manual.html#spatial-analysis-methods).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import necessary geometric objects from shapely module \n",
    "from shapely.geometry import Point, LineString, Polygon "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "note that the coordinates are at Geographic Coordinate system. We recommend you learning about [coordinate systems](http://resources.esri.com/help/9.3/arcgisengine/dotnet/89b720a5-7339-44b0-8b58-0f5bf2843393.htm). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create three Point objects with coordinates by throwing in a pair of x,y "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "point1 = Point(2.2, 4.2) \n",
    "point2 = Point(7.2, -25.1) \n",
    "point3 = Point(9.26, -2.456) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the [Shapely documentation](https://shapely.readthedocs.io/en/stable/manual.html#coordinate-systems):\n",
    "#### Shapely does not support coordinate system transformations. All operations on two or more features presume that the features exist in the same Cartesian plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is a 3D point\n",
    "point3D = Point(9.26, -2.456, 0.57) \n",
    "\n",
    "# What is the type of the point? \n",
    "point_type = type(point1)\n",
    "\n",
    "print (point_type)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see above, our point1 variable is now a ``Shapely`` ``Point`` -object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let's now check the other ``Point`` objects "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "print(point1) \n",
    "print(point3D) \n",
    "print(type(point1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### ``Point`` -object has some built-in attributes and some useful functionalities that can be accessed, e.g.,\n",
    "    - to extract the coordinates of Points,\n",
    "    - to calculate the Euclidian distance between them.\n",
    "    \n",
    "Let's get the coordinates of one of our points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "point_coords = point1.coords \n",
    "\n",
    "print (\"The type is\", type(point_coords))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get the x and y coordinates together or seperately in the next syntax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xy = point_coords.xy \n",
    "print(\"xy variable:\\n\", xy, \"\\n\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also get `x` and `y` separately:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = point1.x \n",
    "y = point1.y \n",
    "\n",
    "print(\"x variable:\\n\", x, \"\\n\") \n",
    "print(\"y variable:\\n\", y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "As you can see, the ``xy`` -variable contains a tuple where x and y coordinates are stored inside numpy arrays. Using the attributes point1.x and point1.y, it is possible to get the coordinates directly as plain decimal numbers.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Now, we are going to calculate the distance between the two points. Since Shapely only knows cartesian coordinates, we assume that we have already projected our coordinates (remember the [North Korea example](https://www.independent.com/2017/06/14/map-rap/)!).\n",
    "\n",
    "The units of the distance measurement would then be the same as the units of the coordinate reference system, e.g. meters for UTM coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "point_dist = point1.distance(point2) \n",
    "\n",
    "print(\"Distance between the points is {0:.2f}\".format(point_dist))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## Creating a ``LineString``\n",
    "\n",
    "``LineString`` -object has many useful built-in functions e.g.,:  \n",
    "- extracting the coordinates or the length of a ``LineString`` (line),\n",
    "- calculate the ``centroid`` of the line, \n",
    "- create points along the line at specific distance, \n",
    "- simplify geometry. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's Create a LineString from our Point objects :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "line1 = LineString([point1, point2, point3])\n",
    "print(\"line1 variable: \\n\", line1, \"\\n\") \n",
    "\n",
    "# you can print the type of line to see what kind of Shapely object it is\n",
    "print(\"type of the line: \\n\", type(line1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to use coordinate tuples having the same outcome "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line2 = LineString([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])\n",
    "print(\"line2 variable: \\n\", line2, \"\\n\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## LineString attributes and functions\n",
    "We can extract the coordinates of a LineString similarly as we did with ``Point``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Get x and y coordinates of the line \n",
    "lxy = line1.xy \n",
    "\n",
    "print(\"The x, y coordinates of the line is\", lxy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the coordinates are again stored as a numpy arrays where first array includes all x-coordinates and the second one all the y-coordinates, respectively.\n",
    "\n",
    "We can extract x and y coordinates seperately:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "line_x = lxy[0] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract y coordinates straight from the LineObject by referring to a array at index 1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line_y = line1.xy[1] \n",
    "\n",
    "print(\"line_x:\\n\", line_x, '\\n') \n",
    "print('line_y:\\n', line_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "It is possible to retrieve specific attributes e.g., ``lenght`` of the line and center of the line (``centroid``) straight from the ``LineString`` object itself:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "l_length = line1.length \n",
    "\n",
    "l_centroid = line1.centroid \n",
    "\n",
    "# What type is the centroid? \n",
    "centroid_type = type(l_centroid) \n",
    "\n",
    "# Print the outputs \n",
    "print(\"Length of our line1: {0:.2f}\".format(l_length)) \n",
    "print(\"Centroid of our line1: \", l_centroid) \n",
    "print(\"Type of the centroid:\", centroid_type)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## Create a polygon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Creating a ``Polygon`` -object is similiar to ``Points`` and ``Lines`` creation, but you need at least three coordinate-tuples\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "poly = Polygon([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use our previously created Point objects (same outcome) \n",
    "--> notice that Polygon object requires x,y coordinates as input "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "poly2 = Polygon([[p.x, p.y] for p in [point1, point2, point3]]) \n",
    "\n",
    "# Geometry type can be accessed as a String \n",
    "poly_type = poly.geom_type \n",
    "\n",
    "# Using the Python's type function gives the type in a different format \n",
    "poly_type2 = type(poly) \n",
    "\n",
    "# Let's see how our Polygon looks like \n",
    "print('poly:', poly) \n",
    "print('poly2:', poly2) \n",
    "print(\"Geometry type as text:\", poly_type) \n",
    "print(\"Geometry how Python shows it:\", poly_type2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Notice that Polygon representation has double parentheses around the coordinates like POLYGON ((<values in here>)). This is because a Polygon can also have holes inside of it.\n",
    "    \n",
    "If you want to know more about it, here is the help of ``Polygon``. Notice the ``exterior`` and ``interior`` attributes. \n",
    "A Polygon can be constructed using exterior and interior coordinates i.e., interior coordinates == hole inside the Polygon:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(Polygon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### Let’s see how we can create a Polygon with a hole inside.\n",
    "\n",
    "Let's create a rectangle make a whole in it. First we define our exterior: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "world_exterior = [(-180, 90), (-180, -90), (180, -90), (180, 90)] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note that these are still not degrees, even though the numbers migth suggest that!*\n",
    "\n",
    "Let's create a single big hole where we leave ten units at the boundaries. \n",
    "\n",
    "Notice: there could be multiple holes, therefore we need to provide a list of holes: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hole = [[(-170, 80), (-170, -80), (170, -80), (170, 80)]] \n",
    "\n",
    "# World without a hole \n",
    "world = Polygon(shell=world_exterior) \n",
    "\n",
    "# Now we can construct our Polygon with the hole inside \n",
    "world_has_a_hole = Polygon(shell=world_exterior, holes=hole)\n",
    "\n",
    "#Let’s see what we have now\n",
    "print('world:', world) \n",
    "print('world_has_a_hole:', world_has_a_hole) \n",
    "print('type:', type(world_has_a_hole))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "#### Let's see a few of the available attributes and how to access them:\n",
    "Note that the units are in terms of degrees (Each degree of latitude is approximately 111 kilometers apart. The range varies (due to the earth's slightly ellipsoid shape) from 110.567 km at the equator to 111.699 km at the poles. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "world_centroid = world.centroid \n",
    "world_area = world.area \n",
    "world_bbox = world.bounds \n",
    "\n",
    "world_ext = world.exterior \n",
    "world_ext_length = world_ext.length \n",
    "\n",
    "print(\"Poly centroid: \", world_centroid) \n",
    "print(\"Poly Area: \", world_area) \n",
    "print(\"Poly Bounding Box: \", world_bbox) \n",
    "print(\"Poly Exterior: \", world_ext) \n",
    "print(\"Poly Exterior Length: \", world_ext_length)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can make a quick plot of the objects you created just by typing its name...don't worry, you will learn how to make nice maps!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "world"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "world_has_a_hole"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# Part 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# ![GeoPandas](https://geopandas.org/en/stable/_static/geopandas_logo_web.svg)\n",
    "\n",
    "[Geopandas](https://geopandas.org) is a Python module that is built on top of [Pandas](https://pandas.pydata.org/docs/user_guide/index.html) extending its functionalities, i.e., all functionalities of pandas works also in geopandas. \n",
    "\n",
    "``Geopandas`` has many good built-in spatial processing/analysis features such as overlay analysis, geocoding, spatial aggregation methods.``GeoPandas`` makes working with geospatial data in python easier.\n",
    "\n",
    "Geometric operations are performed by ``shapely``. \n",
    "\n",
    "Geopandas also depends on [fiona](https://fiona.readthedocs.io/en/latest/manual.html#introduction) for file access, and [descartes](https://pypi.org/project/descartes/) + [matplotlib](https://matplotlib.org) for plotting.\n",
    "\n",
    "\n",
    "You will learn \n",
    "- how to read and write spatial data, \n",
    "- how to deal with projections in Python using Geopandas \n",
    "- and some useful data management functionalities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "#### Let’s Download Land Use data for Copenhagen 2012\n",
    "\n",
    "This bit will download a zip file with data for Copenhagen and unzip it into a folder called `Data2`.\n",
    "\n",
    "#### What is *land use* data? \n",
    "Land Use data refers to data describing what each piece of land is being used for i.e., the usage of a land parcel in the form of pixels of polygons.\n",
    "\n",
    "#### What is the *Urban Atlas*? \n",
    "check [this](https://land.copernicus.eu/local/urban-atlas/urban-atlas-2012?tab=metadata) out!\n",
    "\n",
    "#### What is a *shapefile*?\n",
    " Check [this](https://desktop.arcgis.com/en/arcmap/10.3/manage-data/shapefiles/what-is-a-shapefile.htm) out!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from io import BytesIO\n",
    "from urllib.request import urlopen\n",
    "from zipfile import ZipFile\n",
    "\n",
    "#the data is hosted on the European Environment Agency's website, \n",
    "# so in case it is missing, you can find it in other websites as well\n",
    "zipurl = 'https://www.eea.europa.eu/data-and-maps/data/urban-atlas/denmark/dk001l_kobenhavn.zip/at_download/file'\n",
    "with urlopen(zipurl) as zipresp:\n",
    "    with ZipFile(BytesIO(zipresp.read())) as zfile:\n",
    "        zfile.extractall('Data2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make sure you know where the data is stored after downloading and unzipping\n",
    "import os\n",
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's open our fresh new shapefile using [geopandas](https://geopandas.org):\n",
    "Geopandas can read vector data formats e.g., Shapefile, GeoJSON, KML. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Import necessary modules \n",
    "import geopandas as gpd \n",
    "import os\n",
    "\n",
    "#Set filepath \n",
    "fp = \"Data2/dk001l_kobenhavn/dk001l_kobenhavn.shp\" \n",
    "\n",
    "# Read file using gpd.read_file() \n",
    "data = gpd.read_file(fp)\n",
    "\n",
    "# you can see what kind of object we have now, which is a GeoPandas' GeoDataFrame\n",
    "print (\"The type of data is\", type(data))\n",
    "\n",
    "# let's preview our data but only the first 20 rows\n",
    "print(data.head(10)) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "what you see above is land use data and the associated attributes such as `Cities`, `LUZ_OR_CIT`, `city code`, `land use ITEM`, `Production date`, `shape length`, `shape area`, and `Geometry`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## Let's plot our data\n",
    "Creating a simple map from a GeoDataFrame ==> use ``.plot()`` -function that creates a map based on data geometries. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# a quicklook of the geodataframe\n",
    "%matplotlib inline \n",
    "\n",
    "data.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "if you'd like to make a bigger plot with more element, try ``help(data.plot)``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "help(data.plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "# a bit more decorated\n",
    "data.plot(figsize=(60, 25), column='ITEM', legend=True, categorical=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## How to write it to a new shapefile\n",
    "\n",
    "Let’s export first 50 rows of the data into a new Shapefile by first selecting the data and then with ``gpd.to_file()`` -function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Set the output path for the data \n",
    "outfp = \"Data2/dk001l_kobenhavn/dk001l_kobenhavn_subset.shp\" \n",
    "\n",
    "# Select first 50 rows \n",
    "selection = data[0:50] \n",
    "\n",
    "# Write those rows into a new Shapefile (the default output file format is Shapefile) \n",
    "selection.to_file(outfp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## 🏋 Exercise\n",
    "\n",
    "Plot the subset data! Add your code in the cell below. Feel free to become fancy and add a lot more elements to the plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s print the first 5 rows of the data. You can specify a column like ‘ITEM’ and only show that attributes of that column:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "print(data['ITEM'].head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Since the spatial data is stored as Shapely objects, it is possible to use all of the functionalities of ``Shapely`` module.\n",
    "Let’s prove that this is really the case by iterating over a sample of the data, and printing the area of first ten polygons.\n",
    "We can iterate over the rows by using the ``iterrows()`` -function.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Make a selection that contains only the first five rows \n",
    "selection = data[0:10] \n",
    "# Iterate over rows and print the area of a Polygon \n",
    "for index, row in selection.iterrows(): \n",
    "    # Get the area of the polygon \n",
    "    poly_area = row['geometry'].area \n",
    "    # Print information for the user \n",
    "    print(\"Polygon area in m2 at index {index} is: {area:.3f}\".format(index=index, area=poly_area))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s now create a new column into our GeoDataFrame where we calculate/store the areas of polygons into that column. No need to iterate like before.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Create a new column called 'area' and assign the area of the Polygons into it \n",
    "data['area'] = data.area \n",
    "# Print first 5 rows of the area column \n",
    "print(data['area'].head(5))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## Let’s get the min, max and mean of areas.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Maximum area \n",
    "max_area = data['area'].max() \n",
    "# Minimum area \n",
    "min_area = data['area'].min() \n",
    "# Mean area \n",
    "mean_area = data['area'].mean() \n",
    "\n",
    "print(\"Max area: {max} m2 \\nMin area: {min} m2 \\nMean area: {mean} m2\".format(max=round(max_area, 2), min=round(min_area, 2), mean=round(mean_area, 2)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# Creating GeoDataFrame from freshly baked geometries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Since geopandas takes advantage of Shapely geometric objects, it is possible to create a Shapefile from scratch by passing Shapely’s geometric objects into the GeoDataFrame. This is useful as it makes it easy to convert e.g. a text file that contains coordinates into a Shapefile. Next we will see how to create a Shapefile from scratch.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s create an empty GeoDataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Create an empty geopandas GeoDataFrame \n",
    "newdata = gpd.GeoDataFrame() \n",
    "\n",
    "# Let's see what we have at the moment \n",
    "print(newdata)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s create a new column called geometry that will contain our Shapely objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Create a new column called 'geometry' to the GeoDataFrame \n",
    "newdata['geometry'] = None \n",
    "\n",
    "# Let's again see what's inside \n",
    "print(newdata)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Now we have a geometry column in our GeoDataFrame but we don’t have any data stored yet.\n",
    "Let’s create a Shapely Polygon repsenting our current location (AAU) that we can later insert to our GeoDataFrame:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Coordinates of the AAU campus in Decimal Degrees (but still on a flat Earth!)\n",
    "coordinates = [(12.543269, 55.650061), (12.543806, 55.650511), (12.542054, 55.651096), (12.541585, 55.650672)] \n",
    "\n",
    "# Create a Shapely polygon from the coordinate-tuple list \n",
    "poly = Polygon(coordinates) \n",
    "\n",
    "# Let's see what we have \n",
    "print(poly)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Okay, now we have an appropriate ``Polygon`` -object.\n",
    "Let’s insert the polygon into our `geometry` column of our GeoDataFrame at position 0:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Insert the polygon into 'geometry' -column at index 0 \n",
    "newdata.loc[0, 'geometry'] = poly \n",
    "\n",
    "# Let's see what we have now \n",
    "print(newdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Now, we have a GeoDataFrame with a Polygon that we could already export to a Shapefile. However, typically you might want to include some useful information with your geometry.\n",
    "Hence, let’s add another column to our GeoDataFrame called location with text AAU that describes the location of the feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Add a new column and insert data \n",
    "newdata.loc[0, 'location'] = \"Aalborg University\"\n",
    "\n",
    "# Let's check the data \n",
    "print(newdata)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Before exporting, it is necessary to specify the coordinate reference system (projection) for the GeoDataFrame. \n",
    "GeoDataFrame has an attribute called .crs that shows the coordinate system of the data.\n",
    "\n",
    "### 🔥 If you have never worked with GIS data before, make sure to read this [introduction to coordinate reference systems](https://docs.qgis.org/3.10/en/docs/gentle_gis_introduction/coordinate_reference_systems.html)!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "print(newdata.crs) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "The CRS is not specified here, which will be a problem once we want to combine the data with other data – we simply don't know which area on the Earth it refers to, so let’s add a `crs` for our GeoDataFrame.`fiona` has a function called ``from_epsg()`` for passing the coordinate reference system information for the GeoDataFrame. Use WGS84 (epsg code: 4326):\n",
    "\n",
    "PS. The [EPSG database](https://epsg.io) assigns standardised codes to commonly used CRSs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Import specific function 'from_epsg' from fiona module \n",
    "from fiona.crs import from_epsg \n",
    "\n",
    "# Set the GeoDataFrame's coordinate system to WGS84: epsg code 4326 \n",
    "newdata.crs = from_epsg(4326) \n",
    "\n",
    "# Let's see how the crs definition looks like \n",
    "print(newdata.crs)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Now, we can export the GeoDataFrame using ``.to_file()`` -function.\n",
    "We only need to provide the output path for the Shapefile. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Determine the output path for the Shapefile \n",
    "outfp = \"AAU.shp\" \n",
    "\n",
    "# Write the data into that Shapefile \n",
    "newdata.to_file(outfp)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# Practical example: Saving multiple Shapefiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "A very useful function that can be used in Pandas/Geopandas is ``.groupby()``, which allows you to group data based on values on selected column(s).\n",
    "Next we will take a practical example by automating the file export task. \n",
    "\n",
    "We will group individual land use types in our ”dk001l_kobenhavn.shp” and export those into separate Shapefiles.\n",
    "Let’s start from scratch and read the Shapefile into GeoDataFrame:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Read land use data \n",
    "fp = \"Data2/dk001l_kobenhavn/dk001l_kobenhavn.shp\" \n",
    "data = gpd.read_file(fp) \n",
    "\n",
    "# Print columns \n",
    "print(data.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "The `TYPE` column in the data contains information about different land use types. With ``.unique()`` -function we can quickly see all different names in that column:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Print all unique land use types in ‘ITEM' column \n",
    "print(data['ITEM'].unique())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Now we can use that information to group our data and save all individual land use types as separate Shapefiles:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Group the data by column ‘ITEM’ \n",
    "grouped = data.groupby('CODE') \n",
    "\n",
    "# Let's see what we have \n",
    "grouped\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "As we can see,``groupby`` -function gives us an object called `DataFrameGroupBy`, which is similar to list of keys and values (like in a dictionary) that we can iterate over.  \n",
    "\n",
    "Let’s iterate over the groups and see what our variables key and values contain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Iterate over the group object \n",
    "for key, values in grouped: \n",
    "    individual_LU = values \n",
    "    \n",
    "# Let's see what is the LAST item and key that we iterated \n",
    "print('number of records:', key)\n",
    "print(individual_LU)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Finally, let’s  export all individual LU types into separate Shapefiles:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Determine output directory \n",
    "out_directory = \"dk001l_kobenhavn\" \n",
    "\n",
    "# Create a new folder called 'Results’ \n",
    "result_folder = os.path.join(out_directory, 'Results') \n",
    "\n",
    "# Check if the folder exists already \n",
    "if not os.path.exists(result_folder): \n",
    "# If it does not exist, create one \n",
    "    os.makedirs(result_folder) \n",
    "# Iterate over the groups \n",
    "\n",
    "for key, values in grouped: \n",
    "# Format the filename (replace spaces with underscores using 'replace()'-function) \n",
    "    output_name = \"%s.shp\" % key.replace(\" \", \"_\") \n",
    "# Print some information for the user \n",
    "    print(\"Processing: %s\" % key) \n",
    "# Create an output path \n",
    "    outpath = os.path.join(result_folder, output_name) \n",
    "# Export the data \n",
    "    values.to_file(outpath)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# Map projections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Shapefiles contain information about the coordinate reference system that is stored in `.prj` -file (if the data has been appropriately produced). When reading the data into GeoDataFrame with Geopandas this information is automatically stored into `.crs` of the GeoDataFrame.\n",
    "Let’s read Europe_borders.shp file and checking the crs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "#let's download a file containting European borders\n",
    "from io import BytesIO\n",
    "from urllib.request import urlopen\n",
    "from zipfile import ZipFile\n",
    "\n",
    "#the data is from EEA \n",
    "zipurl ='https://github.com/jimjoker/Geospatial-Analysis-Notebooks/raw/master/Data2.zip'\n",
    "with urlopen(zipurl) as zipresp:\n",
    "    with ZipFile(BytesIO(zipresp.read())) as zfile:\n",
    "        zfile.extractall('Data')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "# Read the file \n",
    "fp = \"Data/Data2/Europe_borders.shp\" \n",
    "\n",
    "\n",
    "data = gpd.read_file(fp) \n",
    "data.plot()\n",
    "\n",
    "# Check the details of the coordinate reference system. It's \"WGS84\", EPSG:4326\n",
    "print(\"projection is:\", data.crs)\n",
    "\n",
    "data['geometry'].head()\n",
    "#the crs is a Python dictionary with a key init that has a value epsg:4326. This is a very typical way how CRS is stored in GeoDataFrames. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s continue by checking the values in our ``geometry`` -column to verify that the CRS of our GeoDataFrame is correct: (note that the units are not decimal degress, it is metric)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "data['geometry'].head() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "WGS84 projection is not really a good one for representing European borders, so let’s convert those geometries  to Lambert Azimuthal Equal Area projection (EPSG: 3035), which is the recommended projection by European Comission, to see the diffrences of the two projections.\n",
    "\n",
    "Let’s re-project our data into EPSG 3035 using ``epsg`` -parameter:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Let's make a copy of our data \n",
    "orig = data.copy() \n",
    "\n",
    "# Reproject the data \n",
    "data = data.to_crs(\"EPSG:3035\") \n",
    "\n",
    "\n",
    "# Check the new geometry values \n",
    "print(data['geometry'].head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s compare the datasets by making maps out of them.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline \n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "# Make subplots that are next to each other \n",
    "fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(15, 8)) \n",
    "\n",
    "# Plot the data in ETRS-LAEA projection  \n",
    "orig.plot(ax=ax1, facecolor='gray')\n",
    "\n",
    "# Add title \n",
    "ax1.set_title(\"WGS84\")\n",
    "\n",
    "# Plot the one with WGS84 CRS \n",
    "data.plot(ax=ax2, facecolor='black')\n",
    "\n",
    "# Add title \n",
    "ax2.set_title(\"ETRS Lambert Azimuthal Equal Area projection\")\n",
    "\n",
    "# Remove empty white space around the plot \n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "You can visually comapre the two maps and see how distances and areas get distorted across the entire area. \n",
    "\n",
    "Let's export the re-projected data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Ouput filepath \n",
    "outfp = \"Europe_borders_epsg3035.shp\" \n",
    "# Save to disk \n",
    "data.to_file(outfp)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# Calculating distances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "The aim is to find the Euclidean distances from the centroids of all European countries to Copenhagen using a metric projection (`Azimuthal Equidistan` -projection) in meters. \n",
    "Let’s first create a GeoDataFrame that contains a single Point representing the location of CPH:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Import necessary modules \n",
    "from shapely.geometry import Point \n",
    "import pycrs \n",
    "# Create the point representing CPH (in WGS84) \n",
    "cph_lon = 12.5683 \n",
    "cph_lat = 55.6761 \n",
    "# Create GeoDataFrame \n",
    "copenhagen = gpd.GeoDataFrame([[Point(cph_lon, cph_lat)]], geometry='geometry', crs='epsg:4326', columns=['geometry']) \n",
    "# Print \n",
    "print(copenhagen)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Import pyproj \n",
    "import pyproj \n",
    "\n",
    "# Define the projection using the coordinates of our CPH point (cph_lat, cph_lon) as the center point \n",
    "# The .srs here returns the Proj4-string presentation of the projection \n",
    "aeqd = pyproj.Proj(proj='aeqd', ellps='WGS84', datum='WGS84', lat_0=cph_lat, lon_0=cph_lon).srs\n",
    "\n",
    "# Reproject to aeqd projection using Proj4-string \n",
    "copenhagen = copenhagen.to_crs(crs=aeqd) \n",
    "\n",
    "# Print the data \n",
    "print(copenhagen) \n",
    "# Print the crs \n",
    "print(\"\\nCRS:\\n\", copenhagen.crs)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Next we want to transform the Europe_borders.shp data into the desired projection.\n",
    "\n",
    "Let’s create a new copy of our GeoDataFrame into a variable called europe_borders_aeqd:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Create a copy \n",
    "europe_borders_aeqd = data.copy() \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s now reproject our Europe borders data into the Azimuthal Equidistant projection that was centered into Copenhagen:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Reproject to aeqd projection that we defined earlier \n",
    "europe_borders_aeqd = europe_borders_aeqd.to_crs(crs=aeqd) \n",
    "\n",
    "# Print \n",
    "print(europe_borders_aeqd.head(2))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s plot the Europe borders and the location of CPH to get a better understanding how our projection has worked out:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline \n",
    "fig, ax = plt.subplots(figsize=(10,10)) \n",
    "\n",
    "# Plot the country borders \n",
    "europe_borders_aeqd.plot(ax=ax) \n",
    "\n",
    "# Plot the Copenhagen point on top of the borders using the same axis \n",
    "copenhagen.plot(ax=ax, color='red', markersize=15)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Our goal is to calculate the distance from all countries to Copenhagen. So, we need to calculate the centroids for all the Polygons representing the boundaries of European countries.\n",
    "This can be done easily in Geopandas by using the ``centroid`` attribute:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "europe_borders_aeqd['centroid'] = europe_borders_aeqd.centroid \n",
    "print(europe_borders_aeqd.head(5))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "First, we will create a dedicated function for calculating the distances called ``calculate_distance()``:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def calculate_distance(row, dest_geom, src_col='geometry', target_col='distance'):\n",
    "    \"\"\"\n",
    "    Calculates the distance between Point geometries.\n",
    "    Parameters\n",
    "    ----------\n",
    "    dest_geom : shapely.Point\n",
    "       A single Shapely Point geometry to which the distances will be calculated to.\n",
    "    src_col : str\n",
    "       A name of the column that has the Shapely Point objects from where the distances will be calculated from.\n",
    "    target_col : str\n",
    "       A name of the target column where the result will be stored.\n",
    "    Returns\n",
    "    -------\n",
    "    Distance in kilometers that will be stored in 'target_col'.\n",
    "    \"\"\"\n",
    "    # Calculate the distances\n",
    "    dist = row[src_col].distance(dest_geom)\n",
    "\n",
    "    # Convert into kilometers\n",
    "    dist_km = dist / 1000\n",
    "\n",
    "    # Assign the distance to the original data\n",
    "    row[target_col] = dist_km\n",
    "    return row\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Before using our function and calculating the distances between Copenhagen and centroids, we need to get the Shapely point geometry from the re-projected Copenhagen center point that we can pass to our function (into the ``dest_geom`` -parameter. We can use the ``loc`` -functionality to retrieve the value from specific index and column:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Retrieve the geometry from Copenhagen GeoDataFrame \n",
    "copenhagen_geom = copenhagen.loc[0, 'geometry'] \n",
    "\n",
    "print(copenhagen_geom)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "When iterating over a GeoDataFrame, apply function is used like this\n",
    "GeoDataFrame.apply(name_of_your_function, param1, param2, param3,  axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Calculate the distances using our custom function called 'calculate_distance’ \n",
    "europe_borders_aeqd = europe_borders_aeqd.apply(calculate_distance, \n",
    "    dest_geom=copenhagen_geom, src_col='centroid', target_col='dist_to_cph', axis=1) \n",
    "\n",
    "print(europe_borders_aeqd.head(10))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "data.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s check what is the longest and mean distance to Copenhagen from the centroids of other European countries:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Calculate the maximum and average distance \n",
    "max_dist = europe_borders_aeqd['dist_to_cph'].max() \n",
    "mean_dist = europe_borders_aeqd['dist_to_cph'].mean() \n",
    "\n",
    "print(\"Maximum distance to Copenhagen is %.0f km, and the mean distance is %.0f km.\" % (max_dist, mean_dist))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# Geometric operations: Overlay analysis"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAABJQAAAEaCAYAAABHKA1KAAAAAXNSR0IArs4c6QAAAHhlWElmTU0AKgAAAAgABAEaAAUAAAABAAAAPgEbAAUAAAABAAAARgEoAAMAAAABAAIAAIdpAAQAAAABAAAATgAAAAAAAACQAAAAAQAAAJAAAAABAAOgAQADAAAAAQABAACgAgAEAAAAAQAABJSgAwAEAAAAAQAAARoAAAAApRg5dgAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAQABJREFUeAHsvQd8VNeZ9//TdGlGvfeGRBEC0auxKcYGF5zYseONnWw2m934n942Wb/7ZpPN55N9N+8m2cQpzptC3Bs2NjYGbGM6GARCCISEUO+9lxlN0f88B8sVsDSadu88J5Glmbn3lO+53LnnOc/ze0ImRAEXJsAEmAATYAJMgAkwASbABJgAE2ACTIAJMAEmMEUCmikex4cxASbABJgAE2ACTIAJMAEmwASYABNgAkyACTABSUDHHJgAE7g6gfHxcdAPlZCQEBQXF6O8vBwGg+FDJ7hcLoSGhmL58uVISkqCXq+XnxuNRuh0/E/sQ7D4BRNgAkyACTABJsAEmAATYAJMgAmoggCvdlUxjTyImRIYHR1FU1OTNBxpNBr09fWhubkZ7e3toNf0c/78eVRXV79nMJps0+l0wmKxyJd2u10aoUwmExITEzFr1ixZJxmdyLgUERGB+Pj4yVP5NxNgAkyACTABJsAEmAATYAJMgAkwAUUSCGENJUXOG3d6GgTI4EM/VMjTyGq1oqenR/5NEmJkDKqsrMSTTz4pjyHDz9jYGBoaGlBVVQWtVisNSg8++CC2bt0KMj59sFCdVH9HRwd++9vfoqamBrGxscjOzkZubi4cDoc0MoWHh2Pu3LnYtGmTNC7R++TZFBUVJeunOqmuyfY+2Ab/zQSYABNgAkyACTABJsAEmAATYAJMIJAIsEEpkGaD++JxAmQwunDhAk6dOiUNO2azGWfOnMGvf/1r6WlEIW2LFi3CunXrkJGRIQ1D5I10ww03YMGCBe8Zj8jQQyFsFO52LR178kIiQ9SkNxIZrg4dOoSWlhbZts1mk4arHTt2yNdkUCoqKsLnP/95UL/ovOTkZNkf6gsXJsAEmAATYAJMgAkwASbABJgAE2ACgUqADUqBOjPcL7cJkMGHPI6eeOIJXLp0CUuXLpXha5OeShRyRsYiOo6MODExMUhLS5O/J41FcXFxoLC1mZbe3l5pZKJ6yYBEr8nridqlQm1Qvw4cOCA/p76QIWpwcBC33nor1q5di5SUlJl2g89nAkyACTABJsAEmIBPCLicLowIb276jZAPNCmehcgLO0xsomk0H/zgA8fwn0yACTABJqAoAmxQUtR0cWevR2B4eFgakMgbicLVSP8oISEBGzZseE9cmzySyPuHBLQDpZBW09mzZ6WRiTygWltb5WsS946OjpZhcuRFRXpMXJgAE2ACTIAJMAEmEEgE6mqaxPNXFfRCMkCr1aH84nnU1FZjWBiVPmg4cgoDk0UYk+bOmYOsjBwZ/i9sTEJzMgHzF84JpCFxX5gAE2ACTGCKBNigNEVQfFjgEhgaGkJtbS0qKipkFrZJDSPSOyK9osmsa4E7go/3jATCd+/ejdOnT0tdJTImLVy4EHPEQ1hmZqZ87+Nn8TtMgAkwASbABJgAE/AcgeGhEQwPj8gKSRJgeGhYbNh1CC3KXvleWNQEKqrP4p2Tx4UsgBGaEC2aOsvR0lODEeu4eP1+X8hhKSLciPlzCxDqSsOEQ4sQDZCZkY2b1t6K4R6heelwQaPVIDomWmyk5bznTW42WxAReSUByvs18l9MgAkwASbgbwJsUPL3DHD7bhMgbyPySnr77bfx6quvyixt5Mlz3333YfHixXLny+3KA+TErq4uHDx4EK+88ooM47v55ptB4uA5OTnSUEau41yYABNgAkyACTABJuAuAQrLpx8KwadnK/qbEpR0dnbgXFmp8D6qlFUbjAY0NjXg8JGDKD55DhAGooVbdSi4wYLwaIMIcRN1OCZQsCYKOQstMOm0EA5IHypkXxqHE4ee6cRgjx06QwiGeu0oO9SHisOAQ+Q90RiAwvkF+PvPfxFOIRFAGpTz5s3FyuWrEB4eKftHxi3qI3lEkRcUaV1yYQJMgAkwAd8TYIOS75lzix4gMDAwIDWSfv7zn0s9pG3btkkhbQpno7AxeshQQ5nUXiIvrOrqaum19Jvf/AYbN27Et7/9baxcuVLqEahhrDwGJsAEmAATYAJMwPcEyNuorbUdpWfL8Mtf/Epu1rkcQFjsOGJyRmBOcAkjDqDTh2DOKgsK1oVhdPiK4SksQguTRSueRYRBRxxDBiSdXgOtOPZaNh6qy24VdbrEH+IgCoUbG3LCOiIsVFSB+GmrsePUrmFhLKKXEzCbQjHea8bFQ+L1hBP5s/Nx8+ZN2LBxvdhkyxTeSxG+B8ctMgEmwASYANigxBeBYgiQkDWJbO/cuRPFxcVSrLqgoADz589Hfn6+iMFPVPUOld1ul9pQpaWlOH/+POrq6qSX0pYtW7BmzRqZIU4xk8kdZQJMwC0Chw4cQ1NDM2zCi0DzgdWaXYj+k+fi0qWLEBUb7lbdfBITYALBRaCvbwB7Xn0TR0++iVF9LWLzhjEW0gq6n7iEp5E5SoP4TC0iEvQSDN1yIuP1iE2kpCXC0iOtPx/3DHIJi9GEsA1dr2iEAerjZ5I1id6dwPCgA+31I1ee6+gtYXwa7HKhocwmD4mOjYC1Jwx9NVGIRD7W33Ar1qxdheSUxOs1y58xASbABJiAhwmwQcnDQLk6zxMYFaKOZEgiPSHSSqIwsMjISCm2TZ46nsjG5vlee7dGymK3f/9+GQZH3lixsbFSvJtC/bKzs73bONfOBJiAVwg4xl0iI2Ur+vr6pLHIZnXgyJGjUqjfoDeK6BIHhnQVsDp7RRiIQyy73l+OUahKaloK4izZ6KoKFQsuF6KjolEwfx6yc7OEB4BDholERUUiPTNVeBPQYpALE2ACwUigq6MHBw8cwTunD6FrvAx2cy1iZo0ge5ERablR0ogjw96Ep5ExVAs9CR29W5zC2OMQoW2T9mwyAUm70uQB7v4Wt7P372jCXCXC2HQfeIfaGXc5pSeTbFvjQmvtEKpOjQqjkgWhtllIsSzGwrkrsHzFMmTlpLvbEz6PCTABJsAEpkGADUrTgMWH+pYA6SPV19fj4sWLuHDhAkhsOyoqCrfccgs2b94clIakj84A8Xn99delsY30lMhTa1K8m8L/uDABJhAYBGi3fmRkVBp2aCVms1rR3dUjjEf9UoBWI3b7B22tuFRTjob6RuF9qINJRHC8/OoL4t7XASFFAoOwE930gAVJWSahOyK0SShc5N1C+iFCvQR1ZcM48LgVQhcXqanRWLJoFXLS54rwFJswKLkQH5eEvOxCJMdliEWhiGkRJTMzA+ERFqlVQl5PZouZDU7vcuVfTEAtBCisrLOjC3W1DTh55jAOFe9EXddJrLk3Cos2xor7ilkYcDTCbP3+fYXGTvcZClELhEKGpJAPqHyTwYkMXK11wyh+oxO1x7SICpmHpYUbsHzhesyekycM7clXvJwCYQDcBybABHxOgNaTpA1Hz0n0d1xcnFxDsu6a56aCDUqeY8k1eYiAQ7haW8Vi66233sKLL74oQ7uWLVsmxbaXLFmiyKxtHkJzzWo6OztxUIh3v/DCCzLT3V133YUvfelLIKMSCVeyePc10fEHTMBrBCZFbumhhYRsKyougvTfNBotWlqbsW/fHhw7fhRGo0mElwDL7tIhLtuFUaElojdqsHRLNDLnhsNk1sIlFnXkkUQCth9Km/SR3ktRXOHpREYrnYhSqT47hLLDfRjosAvDVQgGuhxoPBuCzhqheaK/IsT7/e/9ixC8nS+yN42IvhhQVLQIiQlJCA0Nk4Ylfuj6CGR+yQQUQoDuQWQMslrHxAZdHV7bvQu///0fgbgGPPCjdCy9NVbqHU0lRC2Qh0wOVFpaLA44UH6sD/uf6MKFvTp87atfw4MPfAF5s/LkPZE9MwN5FrlvTMA9AiSJQl7atNahex4Zj0gmxGg0yvfLysrQ398vX5eUlIASHNH6iLKA0xqJnnHoXPqbCtVH0R/87DP1+WCD0tRZ8ZE+IEDi0zt27MD/+l//S3rafPrTn8LatWuRnp4hbwT0j5/LxwnQDZRunoODg1K8+7XXXsMf/vAH6cn19a9/HcuXL1eNUPnHR8/vMIHAJNDd3YNjR07grQN70dRzGuP6dgwOjUptkbh0PWYtNSApTwenXfRf/BtOzDXBHKkVD0BXwknIqETituJZx+3isLtgt5H4rahC1GMTorf9nXaM9DulgYnqHmjVourkGJqqRmAKNSArNQcdlSZsufluPPCF+0RIbYzb7fOJTIAJ+I/AyPAIjh87hV2v7sTxs68ga6UTi7aYkLnAhNBwHQziHkOFjE5KL3Qvo3E4xP1udMiBrgY73tk5jJJdTqwo3IJ/+sqXxLPQEoSGCVdPLkyACaiCABmSSA6FnBDI+aC1tVW+pogWkgEhzd2kpCSxUVYkf5PDAkmG0HozIiIC0dHRInNkOM6cOSMlQ8iY1NjYiE2bNiEtLY2NSlO8StigNEVQfJh3CXR0dEhvpN27d0trMlmP6UYwZ84cFpueJnry8CLB7lOnTklvJbq50s1y69at0sDE3krTBMqHM4FpEmhubMOhQ4dx5MSb6Js4j+i8ASTPcyA8nnbOnGLVA7GY0yJGCNuGx1xJtU1N6A0i8O0D1iOKaJv0MJhWF8gA9e4CkXbuqU56iwrZlZwOl0zvTa+pucFeG/o6bRgbFkYmsUOn1xjRVG5H9TED9P35uHHVLbjzrtuQLTIpvVcRncyFCTCBgCTQL8S239x7EMeK96NtpBSmlHbE5Y8gY34o0maZERaql6Ft5NGoqiLuZ6S9RGYyu7jPtdYMo/rMMNoumODoTENy2EJsXL8VK1ctQ2JSvKqGzoNhAmokQLq5VVVV0rBDBiCbzSbXNrSJPmvWLGRlZcnPyPubkjNRSNvIyAgsFos0FvX29sq/Y2Ji3ttYb29vl/WQk8JkZvC2tja5VqJnIFo30bo0OztbeHbWg6JASKs2PT1dHk91kpGK2udyhQAblPhK8DsB+of56KOPCvHZI/If6ooVK/CNb3xD3gD83jmFd4C0p8hqT8al7u5u6flFGeHIlZMLE2ACniVQXVWH40dP4WLNSXTaSmE11SJ+zijyl4UjdVY4wrT6STuP+E2KR+JHrOeksUf8h8La3jvAs117vzbRDhmRpCv3u2tJkiTRil4IXyjZL5HfCUOD4yh5uxOXTzgQZstHjHY+FuSvEjv8qzCvIF8Yv/ge8j5U/osJBAaB7q5eHDl0DGfOHUPz4FmMGS9Lse285WHInhcFo9BeE/6K0gsyMHrsvV6QMd0g/mMX7pm15/tFKNyQFO+OcM1BRswSLJy3EkuWLBZJClK81wmumQkwgSkTIG8jMvY0NzfLdUpCQoLcICevorCwMLkupE02MvjQseRBlJqaKg1B9LmnChmtKioqpG5vS0sLenp6hARAKMgoRUYoSpxC7cXHx8u+UHQIeUTl5ubK9zzVDyXVwwYlJc2WCvtKMa3bt2/Hs88+K77Yl+Chhx7C7NmzpWFJhcP125DeeOMNPPzww1K0m7SVbrrpJtZV8ttscMNqIuBwiAeg1g5UXb6MEyVv4OCxPbBbqrHwVj2W3JyI+MQwaSNyCDdqJYWVkN5SqDAzDTvsKNnfhTf+1oIIlxC7nbcRK5esx9y8hUhKTpBi3mqaTx4LE1AaAbqvdHV2o762ESdE1raDJ19CU/9pLLnDjEUbYpCSaxZeh1ppWJGhr0ob4Az7SwZ0vfA6IAN+c/UwTr7ejoZiI+INhVhWuBFLF6xDXn6uFO+eYVN8OhNgAm4QIM0jCkEjww0ZcMgryGw2y+zVtBlOhhqKtCDvIdJF8kd2b/KIIs8n2owjYxZ5RFFfL4tnP/pNxiZaxxYWFsr1FfUzmAoblIJptgNsrPSPce/evVI8+stf/jK++c1vIisrK8B6qY7u0M3v3LlzuO2226Se0g9+8AOsXr1aHYPjUTABHxO4EoZGu/wu4Xp9Hn/d/lf87g+/xcZ/NGPTg4mYtSgCRpNYwAiPIyUZka6FkYxL5MFUXTaEt59pFpmUDMhLWI8HHngA6264EbExsdKziQUsr0WQ32cCnidA9yHaSW9pbcGuXS/jN//zeziiavH5/0hH0cYYKeYv71UU58pFEpgU7+7vHsf5w704/Gwf6o5G4B/+4Yvi5x+Rm5MrBXknxXkZGxNgAt4lMDo6KjWLLl26JCMp5s+fj4KCAkVFqdB9mELkyHOKvKko9I7Ws+TVFCyFDUrBMtMBNk4yJv31r3/Fv/3bv+EnP/mJzOBGsalcvEeAtJVOnz6Nf/3Xf5U3vB//+Me45ZZbvNcg18wEVEpgSLg3k9Dt//vjn1HbcwJzNwrR17vMiE3RIyxCJzKxkYIH6R+pA0CI2OGnuDynfULoLNnR02RHzRmRTentceEakYfbNt+Dez/7KaSkJqtjwDwKJhDgBGgRdub0Wbz88kvYd+R55Kx0YdmdochaKMS2LTqZJZL+3cpbkEruQ56aEuIiHEaFeLcTIyIrXGeDA2deG8HRp6xYt+xOfOM7X8GixQulJ4Sn2uR6mAAT+DAB8kqiTahjx46B9IsWLFiArKys9zSNlLZBRWLe9EMeVSdPnpReVGQYo5A8pY3lwzM1tVdsUJoaJz7KgwTIpZEykFEmMtJL+t73vieFzjzYBFd1HQK7du3C448/DtKuIuMSZTIIhpvddZDwR0xgSgSam9pw4vg7OHh0L7rtF2BMbUVygR3ZC8OQmmuRGkQkgyS1kKZUo7IOChFiS8JZCQ4hdtvXZUNz5QjaKkPQVREO7UAOFs5aj8994V6kZ7AmibJmlnurJAJ2+zh+/cs/YOdrT0OX0IrcG8ZQcEMEsgrCYTbpYRdmJJ/osSkJ2kf6SkYlEu8WvzBudwrx7hFcOjWIhlMGaHvzsXz+VnzjWw9JnZSPnMovmQATmAEBMriQEYm8esjgQqFtpEtE2dY8qYM0gy7O6FTavCcRbxoThe5RxjjSraVwONJ8UmvRCi+FH6t1cDyuwCNAmkl/+tOfpFA0ZXD74Q9/iIyMjMDrqIp7RFkJ6AZ+/vx5HD58GDk5OdKCzi7eKp50HtqMCNTVNIo0s29j7/4XcabmVfQbTyJuQRuW3GZG4doYETtvkuFttOutFq+kqwIjY5n4CRFWJUu4HknZYUjI0cBl6kb3UD2a2+pQfq4GHS39QvzXhKjoKKElcMVb66r18ZtMgAlMiwBpeDy2/Wm88OqfMBZZijV/Z8Kau5KQJgzaWp0QoBYq/6q+B02L1vUPJk50P6N7VFy8CbkLRahy9Chqampw5mgNQpyhmJU7S4StmK5fEX/KBJjAJxKorq6WhhbSSiIZDtJBItFtWgNGRkZKA8wnVqKAA2gtRXpPtM4i3SUKgaOxksGMGIyNjclQODI4qamwQUlNsxngY6EwtyeeeALPPfccKEb2u9/9rhRaC/Buq657dLMjgTu6kb/00ktSTG4yHSYblVQ33TwgNwmQ63JLc6sIEy3Bmwdfxr6jT6Fm8A3EFTXjps9FY/FNCYiKvWJIctCqJJjKuwsxGnKoWY+c/CjMWSs8tCJ7sW/P26ituYzenn5YR+3Ca0sPk1iQqe3hKZimm8caGAQok9Dru/fif37/M4TPqcc930/Fmk1J0ISRSCwbktydJTIsOek/wm0pOd2C+CwD6muacfD1MsTHpCAnO0eEvxncrZ7PYwJBS4AMKhSeS0aU0tJSqTNE+kJLly6Vm9lkZFHzuoMyv2VnZ8sfg8GAEydOyAxx5LE0KdqtlqzbHPIWtP/MfTtwskhTprEHH3wQJMD9rW99S/4D820vuLUPEpgU6t6yZQtWrlwpw9/oNxcmEMwESMSWQkBLz5Xib9sfwyuvvoCc9QNYc0805q2MQniUXjVi256eZwoh0Qp2lacH8PqfGtF8OhxrFm3D/Z+9D8uWrZC7kJ5uk+tjAsFAgBZl5CX5ta//f8jfNIg7v5mA/MJIWMVOPxfPEaBQOLr/97bZ8Oi3L2OkJgNf//IP8eV//pLnGuGamEAQECDPHNJGampqkh5JmZmZSElJ8UuGtkDBTWF+VVVVaGhokJng8vLypOTLpHEpUPrpTj/YoOQONT5nWgRoV2379u340Y9+JEPc7r//fmRlZU2rDj7YOwQo1pcs5v/+7/8udFEcUiR98+bN3mmMa2UCCiDQ3d2Fh/7pm9i9byeKbjbjvh/FIyZNi9BwIXQrxLYpS1Awpt6eytRNLsbsQuzWZnWiu8mBS8fGcXbvKCIdC/D1h76LLbfx/WUqLPkYJjBJwCVuOH/6w3b8x/95GEvvcWLLQwlIyRVhWML4wfeiSUqe+033MQpfHu4bx9/+tR49ZWm4984v4Yf/9j2BXHzIhQkwgesSoGxnlASIPJOKiopkWBt56Gi12uueFwwfktcWbVxS+DJltjt79qzUE160aJGi9Ww55C0Yrl4/jpGEn//4xz/ixRdfBBkqHnroITYm+XE+Pto0uZrSrkFERISM7T1w4IC88dN7LNT9UVr8Wu0E6usb8J2vP4y6kX1Y/Gk7bhShbXNWiFh4ix5avYYzJk3hAqDIEZ1gFRaqQ3iCDhFJLsRm2XH25GVUn+8UosHRmD0nbwo18SFMgAmQJ/Ejv/x/eHLHb5G8vA23fzUZqflh0ApvQDJ6cPEOAfFoJNKW64VGnF6E6XTi7LE6DPe4ULRoAYfvegc516pwAiREXVZWBkq8RKFsJKtBHknJyclSbFvNoW3TmToyqtFPaGio/ImJiZHe2+Xl5dIIR7zI+Ka0NRgblKZzFfCx0yJAnkkkwL1v3z7k5+fjBz/4AYe5TYug7w4mt0u6udGXAQl1k3A3fRHwF4Dv5oBb8i+ByspL+P3vf4/DZc9g7ed1uP2fM5A52yKNSKoX2/YwejIqkSYJ7fRHRhqQPiscMNlw6UI9ai50IjoiETm5WYp7YPIwJq6OCVyXAO3uP/n4s3j6pUehy6jEp7+TjuxC8W9JFDYmXRedRz6ke1hiUhjMCUBneyeKD1+GPsQstF9ypS6cRxrhSpiACgiQ80BdXR1I3oR0g8goQkYl+pv1E689wZTVjgxvtP4aHh4GJa6i37T2IqMT6SspxbDEBqVrzzN/MgMCZEx6+umn8dRTT4GyuZEANxmVuAQmAbp5kRGJLOXPPvssurq6EB8fz9nfAnO6uFceJnCpsgpPPS0SBrz2O2z5mgVr74kXD0JG2ElsW/yfi5sEBDtC6BIQcwsi4dSN4kJJNS4UN2NWTj4ShDgnZ4Fzky2fpmoCQ0PD2Lf3Tfz3r38Ky+xa4ZmUhLmLo0UWN0olqeqhB9TgnAJ2YmoYYjK0qL3UhCP7ypEQm4rsrGwYWKg7oOaKO+N7AhS+RaFbpJVE3pSUNZp+yBjCZeoEyHBEYuUk4E1hcLQmm8wIRwYnJRQ2KClhlhTWR7KuvvXWW/inf/onbNu2Dd///vfZmKSAOaQbGN3MVq9ejR//+Mcy+xtlg0tNTVVA77mLTMA9AiQY+fvf/wFPv/IIbv92JG66PwHmSL1Mv+1ejXzW1QiIHFTInCvCB+NceOfgRZw8UIlFRUuQlJR4tcP5PSYQtARIzPbwocP45re/hqSl7UKAOxEFS6NhYwFuv1wTdO+KSjBh1iILzp6owcG9ZxERFodFi4v80h9ulAkEAgEyIJFW0pEjR+QG9IIFCxAefsWDMhD6p8Q+kGGJIkbIc+nRRx9FZWUlCgsLRfitJeCHwwalgJ8iZXWQPJMee+wxmcWNQtwooxtZq7kogwAZlcgavnz5cpmVb//+/dKgxHOojPnjXk6PALlpf+OhH+BA6ZO4/TuRWHdfAkxhQiuJPQCmB3KKR1MIXGyqCenzTGKn/yJOHa5EalImsnMyp1gDH8YE1E2AjBePb38a3/vhN1B4+yg+/b1kZMwzCy8/4ZjE9yW/TT5JcZssOixcH4WqC004vPc8BrrHseaGVX7rEzfMBPxJwCXibouLi7Fw4UKpjUvrBy6eIWAymVBQUCC9vy5cuIClS5d6pmIv1sIGJS/CDbaq+/r68Je//AUvvPACbrzxRnzta19jzSQFXgTkqkqi3GQRp/SWhw4dkl8WGRkZionlVSB27rKPCZCL9kNf+jZqh/dhzQNarBNhbuHCM4kWbbxw885kEFfKlBcZr0dirhZnjlah6lw7ws2xyMvP9U6jXCsTUAgByvzzyK/+iMeeewTxS1pxhwhzS59jFuEPJMDN1iR/TyOtly3hQqg7W4+21naUHKvH2OAEFiwsFDoxOn93j9tnAj4hQGFuHR0dUjMpKytL6q2SZw0XzxKIjIyUG/rE+/XXX5ecKYFSoBY2KAXqzCisXxRDSwLce/bskR5J5J1E4VJclEmARODI7ZIyDZSWluLo0aPydVJSEgt1K3NKudcfIFB9uQa/e+T32H/2Maz5nBabPpeCqGihmeTkRdsHMHnlTzIqURa49GzhdaG3ory0BjXlnUiMTUNmFhutvQKdKw14AhTm9uxTO/DEC7/DRHI57v5uOnKLIoS/EgtwB9LkucQNLDk5DGFxE8Ko1InTR6ph0kXKzVOTyRhIXeW+MAGvEKivr0djY6MU3SaZDPZM8gpmWSmFvpG4+alTp1BRUYHY2Fj5470W3a+ZDUrus+Mz3yVAcbQk5Lx9+3ZpTCLNpLlz5zIfhRMgTyUSUqeb2ZNPPimzD9BuBBmVuDABpRIgY9LTzzyJPz/1f3H3w7G48TOJiIw2YJzEbrn4hAAZlUjsdnZhJBxCqPv86VpUlLQgb9ZsocUQx4KePpkFbiRQCAwPj2D/W2/j//z3TxCWV4M7v56MeUuj4CADN1mUuAQUAbp3JaWbpVB3dUUDjr1ZgYS4VOnJTZtwXJiAGgnQWo8M37W1tTJaoaiINcR8Mc+UAY7WYi+++KLMmEcZ9MjQFGiFDUqBNiMK7A/pJt1888249dZb8fDDD2P27NkKHAV3+WoEaOeBPM1IqJvmlrLAkb4SpwG9Gi1+L9AJNDY2SaHDp17+De6QmknxUoBbLtwCvfMq7B8tzDLnCaHuWCdOHarA2eP1KJg3X3gAsNFahdPNQ7oKAVqgHT16DN/5zjcQu1CEuX0zAYXLYoQANxu4r4IrYN4iravoBCNyhVD3mePVOLS3FNHhiVhYtCBg+sgdYQKeIkDGpPb2dhnqRhvLnLXbU2SnVo/ZbJZrL/JSqqurk2uxQAt/Y4PS1OaSj7oGARK1feqpp6Qr3he/+EWsXbuWdXauwUqpb08KdcfFxcl0luPj45g/f75Sh8P9DlIC3d3d+N8P/xQHzz6Fm/7BgA0PJsJk1rJekp+vB9KHiU01IilXL4S6L0hdpYzUHBH+lu7nnnHzTMDbBCbwzJM78C//+m3kbujDp76bhOxCCwtwexu7h+oXygAIDddi/g0RuFTWKDyVyjEy4MLKVcs91AJXwwT8T8DhcKChoUGGXZHDgFLS2PufnGd7QEal9PR0nDx5UnorkVB3IGXVY4OSZ+c76Go7d+4cfvazn0lD0l133QUyOnBRHwEKf6Mb2cGDB0FGxDVr1kh9JfWNlEekRgJdXV34+j9/Hxe7dmPpZ1y48bPxiBBhbhR6RT9c/EeA8BsMWkQm6BErbEinjwih7vPtiApPQO6sbP91jFtmAl4m8Ltf/1mE3v4KMQsacftXE5FZYIFOGwInC3B7mbznqteI+YoQyRzihVB3c2MbSo7UYXxUI4S654OFij3HmWvyDwHyoKQoFLqWSb+Hwq1oPcDFPwQo1I0SJtFGP3mM0ZqbXgdCYYNSIMyCQvtAhgUyMJBo8ze+IVLcFhbyF6hC53Iq3aabFrlatra2ypA3dnmdCjU+xt8Eamvq8T+//A32lzyOZfe4sFEIcMfFhQoBbg4p8ffcTLZPRj29UQh1z7JgwmBFWXE1qs93Ii0lSxiyU9nrdRIU/1YFAQofeXHHTmx/9jcYjynF3d9JR/6SSHGdC20xvi0pa47FvYumLEUIdYfGutDc1C68LKsRaY4TXpZZMBpZU0lZE8q9nSRA9ynSS+rp6ZGZn+Pj49mYNAnHj7/JQywtLQ1NTU04ceIEyHMpELzG2KDkx4tC6U2/8847MvsXeavce++97LGi9AmdQv/JOk5GpQsXLmDz5s2c3WEKzPgQ/xGor2vA008/hV8++lPc928x2PB3yYiONbIAt/+m5JotS28xoUsyrzAaNs0ISk5cRvmZRsybUyB3RnlX9Jro+AMFEZgQFzplSHroy9+AJq0cn/mXNMxbHi0FuNlbUkET+ZGukh5ccoYQ6k4PQVV5A0qO1mHp4lXSo+Mjh/JLJhDwBFwulzRY0L2KvnszMzMDvs/B1EEyIlHCJJKcIW+l1NRUv3sqaYJpAnisniMwMjIiPZMolOSOO+7gHWTPoQ3omubMmSPT45KrZXl5Oex2e0D3lzsXvARaWlrwhz88ikef+E/8w89TsPbueFii9WxMCuBLghbUoy4HNnw2FXd+LxJ1gwfxk5/+WBiwywO419w1JjB1AvTstHfvXjQ2N6BwfTiWrIhjb8mp4wvYI+ne5RCL8JwFkXjg3zNhtVxE6fl30NPdF7B95o4xgasRIGNSR0cHzp49Kz1fVqxYcbXD+D0/E0hMTMR//dd/SUMSzdXo6Khfe8QeSn7Fr9zGSRTs9OnTMoX8pz/9afZUUe5UTqvntFNBO6ydnZ1SoHvRokWglJZcmEAgEaCHof/49//EgbPP4IYv6LDxgQQhnqpjvaRAmqTr9CWEhLqThVB3jg5H37ogPJWakZqUgYxMFuq+Djb+SAEELl+uwS/+7/8g79YurPxUOKITjXCya5ICZm5qXdQKTSWdMQRjY1bUFI8jMSoXObns3TE1enyUvwmQMamtrQ2HDx/GsmXLkJeXxw4D/p6Ua7QfImKkKWqkoKAAZWVleOGFF2TCJPJe8kdhDyV/UFdBm3v27AHttN1yyy0qGA0PYToE5s6dC8ousHPnTpCHGhcmECgEyNhJD0P/++H/wNn6V7DoLjtuuj8eZmFMosLrtkCZqev3wyVEiUPNOsxeEY7bvxmJ2oED+PXvfoE39u6//on8KRMIYAL9fQMoLTuNQf05LL8rDGn5Ztj5phTAMzb9rrnEfOpNGhRtjMbF5oMoKTsx/Ur4DCbgBwKkmUSSFmScoM1i0ukhowWXwCVA80M/S5YskeG1//mf/4n6+nq/dJgNSn7BrtxGacFGNxy6YJOSkrB8OadHVe5sutdz8kiaN28eYmJipFGppqbGvYr4LCbgYQJjY2NyZ+3NozuQvrIfG4RnUkJCKOuTeJiz16sT4SOU6cokjEorb0/Ams/pUdW7D7v27IDD4fR689wAE/AGgZrqOhw48joKNjmROd8Ek0kHMp5yUQ+BCaHQTR6WSdlhiJszhLqOEpQUl6lngDwSVRIgz6SGhgaQVABlDps9ezYnWVLQTGeJBAC33nqr1FU6fvy4XKPTnPqysEHJl7RV0JbD4RAit0+D1P5JjJuFUlUwqW4MgXYuvvKVr+Avf/kLXnvtNeHePeZGLXwKE/AsgY72Ljz//PPIWe3E2s/EIyXRDKvYdeOiQAJine1yTkivsjs/n4nF24D20XN45/gpsQj37YOSAulxlwOMgMPuxPmLJTh3eT9W3ilSPUcZRKgbX8cBNk0e6w6Fvm24PwO9jgrsfGmXx+rlipiApwmQZ1J1tcisKn4ohIpC3bgoj0Bubi5+9KMfobu7W+pfkfSDLwsblHxJWwVtDQ0N4Wc/+5mM2dy4caMKRsRDcIcAxeh+5jOfkdpZlO3v4sWL7lTD5zABjxEgY/f58+fx6mu7sPIuC1JzLRgDG5M8BthPFVFE0LhIzL14QyJiZg3i2Wee87v4pJ9QcLMKJnCu9ALOXz6O+PlDSMoNhU6vAduTFDyh1+s6eVgKY/iiJXEYwGXsfWsP2tvapf7k9U7jz5iArwnQ5gwZHihRgMFgkHIWvu4Dt+c5AhRFdP/990tJmkuXLnmu4inUxAalKUDiQ64Q6O3txY4dO0Ai3IsXL2YsQU6AvNO++93vor+/H/v27Q1yGjx8fxOovHgZR97Zh5v+PgxZhaEwhtKCjcNJ/D0vnmifwt9y50YiIqsfh06/Kl3zaVeVCxNQCoGnHnsB+448h9u/kiFCooSeG9+blDJ1bvfTLgzhq7fFIueGXjz6+7+yHo3bJPlEbxGg79GmpiZpSCJ9VC7KJkB6ShRBFBkZifHxcVitVp8NiA1KPkOt/IbIHZLCSe6++27k5+crf0A8ghkTuO222xAeHi5E/M6jqsq31vAZd54rUBWBd04W47X9O3D7Q4mITTGCTEnk3cJF+QRoHrXQIG2uAdmrxrBr1yvo7eF03Mqf2eAYwaXKarSPlSJ9kRUZsy1y0HxvUv/ck+D6vJXRyFg6htffehaHDh+CzWZT/8B5hIohQMmV3njjDWmAoDT0XNRBYPXq1VID6+DBgz4bEBuUfIZa2Q319fXhwoULsNvt0pIdERGh7AFx7z1CIDU1FStWrMDg4KBIWbnDI3VyJUxgugQa6ppxsfokDMntyC0Kh8Gk5XCS6UIM8ONtYrc/dZYFK7ZZsPP1Z4ToZF2A95i7xwSuEPjbn5+CJq4RN96XLMyiXIKFAAmuh5uMSJ2rhTapAY8/9hiaGluCZfg8zgAnQNEFxcXFMpsbCXFzUQ+B2NhYUAKlepFAizyVfFH4u80XlFXQBhmTysvLsXXrVtCFyoUJTBKgzAL0ZfTiiy+CNLa4MAFfE3jl5dfQOlKMe76ZLrODceYkX8+A99sjge7wMAPSC/XoC7mAk2eOoaujx/sNcwtMwE0CpGfR0tKKA0eFPklil/BWiYWDQ93cpKnM00j/LTHLjDWficKbB18TepPlQl+JxdiVOZvq6nVnZ6dc133xi18Eeyepa25pNLQuy8zMRFtbm08SmbBBSX3XkMdHRBm8zp49Ky/KO++8Uwq3ebwRrlCxBCi96KpVq6DX6/Hkk09ieHhYsWPhjiuPQHdXD/a8uQsdtnNYKYSbKZSEw0mUN49T6bFdBDKGWfR48EfZePn1J3D82MmpnMbHMAG/EBgZGcXrr7+OxEVdyCrSwhAiPCc5DNcvc+GvRmlzIyxCh9wloTBn9uD4qQOorqrxV3e4XSYgCVC0SVRUFG6//XafGBsYu+8J5OXlYeXKldJLiZLWeLuwQcnbhFVQf2lpqdhla0FOTg7mzJkDEmPmwgQmCeh0OnnTKiwsxCOPPCJTVk5+xr+ZgLcJPPXE84ia3Y5V22JETjderXmbtz/rp8WZ3qhB/tJwjIdfxoHju3DhfIU/u8RtM4FrEmgVz0279+zErNUhyCyIgIPvT9dkpdYPyICo0YYgJtmIdffEY/uTf8Bbbx5U63B5XAohcObMGRw9ehQkW6HRsClAIdM27W7Sep08ZRsaGrwu0M1X0bSnJ/hO2L17t8zktWXLFs5SEXzTP6URz5s3Dxs2bBDC3FU4cuQIeylNiRofNCMC4kGd0t0eOPYqYvL7UbAqBuTBwkXFBMT0iiQmCIvUYtHWMBRffAPPP/2SigfMQ1MqgYH+IZSWlaDLfhYZhRrEJJqEJwDfn5Q6nzPpNy3ojGHinrUxGroYK06XHkVjY/NMquRzmYDbBEZHR6U8BUUVkM4OF/USMJlMUiOrpKQElKndm4UNSt6kq/C66UuwubkZly9flmkIyXWOCxO4GgGj0QjKKkAZACmrQGVl5dUO4/eYgMcIjI6O4dDBw7CGVSEx34noWLFgEzo7XNRNQM6w+M+iDbFwhjeh+NxhsThrkrtw6h45j05JBOrq6nHw2F7krB5HbKpOZikkbxUuwUdgQkgmhYjVVlKOCatut6Cm9TRefGFn8IHgEQcEgYqKCumVRHIVXNRNwGAwICUlRc63t8Pe2KCk7mtpRqOji++ZZ56RYm3r1q3jULcZ0VT/yST+9vOf/1xqbZ07dw60C8KFCXiLAAlKvvr6S8hb40RaPu2yCdcVLuonIBbltDCPSzFiwXoR+hZah+eeeYENSuqfecWM0GF34sLFUpRcegsr7ogXGjp6OMiqwCWoCWj1Gqy5KwnDmjrsef0Nr3sMBDVsHvxVCbhcLqmpQ5vAJGPCRf0EyAtt2bJlMtKIvPq9Vdig5C2yKqiXUsH/9Kc/RUFBATZv3qyCEfEQvEkgRMSipKWlSZ2turo6mT3Cm+1x3cFLYGzUivLKc6juOoicxQZEJ4bygi3ILocQYUBcenMiENOEX/3yN7DZbEFGgIcbqATKL1SirOo4Imf1IDXfJHS/hI4F25MCdbp80y9hBCcP2vlLIpFaoEN1fSX27NkLb3sN+GZw3IoSCDidTlRXVyMsLEw6CpD+KRf1E6C1WUJCAl544QXs3bvXawNmg5LX0Cq7YrJi9/T0yFC3yMhIZQ+Ge+8zAiQAd/PNN6OpqQmnTp3yWbvcUHARqL5chwOH9yF3tQtxWXroDRpesAXXJQCncFNKyglFwQ0RMCUMiAyTT2BwcCjIKPBwA5HAjud2Yf+xF7HlS0LwVggyT7B2UiBOk1/65BStLt8ajexVI3juqZ3o6ur2Sz+40eAjQAal06dPS+MChUFxCR4CZERcunSp1M46fPiwVwbOBiWvYFV+pWRMevPNN6XQ8qxZs5Q/IB6BTwiQ7lZ0dLT0FhgeHvZJm9xI8BEou1CCoyWvYPEtETBH6liKO/guAWlA1AthkrzFkSi4UYe//e1vaBJaSlyYgD8J1FTXo7G3GPHzBpCzIEJ2hbWT/DkjgdU2JY6YszwKS+7UohdncabkNEZHxgKrk9wb1REgY9LQ0BBINyk7Oxtms1l1Y+QBXZ/AmjVrsGTJElC4ozcKG5S8QVUFdZIxoLy8HHl5edKarYIh8RB8RIBCT8jFkt1pfQQ8yJppqGtGRe0paBNakDnfDJ3QpWAx7iC7CN4dLqVhT8kNQ9GtZpRVnEJx8SkMD7EhOzivBv+PmjZUHt/+FByR1VgrUsTrOB23/yclwHpA31URZgPyV4Qhe80o3njrNbS3twdYL7k7aiPQ398PEuOmNV1MTIzahsfjmQKBuLg4KUtCa7SWlhaRddSzcdhsUJrCJATjIRS6ROkG6cLjGO9gvALcH7PFYgFlFhgfH3e/Ej6TCVyDwJ7db6C24zQ2PpAKvYm/wq6BKSjepsWZOVSHvKUWpBQAb769R2yEVATF2HmQgUeAElHs2vUaXJY2FK6Nh8PDD+yBN2LukTsExuGSun8LN4TjROnrqLpcIZ6zKRiOCxPwDoGBgQFQshxez3mHr1JqpXVZTU0NduzYgbExz3pG8tO4Uq4C7icTUAAB8kwil8qMjAyQMHdxcbECes1dVAoBehiiBVt1+0msuSVJhj1xOIlSZs87/SQvJZNFi21fT8bZi8fwzvFTbMz2Dmqu9ToEyDuJEpmMOvuhCbXDKLyThA4zFybwMQJkCDeZtEjKMyAkrh3nK0vQ2sxeSh8DxW94jIBG3I9IRyc8PNxjdXJFyiOQm5uLoqIiKc7d29vr0QGwQcmjONVRGbnBVVZW4ujRo9i0aRPS09PVMTAehc8IkIcbxWzb7XaftckNqZ/A49ufQcTsJqy9N4YXa+qf7imN0CkEj42hWizdEoO5W4ZR130KJcXnp3QuH8QEPEVgSAjCb9/+Nyy41SG8kyLF/SkEfJPyFF311UNZKikiksyOE+KZmwySXJiANwg0NDSAjAd333039Hq9N5rgOhVEgIxK3/72t2UIpCeNSmxQUtBF4KuukhdAd3c36uvrpacJ34B8RV497ZAhiQxK5LHEhQl4gkBnZyeef+kJRM3uwrJbEkDiplyYAF0GdJsxh+uw/LZIdFhLcfjY20JXy7P6AEyaCVyPgEN837W2tCI8fgKWaD3fna4Hiz+TBCgictwqDOJGk8hUygt9viy8Q6C5uRn0/BQRcSVJgHda4VqVQoCug9WrV6Ojo0Ou9T3VbzYoeYqkiuqhnRIyCJAeAHmacGEC0yUQHx8vxdxZmHu65Pj4qxGw2x1oam5CQ0cZzImjSIgMYyHuq4EK0vdoc98hwkjyFkXDaWnGxdqTaGho4l3/IL0e/DFsrUaL0NBQOGziWhyX/kn+6Aa3qRACwidJJJQIQe5CC2rqK9Hc1KKQnnM3lUSAkiuNjIwgKSlJSd3mvnqZAGX5i4qKkgYlyuruicIGJU9QVGEd5JVE4spcmIA7BBYvXoysrCy0tbW5czqfwwQ+RMBmteFy1WWkzNEjMs4IJ+//f4gPv7hCwGIS31uJdow42lBf18BYmIDPCIyP2zA6NgKDKUR4m4TwHcpn5JXZkEtYwSlU96b7k/DW4V0oPnlGmQPhXgc0gZMnT8o08aSbw4UJTBIgx5EFCxZIx5Gurq7Jt2f0mw1KM8KnzpMpVCkyMhI33HCDzNalzlHyqLxJgG5UlJbyF7/4hTeb4bqDhEBfXz+eeeYZzN9gQPpsiwh343CmIJn6aQ1TBkG6hDpJiEZ41/LjzbTg8cFuE6Bnptq6Wpw4cwAZi4DkHDOc5DbHhQlci8C7obrGUA36B/qkF8m1DuX3mcB0CdA9ibJ5JSQkIDExcbqn8/EqJ0Ai7ZmZmUhNTZVi7Z4YLj9xeYKiyuooLS0FWbW3bdsmXbhVNjwejg8IUKgkpaekrDdcmMBMCdBuCmm7iagSIWTKulwz5anm8+12EUzi0sJoMql5mDy2ACLgFHpdA4MDaGythzkWCDPq4RJi8VyYwCcRoKvEKe5ZpKfEhQl4koDVagXJTrD0hCepqquuS5cuyYzcnhgVG5Q8QVFldZAgN6V8T0tLY4OSyubWV8OhTIG0Q0K/uTCBmRKgEFyd1iC1SSjlMhcmcDUCDuG5tnBdHPTxbXhs+xOsoXQ1SPye5wkIgzd9342N2RCinRCZu9jo7XnIKqxRZBMwmDSIz9IhKt6gwgHykPxJgLxQKDEOJ8fx5ywEdtvkwUab/42NjTPuqG6mNdDOMaWd6+/vFzvHbJ+aKU9/nk9zGR4eLsRMG9De1oGa6nqkJFVLo5LT5fRn1/ze9oTYbQwzhyIiMpyNbFOYDdoRoUwCdD2RICAJwHFhAu4SaO9oxYizFVFJGpjMWtYncRekys+jMKNZRZEo29+MPU/tY4OSyuc7UIZHIZYGkaUrIlYnPShFIvhA6Rr3I0AJiFvVleyUkTpkLTTAig6MDtsQZjEGaI+5W0ojQBu6FO5GEiZcmMDVCOTm5kp5kqGhoat9PK33ZmxQqqyslOFRZN0yGNjCPi36AXgweQKcPPkOGpvrcPj4PjR3VIPemwhyTxPKepebnY8Vy9cgPz8vAGcusLoUFhaGWbNmge4PFEK5bNkyvj8E1hQppje9Pf0oKz+DLttZZCzIQWSUCXZOCa+Y+fNlR2kn1m4TIW9Orcd0AXzZf25LmQRs1nFYx4cxe7kJoRYNK7wpcxp93mtxu5IbJIk5elQ3nsfxYyew6ZabfN4PblB9BMiY1N7ejqVLl7JBSX3T67ER0eY/GRyNxpkbsmdsUPrxj3+M6upqFBYWSi8E8nLhokwC9DDe29uH8orz6LNXo9PUjc5mdt2m2RwbG0XH0FYkxWeyQWkKlzddSxkZGYiNjUVJSYnMJsAG5ymA40M+RoC+U0ijhIInQ4R+Et+RPoaI37gaAX4UuRoVfs8LBGpr67HvrdewaluM8FIycBZKLzBWc5UUxm2zjYvQE7uah8lj8yEB0pw8c+YMcnLEJhx7KPmQvLKaouiypqYmGVFCMjczKTM2KLW2tooFdj6+853vIDo6mjVTZjIbfj6XhJQvVV7GX5/8LSr6LuHeh5OFEC4JnPq5Y/5s/t3Vq802gtEzwi15xOHP3iiq7ckQWNZRUtS0BWRnyajksAsLARsJAnJ+AqVTFJpsMGmRkqfH0Ep6vOELJlDmRs39GBsdE1nearBkwxVNHL1QUdIILaVgL07x74817z75KnA6JkDXDGem/GRWfMTUCYyNjfGafOq4gvJIyvRGmskVFRVYvHjxjBjM2KBEi8WYmBjk5eWxtsyMpiIwTraOOJG/MAHaCQtS08zCK4DVAGhmxsWXfUutARNjgTFPSujFpCg3uVKyKKASZiww+6jVCJfcaDOyC8OgN2rYRBCY0xQQvSIHab1Wg+gUHYwxQzh48CBWr1qL0LDAy/g2MjKKqspqsWEjOk2xL1wUSYD0AcvKzqOlqQ1RZX3Qho4jLFwnvCqD16CkEZpSRpHpLjrJiKhEAxuVrnNl01WSmBWKrtZBoWXSdp0j+SMmMHUCtAm3fPly9k6aOrKgPJJC3lJTU4UxW6RQnmGZsUGJFoq0cKT0hKGhoTPsDp/ubwLjNrv4cjNBE2uBdULsL9G3XfA+F703HeMirau4zLlMg0BKSgqys7Olxtqdd94Ji8UyjbP5UCZwhUBrazsqq8pQeKOI8xaC3GTk5sIErkVA+IZAbwqBzTUsw22LipYEnEGpra0dBw/tx8mSg9AZQqSQ87XGw+8HMAFhB9RAh8sVTejrHkZXeQwmBgxyToPZs5s0pRyaQSzYGIrlm2PFhhzfs692Fcvna2FQnr00CpffbkTxyTP44j9+jjfgrgaL35syAfJMIv0kih7i5+4pYwvaA0knOSA0lIJ2BlQ68P6+QUzorIhLFTu6bExS6Sz7ZljktbhmzRr84Ac/wFe/+lWQgYkLE5gugbJzZXjymcfw1e1RMIpwJipiAzw4y7vrMrkQCU4CUxi18KqlMO0JkRHQJDZHAtD75+yZUvzqkf+CK+E84tKNwgAhPO+CObR8CrMaiIdotCEY6LKitV4LvSYHCxM/BY0QhHcNB+dk0gazS1zI5K1VWvUmwuIGsXpznDAocbkWAbqXZ2ZGoGfoLNoulMoNevIa4MIE3CUwPDwsQ5iSkpJAP1yYwPUIUGb3S5cuYfXq1dc77BM/47vWJyIKrgOefvppNGr2YlvRzBXfg4scj/ZqBEgYkMqkntLVjuH3mMD1CNDDtU4nMk2KgzTCnqQhWe5gjRB6d9wusavNRqVrXzXk7k8/8r4TgAYlm9UhvO0msOkrWZi9MhzGUGGECE4bxLUnUQGf6PQhaLw0iFPPT8Bevgw/+emPYDKa4HAGp9biZMTC2ZJz+I9fVgkO5xUwi/7vInndhohdEp1u5mEn/h8N98DfBOjfIT9z+3sWlNO+p64XNigpZ8590tOh0R6MGQbEF1uyT9rjRtRNwJ+C3CMjI3j26edQc7lBZlDhL1hlXWtkMDGbQ3FGZAmsv9SF537Wh6gEnTAqhQStMYUWsAmZJqz6VCwS0oxiwfauy5KyptarvXWIxVlyjhmz19iw8687sWnjLSJhSKRX25xu5WTsQogQ4jUJrRlhTDIaOJRzugwD4XitsGwbxBzq9EJ8Wqt7T/ZBr9EHQvf80gfaAJCegRpyI+X701QmgSiRMLfLGSKevXlZNhVmfMy1CYSHhyMrK0uGMXlCG+faLfEnaiBQVFSE5ORk1NbWCm/JTLf1lPjOpYarwYNjSJ4t9CfC9QiPFR4B/CzgQbLBWRWlK12yZAnoC86XxWq1ofRcKZ577U8wxPUIQ4RePKz5sgfc1kwJhGhCMNjhQMNgL6IzRPKH5CvaJFRvsN2bKMTPbptAb5sdzVVWFKyNQEqa8ITgBdvHLjPy3oqINSI+245Tp05icHDwY8f4+w3aEcSECA8SBkFaSDoNnA3L33PiVvvkMSkM3AhxYnR0VGTFdQS9QYCMpU4H6W/yA+RUrykiRRqBrZp+7HplN267/Va3F3VTbfN6xw0ODqGvt1+YS4M1tvx6dAL7M/Jy6+3vQUNjPdauWesRbZzAHjH3bqYEaJ1GG/CXL19Genq62/ceNijNdCZUdn58pg76GD3MBh2sTvbBV9n0+nw4CQkJ2Lp1KyorKxEbGyszQvqiE5QkoLamGs19ZbjpLgPmrhDCoGN8PfuCvafa0ApvnMunhzDkHEaCJhQbvy/LC8EAAEAASURBVJAgqhaLcGEweDfyy1NNBXw9Gl0IRgccgscwjr/ci3ErMQg2ClOfJp0w2OgMwNDgiDDa8L/7qZPjI6dDgDyU2mpHYe2NxIIF81lMeTrw+Nj3CJBBKX9ZOKytwzjw9kFs2bLZ7UXde5W6+QcZ4Pfs3Y0LF8/CHG7iLQs3OfrjNHoioOySo8N2hBos2HzzZr4n+WMiFNbm+Pg4BgYG0N/fL56v3X9eYoOSwibe2911iF1w2gkX+0veborrDwIClDmAvJN2796NtLQ0nxmUaLGtN+gRm2hGTqEZc+fFwCquai7KIaATc+gc06KrVihMCCPSrEUWaUIJxjsTLVyHhx0YHXSI65oeG4VGkHKm0uc9nWSjFUalYNWz8Tn0IGxQL/5dntrTDfulWfj7Xz/IuiVBeA14asjjVpFJ2B5yJWzSj3sFJ98pxl+efATn6k4gt1DDnt2emmAf1EMak10tLujGYnD7ui9JgXfK4MWFCVyPAHlMk0epzWa73mGf+BkblD4RUXAdMNrvxJCIDRoXu7q8Ax5cc++N0VJM7q9//Ws8/vjjmD9/vjeauH6d4iZJoW5kIKWdGy7KIRAiHo6k+LQwJtGmifRMEl989DvoimBB1+/7m0d+XHEEOPzJq8Nk0SB3uQ7N7bXCU2kBwiPMAd5z7p4SCThECjMyBGhZUFmJ0xc4fRY3ritRgv69t3d19CA5T4s598Zi4/0ZwhvWfY+FwIEbHD3RG4VXd0kfyveEwtUfHGPmUc6cAGnemc1m2O32GYUqs0Fp5nOhqhqqS6xoN41i5V12RMYZ3v2CU9UQeTA+JEChZ21tbYiKioL/dkoml5g+HDg3NWMCZNAe7ndIQ0p0krgXyRp5LmcMVu0V0MJMjDE8WodFmyJwpvQ45s5ahAVF89Q+ch6fvwgIPSzegPMXfOW3S1neEjKNyBbe1LoqcSWRxpqfity4ERsYxlANQrVaaM3+64ufECi2Wb3QvDKZtVe8mPlRSbHz6I+OR0QIXc6UlBl52bJByR8zF8BtZqXNw4RI9WodtSISIl6ACxOYAQF6MJpJTO4MmuZTFU6AQt4uHBxBd5MTN94vsnTxA5LCZ9S33dfqNYiM16OvohejI6O+bZxbCy4CImMfB6EG15R7crTkfZqcakFL7BAOnHgLn6t/ADk5OTNa3LnbP2nLEt+1Tvbsdheh384L0YrIEqGjHqKdACWmYV18v02F4hom76Tm5uYZGbNZwl9x0+7dDq9ZswqLli6AS+QP4sIEZkpgbGxMpOyODvrMNzPlGIzn05dTXakNA+0usXMbJhHwA1IwXgnTHzPtqZMmAOkB6nV6DkeaPkI+gwkwAR8RICOAzeFCX9court6MGYdk7qBPmqem1EJAdJa7Giwor9Fj3kFc4SwOy/xVTK1Xh8GGZR6enqkOLe7jfHV5i45lZ6Xni6Ek6PiMTJop4RKXJiA2wQo1K2xsRHr1q2TIW9uV8QnBi0B0gSgDGf2cXZPCtqLwI2B09Wi0YTAEqkVniMUNskbJG5g5FOmSkBccPy4NFVYfNxHCeiERanqTB96LoXj9pvvwez82WwE/ygkfv2JBAwi5O3Mm32oPWrBA5//Oz/KTHxiV/mAACNAyZMKCwvR0tKCkZERt3rHBiW3sKn3JMqKU1U8iCPPd0EjfF/9GMqtXshBMjIS5D5//jwblIJkvnmYTCBQCFAISViEDusfTET3aB0amxoDpWvcD5URoAyUE1K3mE1KKptanw6HEi5MyES0/tVQ8umguTGPE5gQeycTlCRA6F9xYQJTJUCRJLT5r9FoxAacexmxWUNpqrSD5DhykbQJ4+RAH+/oBsmUe22YdFNyOBwwGMS+ibhJcWECTIAJ+IIAbYTQQn9EfI91tHVgoL9fZgckryUuTMBTBBzC/y1vaShCYs1w2NmL0lNcg7EerfDEhUaE6VpFdAAXJuAuAboNvXsr8qe4u7vd5/P8S+D06dMggW76mW7hVd50ian8eAp5mze3ECmJmWKk/ICk8un26vBoh4S+0EZHR1mY26ukVVz5Bx6OVDxKHpqHCYQIw5FtzIWKdwYRF5GBpGTKXsLGJA9jDvrqnOIZKasgAiHhXXjt1d0zSrkc9DCDGADdmXparRgfNGD27PwZCeMGMUYe+kcIkI4gFyYwVQLkAEA6SuPj41M95UPHsUHpQzj4RZjFhPy82SJ9YBrDYAIzImC1WqXbbX5+Psdyz4hkEJ8sdmxDhOc2O7gF8TXgxtDpwcYx7kJL1ZjYHMlGqjAocWECniZAy7XwKBP6hjrx9v63eePE04CDpD4yKNWVDWOg0Yw1a1exR3eQzDsPkwkEGgGKJnHXs41D3gJtNgOgP+M2u9jddYhwJUpByTbHAJgSxXXBJQQBSNyNLN533XWX4vrPHfY/AVqsRSRqYR1xYnxsAqZw//eJe6AMAnTtUNibzhAiBbkdbmoCKGO03Et/ESBDgEuEVpI3bqhZhL6x6KS/pkLR7dJ11Fo1DmdTGGbPzRf6pfzcregJDZjO05XFhQlMjQB9f+n1ere/x/iuNTXOQXXUaP8EGsvsuPTOkNzl5VCBoJp+jwx2eHgY5eXloHhcLkzAHQJ2uLDu3mjMWR6G8qODsgper7lDMkjPERcLZQgk4zbpKXFhAp4mQFeVMVQjswmGmkWIN4dVehpx0NRH320hIRNwOaXCe9CMmwfqWQIuoadMP2Qc4Oclz7JVe20mk0mE3M6GxWJxa6hsUHILm7pPmnBpMD6qwfCAkJyk7zY2cqt7wr0wOordJu8kd2NxvdAlrlJhBGixFmrRYFTch2pLR648HPG9SGGz6Ifu0sJMNNvbOo5Xf9WDrRs/g6IlhX7oCDepdgLknRSXFoqEfCeGcBlVl6rE9557GXLUzorHd20CIsgEOr0WGpEUhwsTcJcAbcIV3GjG0m1GdLZ3cwiuuyCD8DzSTtq7dy8WLVqE+Ph4twjw3cstbOo+iTSU1t+0XghMkjWJd3bVPdveGR1ZuikMwGazeacBrjUoCFDmpNEhJ0aEUYmK2HMLinHzIN0nMHmFUKhke7UD6anZMBoN7lfIZzKBaxAgzdtQ6KALHUf3QBMaGht4E+UarPjtqxMgL5KW9iHkZBRgw4b10OlYieTqpPjdTyJApuyU3FCEJQ3jwNsH2bj9ScD48/cIUFRJRUUFIiMj3dZwY4PSezj5j0kCubNysWrVarGIE/ZuDhWYxMK/p0GgpqYGdrsdCxYsmMZZfCgTeJ8AmbLNUTpExOmlFg7ZtyfYwP0+IP7rugQmZUhs42zUvi4o/nBGBOg+NW5zYnholDOazohk8J0sw9xEmGRlcR/iw3Oxdt3q4IPAI/YoATJIDg4OoqTkrNAPZG9Jj8JVcWUUVUJrNooscbewQcldcio+z+lyoK9rFNUlIyAPAb5IVDzZXhra888/D7J4f+tb3/JSC1yt2gmQMTtjXhjS54aKBZtLiHNzCK7a59yT4wsRFiWDUQ8tJ5bwJFau6yoESGdSp9MIr2726L4KHn7rOgTIo7Kt2oqu5hGxoGPj93VQ8UdTIEAbKTq9BqZQI2soTYEXH3KFAEWUREdHS+0td5mwrcBdcio+T6cH+jrseOuxYdhGXcKgNBlEoOJB89A8SoAerEkMlwsTcJcArc3CxJfc6JALF4+PoL3WBqdjYkZfeO72hc9TFoFxpx2hRjPuv/+ziE9wTw9AWSPm3vqLgJB8R1iEDjFpOrS1tgrj97i/usLtKpQA2bzJW4ntkQqdwADptlPci6ISDUibHQZzhNDk4s2UAJmZwO4GeSUZDAYRmbRKZnlzt7dsUHKXnIrPM5qMyEjPxPJFS2HQG4TMGxsGVDzdXhkaZZggoxK73HoFb9BUSppJFHVrt1H2G/GH+D+bt4Nm+t0aqE48RLdUjaL41UHMmzcf4eHuZSxxq3E+KegICL9JJGWbMecGAx5/4jF0dXUHHQMesPsExFealJbQiM0TnZ71k9wnyWfKJAEJYdBGDODNI7tQdfmSeAbn9RtfGdcnQAalvr4+tLS0zMgRgA1K1+cctJ+mpCSjaNEChJpJzJRvSEF7IbgxcLo5jY2NISwsDLGxsW7UwKcwgfcJuIRXklMYkwwm4SspvrFYR+l9NvzXxwnQQw1lBuyos8Es7kEscvtxRvyO5wiQQSAsXIfIRA0uX74Mq9Xqucq5JtUSkB5J4tF6qMeOqhPjiA2dhUWLF6p2vDwwHxAQNyOjSBLg0ljR0FwrjQTOGWji+KDH3EQAEKDNf1qzdXV1sUEpAOZDdV0IDRWZAsLC0dU4DqsQnJwUOFXdQHlAHiVAN6b6+nr5RUbGpKysLI/Wz5UFFwEKJ4mI1SE+zYC2Givs4xzyFlxXgHuj1epDoDeFSJFJl8xW6l49fBYTmAoBMng7xL3JYDBBp2Uvk6kw42MoxG0Cw30OVBePI8aUjYRE3oDj62JmBGjDbcIVIqRKRDITIdBN0QJcmMAnEaA1f1JS0ozCJNlD6ZMoB+nn4nsOVqFdcuKlHnQ3j0PHFqUgvRKmP+xWoSMxNDQkQk3CQTcpLkzAXQIizyRyi8xYcVsU9j/ZiZF+B7Qc9OYuzqA4jzxGSGvL5QiB0WiEhr+7gmLe/TlIWrNpdEKY2zCBcRZW9udUKK5tjRYQyhJwTrifXUlxg+YOe40Aff9ptCEwhIZIPTfWMvUaatVUTN5Jer0eK1euZA0l1cxqAA2ELq74uHixg6KR8d0B1DXuSoATICFA3hUJ8ElSSPfo4UhnIG8TjTAmOYWOEmsoKWTq/NJNsh1Z7U40X7ShuzIMN2/ejIiIcL/0hRsNDgKkW2IM1SJtjgH56604/M4baGpsDY7B8yhnSCAEWrH4dwgPNxeneJ8hSz6dCJAwd3yaCcvvtGD/gTfR1dnDYJjAdQmUl5fj8OHDiIyMnNHajT2Uros5eD+kcKVtd90FvVEjvuzs7BMQvJfCtEdOxkgqLMg9bXR8wkcIkLM2eUtOiEUbJQ2k3Vx24P4IJH75HgEScXfYXRjsdmGow4j0tHRxzfBjznuA+A/PExD3J8qEawzTIDRqAr19XRgZHvF8O1yjqghoNORB4sTpN3uxYd2tWLykSFXj48H4hwDJBJgjdUidbcD5C+cwODjon45wq4ohQLp/IyMz/87iJy3FTLlvO6o3aJGbMwvWhgxY+wziFsXC3L6dAWW2Rp5Jzz33HOLi4rB8+XJlDoJ7HTAEyEOJjNrmaJ3YyQVGhN6EXWjisCxAwExRwHWErg2HfQK2UdK0oaQSXJiAdwmQIZNuSuRkQlmVSBuHCxO4HgHaGHGMu1D5zgDmz1mEnNzs6x3OnzGBKRGgO49WpxFJTEKEMWlI3o+mdCIfFLQEaL2Wnp4+4/GzQWnGCNVbAQm7JWtvxEBdJLq7RqRrrnpHyyPzBIH29nbs3bsXc+fOxY033uiJKrmOICbgEguz0HCtcOE2SoNSV6MNo0OUJIAex7kwgY8QEJfF+JhTCpLGxEbPKGPJR2rml0zgmgRICJcMmZZIrdDt0s9I2PSajfAHqiJARkiyO1KWN5MhDCaTUVXj48H4jwA9HZGOkkYnkgU47P7rCLcc8AQoK7fFYkFmZuaM+8oGpRkjVG8FRqMBf/eFuzFQHYPyIwMshqveqfbIyHp6evDCCy9g/fr1WLZsGT9Ue4RqcFdCCbpMYrctMkGPECF629Mmsk4Ok8GACxP4OAG6LlrrRqCxRWPtDatnpAfw8dr5HSZwdQLkv02elHOWhWNCMw6blRdxVyfF704SkNm4yJ1EbI7QL/rhwgRmSoCuI7oXRSZoYdd0o6unjb2UZgpVxec3NTWhs7PTI+s1fi5X8YUy06GRuPKixQsQ6VyAnqpI9PWOcajJTKGq+Pzm5mb8+c9/xi233II5c+aoeKQ8NF8SkA/a4j8y+zttvXG8my/xK6otnTA1Xiruw1h7FDbfuokNSoqaPeV2ljTeSJh77upIPPn0Yzh29B3lDoZ77nUC9BXmFO5JYwPC29ZmQkxULMLCwrzeLjegfgIuIfBuCtMhJd+AYU01LtVcQDcLc6t/4t0cIWUBtNls8sfNKt47jQ1K76HgP65FYP3a2xCrnY+SA23Qygxe1zqS3w9WAr29vSgrK0NSUpIMdyMXSi5MwFME6AHcJERvDUbhxi2+tXg311Nk1VeP3SY0bJxCIDk0VH2D4xEFJAEKXaJ7lDlKLOQWD6GusxjVVXUB2VfulP8JkCD3yIAdjeXjiEQBCucuQVx8tP87xj1QBQFKEqDVhwjlWzvsjnHhoSTE3bgwgasQiImJQUZGBuLj46/y6fTeYoPS9HgF5dGrV69GbsIqtJ4Lw+iIXcZ9ByUIHvQ1CZw/fx7FxcW44447kJiYeM3j+AMmMF0ClLWEdJQWboxCd/M4hnod0JEILhcm8BEC9ECjFboRLjhEBqXxj3zKL5mAdwmQoWDtPXFoHTqHN/bu925jXLtiCRiEJ2VH4wiOvziATTd8Crm5eYodC3c88AjIcEoRh+u0h0Cn1YvvRJHRhAsTuAqBrq4ujI2NITw8/CqfTu8tNihNj1dQHh1q1qMgdzWSdCtRWzYAp8MlvAR4QReUF8NVBk3pJsmg1NLSgq1bt0qBt6scxm8xAbcIOEU4SWiEFos3R6Px4hjaa6xXck7yLcgtnmo+aczlgLVfj1BNDBISZ77jpmZWPDbPEiAvpQlxT1qwLg42Yx1OnzuM1tZ2zzbCtSmeQIhYdVknHGipcKDrfCy2bduG6JgoxY+LBxA4BMiDm0S50+aEwe4cFY4A1sDpHPckoAiUl5eD5Eo84dHNBqWAmtrA7UxhQRHWLbsLZ/YOYtzKabsDd6Z83zMyJjU0NEi3yZycHI+Iu/l+FNxiIBPQiHgSCnmjNMut1WPoarOxUTuQJ8zHfaNwI/qpPjuAwZoE5CWuw/zCeayh5ON5CPrmxEou0mxESqENw4YKHD50hDMNBv1F8WEAOmFRuniyF21lFtx+833COymXv8s+jIhfzZCAS2zCGUI1uO0rKWjtr0TNZQ6/nSFS1Z4+ODgo12ye0HBjg5JqLxPPDiwiOgwZCQXQNq+GbVAL14SIyWUPAc9CVmhte/fuxdDQED772c8qdATc7UAmIPVJhMe2JVaHpVui0dFgw4XDAxB53wK529w3HxIIER6zWvGz6w/1yI1fg3vu2+bD1rkpJvA+Abvwn5y/NhbRs/rw2usvoa+v9/0P+a+gJ2AS4W5Hd3bi8sFw/OM//z20Wg5HCvqLwgsAaIMlzKJHaeUBlJYfh33c4YVWuEqlE8jOzvaYTAkblJR+Nfiw/xnpWbjlxk+j7O1BdLeOQk/quFyClsCEWOnX1taivr4ecXFxWLFiRdCy4IF7mYDY+ScX7iW3RmOgy4nyI8OwOVho0svUFVE9hZCQ51pDxSjQm4nC/FWYMy9fEX3nTqqPAHkHJGdZED97HC3DJSgtLYV1jENO1DfT0xyRWODTd9jFil4YxzKxZslGZGdnsRflNDHy4Z9M4EqSgBBExhuQvMCKhu4SvHP89CefyEcEDQFav124cAHR0dFIT0/3yLjZIuARjMFRSWxcNNbdcCPqD8ehu4G8A4TqG5egJRAitkAeeeQRme5248aNQcuBB+59AvSARKnd4hKNyCwQugDjTlScGJJhTrQTxyV4CVAIiXXUgYPPt2Fh3kYULVzGISTBezn4feR0r9JBi7S5JqQuGsZTTz6N3t5+v/eLO+B/AhMipfvhF1uQEFqAT336DjYm+X9K1NsD8VykN2iwbGsMmofPYOfLO9U7Vh7ZtAmQEPfZs2dhNBoRGxs77fOvdgIblK5Ghd+7KgHaCU5PzUCqfj0Gai3o6R4R7rq8mrsqrCB4k7IDvPrqq5g3bx7Wr18fBCPmIfqbAGV8W/2paMSlGLD/sW447GL1xrcgf0+L39onY6JL/K+/w4GL+yewdvkmzJnNGZP8NiHcsCQwLiQBEjPMWLDJjIOnXsblmktMJogJ0H1qQuy/drdYUVesQUZsEVatYY/uIL4kvD908WjkFAbM/MIYhMS0oLTqAFqaW0WWbvEBl6Am4HQ6QdpJJpPJo0mU2KAU1JfV9AdP6Scf/OK9GKqLRcXxQbEPx6u56VNU/hm9vb148cUXsXnzZixdulT5A+IRKIKAQzwgzZkTLly5TSjZO4LuJhtcjgmQsZtL8BGgbKN9nVbUlziQn7AeBXOKYA4PCz4QPOKAIkDGA6Nei5R8A+ZssGP/oVdRefFyQPWRO+M7AqTxZrc5UXawH7mxa1E0bxV0ev7S8t0MBHdLy25JQPoSKx5/7CnYbLbghsGjB3kndXZ24vbbb0dKSorHiPAdzWMog6MijdBNWrZiMSz2+ei6ZBGCk1YZdhIco+dRThKoqqrC9u3bsWnTJumhNPk+/2YC3iZgJdHbmyy4+ctReOU3bRgZcILCnrgEHwHa0GioGMTZfaO49577kJaWGnwQeMQBScAhrErmCD22fiUVf338URw+eDwg+8md8i6BK95JExjsdKBk9wjWrbwVRYsXeLdRrp0JvEvA7nKhaHUc4ub04b//6zcYGBhkNkFOYGBgAJRMySWuDU8Wfgr3JM0gqmvD2jsRHTIfJQdbRXYdDRuVgmjuL1++jF27diEiIgLz58+Xv4No+DxUPxOwC5ft1DwTCtdFoPKdYbTXWGF3uNhLyc/z4uvmSeB2aGQcbZUajDdl4aab1iMiMtzX3eD2mMBVCZCXklYXguScUKQvH0Vt+2k01DVf9Vh+U70EyItydMiOhnIbdIP5WLxgNZJTEtQ7YB5ZQBGgCDet+F96gQHzt47hpZ070NvDmm4BNUk+7Ex5eTmOHz+OO+64Q4a8ebJpNih5kmYQ1bV69WqRnnk12srMIoOJQ8TlBtHgg3yoR44cwRtvvIEtW7YgKSkpyGnw8H1NYEJkUTIZtFcWakL49sKRQXQ126BnLyVfT4Vf26Nk29Vn+zDYEIV1y2/je5FfZ4MbvxYBvSEEG/4uEc39pdi7+61rHcbvq5SAARr0tI3h9GtDWL/6TuRk56p0pDysQCUwLry6k3LCsOEL0djx8pOorq4O1K5yv7xMQKvVIiYmBnPnzhUayPQU5bnCBiXPsQyqmsIsBszLWY0k3SrUXRiEU3oIsJ6S2i8C8k46evQohoeHpYXbYrGofcg8vgAk4BTi3JZYHdbdF4dLp4bQVDFGSeBYoDsA58pbXSKvtAsHhxHSk/H/s3ceAFYVVx//s72zsCyw9M7Sq0uTjoCFZg8aY0sUjX4xxRS/L81omomJRo0llthjxYIiKmBFpYMivffO9vq+8x+85LkuiOy+5b27/9HHe3vfvXNnfnPnvDNnzpzBRd85T8FGQwVa+Z4wAU60BUwt6j2iEfKj12DegjkWu2L3CeenCyOLAL0oD+YXY9VHhdi+OA1TpkwGd0tWEoHaJMDg3KlpcWjbJwE7ChdiyfJPsH/fwdosgu4VBgRyc3NdzKThw4eHpDQyKIUEa93ItEe3Xhjc+0zMn3HAAg5WwDx7lXxOYMaMGVi3bp0L5taxY0fbnlsixOdNHpbVqzAvpfiEKHTom4Lo2HrYuqoAu20HnWgJobBsrxotlP3OMC7J1rV52LMyFZnxvdA5u6O24K5RyMqsxgiYUSktOQHNe5egIGEF3p37fo3HrqixsiqjGiVA76T3nt+JeU8EMGLAJHTunF3jXgE1WmBl5lsCnHBLTI3C8G/Xx8dL38SihUt8W1dVrGoCXOr2+eefIyYmpuoTqnlUo8FqAqzLl6c3TEGrzO4oX5uDssIY8xCo2QBfdZltONZ9y5YtmDNnDpKSknD++eeHYxFVpjpCgDP/9cxbN7VhDPqObYAd64rw6fuHEGsKvJK/CTAmSawZsl/710Z0bHoqzv/W2f6usGoX8QRKTTfqcWoG6rfbj1dmvIADB/ZHfJ1UgWMToNH7UF4p5r++F4Vbm2LatVeFbCB37JLoWxGAGbEDiI2LQt8xGdi47yMsWvohCgqKhKYOEAiYwrxnzx5r74KQ1lbad0jx+j/z1q3aYOywKVj4xkHs21mAuGg9Un5t9VdeecVtNdmrVy/k5OT4tZqqV6QQMKNSlAW97Tk8DSVFAaxdXIC8wtJIKb3KeQIEom0JSUlROT6yQdrWBQ3Qs+MwDBjY7wRy0iUiUHsEOJhr1j4FmZ2KsPHAx1iyZCmKirR9d+21QC3fid76ZlFa+OZuxBa2wMjhY9ChQ3t5UdZyM+h2/yXgJuFsMiazZQKyehVg/a4F+OiDBf89QZ98S6CeyaLY2Fi0bds2pLEmNfr37SNUOxXLbJyBUSNHYfVbqVj9SSEKSmwnCxmVagd+Ld6FMZNmzZqFrl27YuLEibV4Z91KBKomQAWJgZMat0lAy+wkFOWX47MPDrnlUJwdVvIfATbrwV0lePW+bejWciT69xloAzf/1VM18hcByqpY22mpeZc4NOlxCM8++wwO7NdOS/5q5aDaWHuXFlXgrce32ZLc7jj3gkkKDxCERx9PEgH7reSkTM+Radhf8Rnmvvum81w6SaXRbWuJQHFxsQtV0qhRI2RlZYXsrjIohQxt3ciYGyu1atUW2RlnYceCTKwxo1JxSalt4S0t3y9PwP79+/H2228jPj4e55xzDgYNGuSXqqkePiDAH7Hep6UhJT0aHzy/D+WltDQp+Y0AA9wWFpZbAHbz8pgfgzPGTUJ2105+q6bq41MCpQHbaal1CnqNScLbH7yIteu005Ifm5pLcissCPL6pbnYuDAa7bL6YsipZvhWEoGTTcBUo3LzlmzXMx3xWbuwcuv7WPHZ54rpdrLbJYT3pzPAxx9/jLVr1zovJXoqhSrJoBQqsnUo39jYGPzu979C25hzsGx6Kj7/+JB5DmhQ54dH4NChQ5g7dy7uuece3HjjjRg3bpwfqqU6+IhAmSlIHbunon5mgnkoFWLvNtskt9ytOPBRLVWVGHNF2rmhAEvfKsKgPqPRs0dvJCcnCYwIRAQBLntLiI9Bs85xyOpzEHPfn4m1azZERNlVyOMnwLlUesu+dv929O8xFCPHDDXvpJrdnvv4S6MzRaASARuaJdiYrUNOAuKytmP6iy+hrLSs0kn60y8Etm7dioceegjZ2dnIyMgIabVkUAop3rqTeVxcLKZNuwYD21yKhU+kYtWCfVb5ADirrBSZBBjA7bnnnnOvm2++GV26dFEMgMhsSn+X2hQkmq+zBySj39gUzH5yFwoOlWnHN5+1epwZlNYuzsPb9xfYpgAXoFXrFj6roarjdwLlNtGW0iAOE69vjjvv+Queevw5lJVpMOeXdqe+W1JSjo3L8/HJSwU4dcA49Onb0y/VUz18QqDMNKZOfRsgvf1+/Gf6g9i7f69PaqZqBBNgMG4uc7v++usthluHkG8KIINSMH19rhaBhg3TcfbECzCu3zTMvTsR65bmodyUJcVUqhbWk3Ix19w+8MAD+PTTT/Htb38bvXv3dkveTkphdFMR+BoCZfbD2a53EnqOSsfM+/Zj79bDwbm5JFcp8gnwN2Tt2kPYtMi80dr0wODBg5GSkhL5FVMN6hQBeinF2EYCLToloY+FInzprQdx2x/uDPnuO3UK8kmqLOP20Q9p6+o8fPJCAa669AZMnDAJTZo2OUkl0m1FoGoCXJKZnBiLFl1jkZF9wEJavGkx3Q5WfbKORiyB9957D//+978tLE0rJCQkhLweUrdDjrhu3aBZi6aYePp5GNLxcix7MRkrP863uBcK1B1JT0F+fr4TQps3b8bAgQMxduzYkFu2I4mPyhp+BKggJZkbd8vsRDRuHYePXt6HbeuKEGMWJQXoDr/2+kYlsoEal7u9+/x2rP8oHudeOAnNmzeXt+Q3gqiTw4UAvSlj46Mw8dqWSM3egBdnPoiHHvg3cnNzw6WIKscJEIizZW1rlu3HstfroX3KBFz1vWvQrVuXE8hJl4hA6AlUmJdSi04pyJmcgnv/eR+2bt4e+pvqDrVGoKSkBCtWrMD69euRlFQ7oQFkUKq15q07N2rSrBFu/NmP0eDAeKx6IxXrFheaG3CZBequOwwitaYHDhzAnDlz8Oabb2LEiBE499xzI7UqKnddImBGh2JUIC0jBuf/rAkWv7UfH764F3u2F2uDgAh+DmgMtP+xfVueBeKOReOYvphy9kTExcVFcK1U9DpNgBYlS+06peKsa5ohvfsm3PvwX/DS9Je189thNBH1L2UUl7pt3XQQS16tQMyWgbj8whvdEhMnvCKqNipsXSFQZpNw6enx6HJqIhZ8/h7mffIBCguK6kr1fV9Pri6hnjRmzJha8U4iUA3xff9YnZwKcqeLW/74G7SNPhvLLVD3yvlypzw5LXH8d+UM6TvvvIM777wTP/3pT3HGGWcc/8U6UwROMgEuJ4lLikLnAWlo2ysJ7z+3BzMf2IGiXIvQrRR5BGygxr0dCq393nxoBwb3mIybb/0lOnbUzm6R15gqcWUCheXlyO5dH5N+kIVWQ/fghz/+H8ya9SboIawUOQQoo0ps98nHbl6H9ENjceP1v7fdJztETgVU0jpLoNy8lBJTYnD5H5vikWdux2OPPmm/uV9YvOsslciveGlpKTZu3OiM2hMn2trqWkoyKNUS6Lp4m4SEeAvUfS36NrsIHz+WZAFVDxgGBeoOx2ehsLAQzz//PJ555hn89re/NVftblpSEo4NpTIdncAXelBiWjTO+XEL9BheH4tmHcATN292SlKUfu2Ozi4Mv4mxSYn8gyV46R8b0abehThv/DXoP6B3GJZURRKBEyNQYoO35h2TMOn6phh6RQDX/s+VuP+eh1Fuxial8CZAz6Ro80w6uKcY91y/Aqf3+KnFTboRrdtqs4DwbjmVziPgJuFs+e3gSZlIaL0RDz91B+782z+9r/UegQS4yuQ3v/mNMyiFele3ynikYlcmor9rlEBGo4Y4d/JFGNtnGmbfFe8CdXNXEwXqrlHM1cqMAbgffPBBLF26FFOnTkXfvn0VgLtaRHXxSSNgRiUurW3UPB5jL2uC7mZUWv7eQTz1uy0ozCtHnAV31jKEk9Y6x3dja8Ok6Bjs2HwIb/xrN1pWnIPJY69Er559ju96nSUCEUIgYF6V9OZu0jYRY76Tie4TCvHkS/fgr3+6C0VFWn4Srs1IYxIN3ltWHcKse/djUOurMOG0b6NTJ3kmhWubqVxfJeCckUwlSm0Yi9OnZaJhr414wWK63XXHfSgtObyxyVev0pFwJUDPpN27d6N+/fo45ZRT0K5du1otqj1KSiIQWgItWmVh0hkXIKfVJVj+UgpWf2KBuosYqNt+lZVOKoHCwgI8/vjjWLduHXJycnD66acrAPdJbRHdvNoEzCDBmbdWNvM/8uJM56nEIN1vPbIbu7YVI9ZclTggUApPAgkxFtx27R4serUMaQdGYfLoaejbpx9i4qLkjh+eTaZSVYNAucUyYWppO7+d98O2SGq/Fs+99i888tDjyMvLq0bOujRUBPgbsu6zA1gyI4DMonGYeva1yO6cHarbKV8RCB0BEz/lpi916d0Awy+qj5TOa/HUS/fi6aefQe4hyZ/Qga/ZnBl8m7FvuWSaGyn16dOn1h0DZFCq2TZVbkchkNU8Ez/9+Y9Rf99YrJyVivVLLFB3abkCdR+FV20cPnTokIuZNHPmTAwfPhwXXHBBbdxW9xCBkBPgzFtxeQXadU7B2Tc0R6f+KebtshMfPL8X+/eUKFB3yFvgBG5gRj4a+vYdLMQn04tQvKIvLp30Cwwc3N+C3h7Or54sgScAVpeEOwHKKwbJbZ+dgonXNUea7f52131/xquvzMDBg4fCvfh1p3wmnxiAe9vmXCw2g3e9jTm48ls/R5eunaXL1p2nwH81NflTUlGBrqdkYMwVGYhrswp33PN7zHpzFgoKCv1XX5/VaO/evZg9ezbee+89pKamolevXkhMTKz1WsqgVOvI6+4NY2Nj8bvf/xqtApOwbHoyVi88aKtP5CpwMp4IznwyAPdf/3o7fvKTn+Css846GcXQPUUgpASKLBZJ/cwYXHlbW3QemIK3H9uF1+/bgeJ8xSgJKfgTyZxKbWEFnv3zemTmnY5rLvkN+ub0OJGcdI0IRCSBQjOCd+mbboG6m6LZwJ244Uc/wJzZc8AYh0phQICGv9IKPHnLWiTvHokfTrsFXXtok4AwaBkVoZoEnFHbjErteqRh6i9bIbbVGtzyp5vwwvPTq5mzLg81AYYsWbVqFSZMmID27duH+nZHzV8GpaOi0RehIJCUlIjvX3s9emVOxYcPJWD98gNuVpqzPkq1Q6CkpATTp0/HE088gV//+tfo3r27xXGQKKgd+rpLbROgohSfHIULb2qJHiPSsHDWfjx962aUlwUOzypL9NR2k3z5fsafwW3zDpTgxTs2oFPCRTh3/DT06df9y+fpLxGoAwRKTWC1zE7GlB9lYcBFpfjetEtx/z8ftmW8FXWg9uFZRTpGUkbl7i3BnVcvw5jsH+OaK25Ch05tw7PAKpUInAAB6krc5S2jWTyu/Es7pGZvwh3//D3uuP2eE8hNl4SaANvquuuuQ+fOnfHzn//cxb8N9T2Plb9Gkceio+9CQoCBus8/+9s4zQJ1v3VnrAXqzoUCdYcE9VcypTGJ1uyFCxe6ANz9+/dHQkLCV87TARHwC4HDShKQ2SIeYy5tgq6n1sen7x/CM3/YgqK8ChdTyS91jcR6xNt6tl22hOTdJ/chq3gSzh77XfTu3TcSq6Iyi0C1CTD+GwN1Z7VPwGmXZSD79Dw89vxduP3Pd4O/30q1S8AZk6w9tq7Nxaz792Ngq+9ZAO5LbBDXsXYLoruJQC0QoL7EJeZNmidiyg+bokHPjfjPS/fjnn88AAZ9VgoPAjt27MBtt90G7uTWtWtXF4g7Li7upBZOBqWTir/u3rxl62aYcta30KfpRVj+cjLWLiiwQN0lhwN1y2MgJA8G3eafeuoprF69Gv369XPL3LgMUUkE6gIBBp5sazP/I76Vic4D0jDvpX2Y/fgu7Ntdgnju/qZU6wQSoqOxcf0+LJlZividQw4H4O7bH7EKwF3rbaEbhg8BF6jbBnatuyTj/B+3RXzbNXj21Qfw6CNPKlB3LTdTjHlvb1x5EEtnVCA9dzQuPvc6dO3StZZLoduJQO0RCJgzZKktv+3WsyFGX9rAAnWvc0bt5559AYcO5tZeQXSnKgns2rULH330kXPEmDx5Mlq0aFHlebV9UFp0bRPX/Y4QYKDun/3iRqTuOQ2fz0rG+qUFKC4pc98r9uoRTNX+wFnNrVu34t1338Vrr72GYcOGOe+kamesDEQggghw5q2UMQK6JmPCtVno0DcZr92/E/Om77XAt6Uu2GoEVSeyi2qTBlzmzADci14rQuGKXjh39I8xZOhARMccnlFQAO7IbmKVvnoEKK9coO6uqTj7By3coO6Oe/6A11+fCW6ooRRiAl/IqB1bLQD3jFIE1vfHlRf+HF27dVEA7hCjV/bhQSC/vAzdcxrhjO9nIr7dGtx57x/w5ltvIi9Xu7+drBbavXs35syZgw0bNmDatGnmzd07bFaZyKB0sp4K3dcRiI+Px823/ArNSydg+fQkrPz4IAoOllu8AAGqCQK5ublO+Nxyyy24/fbbccMNNygAd02AVR4RSYCDtBKbeWvcKh7T7miH9r2T8PJdO/DK3TtQViKhU1uNSpNRwHa1evLWtUjYOgbTpt6MoSMG1NbtdR8RiBgC3K2SgbonXJ+FJv132CYaN9hv+lwF6g51C9pvRUVZBZ7+4wbEbRuG6674LXr07hLquyp/EQgrAqUV5WhtRu1v/aIlYluvxG1//y1emv5KWJWxrhSGSw6ffPJJLF261MVLSk9PD6uqy6AUVs1R9wpDT6SUlBRc9b1puHDYLdjz1qm4fdpSzHtlG0qLS5FoS1EYDFHpmxN47rnncN555zlvpGXLluGmm25Cz549jecXe3B/8yx1hQj4ggDjlMQkROGSm1uj18g0zH9tH56yQN1lpYcDdctDMnTNHGMyPd88wh67ZSV6pl+KiyZdh959tJtb6Igr50gnYP4xaN3NAnX/uCl6Ti7A96+/Cv9+6AkLoCsjeE23LeNXxToZVYK/XbMEQ1t/H9dc9gt06ard3GqatfILfwKchKM+1MhiUH7r/1ohsf1a/OP+P+KuO+4N/8L7rIT0TurSpQu4zG3QoEFhV7uYsCuRClQnCTRr3gzJSfWR1bgVTuk5Ais2fIDpt3yA+u12ovvIFLRr1wAMB2eT2gjYYFCpagJc2jbH3CGfffZZF2shLS0NV199NUaMGIGcnByc7KBtVZdaR0WgdglQSbIxGho3S8C4K5ogJj4Ky989hPjEbTjz6qZIqR/jlpvUbqn8f7c4G6ht33IIC2fko1nJZEw6/XL0sQDc9TS15f/GVw1PmECFKT6cWGveIRFjLs/Agf0b8dCTf0dBfjGuv+FqWyaqDnTCcHmhDZi5xDbBBFFhRRlWLjiA5W8W4ZSsKzHxtEst6K08k6rFVxdHNAHarbkUvVmHJJx1fSZeu3sdnnj+PsTFx+G7V10W0XWLhMLv3bsXDzzwgNuRu1u3bmjevDliYsLPfBN+JYqE1lUZQ0KgfoNk9G7Qy2are+Hd2T0x+9022LtxCTbO2Y79G/Yjo105MrISEB8fy7GgLYuzf/lBye3+MnfuXLzzzjtYsWIF1q1b54xIZ5xxBk499VQkJyeLkgiIQCUCJYyp1C0FFVO5vCGAj17eh+S0GAw5JwONGseBy02Uqk+As/4xNmDbuukgVswuBzbm4Nxx30ffPv00GK4+XuVQBwgwUDcNr+16pOCcH7a25aIr8MLMf6FB40T07NbPedAkJmnH1m/yKNDzop7JJprjSm1522fLD+KzDw9i1+cJaJkwApdc/AOLmaQA3N+Eqc71JwE3CWdjru69M1B6RQVm/WsN/vXo39GgYX2MHz8OKakaY4Si5fPy8vDBBx9g8+bNGDhwIJo2bYpw3UxJBqVQPAHKs9oEho7MAV/LlqzAzDdew+J5s9Bq2E4UdClGw+YBJKdHu52AOLPkBF217xiZGXDntj179mD58uV44oknsGbNGgwYMAA//elP0aNHj7AJ1haZdFVq3xMwg7QL1N0jGQlXZ7n4ba/dvwOxCfUw9LxGSEyJRoUJGM+jyfc8ariCHADXs/9KLT7VgUPFWPJGEcrX9MXZI3+IU4eFn8t2DVdf2YlAjRKgp0CZzaJ16pmKc25ojXkvbsJDL/4O3RcPw5jhE9Gze2+bvW6GhMT4Gr2vHzNzXpEm/4vyy3FwZzG2rSnCu08fwuoPgD6dB+DKX92E7j1kTPJj26tOJ0aAelCBBeruOyATSamxeOrm1fjT7b+2lQ+xGDNmNJKSk04sY11VJYGioiK8//77WLVqFa644goXgDucNyuRQanKZtTBcCHQo1cXZHftYNbZyZg75x3Mfuhx1Gu+El1HxKKNeRYkJNfNeEAV5llBQ9KCBQvw2GOPuSBtU6ZMwS9+8Qu0adMGiYmJ4dKEKocIhDUBKkllNvPWpG08Lr21Ne790ToL0r3djEtlmPyD5jImVaP1PGPShqX5mPnoBnRImIyLJvzIdnNTAO5qYNWldZwAPSc75aShdfcUbFtdhAUzZ+C6nzyN7u2H4aqrrsaUKZPNg7sCUbblvdJRCJjczz9QhlXzczH7sd344JkC5PQ9FT+65jzQs7t9+/ZHuVCHRaBuEyiwZaGtuqTg0t+3wSM3rcWf//ZbC7FRgAu/dZ5kTg09Gtyd+7777gM3VjrrrLPQq1evGso5dNnIoBQ6tsq5hgjQva916zY455xMnDZ+FJ546mE8fdN9iG+xDqMuaIb+4zKREhOHQtjucAyy5PPEOEkvvfQSXn31Vbd9MANv33jjjcaotQtwHo5ra33eJKpepBOg2DBvx+QGMbjiD23xzn9sgPHiAbz3/EGMs5glfcemI6NZvC2NqGfGp4o67RX5dU0dZbFeEozU/twSLJmzD3Oe2oEFr5fhyksux2UXXo0+vfp/XRb6XgRE4DgIxCdGoUV2gnltZyLnjHSsnb8at/3rOvzyZ7/HD350Dc6ccDqaZjU+jpzqxin0SoqxfxIRjVef3owX7tiO3E1JOHP8JFx535no3/8UNG/RDElJmpCrG0+EankiBDgJx+WiGc3jLFB3Szz/l1W475HbsWvbPnx32neQnCJPpRPh6l1z8OABPPjgQy6UyfDhw91mSt534fwug1I4t47KdoRAtAVzTauf6l7nTrgUWZnt8fmGj7H+9YWY98JStM+JwqCzmqJls1TYogoL3m3LVHwW/mThwoWYOXOmW0ublJSEkSNHom3btujXr58zJh2BpQ8iIALfmICnJDXOisep5zayXU0SsH5JIdYvy8faRfnmDZCE7AGpaNElCQmxUW7pSV0wYB8vSAYNNh0Te3cUYeHbtgx37kHsX5eEzMRB+MnVo21nkkm2JKc7aHBSEgERqB4BT7+Ji4+2pe3RaJSZgLTGBYhvaPLq47WYs/IefLj4dfTtOtziKA5Dtx6d6+QOrxz4RlucJP6Xn1+CFUsPYvYz21GyqTWGdJyKvlOGYNCpOejWPdtNyFWvVXS1CNQNApQ//M1vmZ2M0y5viNUfbsOcz+7FJ1e/h3GjJ5gxexwyGjWoGzBqqJafffYZPvnkE7dBQOPGjd2Obp07d3Z/19AtQpqNDEohxavMQ0GgXcdWaNfxIqxfOxQzX3sTc96bhe3vrcJ7B/ejSdcD6NA7HemZ8Yi2QZ8XtJuDxUhMpaWl+PTTT92SNsZHYrDtjIwMDB48GKNGjUJ6enokVktlFoGwJEA5weUkjZsnoOk5CegxrASL3jqAlfNyzbiUj/07StHsswLbwjvJdjxJsFgl0U7EMGBuXU1UKpkKC8qwdXUe5j2fi43zYxFX1As5HU7BaaedZsrleNtMQXFd6uozonqHjgA3Jzk8dxZARpNEjJzUEv3GFOOjGevw6TvLcWDZSqzftgxdPhmIgYP7o137tnUitiK9kaLMmsQYePv3FJlsKsCWT8uxb00qAmuGYGS/sZgw8SzbBKZH6BpHOYuAjwm4cZXJn26DGqJp20KbSNqChTNXYtairVi9cRl6dx2Avv16o237Vj6mUL2qcWkbx3bFxcVYu3Ytdu3a5RwFGIC7RYsW1cu8lq+WQamWget2NUeAQurq71+O86dOwYxX3sCzzz+JBa9+iMEXFKH9KeXIah9vwbttRzjTtjhL5abPefsvxn7hamQKWMEoXBgj6fPPP8ezzz6LlStXOrfHadOmucBs2rWt5p4j5SQClQmUfrG7W/3MWIy6sLEtq03H0rmHsNiMS+sW52HT8iR0G5pmy02SkN7EjCcJUSZnvgjeXTkzn/7NARtn/Qvzy3BgV4nFcinB6vcDWP5sEwweOAznnD8Zw4adqnhuPm1/VSv8CJSZYbvMPLQTkk1undcKOeNLsfz9nXjvqcfw8vTXMGz+KJx5xkQXjyMrK8uMvHHhV4lqloi6HgPX0ns015bd7t1ajI3LyrDGgm3nb2iOnu1H4v/++B1kd+lczTvpchEQAY6jOKHWsGkizpjaFqdNLcf701fhkVvn4o05XXD2hAsxYuhYtG3dHvXT6/tS5pzIU8AYd4cOHcK2bdswe/ZsF3uKcdumTp1qGys0P5EsT/o1Miid9CZQAapLoGHDBrj4kgtw1sRxeOedd/CPO+/C8reWYPDUWPQc2QBlJQGkWGyU6JjDM+m8XzgZk2hAomChpZqKECP7r1ixAo8++ijmzZtnATan4K677kK7du3cLH84R/mvblvqehEIJwJUlMrNAp2YFoOBkxqi16j62PRpAT6ZsR/3/2QT2vZIwJjvZKLPGHPtrktOShSlVl/u3rbk7YN49z+7ULApC6MGnI/rHpzijN/aGCCcnmSVpS4RoDHFIr0hPinajEpNbQOTYnz44j68eu/DeO6FpzFk4Ei3a9DEiRMjZjnFN2o/k037tpfgk9f24aPp+9Ggoi/GjzoHk64/Cx07dvpGWelkERCBrydAmcNwI0wDz2yKIZOaWsD7PLz12F9w33f/ijEDL8D5512AnAGnIDU19esz9PkZmzdvxpw5c5zjwLBhw9C1a1dEuqOADEo+f2jrUvW4/GvM6NPQq2cfLFsxH+/Om4m//H0mVq3bgDGXxSGtUYy5QEchtUE82vRIRYdeKRaasZ69zNhkLw4cy75wHic3jg9DsZSFlunp06c710YOumhMevnll50RKT8/3wnb7OxsnH766bj++uud+2P9+vWhYNt16WlWXcOKAGe97b/ElGi0652Cxm0SMHhKBl65ZzvuuGIjMrJ2YcpPGiFnQgMkJxyWJeU+81jizH+UxSJJMmm5J78Ibz66FTPu3Ye+nYbivJHXIef6U7+YhUyzbYT95/kQVs+jCiMCx0mAekzDrFiMvKgReo9JxWfv5WHZmx/g+h/Ow+9+fZsF774Wp58xFg0bNjzOHMPvNHpLxpiASjDZtGzZPrz2r43Y+EkcBvc6EzdcfDr69MxB82YtkZqWHH6FV4lEwGcE6tnm29SXWvdIwtk/bIbtk4uxedHb+Osj7yHrP/0wauhZboOlzMxGPqv5savDsd+iRYuwePFiF6+Ny9poXKPsTUhIOPbFEfCtDEoR0Egq4vETSEpOROvklsjMzERWRjt0az0Se/fvQpueCUhOjTX3wu1Yu2EVtszcieIViTZACmDRwmXYumsNOvRLQfveyYhPjEO9gAUBz4hH8zZf7PYRHbB4KTFIjjm2QpKXl4cdO3agoKDABcCk99GCBQvAgNrbt293QqRnz57OzbGwsNDNDtJQxEj+o0ePNs+pgBuMce3sKaecgjZt2hx/5XWmCIhAaAjYqIyxOExHQpztrJRpMiIzKw7jrmiCRs0TbUlFIeY+vRvvPnMQfU5LNU+mVGS1S0CsbdtdakoEL43I9IURiYZ320MTW9fm4r2XtmPL4ji0bzQI351yCvr3HmhxSHqhecumEVlFFVoE/EyAngM0BNdvGGsDlzg0aBKH5p1zsfzdfVi/aCEeffE2fDTvY0w4ayKGjRiM+ITwMwYHB/I3kWRmI04EcshqctleeUUlWLXsEN54bCsSi9qjfer5GH9uDvr1yUH3nl3RoGF9Pzex6iYCYUWAYUbKbEqeGwY0bZVsu8El2CYn+WjY7gCWz3wb//z3Z3jr7Vk27hmJ0WNGuJ0Vw6oCNVwYjvX27t3rnAjoNNCxY0dwyTF35vbT5JsMSjX84Ci78CCQlJyAfjm93Cu4RAf25uPTRp9i08ZNqFdWzxlw1u95H7s2LUFSszSktUzFjp3bsXXLDiQ3LEfxYFOuTGvJNSGwdZ7NeuW9iejEEnTo0AFpaWkuSDZdF2kU4lK0rVu3YufOnUcMSvRAojDhOtkNGzagQYMGoPdRjx49EBsbCwbdTklJwYgRI9x3wWXVZxEQgTAjQMOSW05iYsHkQo+B9c1jKRlrFuTh41f2Y9mcAnz8cil2rCtA254p6Ng/Gc06cpNq83akVSnIsBR2RiarDwdoXuIglMUtzCvFni2F2PhpEfavTcLehV3QrF5PTBhyNk4/8zRTGmO9S/QuAiIQhgQoaw7HVwogzYxKfYc1QnZOA6xbcgjzZnyO5St3oPC1zdi4bRU6tulu55aanKtAs+bN7JVlSzESnK4UFR1testxDBuCBcnX8fhC7lDW8DLqUSZ67PPhTCwynelQNjw1j0/KpOLCcuzcVODitjEgOQevBXujsWtFMspXDkLPvmMwYcKZ6D+gjzv/626v70VABEJDwFt6y37bqkOavVItru0eW377ORavWIONzy7G8hVZm3r7AABAAElEQVSLMGjAUNutui9atY6sINTHokbngNzcXDcmpKOBN0Zs2bIlWrVq5csdN4/jl+FYyPSdCEQWgfSMZAwZloMhyDlS8CnnnWXLznJtpBhlShQwd+5cFB+cj7xde7G76LChaNGy+Vi1dA3K8rdg6WfzccEFF7gI/DNmzHBxmzx3RS5NowsjhUeUeSc0atQIF198MShEysrKnBChUUlJBEQgsglwkMYd4eJsy+5uQ+qjo3k4bp1aiHef3oPFb+daAO9c9DZPpdMua4xGreJNJrgo1q7SNEbxdcSCY3nVtoHpS/e3onBgZuOzw+Wzf4uKynBwTym2rSrBuo8qsO79ZLRvMgTXX/ItnDZupC8VosO1178i4F8CXly42PgoZA9IR/s+9bF1TR5mPzUD9z/7DlIq2iNQr9wmu0osqP5wt1wsMS7d5FOFm/xq07qt2122vMyUJZNhsTGxZnCyhbAxNJsfTp4xiMYhZxbiP19YjNzf3olB7+4UO6ekuAzFBRUoL6MBiQZ8YP3iQpQWB2znXjMeHQQ2LCrH1s8sREFpBaKjopESZ1tsNx+BX/79KnTq3DEoV30UARE42QSoW5iJ2hWj//BMdB+SgVULDmHuUyvw8EvzLXj3dEwcfx4mTTjX7ULZoEHk7l7NZW30Qjp48KB75y7dTL1793aOCO4Pn/4jg5JPG1bVOn4CNAYd3tL6sPJz3oWTcd4Fk03/oUIT5WbMbrvtNqRiMTq074BBgwa55WsffvihC6Q2ZswYZyyiRZqBs7lcjcHVKFicQkWlyl78nu9KIiAC/iHAWTim6NgotOqajAv/NxFjLrWAsDMO4PUHdmKLbVc9cmomOvSloflw/49LjkJiqkVus2soF052olzibm0FB230ZsXhQG714lwsmpmL3DWNLU7SePz+f8+z3aF6R3zgyJPNWvcXgXAgQO8eJrPHoGWnFFzyy87YvrYI65dvcxuYUFWJjnoJD73+OGY/sdcF+M5q0hg5vYeiSaNWKCwotmuj0LpVa4wbP850n7bO45p6T1l52bGr6KxH5rVpRqOSogpnNOL9WKYtnxeYQf6AeUUWIzE5BnHR8Zj/vG1WkmdeVhWltgNSM5x99jmYesVhL+9o85pq1sx2b7NQAtKvjo1d34rAySbASTj29c7909CpbxpOv8K8ul/dh4efuBV3/f0+3PrHW3DhhReAk/ORmLgahWPD3bt3Wz0uxLnnnhsWOl5tsJRBqTYo6x5hTyBYEWHg7q+memjUsJHbbnfUqFGgCyN3ZaMhii9vUMj1sPROYqKiE5yC7xF8XJ9FQAT8QYCKEmf+G7eOx8iLM22XyfpY/ckhzHl8Lx7+2TYU59VDjK2ibdcv1gxMyUi1uCYcRKVlxGLoeY0sRlu0M2RzqMeXmXdcYryQmFhv3v+Lg0d5YxnoYm6m8C+dQal2+FUPeeUlWLP4IJbNPeAGj7EJUVj+Th6Wzy1xN+VYc2BODk4bcQUGXTgCHdt3QXoDbfn7JaD6QwR8QoAqTz2TGU3axiO9Kb2yD1eMek2rHvEYc0W6fQ8LHBuL3euWYfvqRcjbX+bOs1X+ePiSe7Bvo0mcmArTf6IxdvwYbFq/B/EVBZj96lZktU9Eh+w0rFmRiw3L81BqRiTa0betzcf8GQXYuf6wYctUKVxz9TR0T2yA/KZF5gmVZjP7ffDLi9uZccs2OzBjVVwcY0E1QFJS4mEDkhWW4QOkX/nkYVQ16gQByhPGRmveKRHjsprYhgEmH+YX4uX5f7Adrh/B+NFTMGBwPzNcFzoeI0ePtKD6X8S0DSNC9EaaOXOmGw/yMzeH6tKlC0aOHOk+s6h1RTbJoBRGD6aKEr4EKszvmi7dSUlJTnnRsrXwbSuVTAROFgEOkviKMk+k9IYxNvCJNWNRNFIz4rBrQ4kt2zBDs32X0cJifmwswo71RU7ZOLCrFE/fugW7N5Qjf38FYpPsnFZR6DfWjDgWBNzGUXbe4Wuj7Vc7MTXK7bZWbpsF8L9YMxV5xiJu3btza77FhSt0u1Tyfgf3FOHzjw9g+yqLC1dUjuzO2Wid1Q/x28pdvuWBMoywQJHn9+1k5TdvSrth+w62qUH3LmjTtpVT/E4WU91XBEQgtASc3DJPy2gzWqeYwSY4pdW35bqtD8dT47/1mxQio32Fk2U8L/9QKdIt2O7+HaXO6B2oqIeNa1/D3j17kFGvFfZ9mI3ilcDS6ZsQVdQU5flJbjktJ+6apzRBjwkdER+VbrKq1OlWo2zgmFY/xf4ud5N1TbOaIinZLE1KIiACviHAZXDUXShzGtgGSOn2ato2ESvb2mYBC9dgxZ7Hseb11+2Mw9NqyzbPQMm++ti/s9yW3iahf/9+6NylAxKTEtwKEU7mt27TwhnGaxLS/v373UZLNK5z5QmXsm3atAnFxcVo0qQJGjdu7F6e0SgjI8MF22b83LqWZFCqay2u+p4wAQoUunMriYAIiMCxCByOGWCWJVOHGjWLR9YkBub+sofRcvMQ2rGh2ClA5RYLZNUnedi3vcQMSgHE2A6yRWZ84gCOChdlD68uLghg79YyvP/iLuy3DQbSM2ORlBaLDZ/moTCXsslkVEkMivcmIX93mi0lqed2Wom1Bbvp+zKQa4O+qHxThFrY7mxZg5HUwYLt2n8lxSXo268XWplCpiQCIlA3CVBulZs8CE78uzToQEp6Auqnm3dQ0LG+Qw7LKnpFcnenF+5ej3WfRaFby+EY2+tc5BXkYemyhWjVog0yOzSxpXEWT9LW1bZp29o2TulthqPw21kuqHr6KAIiECICXnwlTphxw4CBY7PQa3gplr23G5s+X+c8vnnrdXnLzJBTgZ0bSk3+pKCsySAsXt8S9cpjbcOTMiTEJaFT255okNrYJsmizZsyHo2bZKJFq2bOm/FoxS8qKnKbKNEgxN25Dx06ZAbvIrs+we0Wzs2UuEMbDUSMg0uDEs/jChVuuETvyD59+riVKYdDpxztTv4/LoOS/9tYNRQBERABEThJBA7vruQtXvtvIdr1TkGH3qlfHAhgyNkZFlPksMG6rDSA3ZtL8I4F+C7Kt2UeifXM2FSGxbPycHBbNGb9ox42fxKNHsPjUb9pFGY/Wg+7N9EwVI70+mno0bkv2rfJtt2SbFemilj06NkdkyafZUG3bXhoxilvWe5/S6NPIiACIvD1BLydm75ypg0IOSjk99kDG2Dju7EYkjMI3PSEqbBwkhuUfeU6HRABEajzBEwtcR7VNGDHxEWj/+gmyBnd9AiXEvNYzNtfars8VtiSfjNyly7Dy3e/7mKtxSSY4CmLxRvvNrOZuFQL1h8APYV6dOuDnL5D3c7a3CCJK0wYF5cGIybqQTQWcaduHt+3bx+2bNmCAwcOuOVqAwYMcAampk2bOuMSvaC4Q9vQoUOPlEsf/ktABqX/stAnERABERABEagVAt5uS8E3q8do2JZibYUH446c/7PDHkMxcabsbC9Cux4HMP/fDXHDdb/A4CE5WL9hPRYvXoR//Ha4Wx5SZjsvcZYsMzPTlts1dEqS826ykR5jFkTbf0oiIAIiUOMEnGPTYe+mchvQUb6Vlv7Xt6kuLgGpccbKUATqAAF6LdHTkR7X/022HLdBHFK+2ACO31z2u87uDMaMZFy2bWuKkLcvz3lPFuQewmcfLMVfp/4dZ50+GT169EB2l2y3PI1Bs+lhNHjwYAvo38x9pkGpbdu2trPlMKSlmXe3rUbxlrF5Zaj8t3dc74cJyKCkJ0EEREAEREAEwowAbUv1Yg4bmLhcLiE5Go2axyM11Tyb2ne0XZXao2kT292oex9zAU9zxiMGcKLSE7yFd5hVS8URAREQAREQAREQgW9G4AsvSF5kH51+xHem6JQotMhOQIXtHMlkMfyR3rwePp6xF8NGDsG4MWeiqQX/pn7UubMZokxXSklJsSD/FnupdWt3nMvXvGVrNDApfTMCMih9M146WwREQAREQARCToAu4AGb5WeylWtOAWKwb27VzZ2OmJKSE9zL/aF/REAEREAEREAERMCPBKgTBdXL04+8Q3HxpihZ/EmmBNumJCMrAakNYtC0aRNbqtbaAngfDu5PDySlmicgE1zNM1WOIiACIiACIlCjBJyBiQqVfdDmADWKVpmJgAiIgAiIgAhEMIGKCtONbBKOLy6ZYywlfi61YNpc4qYUWgIyKIWWr3IXAREQAREQAREQAREQAREQAREQAREQAd8RkEHJd02qComACIiACIiACIiACIiACIiACIiACIhAaAnIoBRavspdBERABERABERABERABERABERABERABHxHQAYl3zWpKiQCIiACIiACIiACIiACIiACIiACIiACoSUgg1Jo+Sp3ERABERABERABERABERABERABERABEfAdARmUfNekqpAIiIAIiIAIiIAIiIAIiIAIiIAIiIAIhJaADEqh5avcRUAEREAEREAEREAEREAEREAEREAERMB3BGRQ8l2TqkIiIAIiIAIiIAIiIAIiIAIiIAIiIAIiEFoCMiiFlq9yFwEREAEREAEREAEREAEREAEREAEREAHfEZBByXdNqgqJgAiIgAiIgAiIgAiIgAiIgAiIgAiIQGgJyKAUWr7KXQREQAREQAREQAREQAREQAREQAREQAR8R0AGJd81qSokAiIgAiIgAiIgAiIgAiIgAiIgAiIgAqElIINSaPkqdxEQAREQAREQAREQAREQAREQAREQARHwHQEZlHzXpKqQCIiACIiACIiACIiACIiACIiACIiACISWgAxKoeWr3EVABERABERABERABERABERABERABETAdwRkUPJdk6pCIiACIiACIiACIiACIiACIiACIiACIhBaAjIohZavchcBERABERABERABERABERABERABERAB3xGQQcl3TaoKiYAIiIAIiIAIiIAIiIAIiIAIiIAIiEBoCcigFFq+yl0EREAEREAEREAEREAEREAEREAEREAEfEdABiXfNakqJAIiIAIiIAIiIAIiIAIiIAIiIAIiIAKhJSCDUmj5KncREAEREAEREAEREAEREAEREAEREAER8B0BGZR816SqkAiIgAiIgAiIgAiIgAiIgAiIgAiIgAiEloAMSqHlq9xFQAREQAREQAREQAREQAREQAREQAREwHcEZFDyXZOqQiIgAiIgAiIgAiIgAiIgAiIgAiIgAiIQWgIyKIWWr3IXAREQAREQAREQAREQAREQAREQAREQAd8RkEHJd02qComACIiACIiACIiACIiACIiACIiACIhAaAnIoBRavspdBERABERABERABERABERABERABERABHxHQAYl3zWpKiQCIiACIiACIiACIiACIiACIiACIiACoSUgg1Jo+Sp3ERABERABERABERABERABERABERABEfAdARmUfNekqpAIiIAIiIAIiIAIiIAIiIAIiIAIiIAIhJaADEqh5avcRUAEREAEREAEREAEREAEREAEREAERMB3BGRQ8l2TqkIiIAIiIAIiIAIiIAIiIAIiIAIiIAIiEFoCMiiFlq9yFwEREAEREAEREAEREAEREAEREAEREAHfEZBByXdNqgqJgAiIgAiIgAiIgAiIgAiIgAiIgAiIQGgJyKAUWr7KXQREQAREQAREQAREQAREQAREQAREQAR8R0AGJd81qSokAiIgAiIgAiIgAiIgAiIgAiIgAiIgAqElIINSaPkqdxEQAREQAREQAREQAREQAREQAREQARHwHQEZlHzXpKqQCIiACIiACIiACIiACIiACIiACIiACISWgAxKoeWr3EVABERABERABERABERABERABERABETAdwRkUPJdk6pCIiACIiACIiACIiACIiACIiACIiACIhBaAjIohZavchcBERABERABERABERABERABERABERAB3xGQQcl3TaoKiYAIiIAIiIAIiIAIiIAIiIAIiIAIiEBoCcigFFq+yl0EREAEREAEREAEREAEREAEREAEREAEfEdABiXfNakqJAIiIAIiIAIiIAIiIAIiIAIiIAIiIAKhJSCDUmj5KncREAEREAEREAEREAEREAEREAEREAER8B0BGZR816SqkAiIgAiIgAiIgAiIgAiIgAiIgAiIgAiEloAMSqHlq9xFQAREQAREQAREQAREQAREQAREQAREwHcEZFDyXZOqQiIgAiIgAiIgAiIgAiIgAiIgAiIgAiIQWgIyKIWWr3IXAREQAREQAREQAREQAREQAREQAREQAd8RkEHJd02qComACIiACIiACIiACIiACIiACIiACIhAaAnIoBRavspdBERABERABERABERABERABERABERABHxHQAYl3zWpKiQCIiACIiACIiACIiACIiACIiACIiACoSUgg1Jo+Sp3ERABERABERABERABERABERABERABEfAdARmUfNekqpAIiIAIiIAIiIAIiIAIiIAIiIAIiIAIhJaADEqh5avcRUAEREAEREAEREAEREAEREAEREAERMB3BGRQ8l2TqkIiIAIiIAIiIAIiIAIiIAIiIAIiIAIiEFoCMiiFlq9yFwEREAEREAEREAEREAEREAEREAEREAHfEZBByXdNqgqJgAiIgAiIgAiIgAiIgAiIgAiIgAiIQGgJyKAUWr7KXQREQAREQAREQAREQAREQAREQAREQAR8R0AGJd81qSokAiIgAiIgAiIgAiIgAiIgAiIgAiIgAqElIINSaPkqdxEQAREQAREQAREQAREQAREQAREQARHwHQEZlHzXpKqQCIiACIiACIiACIiACIiACIiACIiACISWgAxKoeWr3EVABERABERABERABERABERABERABETAdwRkUPJdk6pCIiACIiACIiACIiACIiACIiACIiACIhBaAjIohZavchcBERABERABERABERABERABERABERAB3xGQQcl3TaoKiYAIiIAIiIAIiIAIiIAIiIAIiIAIiEBoCcigFFq+yl0EREAEREAEREAEREAEREAEREAEREAEfEdABiXfNakqJAIiIAIiIAIiIAIiIAIiIAIiIAIiIAKhJSCDUmj5KncREAEREAEREAEREAEREAEREAEREAER8B0BGZR816SqkAiIgAiIgAiIgAiIgAiIgAiIgAiIgAiEloAMSqHlq9xFQAREQAREQAREQAREQAREQAREQAREwHcEZFDyXZOqQiIgAiIgAiIgAiIgAiIgAiIgAiIgAiIQWgIyKIWWr3IXAREQAREQAREQAREQAREQAREQAREQAd8RkEHJd02qComACIiACIiACIiACIiACIiACIiACIhAaAnIoBRavspdBERABERABERABERABERABERABERABHxHQAYl3zWpKiQCIiACIiACIiACIiACIiACIiACIiACoSUgg1Jo+Sp3ERABERABERABERABERABERABERABEfAdARmUfNekqpAIiIAIiIAIiIAIiIAIiIAIiIAIiIAIhJaADEqh5avcRUAEREAEREAEREAEREAEREAEREAERMB3BGRQ8l2TqkIiIAIiIAIiIAIiIAIiIAIiIAIiIAIiEFoCMiiFlq9yFwEREAEREAEREAEREAEREAEREAEREAHfEZBByXdNqgqJgAiIgAiIgAiIgAiIgAiIgAiIgAiIQGgJyKAUWr7KXQREQAREQAREQAREQAREQAREQAREQAR8R0AGJd81qSokAiIgAiIgAiIgAiIgAiIgAiIgAiIgAqElIINSaPkqdxEQAREQAREQAREQAREQAREQAREQARHwHQEZlHzXpKqQCIiACIiACIiACIiACIiACIiACIiACISWgAxKoeWr3EVABERABERABERABERABERABERABETAdwRkUPJdk6pCIiACIiACIiACIiACIiACIiACIiACIhBaAtU2KJWXlyMqKgpJSUmhLalyF4GTSCAuLs4954FA4CSWQrcWAREQAREQAREQAREQAREQAREQgfAgEFPdYtCQlJeXh1WrVqFhw4aoqKiobpa6XgTChkC9evWcIWnHjh0oKioC/1aKMAJqswhrMBW3KgJ8jCV9qiITocfYmGzTL9pVYipC21HFPkKgXlQ91LNJN++ZPvKFPkQWgS/kEgstuRRZTafSHibgfl7dPyJSWwSqbVAaMmQInn76aYwYMcJ5KcmDo7aaTvepLQI0Im3ZsgXjx48HPZWUIpOAFKPIbDeVWgT8RsDzc6W+e8RK6P7wW01VnzpJQM+yL5r9S/LJFzVSJeoKATNt15Wqhk09q21QuvzyyzF48GDs2rULMTHVzi5swKggIhBMoKysDC1atECXLl2CD+tzmBLggK2srBz5+SU4sLsI+w4UobigPExLq2KJwLEJREXXQ97BEhzaV4yKcnkBH5tW+H9bUVGGfbtzsfDtvdi9PQ+x8VEIqFnDv+FUwmMS4KQNFylsWZWHvTvSzViqh/qYwMLwy7LyUuzafAhbi3YjPikGZSWe+TsMC6siicBRCMTG18OuTUXYubkcZaWm+9fTc3wUVDV2uNoWoDZt2qB58+YoKSnRcqAaaxZlFI4EoqOjZTQNx4apokyJCfFo07oNYgua4427tmHeU9udolvFqTokAmFPgAO10uJy5O4NILXUluHGaPYt7BvtGAXs1j0b5066CEvWvo1d+6KtPS2cpfTdYxDTVxFDwJa8leSWonWThmjWomnEFFsFPUxg8JCBWLVhNBatjMH2OWkydOvBiEgC9ewntSi/DD06lCCjYSNw/KYUWgK23FlRhkOLWLmLgAicDAJ79+zFgw8+gtWfr0dxUYmLhXUyyqF7ikBNEAiYxSEuPgbtO7XCRRdNdR6TNZGv8qh9AtzMZNOmzZj91rsIlJslSetxa78RdMeQEahn3gAp9ZPQs1d3ZHfODtl9lHFoCMz/ZCGWLlmOKFTb5yA0BVSuInA8BEwOxSVEYcCAHLRs1VIhS46HWTXOkUGpGvB0qQiIgAiIgAicDAKcC9ImASeDvO4pAiIgAiIgAiIgAiLgETCnMCUREAEREAEREIFIIiBjUiS1lsoqAiIgAiIgAiIgAv4kIIOSP9tVtRIBERABERABERABERABERABERABERCBkBHwrUFJoaFq9pmpimdVx2r2rspNBETgRAgcb9883vNOpAzeNbVxD+9eeheBcCNQ15//SKh/TZexqvyqOhZuz6rKIwKRQKCm+1JN5xcJDFVGEahpAhEbQ2n+/PnYtm0bGNwyKsq23P0itnhaWhratWsH7j6nVLMEcnNzsWzZMiQnJ6N9+/ZISUmp2RsoNxEQgRojcOjQIaxcuRI7duxAWVnZETnJpVLsuzk5OUhNTXX3y8vLw9y5c5Gfn+92MuQ5vCYhIQGnnnoqGjRoUK1y7d69G6tWrUJ6ejo6dOiA+Pj4auWni0XgZBPYu3cv1qxZgz179rhYVtRDGjdu7HSPmJgYLF++3OkiTZtqp6vNmzdbEPJNyMjIQMeOHUO64w51wo0bN4L3PHjwoGubzMxMJ+/YRkdLlJefffYZKKuYBxPlIF9xcXHu3ZOjlGPNmjVDy5Ytj7rZA58Lyt/69eu7Olcl87g7Mrls2LDB3aO0tNQ9P21Mf+WuRPv373dlKi4udvenPO7WrZuT2yyXkgjUVQKUv+xf7GcVFRWuf5AF+4XX19hf2Y+ov3j99Wi8OJ5cv349GjZs6PorZXjlVFRU5M7ZunXrkf5KedaiRQsnB1gWyhDKD45JqV917doVSUlJR8pXOU/9LQJ+IfDVHhMhNXvyySfxxhtvHFEaWOzmzZujb9++mDp1qvtRrm5VKDwoVKoSLNXNO1yvp+JCpSs2NvYrRaSwve6669C5c2fcdNNNTrH5ykk6IAIiEBYEaEh69NFHMWfOHGcI9grF/tu6dWvcc889RwxKVM5+9atf4dNPPwXlHhOVMipLTzzxRLUNSh9//LGTGQMGDMBvf/tbNGnSxCuO3kUg4gjQYPHmm2/ilVdeccYADlpodOjZsyeGDx/uBji33Xab61OTJ0+OuPoFF5jygPpAdbZdfu2113D33Xdj9OjRuPnmm90AK/geNfmZRpkPPvgAL7/8MnhfToSNGDHCfaZB5mhp+/btePjhh/Hhhx9i6dKl7jQOTrt06eLkH/UiGtz5TkNS9+7dccYZZxwx8FTOd+HChbjxxhvRv39/J/M4oK2ceM7TTz+Njz76yA1K33vvPZx77rm4/vrrnYx87rnn8Mwzzzgd9MCBA+Dr1ltvxWmnnRZShpXLqb9FINwIcDxy//33Y968efj8889d8SiDO3Xq5CaueID9lbKrVatW6NWrF84880w3fqGBp3KaPXs2/vznP2Po0KH43e9+5wzBlc+hbGB/5WQBJwrYX7/73e/ie9/7HhITE/HUU0/h+eefd3rVli1bnMy85ZZbnPypakxVOX/9LQIRTcCsqBGZbGbHlfuhhx6yPXcRaNu2bWDRokXumFmlq10ns3gH7Mc8sGLFimrnFUkZzJw5M2BCs8oik68NRgMTJkwILFmypMpzdFAERCA8CFCG2eDKvc4++2wnJ20A5Apng8QAvw9OZkwOzJo1K9CvXz/Xz20wFigoKAg+5YQ/v/jiiwHzagxcdNFFAZsJPOF8dKEInGwC5vUSOO+889zz/Je//CXAv9l3bPASMCNGYPz48QEbsATMGyfw0ksvneziVuv+rJdN3gXMi6Za+dx5552O11VXXeU4VSuz47jY0wG///3vO7l3+umnBwoLC495ZVXykm3JvCgv+SIPM9QHbOAYMONhwDzSAmY4rFJOzpgxw9X5wgsvDJi31FfubZ4MATN0BXr37h0wA5jL2wz4AfOQCPzv//5v4D//+U/ABqEBG+Q6Zg8++KCry8iRI6vM7ys30AER8DEB8wIK8GUG4wD7BMeBl156qatxcH81w07ggQcecDpNVlZWgP2I31dOHEva6pbA5ZdfHjDPwC99zfuY8df1ed6L4x8z7joZb95Jro/ee++9AZuEC5iRy+lcZpRyZeJ4iecqiYDfCUSsh5Jn7eUSNyYu4fA+V2cmzWVm/9BVmjNDF1xwAbKzs73Dvn43ZcnN6nE528CBA79SVy5Veeyxx9zMWGvzcFASAREIXwKcXfe8Kz3Z6MlNzyU8uPSc3Rtu3hX09KTHkg123JK34HNO9DOXzdGrlMs/uOxFSQQilcALL7zgZqZ79OiByy677Ijewf5DzxH2MS6FW7duHcyIEanVdOXmEg56KHKGn145J5rMoA16J3I5ybG8hE40/8rXeTrgN7lXVfKSbcq8vPx4H3pXjhs3zr2z/enVwPSjH/3IvXv/sL6UeZS9XHJXOVG/pKcD5eygQYPc19S7bNDhlsw8/vjjTn4PGzbM6VxcVjN48GB37jepV+X76m8R8AMBb/lq8NjP6xfB+g31mSlTpsCMSbjkkkuctyD7mBmOvoTBjM7O25DLWSuH8+ASNvZHekXRg4meqDxGvYZ9m3LSJuNdP2V/pt7F0Cv8/pRTTqlyxceXbq4/RMAHBCLWoOSxZ6dm4g8+18vWRLKZKrc2l27tNkNfE1lGRB5r1651ArMq5YcVoJCl+7aSCIhAZBHgYImJitSxEgfDdAdnTCVPOTvW+cfzHe9JI5IMScdDS+eEOwEue+DyqFGjRn1lKSgHMkOGDMF3vvMd/N///Z+bmAr3+hytfIzvw6UkjEXEyaYTTez/XO5V1ZKvE83zeK+jLnciyRusHk1ectDJgePEiRNh3g9uIs48EdwSYcpaXkfjGV9HS4xZx7hNbYLifV5zzTXO6P7222/jrrvucsZ9xuViYuwkLhfmwNibIDha3jouAnWFAMeAX6ffsB+OHTsWZ511ljOQv/rqq275rTcxzv5KQ/HRluLzHuyTHGN615Avjchc+kbDMX8XOOneqFEjh55GYuZHg1awgauutIvqWfcIRPmlykf74T+R+jE4G9euc3axsqX6RPKLhGvokcWZ1507dx63slKTzCOBkcooApFK4Jv0VSpnnoJWE/WtybxqojzKQwSqQ4DGViZOOHEQUXnSiToDPZsZe4ceLpGaOBtPnYCpqpgjx1uvSOz/xysv6dXAmC00vDGmp2fAqqrOlfOkUZLn02vTS/RQ4nNDIxKNTTRceYNRGuQYg4pBfiP5ufLqqncRqGkClftYcP70GqLRl56WtmQNcyy2pJe+rr+yvzEQNx0XPGMuP9NbiXEmGXyb8ZoY/NvzAqehmJMOjFkZ7OHo3dN7P1aZvXP0LgKRQCDiPZSqgswfaRqD+INMQcHOzNl2CgUeo5WZn6kkMZCal3gdd/hgsM377rvP7WbGgI4M8kglkufzOubJczlr5ymTFAoUJMyP5zDZOl0nZHgey8Dd0fgdr+V3nAGkIKLCwBkx3of3Y15efjzOe3B231u+wryZJwWYJ4x4HsvnKR88JziRB1+0tPNcrzzMkzuJUDG2NcCunjzPK4unSHrX8zq+qAR5gtO7j8eddfMS68wyeQKb5WYbMH/em3nwHNaF37FsZMhj/KwkAiJQuwS8vhp8V8pM9lG+PBnAPkr548ks9nMe8+Qfr2ef5veUdbyOfZ6DJH4OTpQdngzg/fniwNzLi/dhXjR881xez+8on7z8mTf/5quqOgTfT59F4EQJcCkDZ6Xff/99/PznP8cPfvADF/CVzzV/z9gPeM7VV1/tjAN8dvlse4Yo/r7xxXP5LLNveN+x//A7r994v5N8tvkddQH+vvL553POY/xN53FPF6Guw+OVf59ZX96P1/P3nIll43nMyzuf59DYwYDWFm/E1ZXBoFkmXsv8WWbm4elSvB/L6B1n/6Nuw7KxbqyH12fJqXL/ZDmYH/OvrKN45fLOYTlYRv7Nd/IOLr+rWC39w3amx9DixYvd88A2Z/mpU1Ul8/gdWbCeZMc6ePKLx1hXtv+uXbtcDZgHJ/noYcG68jsm5sH25rvHkt+xbag38Vxy5zm8h8eI33ltxPaijPX4smxeG3h5MD/m68lrXkv9zasbr2UevI4vloXn89jR9DdeG1x2nsdr+PIS2TA/tvXRyuKdq3cROB4C7Ac0/nAFBncJZyB8epLyeaRe4ckoPtOUUXzu+ByyL3FcyHe+KAt5nM+w11/37dvnisDnlZuh0PjL69nvmLxnnn3MS5RZ7F/sMzzO/spnnuXkcX7P75gnjzMvT84yD68MzJvn8Tr2IV7L/sq/+V11x1u8P+vK+/E+zNu7j1cXsmMZK5eF55GDUh0hYA9dRCcGLrSmClgE/4BtIenqYj/GAQaBtA4YsIc50KZNm8CPf/zjgLk5BizugTtu61pd8DTrLEfqb9tau+BqtvY9YB3nSL4M6MjAigyyaIqhO5+BGf/5z38GzAIdsLgJAXNvdAHhzC3ySH4MRGsxSVw+JggCv/nNbwKzZ892wRa/9a1vuYCODORmwshdY0YdF/Caedra28CVV14Z+OEPfxgwxSxgCsuRfE3wuDwY7I0B4sxKHjDXaxdE3JSSI+d5H0xhcXW1HQ5c/Xm+7TzjymKzkO473tMGaK6sDAzJc5m/7VjgAouec8457juWhaxtZwUv+yPvDLbLAJIsO/Oztf8Bi7l0JMCdCZ2AzeK5wHZsM3MNdQFMyYx8ycjcQ12dbTbgSL76IAIiUD0Cl1qwSvY5yomvS5RNNhP+pdMoG70At6boBGyHowBll+0U5/osj5k7ecCM8V8K9k2Za7Pu7t6mpDnZYJ4dX8qbf1i8mcBPf/rTgC0XcoFqKVMpbynrmCh3LeZIwBQ1lxeDa9punoE5c+Y4WWNKVIByy3aQUgBMR0z/hIqAxb1xzyj7E1989myCxfWZX/7yl4Hp06e7wPP8vWNiQGsGaOU51CuuuOIKt/GFDWIC/J175JFHAl26dA2Y50ngX//6l3t++awzWDPzN2OCe8apT5jnk7sffyfZXxg81nYkc33TJp0CNpAP8LfajF1VVp9l+dvf/uZ+m20JiOubZgT50vnMj7oKf+d5f+pQOTk5AYsbFLAdjQKvv/564Fe/+pXbCIXfUxYwoO27774b+Otf/xqgbsU+z99622XS1curBzcFoD5SObEe1Nl4D+po1B1stzOXpw1k3Ok8h2WfNGmS488y2RITF0SXQa2r0n2oP/HexxOUO7hMFhvJXUcd6FjJBpqB888/353L9mJZKQOpA/G+bHPqhhZOwGXD4NysJ/UrykyeYx4TAW6aYMtnAv/+97+drmpxO9131MmY17e//W2n83llseVyAdtZyn03ZswY9858vQ0PWC7mx2eH96AeartTBd555x3XtuZV5cpgA2svS/cc2e69TgazHagH27LNgO386c6xwbSTt548Z97UP9nOfKapG9oSn8C1114boF5ZVaK+TblO2c3ni2wY+Jy/JZ5uzeu4AYztkHekLLynLfc7ouNXlbeO1U0C7HMcz/A5Z584WuLzywDdHJvwXPYbJspVT9aZh5E7zv7DtHr1ajem4bne+Iibi7DPU/5yfMPnnZtCMU8zoLg+QX2LcpKJ9+VmRxdffLGTa7w/+w1lvcWqdOds2LAhQDnMMRDz8cZeHC/+4he/CNgyO7ehCX97vGQ75wZuuOEG10coy9mX/vCHPxyRNWbccX3NPB5dnmZIc32eefI3gvKHsucnP/mJ22jAyzf43YxETu+iLKZcpn5G2cWNGvidlyh/qR9SVlF3s1hvAW5YUZWu512jd/8RiP61JXuAIzbZThkueDZnDG3nFefJw5lrBpG0juu8jbh9I2d/6OVjRhoXcHrBggWwjuVmuRkXiJZXWlO5Tp2zTlwTy/M5+8gtXBkXga6L9MwxIYM77rgDtoOL2wqbASfPOON0bNu23a2n5xpbulUyiKIpGW4d7VtvveXW3tLKy/JxdojreGkJt07o8mTAa9aB+ZlAcVZrbkFJ13puS8mZMM5WmdKB22+/HbZjCkzRcFvXpqam4E9/+pObyeK9aaVmYrA4XsutNUfYmn8Tgu6d1zMgpBnbYAqaK6cZ55wF3oQhbJcR585JHmTJYJC8P2ctaSlnkDseZ7Ju4Wbk/ud//sfFnmKwO9vNydWNnl4sM/mx3pxpM2HqysU2oFWeFmxurcwX+TAIKC3idO0mbyUREIHqEbBBrptFpwwzg9ExM+NyX3pUMICll0z5cNtks09yGQyDdrOvcnkG+ztn/kwxg+2K6WIIcBkIE2Uo+zsTt9ilvDPjsXvnMcoOymcGt+VnbsFL+Uf5aINTJ8coQyk7OLPIoMc2wHFLTCgjKItMUXPxEehVwbJxtq9Pnz5HZCDvoyQCNUWAzzQ36uCyJM5Wc4k8Z3FtEO36GGNtsC9w1pq6BH+P+dzzXHqy8LfQBtHuej7XKSnJWLhwgfWlruBvKGehGXuDy5uot9B7mNfyt5C/kabQu5lgBnU2Q6/THRggnLqNDWzcttY8Tp2IG2x4ib/3puS7fmiGL/e7P378ONdnn3jiSVdWxvtgWdnn+FtNucHZdhvUOD3CBvZuiRflAD+z7suWLXPHWEaey7/pvUVG1GP44qy9GQlcmcyYcGTmnmVj/2fZeR25mCHBBbJlsGvGdaQsofxh2c1o4sphgzgnCyh3zKjt8qY+whgmwYnL0MiP53EZIuXK8aTjlZeUNTZwhBnhnLy57rrrXB25FIbeC4yTxHbjvcmT8optxGfAdgR03gP87temhlMPpTylPsiZfuqnrA+3MLcBndNN6R1AvZHxucyQBxtQOmb9+/dzbKjv8Rrqanz22EZc3kMdjHopPaf4Tn2Ueizz6Nu3L2xACDMiOr3LjEpOpxw9epTTcfkd+fF54pIeBgmnNwK50ruKujOX/rGM5MvgxTzOdmNf8RK9GPjcmZHRPZc/+9nPHBeyJkN6tNkOo46ZGSydJwQ9AJm3GcRA/ZT9h/Xz4kp5eeu97hKg3sBng8tO+fxQ3hwtUWehjkC5wmeT+gafJY4V+cyakcbpL5RRlMMcE/I5Zt+kvOW9qKvw2WW/oYxjf2D/4rWUw2bUcTKafZn34xiOz7MZpcAYaRwr8jvKN/YJymjKTfYxM5K7fsX+xf5G7yLKNfYRM3I5T1iOybiKhnlStrAs7CNDh54Km5xwHLw+wjqybtSNKH+p11GGmEHJceLvixma3e8MZQH7oJeoX5lB3skAM7g5b1zGjOM9KN8ZK4rXsO9SHvF3zhw3XFkGDRoI2+3OjV1ZP8o+pTpAwDpIRKeqPJSCK8RZI2tGZzGlFdVL9P7hcVqKrdN7h907Zw05K8dZQM6+VE6cKbGOFDDh8KWveK4JBTdTH7xNLC3CvBdn1f/4xz86CzJnhmjF5owUZ21orea13kwWMzblIGDCwl3rzSSxvrSmm0HnSzM6tLxzG1taiG3g5splCl7ADDquLjboC9jOM0fKS+8gWtVpafaSCTA3w8nZoqrSJ5984spCC33wDCgt4Q8//LCboeIMqjczy/Kb4AyYADxirffyvfXWW11enFHlDK0JL/eVCfuAKSUBlsUUFu90vYuACFSDAGUNZdCJeih5t+YMOPsn86K3gidT2NfNiOOOc2a8crJBhvvOFL4AZ+O8RC8NM26777hVtpeYL2e8OfNWWQ6Ym7o73xQtt007ZwCZ6B3BcpkCqG21PZB6DxkBzr7y94q/ZZydpecMf+tMYXfPoQ1UAjbxdMQjh94hnIGmt0zwMz137pyAGRTcTHnlwtILms80vaA8XYHnPPvss+449RR6RdFjmol6AH+feQ31lOBEL23OMlOnCE5PPfWUu4az0J63NL+3CR+XD/tTsId08LXMk/eiV5IZSJyHzEMPPeR0I3rqeIll4XmcnbdBknc4QE8Aeh3xO5u8c95c/JL6AGf76bVlxo+AGeucJzbPo3eOlyg/bPDkrqccqZxC7aHE+1900VR3f5aDupCXzDDvjlMH9eSk9x3fWTfK0ltuuSX4sPvMNvG42mD1yPf0mKCHebCe531JnZDP39///nfvkNPFPO8NytI55tFJGW6Tj84zgZ7mZG0TmQEbILvt1Y9cbB8ok1kvm1QNPuy881k+M5q5z57XGduc7WYDSOeV4V1khjf3fFBm07Ms+Plnf+DzzWeEXkrsR/yd4G9GcOLzQR2ZnktKIuAR4PjNe8a/zkNpz57dzoOGzy693oIT+w2Pm7HXyaXg7/iZ3nfUSTimqpzojcNrOZ4M9pTkmImefpR7Ntn2pct4jLKBXpxeokc2vXuYF72I2E8o0+kpSDlvxmEno+nVx37EFSTBifWnvOb5wYnXMk+bLHD9yisj5Qx/rzjWXbhw4ZFL+D09GumxxXsFy3+Wmf2bHpH0sKI3O1l6HlleJvSqpKfl73//e++Q3n1O4Pima+xJjNRkgw1XdFqZzSXwSDVoNWXiLE3lRMssE6+lBTc4cRbOXLudhZuz5ZyJZOLMk3VIZ3F+5ZVXnIcPPZSYOIvFRCsxZ4c4c8RjtDDTMs1ZSt6Ls/70PrIO62aYaF3mzAzvR2swzzEB5WacaB2nhxRn7mhxZsA4zg7RC4BeW5wB42wVrcn2DLtZLK/OtJrTS4EzkeTCxHryPL4421pVojWdiTNSwYn50aOL1m/OXLIcTCw/PSLo1cTZOG616yVew0QLO63enAlgoscUmZhrvvNecgf1jwiIQK0RoAw4WuLMuZcoTz2Zwr5OucZZMO4KVTl5sqPyccpTeiRSZlB2eIn5cpaNcQ5sEORm2r3vPNnBGTzu2uIl3p/yhN4ZweX0vte7CNQEAcbf4O84Z475otcQn2POevO3lzPVs2bNct661A84G83Zaf5mcwbb846m9wh/05999nk73t3NTlcun9dv6BlCfYO6AhOfdSZez1ltejQx0TOFHjAsixeHh/2ZugA9Xqgn8G/qDtRB+LtLzyDqPPRMpDcM78XEWXIm3sP77A588Q+v92QF82d9qKfQe5llpbeTl7w+6/3tvdM7iZ7Z1GHYlznLzcS86SFJHYqs6XFNTwDG5qGsYZnoFUA9gceYqH/UdiosLDiiL7G+wbqR13ZHKxPrwFSVrPJ0MJ7j6aNkQtlGfZB86RXF5473YTswH3pvUl7Sq54puAzkQ4Yso00wuHeytFARLqA4y89nhh5QvD+ZU2elPKf+ZoavI6y99iR7tpvXBmwrthn1anpXeIllp+cS25vPB58VL9ETj94Z9C6hlx+3X2ef4bNqy+Ocbkr9mvoxn2uWj95ULJ+SCHwTAvn5BUf6G2VlcAruK8HH+Znyj32ML69vBp/j9WE+5+yvfF7ZR9gfqcPwmadXjxmWj8hT5kWPO64eoZc1E8vgyRBe7/WTadOmOVlN/Ya/P+wjlPsbNmxwfYRlokekGbidjOD39Ab0xmJef6Vs54oUlo+J4y3qUexrnpzhcdaXY0d6NNLDnB5KXuJ4jd/TU4qyg96OlCv8zPqzLJTnlDf0GuNvIVfGePXy8tG7/wjUGYNSZeHBB56disutjjfxx50/aOw07FB33323c7f2rmeefFGpYmeiwGAn8hQHdkx2YCae533m33QdpPCgQkchwMEUjTF8cWtKuiLy3hQYTPxxphHHS6wH60MF0bsfl5NRmNGlk4qDl5g/l+wdT2JdmTfrwfpUlXgOjWn8gQ++D8+lcYgChkprcPLKSAbBbcDjnqIgARRMTJ9FoPoEjtaHg3Nmf6RyUlVi/2QedGGmchOcvMHuN+m3NGRzgM3lbMHu1rwHB5c0UHEQFZw82cHBfOVE2UHXdSURCBUBLr3isk0q5l7i5AwHxHxxsobLI8xryS0rp1JNgxInf7gUnBNJnPihMYB97eOPP7JlTkOcMcjLz3v3nvWqfid5Dgfw3mCef3Pg4PVD73eVfYkDGfZLLsMzjyDYbDxPd4n6AL+jIYD6BX/PqZ9UlZgXy8Tv+dkbAJn3y5FJIcqFyrKB51aVWB4OkNiXydBLrJN5Z3l/uneLj+aMSNRraByj0YIGPIYEYDpamd2XIfpnz569xuyQ01k4ERgs+45W5xMtCp8X1pX34CCOhiV+9u5D/Y9ym7tQ8Vzy5Hcc7HL5DvVPr3yUt0zUY8mRcpMsaajhc+Plyfx4DQedNOjwGeb3fEaYaLz0Bq38m8c9/Y1/e4nPDHVRtlGwoZHfU79lYjk5GcvrKfdZnuCy8LlmeajPsizU6b1n3GWgf0Tgawjs2rXbTVRTPnHMFcpEIzuXBDPRIMuJda+/8t3rr3zeOV5kH/X6K8eQXPrqJW/8St2GspzjPPYnGqQq9xH2E47D2EfYV3gvr7+2sfEW//YS71dVf2VfNO8o99viTVZ411BGMLHM/J3j/cz71k2kePXj9179aMRiWVim4HvzHCV/EfC9QYkdhh0u+EfPa8Jv8nBTcbLAsc44wkEQf0jZsahUBndo5s0fTXamyvlT+fO2nPTK4L1zjS6NPFwbyzXv/EFlLAauv6Wg4Zperg3mvfg3lSsqrhSMnvXZu7enUFKI8DsadDzF1Lvf8byTHY1YVCaOVm4vH0/x8P6u/F6ZBfNm8soafD7rqCQCIlBzBLz+5ykWR8uZso2GHc9j8Gjnsd9W7qfePY52zdGOU3YweTLBO8/7u3K+3nHvOu98vlcuU/B3+iwCNUGARhcq9d6gvXKe1DUY54ixLTjD68188znmbzg98ng9vZdoNOVEEvWJYyVe6z33wedV7hv87mh9gDPQHERx5tqCY39JJ2A+nt5ytOuZNw0Q9FjhhBevoV5BWUGPqWNdx2urSl6djkdHYRwmxh2hAYqz8xYY182M0yDNOD8nouNUVaZvcoztSC8eTpwxblJV7fFN8vu6c1lHMmMsFsbGZAwXsvMS24By0dOreC6voZ5I/bOqxDJT7vOZYAwUercH58nv+UxXpdN69wnO92jPAXVR9oVgfTX4Oq/9OJhlHD3Gh2I7H29ZgvPSZxGoTIB9YdWqlS4WEeUV48KGOnn6lgW+xqXmFcg+yOfbS8H91ZOF7Ac0vngrXLxzg9/ZjzgmZMw2er5W7iOeDKgsj6rSmSqf493HGzt6dfCOe++eLOL3jKNJLyoLu3LcZfHy0bt/CPhm5M5OcbSOUdXxqo4FN6v34+b9OFKJtFgJTnnizA+tz1QWaWjhbBoHYd6LxqBgd18vX/6oV/UDzO+pENH12+KDuEDfDzzwgHNBZ3BsKoIWl8DNvFNxoWWY7t6cZaIy592X73SxpOLFxHKxg1PxpUtj5UQLN8vqpcp15nHOxnKm6FiJiigNYlRKGDguOPG+LCPLUVX6unao6hodEwER+GYEOMNFQzT7Jz0mjpZoxKaMoqfFsRL7beW+W/nvY13vfccZb8oOzmAFy0zmxeUsnN2il2ZVqfL9Kv9d1TU6JgLVJcD+wYkeelMcLVFvoELO32tvdpnnckDBmWcaZu666y4XTJYDmzY2c3wi6XieeZ7D/GkQYP/nLDdlQbDewN9oehlTt/CSN8Dh9Z4exCVLFgfSncfffZ5DYwN1Iu8c7/qq3iuXlzoR+zeXWtEwUznRC4cTefTM+fWvf+1mxDmRxuDlNDaQp1fmynlXzuub/H28eTHALz0ouSSEy9CqYvB1eX3d9165ycrzyqTnAz0L+CwGtyN1Py498XQ55s02YttXNSnI6/ls8Hxey2e2cp4cAHN5nWcY9crD98plr/y3dy65eMtmKNcrpw0bNjiPKxpWWQY+Cyx35bJQF6a3O89REgGPgPfcee/e8eB3Pk8MQM2+YzHfXDiR4O+9z8fKwzvneN7pxU35z8Tnmf2P/TC4v1KucZkny+YlfvbGdd4x7536Evsr+yLzZFkr9xHqURwHVtVHjrduNDzR45L3qSqEATdf4e8FDXM0KLEsnhEruH70QOW5RzNKefXSuz8IRLxBiR2UiQ8+O1vlxB8ydtDKllnvb3YwKkbB6f/bu3+cOpIgDOCs9iLcwUvOIZAISFYiQIIbQLIhAQEBMeIKhN6IQyCtQ2t9jbW0/Wvps4vhobVZ2/J7r1oa5jHT0139TdfX1dV/xjX3NbA6OoJpwWbrIAn7GsnLSGNdysVItF7U/gR1PWpkRBIMuVWB04bTSMNtSrD1qWPTtjkDKiOHHEg6X7uDUHzVQOevBiN3DN04lEybzNpcaTOA48mW38NYa5/p4tJBTIyz7L2gvMqtXEKdGqksCbBi5CkbIkse0tF5tT7eutsa5CV4P/W9STekF9zqc/27EWgEvh4BX0/y5TMOZMtqzYKsgTPYtOXb29vpqBZ/Vah6v9TP8Oiq2aBV32saOkk4ivGDOxJ00oz+m71hqUYNST/55V6Vp+aR+31uBL4FAuqfukpXtLe1rZe+JVnaWwNO+/v7O2NT1ifZZs8Y7TWjXFufJUhPIo5/ojfaxNpO1ra4/hY/To3oiWc5tcjCsLdcXruewMFkrxpHdRpEvwwKGZQS6KUOUO7Jw29tf9W/pJ1zyrFs743Csx3YaONjINN+wgX+Zwv5SprBOx2Wu7u7aZuYQWNmF9zj0JFPZEqezrHz5FtxqnH8ll8N+X9VmuIZ0PN1OV+Rw1HHx8fPHDbBA/6r+Gj5nmr+wau+d1jb85IjDVeb4RbbzLPsNXhkGYxryiwfz9alke4Jyscpp056rzrcsf/cV0727MOwFSsWkW/5PpWTzELK77d4lomq5+xQX4CiH4IOK1tW3dJBtTyU3cjGjf0tHkeAMlsumrrteoftQyD6qeTqZRwWL9WLDx/+nkt9Oajt5Wp2z1InU6dX6WuNm3gV9coziSsebjeYbrmbel/5laNF+6GuyzP5OkvvJX3ldNU3pOf2vcs2KOShS3TEBIXI5Hp+L/UVP8hPqPzoN7tMWSxr02Zo1wQ8Q1/1HQ0GmJlpcoKy1O0GOInxkX1/X3ovM8H+szEI/PrHCOtYGo4cHQ6bj1nrSVE1jMhF5wjhLImDVAAACgJJREFUUFhGH4+tZVucMZSUcUJJbJ5JmfaHoUWBKJ3nOUQYM0bGKIKGjwFmdIgic/hIUxxGDS+sxk5nDTnwRDMypEVG6/05mjh4TGXklKKolDuB0iECZZA/I8pIjjQ1tKZ3my7vvrKZ5i1v8RgA8EAw5DUFUj4adLIiAHJJj6zSsxGc8poCz6gzQqWTaXSLQQdPaTIAlIXByOC0IaJ7nrNxbggJNggGLka0GDo+JylPziQEDg+EQ174y48MPOHej/dmyR+S9BsRI0/PhRCDV58bgUbgyxHQAfv48Z+po7iAfuFB+okf6KJNIzmUbbK6nKGEk/AgY16nw/+MCTyAcziozbCUnllFWf6BXz2XjgCJ6TteYnA5/CaTNHGLTgTucP59TBPn3MKV0tG5lQ+uMRKG0/Ep3rFPAY5iXNowHHeRzdGhEfhWCKh/2lZtlvZMe6njbfaOTjl7RF3UJo6vyT7TJfVS/aVv2mqGe3UWkZPxbjYRg5ydog2kN2wI7b49dHT06YtOi7aYHGwA8pHH3hf0gA5qp9kEGTGm537TJ04tdgo94vxKWyuO9JTJ82wcHRadfkveXNdWc67tjkEus5ToJJ1OJwXPkAVnsJekzbENA7+VG470n70mDzaN/5XbfTqunCmvfUH8zyGHj7wDWLG7yMVWgZF8x1eMZrrel3zZf/Ktthe8ySsd+TrsT5nZ1ZaV4TXv7P3oJLGRvGMdRBifnZ09GTAjG2zJyyHGXjM7BweHx5R1fNl38pZZDJkNgLvwnGeVTfBsZrc7O/A1nMTFk94fvNQb8qoPBhbZjxwzuF59VAfgA4O8I9h4996nMsKB/H7D2LNsRDKqE7VeKg9ccTi7j60rP/FsrK2NEEce8icz/BxkZWPrjOoUi2/Teu0He1tdw+vKSDY2srzZwOJ12F4E1F31Qj3CUfRV3cIB6j49zfHu3V9Tl96+/XPy1Pjy45NBqvS3xpczZ33Hteo0XhLkQw/0JfGKPgnOVd9xLv3zLB5la9ATfIm76CrO0V9Tx8V1kFt9loa6jBvwjXg4i27QOXwffc3b1nfEtXgKX0qXvtIlHOA3GX3khP6Ef/E0PSQf3lIWDnFtDAxcTx9aHmSSvnzYiPJxpq/is/0cOF/a8jHwIH/xlM//OBkPdNh8BNbWoXR5ebkzPhM9K63GUaBMOjaUGuHc3NxMQ4DSUgCNNePE6CGPMQVBQBSB4iIQREDpjARRDJ5Yhp9743OlU5HFZQhSJAppN3wdJsRjt/vxOe2Zjkbx/v5+HhpYyqVxlgcFi9xk52Sh4DpiZj7pGOlUMTw4k3y1g+GBXHT2/OaFHp+cnA04BUaCJycnn4gQBhrn/eEwQ25IEVkxchCH5XWZhiyumURwki4HGDmtqTe6yXiSlzIykBAV44Tx4llxTd1XTgYFzMhuXW0MZk4rpAcvI2wMF9ekRXakND7JOd+Z9LwfAbEi5g6NQCPwOgTo6Js3v01+wjs6K3gGH9B1nWE6SN9XLU/FheLRXTxGbxk9+AUXhYvxB/4Vn0EiL51OHaQ8R6/xj4NceBRv4jCjXLiDPJb72tdDBxWHM2zwg84UrtRxwyEc3nhYefAoTmdc4jyyrRrpex2K/VQjsDPbZbNq2APaRqO36q3BLW324+PjrM/2KmKkL4O2TeeBfcAZwSivtoD2VWdIe2uARR1nu2gndQS0k74Gqy2OHrIL/DbLh36IL9A1eshm0Y7ScfaNgaG7MeNHu06XLNfyVbUqB/3RGeCYUDadHXaL5WZxILsmLx19uken2UqcCMqh86Qj9jBmuKS9J5N0xKP/nEtsFLzB9sBJzvK2N8fucFaRmdOCTUFmXET/zWaxXxU8DZKRlWOCLLjKNTaJjh/+gZ8DpywDp4V9LG0GzpkhTga5cJJ3ATdOE/wi3/Pz82cOQ++JDcQmrJynHMrKDksesU1hAm9Yek4+0ol9RH5cS3YdUPapOgIr9UQnj8NmfKp7LufxTsmsnqhnnocJXoY9BxLsE3TCpQsrdRl2nJZk8oUqZfU+YcRhmXpJPrKrw5xuPlTDxlOPtA/Ko/6RXb3n4ISDQU0YKad3a88kG7sL7FDpwU7bQRb1hw14cHDwKV5k7/N2IoBnfVzA9iD4Rz1SZ9TN6Cuu0N/xdTdf47y4uJj8EMTUaTpBT9gQ0blwFNtEWuwO+qM+GyxXt+kG/eXMwgnhNzII+kf0Fedy8OBMNor2gvOJDuM3jh/P6BORA/dIV7+UHJxXbJroK27V73JduQ3Oe87zeP7w8HA6zMjK0ay/i0+jj9LUf+Vs5tRmaymL9oQzTTshP22K5dkc+NocgwLaAdxhtqGZogL5ycP5RKfJAjNlN5Ggbi4+H+g/G4vAL0Ohns71XZOiUmoGhZAiULjMZqF0lFJD6r7/NXJIx3WHxtAzDBn34oiSJmVkTHheHPc1iFFqaUqDDOIKnmeEiCu47vnEqfKJk7TEFcfhujPZElyTbuLLO/GlL6TcNV6ej6wps+viwQIuCeRNmTwj35RZfu6lrPLzvDgJyUfcBHGWuEb2ir98/K9DuiyTPKqcSbvPjUAj8HUI0DF6TMfoMp2lW+HHVfwhB/Git1U/xaefuNL9pBcuxFl5LtzhXngh0le5XPPcMk74Sf7hjnCUPCo/ycM9R3gzefW5Efg/COisp36lzXIW1H/1LfX/pXaL0X03HDoMcB2PZR2NPqjXSVtdpm/yyACZPLXFrjun/VzqoXuCtKIrSTeySr8GadC5pElG+dBL16UTXVTO2ATJK/nFppBfsJHGf8Ujj3jSJov8yFLTIbs0g4f/2RKJH/mSb3Dy/zKkrM7BJvHyHqUr+J98jsRJeuKkzNISyJUyuxYOFlda7ictMitrfdY9cSKHOMs0xFE+caSbdJxrPsE08uYsP/k6KsaRWzzXU7ZgVDFPXRHX+yWTI3JHrsguXrBxTghGZPEMjFO+Wm8Sv8/bh4A6kvqmjgjqSQ7/57o6o36pQy/pqzqZOi1e9CS6Jj/pqcvuS09aVWfkmXvup94nTtXF5CEN6dZ8pONZaUSnXavhS3REnORddRpHBj9nAUbyqnKTy/OwSbzIXfWQ7MnHM8qUeMGgyt6/NxOBtXUobebr6FI1Ao1AI9AINAKNwKYgYGTbLBFGvNnFRoQ5kozoXl1dbUoxuxyNQCPQCDQCjUAjsKUIfB4S2FIAutiNQCPQCDQCjUAj0Ah8DwQsxfe5e0sfLMG07MESraOjo++RXafZCDQCjUAj0Ag0Ao3AD0WgHUo/FO7OrBFoBBqBRqARaAS2BQF7Vtg35/3Y88b+NWYr2efH17o6NAKNQCPQCDQCjUAjsO4IfN5AZ91L0vI3Ao1AI9AINAKNQCPwEyGwt7c3Pytvs9Xr6+u5Yerp6elPJGGL0gg0Ao1AI9AINAKNwOsR+BfEhSZukc/jwQAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "The aim is to make an overlay analysis where we create a new layer based on geometries from a dataset that intersect with geometries of another layer. \n",
    "Here, we are going to clip land use features from `dk001l_kobenhavn.shp` within a box `CPH_Box.shp`\n",
    "\n",
    "Here is a set of functions for overlay analysis.\n",
    "![image.png](attachment:image.png)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import shapely.speedups\n",
    "%matplotlib inline\n",
    "\n",
    "# Let's enable speedups to make queries faster\n",
    "shapely.speedups.enable()\n",
    "\n",
    "# File paths, make sure you know where you data is \n",
    "border_fp = \"Data/Data2/CPH_Box.shp\"\n",
    "grid_fp = \"Data2/dk001l_kobenhavn/dk001l_kobenhavn.shp\"\n",
    "\n",
    "# Read files\n",
    "grid = gpd.read_file(grid_fp)\n",
    "hel = gpd.read_file(border_fp)\n",
    "\n",
    "# Plot the layers\n",
    "ax = grid.plot(facecolor='grey')\n",
    "hel.plot(ax=ax, facecolor='None', edgecolor='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Here the red area is the borders that we want to use for conducting the overlay analysis and select the geometries from the Polygon grid.\n",
    "Always, when conducting GIS operations involving multiple layers, it is required to check that the CRS of the layers match:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ensure that the CRS matches, if not raise an AssertionError \n",
    "assert hel.crs == grid.crs, \"CRS differs between layers!\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "let's reproject the ``hel`` geodataframe into the same projection as ``grid``is in\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hel_new = hel.to_crs(grid.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s do an overlay analysis and create a new layer from polygons of the grid that intersect with our CPH layer. \n",
    "\n",
    "We can use a function called ``overlay()`` to conduct the overlay analysis that takes inputs like\n",
    "- the GeoDataFrame where the selection is taken, \n",
    "- the GeoDataFrame used for making the selection, and \n",
    "- parameter how that can be used to control how the overlay analysis is conducted (possible values are `intersection`, `union`, `symmetric_difference`, `difference`, and `identity`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "intersection = gpd.overlay(grid, hel_new, how='intersection')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s plot it!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "intersection.plot(color=\"b\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let's check whether the attributes are also there!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "print(intersection.head()) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "==>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s save our result grid as a GeoJSON file that is commonly used file format nowadays for storing spatial data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Output filepath \n",
    "outfp = \"CPH_LU_Box.geojson\" \n",
    "\n",
    "# Use ESRI Shapefile driver \n",
    "intersection.to_file(outfp, driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see what data formats `Fiona` supports  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fiona\n",
    "fiona.supported_drivers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "### Notice:\n",
    "if your dataset contains multipart objects, you will have to convert them into single object features.\n",
    "\n",
    "So, in that case you will have to run the following function on your dataset.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "#import fiona\n",
    "#from shapely import geometry\n",
    "#upcast_dispatch = {geometry.Point: geometry.MultiPoint, \n",
    "#                  geometry.LineString: geometry.MultiLineString, \n",
    "#                  geometry.Polygon: geometry.MultiPolygon}\n",
    "\n",
    "#def maybe_cast_to_multigeometry(geom):\n",
    "#    caster = upcast_dispatch.get(type(geom), lambda x: x[0])\n",
    "#    return caster([geom])\n",
    "#intersection2=intersection.geometry.apply(maybe_cast_to_multigeometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What's next? Aggregation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Aggregation means merging Geometries together by some common identifier. Suppose we are interested in studying continents, but we only have country-level data like the country dataset. By aggregation we would convert this into a continent-level dataset.\n",
    "\n",
    "Here, we will aggregate our land use data by type (column `ITEM`), i.e. the grid cells that have the same land type will be dissolved.\n",
    "\n",
    "For doing the aggregation we will use a function called ``dissolve()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Conduct the aggregation \n",
    "dissolved = intersection.dissolve(by=\"ITEM\") \n",
    "# What did we get \n",
    "print(dissolved.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let’s compare the number of cells in the layers before and after the aggregation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "print('Rows in original intersection GeoDataFrame:', len(intersection))\n",
    "print('Rows in dissolved layer:', len(dissolved))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "What happend? ===>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "print(dissolved.columns) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "As we can see, the column that we used for conducting the aggregation (ITEM) can not be found from the columns list anymore. What happened to it?\n",
    "\n",
    "It is now used as index in our dissolved GeoDataFrame.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "print(dissolved.index) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Now, we can f.x. select only such geometries from the layer that are ”Green urban areas” :\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Select only geometries that are Green urban areas\n",
    "sel_green = dissolved.iloc[15] \n",
    "# See the data type \n",
    "print(type(sel_green)) \n",
    "# See the data \n",
    "print(sel_green.head())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select only geometries that are Green urban areas\n",
    "sel_green = dissolved.loc['Green urban areas'] \n",
    "# See the data type \n",
    "print(type(sel_green)) \n",
    "# See the data \n",
    "print(sel_green.head())\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Now, we have a ``Pandas`` Series object containing basically one row from our original aggregated GeoDataFrame. Let’s convert the Series into GeoDataFrame and plot it to see where those green areas are located:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Create a GeoDataFrame \n",
    "geo = gpd.GeoDataFrame([sel_green.values], crs=dissolved.crs, columns=dissolved.columns) \n",
    "\n",
    "# Plot the grid cells that are green\n",
    "ax = dissolved.plot(facecolor='grey') \n",
    "geo.plot(ax=ax, facecolor='green')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 🏋 Exercise\n",
    "\n",
    "Select `Industrial, commercial, public, military and private units` geometries and plot them. Add your code below!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}