Practical-07-Mapping Data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><h1>7SSG2059 Geocomputation 2018/19</h1></center>\n",
    "\n",
    "<h1><center>Practical 7: Making Maps</h1></center>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with GeoPandas and PySAL\n",
    "\n",
    "We've worked -- and will continue to work -- a lot with pandas, but by default pandas doesn't help us much when want to start working with explicitly geographical data. Ways of working _computationally_ with things like distance and location are increasingly important not only to geographers, but also to data scientists, and this is where we start to move away from purely aspatial statistical analysis into the foundations of a more _geographic_ data science.\n",
    "\n",
    "There are a huge number of modules (a.k.a. packages, a.k.a. libraries) in Python designed to help you work with geodata, but we are going to focus on the two most important higher-level libraries (since they also provide 'wrappers' around some of the lower-level libraries):\n",
    "\n",
    "1. [GeoPandas](http://geopandas.org/) -- which offers a pandas-like interface to working with geodata. Think of this as your tool for basic data manipulation and transformation, much like pandas. You will almost certainly want to [bookmark the documentation](http://geopandas.org/data_structures.html#geodataframe).\n",
    "2. [PySAL](http://pysal.readthedocs.io/en/latest/) -- the Python Spatial Analysis Library  provides the spatial analytic functions that we'll need for everything from classification, clustering and point-pattern analysis to autocorrelation-based tools.\n",
    "3. There is some overlap between the two libraries: both can do plotting but, for reasons we'll see later, we'll normally do this in PySAL and both can do classification (remember, we did _quantiles_ with pandas!) but (again), for reasons we'll see later, we'll often do this in PySAL from here on out.\n",
    "\n",
    "PySAL is complicated enough that the best way to understand how it fits together is to use an image:\n",
    "\n",
    "![PySAL Logo](http://darribas.org/gds_scipy16/content/figs/pysal.png)\n",
    "\n",
    "We're going to concentrate primarily on the ESDA (Exploratory _Spatial_ Data Analysis) and Weights components of PySAL in this module, but you should know about the other bits!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Required Preamble\n",
    "\n",
    "It makes life a lot easier if you make all of the library import commands and configuration information (here having to do with `matplotlib`) the first exectuable code in a notebook or script. That way it's easy to see what you need to have installed before you get started!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "mpl.use('TkAgg')\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.cm as cm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings \n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two new packages we're using this week are PySAL and GeoPandas:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pysal as ps\n",
    "import geopandas as gpd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a GeoPandas DataFrame\n",
    "\n",
    "There are two primary ways in which we can create a GeoPandas DataFrame:\n",
    "1. Use an existing Pandas dataframe that contains a `geometry` Series\n",
    "2. Load a file containing some kind geometry (and possibly data)\n",
    "\n",
    "As the LSOA data we have been using already has a `geometry` series (column), we'll use approach 1. first, then look at 2. later.\n",
    "\n",
    "So, first let's read our LSOA Data as a Pandas DataFrame as usual (though you could _also_ read this directly into geopandas): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdf = pd.read_csv( # pdf == pandas dataframe\n",
    "    'https://github.com/kingsgeocomp/geocomputation/blob/master/data/LSOA%20Data.csv.gz?raw=true',\n",
    "    compression='gzip', low_memory=False) # The 'low memory' option means pandas doesn't guess data types\n",
    "\n",
    "print(pdf.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice in our Pandas DataFrame that one of our columns (Series) is named _geometry_. This contains information on the shape of each LSOA; let's take a look at it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pdf.geometry.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So each element of this Series has text indicating the type of shape the geometry applies to (e.g. _POLYGON_) followed by a bunch of numbers. These numbers are truncated in the view above, so let's look in a little more detail at a single entry in the _geometry_ column:    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pdf.geometry.iloc[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we see that what we have is pairs of numbers, each separated by a comma. Each pair of values is a set of corodinates for each the [vertex](https://en.wikipedia.org/wiki/Vertex_(geometry)) of the [polygon](https://en.wikipedia.org/wiki/Polygon). GeoPandas can use this information about the type of shape and the co-ordinates to plot the data as a map (and perform other spatial analysis functions)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now we can create a GeoPandas DataFrame directly from our Pandas DataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf = gpd.GeoDataFrame(pdf) # gdf == geopandas dataframe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to tell GeoPandas which Series contains the spatial information: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf = gdf.set_geometry('geometry')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did you get an error? Can you work out why? [This thread](https://gis.stackexchange.com/questions/267801/csv-to-geodataframe-how-to-have-valid-geometry-objects) might help you solve this problem below..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 1\n",
    "\n",
    "The problem here is that _geometry_ column is a `string` type, but it needs to be a `shapely.geometry` type (details [here](http://toblerity.org/shapely/manual.html) if you want to find out more). We need to use pandas `apply` functionality to apply some function to every Polygon string:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Geopandas geometry type: ' + str(type(gdf.geometry)))\n",
    "\n",
    "from shapely.wkt import loads \n",
    "gdf['geometry'] = gdf['geometry'].apply(???)\n",
    "\n",
    "print('Geopandas geometry type: ' + str(type(gdf.geometry)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What have we done?\n",
    "\n",
    "We started with a regular pandas data frame and now have somehow ended up with a geopandas data frame, but how?\n",
    "- `gpd.GeoDataFrame(pdf)` took a pandas dataframe (pdf) and passed it to the geopandas 'constructor' which creates new geopandas dataframes.\n",
    "- However, that step doesn't mean that we yet have a _geometry_ that can be mapped/analysed, so we still need to tell geopandas where to find the geodata.\n",
    "- And it turns out that in the 'raw' form that we loaded it the `POLYGON(...)` data is treated as a String. This actually makes sense since pandas doesn't know anything about geography; so it gets to this column and 'thinks' \"Ok, well I guess this must be some kind of text...\"\n",
    "- And geopandas responds with \"Woah, I don't map text, I only deal with _geometries_.\"\n",
    "- So we use the `shapely` library to convert the String to a Geometry object and we use `apply` (a pandas function) to apply it to every value in the `geometry` column and overwrite the original string.\n",
    "- The last thing we need to know is tell geopandas that it's now a 'proper' geometry..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf = gdf.set_geometry('geometry')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully there was no error this time!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A Quick Revisiting of Inheritance\n",
    "\n",
    "Lt's check we understand how GeoPandas _inherits_ functionality from Pandas. First, let's check what data type `gdf` is using the [`isinstance`](https://www.w3schools.com/python/ref_func_isinstance.asp) function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if isinstance(gdf, gpd.GeoDataFrame): # Is gdf a GeoDataFrame object?\n",
    "    print(\"\\tI'm a geopandas data frame!\")\n",
    "\n",
    "if isinstance(gdf, pd.DataFrame): # Is gdf *also* a DataFrame object?\n",
    "    print(\"\\tI'm also a pandas data frame!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So this is interesting: `gdf` is _both_ a Pandas DataFrame _and_ a GeoPandas DataFrame. This means we can use both Pandas and GeoPandas methods on `gdf`. Here are some Pandas methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"What is the lsoa11nm column type: \")\n",
    "print('\\tNAME type: ' + str(type(gdf.LSOA11NM)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(gdf.Owned.describe())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.sample(3, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can get a bit confusing, but the next two plots demonstrate something quite interesting (and powerful):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(12, 5))\n",
    "gdf.MedianPrice.plot.hist(ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(12, 10))\n",
    "gdf.plot(column='MedianPrice', ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully driven home the fact that we can use _both_ the new geopandas functionality and the existing pandas functionality with our geopandas dataframe!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making Choropleth Maps!\n",
    "\n",
    "So now we know how to make something much more interesting by telling GeoPandas to colour the polygons according to some value, like one of the variables in the data: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.plot(???)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cool! But what do the colours mean? We can add a legend like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.plot(column='HHOLDS', legend=???)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding things like a plot title or legend labels isn't quite so easy. We'll take a look at how to do this below, but if you take this simple approach for plots in your reports you will need to make sure you are very clear in your figure captions what the figure shows (e.g. including what the units for the numbers are).   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also change the shade of colours used using the `cmap` (colormap) option: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.plot(column='HHOLDS', legend=True, cmap='OrRd')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a full list of colormaps possible, see the [matplotlib website](http://matplotlib.org/users/colormaps.html). There are [numerous resources online](https://tilemill-project.github.io/tilemill/docs/guides/tips-for-color/) to help you think about the most appropriate map colouring for the data you want to visualise. The [colorbrewer](http://colorbrewer2.org) tool is particularly useful and [this overview](https://www.e-education.psu.edu/geog486/node/1867) is also helpful. The viridis colormap is [used by default](http://bids.github.io/colormap/) and one of the best all-round colour schemes for [numerous](https://stats.stackexchange.com/a/223324) [reasons](https://www.youtube.com/watch?v=xAoljeRJ3lU). \n",
    "\n",
    "The four main matplotlib colormaps are:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![colors](https://matplotlib.org/_images/sphx_glr_colormaps_001.png)](https://matplotlib.org/users/colormaps.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But there are many others:\n",
    "\n",
    "[![colors](https://matplotlib.org/_images/sphx_glr_colormaps_006.png)](https://matplotlib.org/users/colormaps.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternative to the continuous colouring approach used above, we can also classify the way colour maps are scaled using the `scheme` option. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.plot(column='HHOLDS', legend=True, scheme='quantiles', k=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `scheme` option uses functionality from the [PySAL](https://pysal.readthedocs.io/en/latest/library/esda/mapclassify.html) package and can be one of the following three types:\n",
    "1. `equal_interval`\n",
    "2. `quantiles`\n",
    "3. `fisher_jenks`\n",
    "\n",
    "In each case, we can provide an additional parameter `k` to the plot method to specify how many classes should be created. For example, above we used `k=3`.\n",
    "\n",
    "You'll note that the legend is overlapping the map, which is not ideal, but we'll deal with that later. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 2\n",
    "\n",
    "1. Create an equal interval map with 6 classes for the _SocialRented_ variable using the _magma_ colormap. Do not include a legend. Set the figure size to 8 by 6, and the dpi to 200. Try again with the dpi set to 100.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Create a map of the _Owned_ variable, classified into quintiles. Change the colormap to inferno and include a legend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Revisiting Classification\n",
    "\n",
    "Classification matters. So before we continue with learning about the various mapping options we need to think more about this. \n",
    "\n",
    "I realise that this topic can seem a little... _dry_, but perhaps this video will give you some sense of how the choices you make about representing your data can significantly alter your understanding of the data.\n",
    "\n",
    "[![Do maps lie?](http://img.youtube.com/vi/G0_MBrJnRq0/0.jpg)](http://www.youtube.com/watch?v=G0_MBrJnRq0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PySAL has a classification 'engine' that we can use to bin data based on attribute values. Although we have just seen that it is possible to do classification _without_ PySAL, and it is also possible to do some mapping _without_ GeoPands, the combination of the two is simplest: PySAL has additional classification methods that are _specific_ to geographic analysis problems, and GeoPandas is just plain easier to work with. \n",
    "\n",
    "Let's see how (spatial) classification in PySAL is equivalent to the (a-spatial) classification we've done previously in Pandas by calculating and mapping quintiles. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, calculate the quintiles using the Pandas `quantiles` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Quintiles = gdf.Owned.quantile([0,0.2,0.4,0.6,0.8,1])\n",
    "print(Quintiles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that quntiles split the data into five parts containing roughly equal proportions of the data. \n",
    "\n",
    "So from the output above you should be able to see, for example, that in the lowest 20% of LSOAs fewer than `199` households are owned by their occupier.\n",
    "\n",
    "We could also [calculate these values using PySAL](https://pysal.readthedocs.io/en/latest/library/esda/mapclassify.html#pysal.esda.mapclassify.Quantiles) (as `GeoPandas` does), but PySAL expects `numpy` arrays instead of a pandas Series. Fortunately, it's easy to convert one to the other since `Pandas`/`GeoPandas` are built on top of `numpy` arrays...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mapclassify as mp\n",
    "\n",
    "owned_q5 = mp.Quantiles(pdf.Owned.values, k=5) # Classify into 5 quantiles using the 'values' because that is a numpy array\n",
    "print(owned_q5) # Print summary metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's more information in PySAL's output and we can see how many LSOAs have fallen into each class - it's roughly the same number in each class, which is good because that's the point of quantiles. \n",
    "\n",
    "To visualise this distribution numerically, let's plot a histogram of these quantiles using a Seaborn `displot` with the `Pandas` DataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.distplot(pdf.Owned, bins=Quintiles, norm_hist=True, kde=False) # Notice norm_hist!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this histogram we forced the y-axis to be the 'frequency density' so that we can see the relative frequency of each of the classes (known as 'bins' in a histogram and shown as the height of the bars). For example, we can see that the bin with highest frequency is the third quintile that contains all the LSOAs with values of 'Owned' between `291` and `366` (where did I get those values from)? \n",
    "\n",
    "You should also be able to see where the first bin stops at around `199` (because remember we saw above that the first qunatile is all values up to `199`). \n",
    "\n",
    "What if we had plotted the absolute frequency (i.e. the count) of LSOAs in each bin on the y-axis instead of the frequency density: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.distplot(pdf.Owned, bins=Quintiles, kde=False) # Notice no norm_hist!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Think about why all the bars are roughly the same height... it's because by definition each quintile should contain roughly the same number of LSAOs (i.e. 20% in each bin). \n",
    "\n",
    "It's difficult to see where one bar stops and another begins, so let's add some lines to help:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What is going on here? Add some notes\n",
    "sns.distplot(pdf.Owned, bins=Quintiles, kde=False)\n",
    "for b in Quintiles:\n",
    "    plt.vlines(b, 0, 1000, color='red', linestyle='--') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 3\n",
    "\n",
    "1. Okay so, we've now throughly examined the _numerical_ distribution of data for the _Owned_ variable, but **where** are the LSOAs in each of those five bins? What is the **spatial** distribution? Where are the LSOAs in the lowest quintile (i.e. LSOAs with _Owned_ values of 5 to 199) vs those in the highest quintile (i.e. LSOAs with _Owned_ values of 452 to 710)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that `GeoPandas` is using the `PySAL` functionality to do its classification. \n",
    "\n",
    "See how the values in the legend match those in the quintiles we calculated above! And so now we can talk about how the numeric distribution we saw in the histograms above has fallen out spatially. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. In the box below _describe_ the spatial distribution of 'Owned' households across London using the quantile maps; you may find it helpful to go back and modify the output so that you can see it in more detail by making it larger..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Your answer here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But how do we know that equal intervals is the _best_ way to create our classes for visualisation?\n",
    "\n",
    "One of the best-known geographical classification is Fisher-Jenks (also sometimes known as 'Natural Breaks'), which groups data into bins based on the sum of squared deviations between classes: in other words, the algorithm iteratively looks for ways to group the data into a specified number of bins such that moving a data point from one group to another would increase the total within-class deviation observed in the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try the same as above but using (five) Fisher Jenks breaks. Pandas cannot do this so we can only [use PySAL here](https://pysal.readthedocs.io/en/latest/library/esda/mapclassify.html#pysal.esda.mapclassify.Fisher_Jenks) (_**this may take some time depending on your laptop!**_):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "owned_fj = mp.Fisher_Jenks(gdf.Owned.values, k=5) # Classify into Natural Breaks\n",
    "print(owned_fj) # Print summary metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully you can see how the five classes created using this method are different from the quantiles methods we used above. Note how they are unequal in size, for example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a histogram to visualise the numeric distribution of the Fisher Jenks classification. First we need to get the values to create the bins in the histogram. Looking at the example in [the documentation for the PySAL `mapclassify.Fisher_Jenks` function](https://pysal.readthedocs.io/en/latest/library/esda/mapclassify.html#pysal.esda.mapclassify.Fisher_Jenks) we can see that we can access the bins direct from the object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(owned_fj.bins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But note that these are the upper limits of the bin, so we need to add the very lowest value. [We can do this](https://stackoverflow.com/a/36998277/10219907) using the `numpy` `insert` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "owned_fjb = np.insert(owned_fj.bins, 0, gdf.Owned.min())  # insert the minimum value into the zeroth position of the array \n",
    "print(owned_fjb)   # print to check what it looks like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.distplot(pdf.Owned, bins=owned_fjb, norm_hist=True, kde=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's create a map using the Fisher Jenks classification (remember the map may take a while to plot because GeoPandas is using the `fisher_jenks` function from PySAL):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.plot(column='Owned', scheme='fisher_jenks', k=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Compare the map you just made to the one we made using quintiles. Do you think their representation of the spatial distribution is very different?\n",
    "\n",
    "> Your answer here.\n",
    "\n",
    "Let's see if a map created using an `equal_interval` classification looks any different. Combine the information above to write all of the code in one cell to output the classification, histogram and map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Allowing for differences in the colour mapping, your map should look like this:\n",
    "\n",
    "![Equal Interval Map](https://github.com/kingsgeocomp/geocomputation/raw/master/img/Equal_Interval_Map.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This visualisation of the spatial distribution of _Owned_ households using the Equal Interval classification looks quite a bit different to the quantile classification. This is mainly because there are fewer LSOAs classified in the highest (and lowest) classes, meaning that the extreme values are more obvious on the map. You should also be able to see the difference between the shape of the histogram for a Equal Interval classification versus the other classifications.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 4\n",
    "\n",
    "For the next two tasks, you have already produced all the information you need (so no more code needed). Hopefully, answering these tasks will reinforce your understanding about how the three classifications we have tried differ.\n",
    "\n",
    "1. How many LSOAs are in the top class for each of the three classifications:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Your answer here:\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. What is the interval for the top class for each of the three classifications:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Your answer here:\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Create three more maps for the _Owned_ variable using the three classifications (Quantiles, Fisher-Jenks, and Equal Interval), but this time use **seven** classes in each.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Now create three maps of _GreenspaceArea_ using each of the three classification methods. For each of these three maps, also output histograms. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qu = np.insert(mp.Quantiles(gdf.GreenspaceArea.values, k=7).bins, 0, gdf.GreenspaceArea.min())\n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Why are the differences between the three maps for _GreenspaceArea_ more obvious than for the _Owned_ variable we looked at above? \n",
    "\n",
    "> Your answer here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Map Modifications\n",
    "\n",
    "GeoPandas uses much of its plotting functionality from matplotlib, and we can use this to modify and hopefully improve how well our maps communicate. Seaborn also uses matplotlib so hopefully some of the code patterns below will look familiar. \n",
    "\n",
    "First, in the classified choropleth maps above when we included a legend we found that it overlapped with the actual map. Legends are still pretty rudimentary in GeoPandas and moving them is not straight-forward (if you can find an easy way, please share!). So currently the easiest way to ensure the legend does not overlap the map itself is to make the size of our plot larger using the `figsize` option when creating a subplot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(1, figsize=(12, 10))   # create figure and axes for Matplotlib \n",
    "gdf.plot(column='Owned', scheme='quantiles', k=5, legend=True, ax=ax1)  #include ax argument!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we could keep the plot smaller but change the axis limits of the plot so that there is more white space on one side of the map in which the legend can be drawn:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(1, figsize=(8, 6))   # create figure and axes for Matplotlib\n",
    "gdf.plot(column='Owned', scheme='quantiles', k=5, legend=True, ax=ax1) #include ax argument!\n",
    "ax1.set_xlim(500000,586000)   # you can play with these values and see what happens!\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the code above, I just used trial-and-error until I found the right values for `set_xlim` given the size of the plot I wanted. These workarounds aren't the most elegant or efficient but they work for making single maps. \n",
    "\n",
    "Using more functionality from matplotlib we can make our map look much cleaner: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(1, figsize=(15, 12))   \n",
    "gdf.plot(column='Owned', scheme='quantiles', k=5, legend=True, ax=ax1, edgecolor='grey', linewidth=0.2)  #change line style\n",
    "ax1.axis('off') #don't plot the axes (bounding box)\n",
    "ax1.set_title('Ownership', fontdict={'fontsize':'20', 'fontweight':'3'})  #provide a title\n",
    "ax1.annotate('Source: London Datastore (2011)',\n",
    "             xy=(0.1, 0.1), xycoords='figure fraction', \n",
    "             horizontalalignment='left', verticalalignment='top', \n",
    "             fontsize=12, color='#555555')  #add source info on the image itself\n",
    "ax1.get_legend().set_title(\"Households\")  #set the legend title"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the code above check you can see where we:\n",
    "1. change the colour and width of LSOA outlines so they are easier to see\n",
    "2. set the figure title to communicate the variable being mapped\n",
    "3. set the legend title to communicate the units of the data being mapped\n",
    "4. add an annotation show the source of the data for the map."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's another trick to modify the legend labels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(1, figsize=(15, 12))   \n",
    "ax1 = gdf.plot(column='Owned', scheme='quantiles', k=5, legend=True, ax=ax1, edgecolor='grey', linewidth=0.2)  #change line style\n",
    "ax1.axis('off')  #don't plot the axes (bounding box)\n",
    "ax1.set_title('Ownership', fontdict={'fontsize': '20', 'fontweight' : '3'})  #provide a title\n",
    "ax1.annotate('Source: London Datastore (2011)',xy=(0.1, 0.1),  xycoords='figure fraction', horizontalalignment='left', verticalalignment='top', fontsize=12, color='#555555')  #add source info on the image itself\n",
    "\n",
    "#edit legend labels\n",
    "leg = ax1.get_legend()\n",
    "leg.get_texts()[0].set_text('5 - 199')\n",
    "leg.get_texts()[1].set_text('199 - 291')\n",
    "leg.get_texts()[2].set_text('291 - 366')\n",
    "leg.get_texts()[3].set_text('366 - 452')\n",
    "leg.get_texts()[4].set_text('452 - 710')\n",
    "leg.set_title(\"Households\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 5\n",
    "\n",
    "Have a think about how you could change the hard-coded legend labels into a loop that dynamically formats the legend as integers/floats with no decimal places:\n",
    "\n",
    "```python\n",
    "leg = ax1.get_legend()\n",
    "leg.get_texts()[0].set_text('5 - 199')\n",
    "leg.get_texts()[1].set_text('199 - 291')\n",
    "leg.get_texts()[2].set_text('291 - 366')\n",
    "leg.get_texts()[3].set_text('366 - 452')\n",
    "leg.get_texts()[4].set_text('452 - 710')\n",
    "```\n",
    "As a hint, part of your code should look like this:\n",
    "```python\n",
    "    bounds = [float(x) for x in i.get_text().split(' - ')]\n",
    "    b_txt  = \"{0:.0f} - {1:.0f}\".format(bounds[0], bounds[1])\n",
    "```\n",
    "If you don't understand what this is doing, I'd suggest adding comments and printing out variables as you go..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above we saw how we can change the axis limits to zoom out to make room for the legend. But we can also use this to zoom in on our map to a specific location we want to show in more detail:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(1, figsize=(8, 6))   \n",
    "gdf.plot(column='Owned', scheme='quantiles', k=5, legend=True, ax=ax1, edgecolor='grey', linewidth=0.2) \n",
    "ax1.set_xlim(535000,550000)   #play with these values\n",
    "ax1.set_ylim(177000,190000)   #play with these values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Saving maps!\n",
    "\n",
    "There are a few  ways you can get your maps out of a notebook and into a Word document (_e.g._ to use in your reports). One might be to right-click (or Cmd-click on Mac), select 'copy image' and then paste directly into Word; however, this doesn't always work very well, so a better way is to save your map directly to a file and then insert it in Word. \n",
    "\n",
    "You could save your map to a file by (again) right-clicking (or Cmd-clicking on Mac) and select 'Save Image As'. Then you could insert the saved file into your Word document (by going to Insert -> Pictures). \n",
    "\n",
    "However, probably the best way to save a map to file is to use code, and note how this is the _same_ code you'd use for a Seaborn figure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(1, figsize=(8, 6))   \n",
    "gdf.plot(column='Owned', scheme='quantiles', k=5, legend=True, ax=ax1, edgecolor='grey', linewidth=0.2) \n",
    "ax1.set_xlim(535000,550000)   \n",
    "ax1.set_ylim(177000,190000)   \n",
    "\n",
    "# And here is the key line of code\n",
    "# Note that you can also save to a PDF \n",
    "# though there's no guarantee that the\n",
    "# output will be vector...\n",
    "fig.savefig('map_export.png', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Go and check your working directory (usually where you have saved this notebook) to check the image was created!\n",
    "\n",
    "You can [read more about the `savefig` method](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.savefig.html) to see more options available, but `dpi` is the argument that will change the size of your output image. Remember, files will be saved to whatever your current working directory is, but you can include a path if you want to save elsewhere. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Map Layers \n",
    "\n",
    "So far we have looked at only a single variable on a map, but it is possible to present data for multiple variables by layering them on top of one another. For Geopackage files Geopandas can read these directly from a URL (**very cool**) with no need to download and unzip the data first. Since a Shapefile is actually made up of at _least_ four separate files, this one needs to _first_ be saved locally and unzipped in before it can be loaded. So for obvious reasons we rather like the new-ish Geopackage format.\n",
    "\n",
    "Here's code to do the download the data and write it to a directory (creating it if it does not exist):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.makedirs('data', exist_ok=True)\n",
    "\n",
    "# Load Water GeoPackage\n",
    "w_path = os.path.join('data','Water.gpkg')\n",
    "if not os.path.exists(w_path):\n",
    "    water = gpd.read_file('https://github.com/kingsgeocomp/geocomputation/raw/master/data/src/Water.gpkg')\n",
    "    water.to_file(w_path, driver='GPKG')\n",
    "    print(\"Downloaded Water.gpkg file.\")\n",
    "else:\n",
    "    water = gpd.read_file(w_path)\n",
    "    print(\"Loaded Water.gpkg file.\")\n",
    "\n",
    "# Boroughs GeoPackage\n",
    "b_path = os.path.join('data','Boroughs.gpkg')\n",
    "if not os.path.exists(b_path):\n",
    "    boroughs = gpd.read_file('https://github.com/kingsgeocomp/geocomputation/raw/master/data/src/Boroughs.gpkg')\n",
    "    boroughs.to_file(b_path, driver='GPKG')\n",
    "    print(\"Downloaded Boroughs.gpkg file.\")\n",
    "else:\n",
    "    boroughs = gpd.read_file(b_path)\n",
    "    print(\"Loaded Boroughs.gpkg file.\")\n",
    "    \n",
    "# Greenspace GeoPackage\n",
    "g_path = os.path.join('data','Greenspace.gpkg')\n",
    "if not os.path.exists(g_path):\n",
    "    greenspace = gpd.read_file('https://github.com/kingsgeocomp/geocomputation/raw/master/data/src/Greenspace.gpkg')\n",
    "    greenspace.to_file(g_path, driver='GPKG')\n",
    "    print(\"Downloaded Greenspace.gpkg file.\")\n",
    "else:\n",
    "    greenspace = gpd.read_file(g_path)\n",
    "    print(\"Loaded Greenspace.gpkg file.\")\n",
    "\n",
    "print(\"Done.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See how this function allows us to download data to a particular directory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "greenspace.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the data read successfully, you should be able to see the top three lines of the data, with a _geometry_ column on the far right. \n",
    "\n",
    "Let's see if we can plot this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "greenspace.plot(color='green')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully you got something that looked a bit like a map (well, maybe just a bunch of green blobs...). If so, then we have loaded our data properly as a GeoPandas DataFrame. \n",
    "\n",
    "So, let's try to plot our greenspace (parks) data on top of the Owned data (we'll also zoom in so we can see the result in a little more detail):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(1, figsize=(12, 10))  #setup the figure\n",
    "\n",
    "base = gdf.plot(ax=ax1, column='Owned', scheme='quantiles', k=5, cmap='OrRd')         #plot the Owned data layer on ax1\n",
    "parks = greenspace.plot(ax=ax1, color=(0.1, 0.5, 0.2), alpha=0.75, edgecolor='green') #plot the parks data later on ax1\n",
    "\n",
    "ax1.set_xlim(535000,550000)   #zoom\n",
    "ax1.set_ylim(177000,190000)   #zoom\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The key here is that both plots use `ax=ax1` ensuring that they are plotted on the same set of axes. The first plot is on the bottom, then additional plots are layered on top. \n",
    "\n",
    "Try plotting but removing `ax=ax1` from one of the plot calls to see what happens when we don't plot on the same axes: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(1, figsize=(12, 10))  #setup the figure\n",
    "\n",
    "base = gdf.plot(ax=ax1, column='Owned', scheme='quantiles', k=5, cmap='OrRd')       #plot the Owned data layer on ax1\n",
    "parks = greenspace.plot(color=(0.1, 0.5, 0.2), alpha=0.75, edgecolor='green')         #plot the parks data layer but NOT on ax1\n",
    "\n",
    "ax1.set_xlim(535000,550000)   #zoom\n",
    "ax1.set_ylim(177000,190000)   #zoom\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can create a map template that we could use across a number of outputs in order to ensure that each one is fairly consistent:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plt_ldn(w=water, b=boroughs):\n",
    "    fig, ax = plt.subplots(1, figsize=(14, 12))\n",
    "    w.plot(ax=ax, color='#79aef5', zorder=2)\n",
    "    b.plot(ax=ax, edgecolor='#cc2d2d', facecolor='None', zorder=3)\n",
    "    ax.set_xlim([502000,563000])\n",
    "    ax.set_ylim([155000,201500])\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['bottom'].set_visible(False)\n",
    "    ax.spines['left'].set_visible(False)\n",
    "    return fig, ax\n",
    "\n",
    "fig, ax = plt_ldn()\n",
    "gdf.plot(column='Owned', scheme='quantiles', k=5, cmap='YlGnBu', legend=True, alpha=0.9, zorder=1, ax=ax) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 6\n",
    "\n",
    "Duplicate and modify the `plt_ldn` function to include greenspace as well, and then see if you can write a _loop_ that iterates over White, Asian, and Black, creating and saving a Fisher Jenks map (with title) of each. You should _also_ implement the tricks we've seen above in terms of adding a title, formatting the legend appropriately (you might even want to implement this as a function), etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merging Spatial and Non-Spatial Data\n",
    "\n",
    "Finally, let's consider the situation where we don't have an existing (single) dataset that contains geometry, but rather we have two files:\n",
    "1. a file containing geometry for spatial elements\n",
    "2. a file containing data for spatial elements but no geometry\n",
    "\n",
    "In this case, as long as we have a common identifier for each row in the two files, we can use the Pandas `merge` function to join them together. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, let's assume the data we have are:\n",
    "1. A Geopackage delineating LSOAS in and around London \n",
    "2. A csv file containing air quality data for each of our London LSOAs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's download the shapefile data (using the `download_gsa_geodata` function we defined above):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lsoas = gpd.read_file('https://github.com/kingsgeocomp/geocomputation/raw/master/data/src/LSOAs.gpkg', low_memory=False)\n",
    "print(lsoas.columns.values)\n",
    "lsoas.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks okay, but note the map - it looks like there might be more LSOAs than in our previous maps? Let's check the dimensions of this DataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(lsoas.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully you can see that there are more rows in this DataFrame than the one we have been using before. \n",
    "\n",
    "Next, let's load the csv data file as a Pandas DataFrame and look at the top few lines:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aq = pd.read_csv(\n",
    "    'https://github.com/kingsgeocomp/geocomputation/blob/master/data/LSOA_AirQuality.csv.gz?raw=true',\n",
    "    compression='gzip', low_memory=False) \n",
    "\n",
    "aq.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And check the dimensions of this DataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(aq.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's certainly more rows in one DataFrame than the other. \n",
    "\n",
    "But look at the column names of our two data files; can you see a common column between them?\n",
    "\n",
    "The _LSOA11CD_ and _lsoa11cd_ columns contain the same codes so we can use these to merge the data into a single DataFrame. We'll do an _inner join_ (see week 7) with the Air Quality data on the left; this is because there are more rows in the shapefile than the Air Quality data, so we only want to retain rows from the shapefile where we have air quality data: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lsoa_aq = pd.merge(aq, lsoas, how=\"inner\", left_on='LSOA11CD', right_on='lsoa11cd')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check what the top of our new DataFrame looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lsoa_aq.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's check the dimensions of our new DataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(lsoa_aq.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we have the same number of rows as the Air Quality dataframe but the number of columns is the total of the two original DataFrames (check you can from the outputs above how we know this). \n",
    "\n",
    "And finally, let's check what _type_ of DataFrame we've created:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(type(lsoa_aq))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've created a Pandas DataFrame. So we need to convert to a GeoPandas DataFrame if we want to make some maps; to do this we need to specify what column contains the geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lsoa_aq = lsoa_aq.set_geometry('geometry')\n",
    "print(type(lsoa_aq))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now let's try to plot it spatially:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lsoa_aq.plot('NO2max')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully that was successful! Does the this map have the same shape as our previous maps?\n",
    "\n",
    "And now we have created this new GeoPandas DataFrame we can start to examine the Air Quality data spatially. For example: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(1, figsize=(15, 12))   \n",
    "ax1 = lsoa_aq.plot(column='NOxmean', \n",
    "                   scheme='quantiles', \n",
    "                   k=7, \n",
    "                   legend=True, \n",
    "                   ax=ax1, \n",
    "                   edgecolor='grey',    \n",
    "                   linewidth=0.2)       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 7\n",
    "\n",
    "1. _Join_ the `lsoa_aq` and `gdf` files so that you can access both demographic and pollution data. I've given you code to get started but you will also need to to do some careful reading of the error messages. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You only want to run this _once_ since it manually \n",
    "# sets the index so that we can use the _join_ syntax\n",
    "lsoa_aq.set_index('LSOA11CD', inplace=True)\n",
    "gdf.set_index('LSOA11CD', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Now create a filtered data frame (let's call it `fdf`) containing only LSOAs with more than 2,000 usual residents, and check that it worked by creating a map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Now map mean PM2.5 using five equal intervals for the filtered data frame. Use a grey background for the rest of London (you might want to check out: [this file](https://github.com/kingsgeocomp/geocomputation/raw/master/data/src/London.gpkg)) and add the roads ([from this file](https://github.com/kingsgeocomp/geocomputation/raw/master/data/src/Roads.gpkg)) in red. You'll need to do some layering for this so I'd suggest copying _some_ of the code out of the `plt_ldn` and modifying it accordingly, but it's probably not worth a whole new function. You will also need to work out how to specify a light grey, since the default is quite dark... suggest googling matplotlib colors.\n",
    "\n",
    "The final map should look something like the one below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![ownership](https://github.com/kingsgeocomp/geocomputation/raw/master/img/PM25_in_Populated_LSOAs.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Create a map of _Mean PM10_ for central London with roads layered on top. Classify the air quality data into five quantiles and zoom in towards central London. Your final map should look something like the map below. _Note_ that you will need to investigate how to do both mathematical symbols and superscripts in matplotlib/python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![MeanPM10 in Central London](https://github.com/kingsgeocomp/geocomputation/raw/master/img/PM10central.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. How well does the (incomplete) road map align with air quality would you say?\n",
    "\n",
    "> Your answer here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### More Fun!\n",
    "\n",
    "You can find some more nice examples and applications from [GeoHackWeek](https://geohackweek.github.io/vector/04-geopandas-intro/)!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting More Help/Applications\n",
    "\n",
    "A great resource for more help and more examples is Dani Arribas-Bel's _Geographic Data Science_ module: he has put all of his [module practicals online](http://darribas.org/gds17/) (as we have too), and you might find that something that he does makes more sense to you than what we've done... check it out!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Credits!\n",
    "\n",
    "#### Contributors:\n",
    "The following individuals have contributed to these teaching materials: Jon Reades (jonathan.reades@kcl.ac.uk), James Millington (james.millington@kcl.ac.uk)\n",
    "\n",
    "#### License\n",
    "These teaching materials are licensed under a mix of [The MIT License](https://opensource.org/licenses/mit-license.php) and the [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 license](https://creativecommons.org/licenses/by-nc-sa/4.0/).\n",
    "\n",
    "#### Acknowledgements:\n",
    "Supported by the [Royal Geographical Society](https://www.rgs.org/HomePage.htm) (with the Institute of British Geographers) with a Ray Y Gildea Jr Award.\n",
    "\n",
    "#### Potential Dependencies:\n",
    "This notebook may depend on the following libraries: pandas, matplotlib, seaborn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Geocomp 2019",
   "language": "python",
   "name": "gsa2019"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}