Practical-2-Mapping and Classification with Geopandas and PySAL.ipynb

{
 "cells": [
  {
   "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 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/user.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": [
    "## A (Semi-Brief) Discourse on Families & Inheritance\n",
    "\n",
    "GeoPandas objects are deliberately designed to resemble Pandas objects. There are two good reasons for this: \n",
    "\n",
    "1. Since Pandas is well-known, this makes it easier to learn how to use GeoPandas.\n",
    "2. GeoPandas _inherits_ functionality from Pandas. \n",
    "\n",
    "The concept of **inheritance** is something we've held off from mentioning until now, but it's definitely worth understanding if you are serious about learning how to code. In effect, geopandas '_imports_' pandas and then _extends_ it so that the more basic class (pandas in this case) learns how to work with geodata... pandas doesn't know how to read shapefiles or make maps, but geopandas does. Same for GeoJSON.\n",
    "\n",
    "### The 'Tree of Life'\n",
    "\n",
    "Here's a simple way to think about inheritance: think of the 'evolutionary trees' you might have seen charting the evolution of organisms over time. At the bottom of the tree is the single-celled animal, and at the other end are humans, whales, wildebeest, etc. We all _inherit_ some basic functionality from that original, simple cell. In between us and that primitive, however, are a whole series of branches: different bits of the tree evolved in different directions and developed different 'functionality'. Some of us have bones. Some have  cartilege. Some are vegetarian, and some are carnivorous. And so on. When you get to the primates we all share certain common 'features' (binocular vision, grasping hands, etc.), but we are _still_ more similar to gorillas than we are to macaques. So gorillas and humans _extend_ the primitive 'primate functionality' with some bonus features (bigger brains, greater strength, etc.) that are useful, while macaques extend it with a slightly different set of features (tails, etc.).\n",
    "\n",
    "![Tree of Life](http://palaeos.com/systematics/tree/images/treeolif.jpg)\n",
    "\n",
    "### The 'Tree of Classes'\n",
    "\n",
    "Inheritance in code works in a similar way: *all* Python classes (lists, pandas, plots, etc.) inherit their most basic functionality from a single primitive 'object' class that itself does very little except to provide a template for what an object should look like. As you move along the inheritance tree you will find more and more complex objects with increasingly advanced features: GeoPandas inherits from Pandas, Bokeh and Seaborn inherit from matplotlib, etc. \n",
    "\n",
    "I can't find an image of Python base class inheritance, but I've found an equally useful example of how _anything_ can be modelled using this 'family tree' approach... consider the following:\n",
    "\n",
    "![Vehicle Inheritance](http://www.mkonar.org/dogus/wiki/lib/exe/fetch.php/python/vehicle.png?w=750&tok=4c7ed7)\n",
    "\n",
    "If we were trying to implement a vehicle registration scheme in Python, we would want to start with the most basic category of all: _vehicle_. The vehicle class itself might not do much, but it gives us a template for _all_ vehicles (e.g. it must be registered, it must have a unique license number, etc.). We then _extend_ the functionality of this 'base class' with three intermediate classes: two-wheeled vehicles, cars, and trucks. These, in turn, lead to eight actual vehicle types. These might have _additional_ functionality: a bus might need have a passenger capacity associated with it, while a covnertible might need to be hard- or soft-top. All of this could be expressed in Python as:\n",
    "\n",
    "```python\n",
    "class vehicle(object): # Inherit from base class\n",
    "    def __init__(self):\n",
    "        ... do something ...\n",
    "\n",
    "class car(vehicle): # Inherit from vehicle\n",
    "    def __init__(self):\n",
    "        ... do other stuff ...\n",
    "\n",
    "class sedan(car): # Inherit from car\n",
    "    def __init__(self):\n",
    "        ... do more stuff ...\n",
    "```\n",
    "\n",
    "This way, when we create a new `sedan`, it automatically 'knows' about vehicles and cars, and can make use of functions like `set_unique_id(<identification>)` even if that function is _only_ specified in the base vehicle class! The thing to remember is that programmers are _lazy_: if they can avoid reinventing the wheel, they will. Object-Oriented Programming using inheritance is a good example of _constructive_ laziness: it saves us having to constantly copy and paste code (for registering a new vehicle or reading in a CSV file) from one class to the next since we can just import it and extent it! \n",
    "\n",
    "Inheritance also means that you can always use an instance of a 'more evolved' class in place of one of its ancestors: simplifying things a _bit_, a sedan can automatically do anything that a car can do and, by extension, anything that a vehicle can do. Similarly, since geopandas inherits from pandas if you need to use a geopandas object _as if_ it's a pandas object then that will work! So everything you learned last term for pandas can still be used in geopandas. Kind of cool, right?\n",
    "\n",
    "Finally, looking back at our example above: what about unicycles? Or tracked vehicles like a tank? This is where _design_ comes into the picture: when we're planning out a family tree for our work we need to be careful about what goes where. And there isn't always a single right answer: perhaps we should distinguish between pedal-powered and motor-powered (in which case unicycles, bicycles and tricycles all belong in the same family)? Or perhaps we need to distinguish between wheeled and tracked (in which case we're missing a pair of classes [wheeled, tracked] between 'vehicle' and 'two-wheel, car, truck')? These choices are tremendously important but often very hard to get right.\n",
    "\n",
    "OK, that's enough programming theory, let's see this in action..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting the Data\n",
    "\n",
    "Two of the data sets that we want this week are downloadable for free:\n",
    "\n",
    "1. Grab the 'simple' world dhapefile from [ThematicMapping.org'](http://thematicmapping.org/downloads/world_borders.php) and save it into a directory called 'shapes' in your notebook directory _before_ unzipping it.\n",
    "\n",
    "2. The London AirBnB _listings_ from [InsideAirBnB](http://insideairbnb.com/get-the-data.html) which you will save into a directory called 'data' in your notebook directory and _not_ unzip.\n",
    "\n",
    "3. The third data set is one that you can download from [here](http://www.reades.com/LDN-LSOAs.zip)\n",
    "\n",
    "That _should_ cover it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading a Shapefile with GeoPandas\n",
    "\n",
    "Similar to what we did last week, I want you to start by creating a map of the world coloured according to population. But this week we're going to do it in GeoPandas instead of Folium.\n",
    "\n",
    "First, you will need to download the zipped World Borders shapefile from [Thematic Mapping](http://thematicmapping.org/downloads/TM_WORLD_BORDERS_SIMPL-0.3.zip) and save this to your 'shapes' or 'data' folder since GeoPandas can't read 'remote' files. You will also need to unzip the file so that Geopandas can read it.\n",
    "\n",
    "_Remember:_ as long as you have saved the file to somewhere accessible to Jupyter (i.e. a directory _beneath_ the current one where you are working) then you can find out the 'path' to a file by switching to the 'Home' tab in your browser and navigating to the file. You can _also_ find out where the current Notebook is running like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!pwd # !cd on Windows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In other words, depending on your setup you might need to tweak the path in the code below to get your own code working..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "\n",
    "# Make numeric display a bit easier\n",
    "pd.set_option('display.float_format', lambda x: '{:,.0f}'.format(x))\n",
    "\n",
    "# Make sure output is into notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gdf = gpd.read_file( \n",
    "    os.path.join('shapes', 'TM_WORLD_BORDERS_SIMPL-0.3', 'TM_WORLD_BORDERS_SIMPL-0.3.shp') \n",
    ")\n",
    "\n",
    "gdf.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far so good! Notice that even the output of the `head` function is the same as from pandas. The only obvious difference so far is that `geometry` column on the right-hand side which doesn't contain numbers or strings, but `POLYGON`s. In a nutshell, aside from the ability to read geodata directly (we just loaded a shapefile using only code!) this geometry column is the other marker of a big step from pandas to geopandas.\n",
    "\n",
    "But let's dig a little deeper:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"Let's look at inheritance: \")\n",
    "\n",
    "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!\")\n",
    "\n",
    "print(\"-\" * 50)\n",
    "\n",
    "print(\"Let's look at some pandas methods we know: \") # See inheritance in action\n",
    "\n",
    "print(\"Using the 'sample()' method: \")\n",
    "print(gdf.sample(3))\n",
    "print(\"-\" * 50)\n",
    "\n",
    "print(\"What the POP2005 column type: \")\n",
    "print('\\tNAME type: ' + str(type(gdf.POP2005)))\n",
    "print(\"-\" * 50)\n",
    "\n",
    "print(\"Describe the column: \")\n",
    "print(gdf.POP2005.describe())\n",
    "print(\"-\" * 50)\n",
    "\n",
    "print(\"Let's look at some new methods from geopandas: \") # Here's stuff pandas can't do\n",
    "\n",
    "print(\"What's the geometry column type: \")\n",
    "print('\\tGeometry type: ' + str(type(gdf.geometry)))\n",
    "print(\"-\" * 50)\n",
    "\n",
    "# The next line is very slow, don't run it:\n",
    "# print(gdf.geometry.describe())\n",
    "# print(\"-\" * 50)\n",
    "\n",
    "print(\"Is there projection information: \")\n",
    "print('\\tCRS: ' + str(gdf.crs))\n",
    "print(\"-\" * 50)\n",
    "\n",
    "print(\"Help with the centroid method of the GeoSeries (geometry) column: \")\n",
    "help(gdf.geometry.centroid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This next command may take a while to complete -- the geopandas developers know that their mapping code is sloooooooooow, and that's one of the reasons we're going to use PySAL for most of our work. But 'watch this space' since they are working on it and there may even be an update this year that addresses the single biggest complaint that people have about geopandas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gdf.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some questions for you:\n",
    "1. What type of data structure is the CRS attribute of the geopandas data frame?\n",
    "2. How would you confirm this?\n",
    "3. Can you figure out how to reproject the data using the geopandas documentation?\n",
    "4. What is the data structure of the geometry column (broadly)?\n",
    "5. Are there any differences between what a similar pandas data frame would give you if you ran any of the non-geographic commands?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with the AirBnB Data\n",
    "\n",
    "That should have given you a sense of some parallels/differences between Folium and Geopandas: Folium is more of a web-based, interactive visualisation library, while Geopandas is really built around large data sets and data analytics. We're now going to move into areas where Folium would have a lot of trouble keeping up: a data set scraped from the AirBnB web site. \n",
    "\n",
    "### Randomness & Reproducibility\n",
    "\n",
    "In fact, it turns out that rapid visualisation of the _full_ AirBnB data set using Geopandas/PySAL is hard: there's simply so much of it that visualisation is slow unless you're in a dedicated environment such as QGIS. So, for the _exploratory_ part of our work we'll want to work with a _sample_ -- but what happens if every time we take a sample we get a _different_ sample? That obviously makes things a bit harder, it would be handy if we could get the _same_ random sample every time _while_ we're doing our testing and development before expanding to the full data set.\n",
    "\n",
    "That's where `random.seed` comes in: by setting a seed we ensure that any process based on a random/random sampling process will be reproducible. In other words, we'll get the _same_ random sample each time. To under why this happens you'd need to read up on pseudo-randomness and computers but that's not really relevant here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import geopandas as gpd\n",
    "\n",
    "import random \n",
    "random.seed(123456789) # For reproducibility\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "df = pd.read_csv(os.path.join('data','airbnb','listings-summary.csv.gz'))\n",
    "\n",
    "print(\"Full data set shape is: \" + ' by '.join(str(i) for i in df.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sample = df.sample(frac=0.1)\n",
    "print(\"Sample data set shape is: \" + ', '.join(str(i) for i in sample.shape))\n",
    "sample.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can _always_ do a 'quick and dirty' scatter plot to see if the data seems vaguely sensible -- it's obviously limited as a geo-visualisation but it can give you an _idea_ of whether or not you've done the right thing with your data. For example..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#help(sns.jointplot)\n",
    "sns.jointplot(x=\"longitude\", y=\"latitude\", data=sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll notice that this 'map' isn't particularly good, but it does tell us that the longitude and latitude values are reasonable: you'd expect to find more AirBnB listings towards the middle of the city and there's a _hint_ of the Thames and the Lee Valley in there (though this is a bit of a stretch). To actually _map_ the data we'll need to be a little more rigorous... \n",
    "\n",
    "Let's step through what's going on below:\n",
    "1. We need to import the `Point` class so that we move from separate x and y columns to a single 'point' that Geopandas can work with.\n",
    "2. We then 'zip' up the x and y (i.e. lat and long) into pairs -- think of this as a simple way to pair _each_ x and y based on their row position and this allows us to move from separate columns to actual points.\n",
    "3. The next step is to tell Geopandas what projection our data is in -- raw lat and long are _usually_ recorded in WGS84 which has the EPSG identifier 4326 (_i.e._ epsg:4326).\n",
    "4. You'll notice that to create a new `GeoDataFrame` we do so _slightly_ differently from how we created a new `DataFrame` last term: we pass in the existing `pandas` data frame (`sample`), the CRS (projection), and finally the `geometry` that we created from the `zip` process.\n",
    "5. The last step is to reproject the geometry into OSGB (Ordnance Survey GB) which has the EPSG identifier 27700. \n",
    "\n",
    "You can see the results of this step in the final step where we print out the first 3 rows of the reprojected data: notice that the point coordinates are no longer in lat/long!\n",
    "\n",
    "Some of this _might_ seem a little tedious, but it's incredibly useful to be able to automate this process: we can reproject a whole series of shapefiles (e.g. every single file in a directory!), we can convert CSV files into shapes that load automatically into QGIS instead of having to do this process manually..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from shapely.geometry import Point\n",
    "# Convert x,y to Points using zip(...)\n",
    "geometry = [Point(xy) for xy in zip(sample.longitude, sample.latitude)]\n",
    "\n",
    "print(\"From the data frame...\")\n",
    "print(sample.head(3)[['longitude','latitude']])\n",
    "print(\"-\" * 50)\n",
    "\n",
    "print(\" \")\n",
    "print(\"From the geometry zip...\")\n",
    "print([\", \".join([str(p.x), str(p.y)]) for p in geometry[0:3]])\n",
    "print(\"-\" * 50)\n",
    "\n",
    "print(\" \")\n",
    "# Initialise to WGS84\n",
    "crs = {'init' :'epsg:4326'}\n",
    "sdf = gpd.GeoDataFrame(sample, crs=crs, geometry=geometry)\n",
    "\n",
    "# Reproject into OSGB\n",
    "sdf = sdf.to_crs({'init' :'epsg:27700'})\n",
    "print(sdf.head(3)[['neighbourhood','geometry']])\n",
    "\n",
    "# And save it as a new shapefile\n",
    "sdf.to_file(os.path.join('shapes','AirBnB-Sample.shp'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with the world map above, plotting the AirBnB sample will take some time... so be patient! You'll notice that the results are now also reported in OSGB units, not lat/long, so this is one way in which GeoPandas is more 'knowledgeable' about geodata than pandas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sdf.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Saving to a Shapefile!\n",
    "\n",
    "_Aaaaaand_ as a final step: let's save our data sample as a brand new shape file that we can open directly in QGIS or ArcGIS. This is going to be rather useful: you could create your entire workflow in an automated, replicable fashion, and then just write the results out to a Shapefile so that you can get your visualisation _just right_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sdf.to_file(os.path.join('data','airbnb','listings-summary.shp'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with PySAL\n",
    "\n",
    "Now that we've prepared the raw data, we're ready to begin a more substantive analysis using PySAL. Of course, at this stage we're still going to take it fairly easy and not jump feet-first into using concepts like spatial autocorrelation or spatially-aware statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import pysal as ps\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "from pylab import *\n",
    "from pysal.contrib import pdio\n",
    "from pysal.contrib.viz import mapping as maps\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_Breaking News!_** Before we get into the guts of PySAL, I wanted to draw your attention to a _brand new_ feature in PySAL 1.13: much tighter integration with/support for pandas. It is now possible to ask PySAL to read a Shapefile or DBF file and return it to you as a pandas dataframe. This expands the range of options for working with geographical data by a significant amount:\n",
    "\n",
    "```python\n",
    "table_path = ps.examples.get_path('south.dbf') # Get a sample data set from PySAL\n",
    "table = ps.pdio.read_files(table_path) # pdio = pandas I/O (input/output)\n",
    "table.head() # It's a pandas dataframe!\n",
    "ps.weights.Queen.from_dataframe(table) # Create weights from dataframe\n",
    "```\n",
    "\n",
    "You can find out more about the _new_ functions available [here](http://nbviewer.jupyter.org/gist/ljwolf/b960c198742eb26c335d9e32bdc7b76c?flush_cache=true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading a Shapefile in PySAL\n",
    "\n",
    "Let's load a shapefile of London LSOAs -- these are the geographical boundaries for the NS-SeC data that you were working with last term. We kept geography out of it since getting to grips with Python is hard work enough, but now we want to begin exploring how adding space to our analysis starts to change _many_ things..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lsoas = ps.open(os.path.join('shapes','LDN-LSOAs.shp'))\n",
    "print(\"LSOAs class is: \" + str(type(lsoas)))\n",
    "print(\"There are \" + str(len(lsoas)) + \" polygons in the LSOA Shapefile.\")\n",
    "print(\"Shapefile headers:\")\n",
    "lsoas.header"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "listings = ps.open(os.path.join('data','airbnb','listings-summary.shp'))\n",
    "print(\"There are \" + str(len(listings)) + \" points in the listings summary Shapefile.\")\n",
    "listings.header # See what we loaded"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Revisiting Classification\n",
    "\n",
    "Like pandas/geopandas, PySAL has a classification 'engine' that we can use to bin data based on attribute values. You might wonder why PySAL has its own process when we can already do this within pandas -- there are two reasons: first, now that we've got a shapefile we don't _really_ need to use pandas/geopandas any more and can do our work solely using PySAL; second, PySAL has additional classification methods that are _specific_ to geographic analysis problems. \n",
    "\n",
    "The best-known of these 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.\n",
    "\n",
    "But first, let's just go with quantiles and see how you can now work only with PySAL:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Notice: listingsdata != listings (DBF != SHP)\n",
    "listingsdata = ps.open(os.path.join('data','airbnb','listings-summary.dbf'))\n",
    "print(\"There are \" + str(len(listingsdata)) + \" rows in the listings summary DBF.\")\n",
    "listingsdata.header # See what we loaded\n",
    "\n",
    "prices = np.array(listingsdata.by_col['price']) # Retrieve pricing data\n",
    "pricesq5 = ps.esda.mapclassify.Quantiles(prices, k=5) # Classify into 5 quantiles\n",
    "print(pricesq5) # Show the classification result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pricesq5.yb # Check the class assignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# dir(maps) # If you want to explore this class\n",
    "# help(maps.map_poly_shp) # And here's how we plot polygons!\n",
    "fig = figure(figsize=(8,6))\n",
    "\n",
    "base = maps.map_poly_shp(lsoas)\n",
    "base.set_facecolor('none')\n",
    "base.set_linewidth(0.5)\n",
    "base.set_edgecolor('0.8')\n",
    "\n",
    "lyr = maps.map_point_shp(listings)\n",
    "lyr = maps.base_choropleth_classif(lyr, pricesq5.yb)\n",
    "lyr.set_alpha(0.5)\n",
    "lyr.set_linewidth(0.)\n",
    "lyr.set_sizes(np.repeat(2.5, len(listingsdata)))\n",
    "\n",
    "ax = maps.setup_ax([base, lyr], [lsoas.bbox, listings.bbox])\n",
    "fig.add_axes(ax)\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's unpick the steps involved in this process:\n",
    "\n",
    "First, we create the figure and specify an output size (this was done in 'inches' [as the screen sees them] but you can specify other units as well such as pixels and cm):\n",
    "```python\n",
    "fig = figure(figsize=(8,6)) # Set the output size in 'inches'\n",
    "```\n",
    "\n",
    "Next, we create a 'base' map that uses the LSOA polygons loaded from the LDN-LSOA Shapefile:\n",
    "```python\n",
    "base = maps.map_poly_shp(lsoas) # Create a polygon map from the LSOA object\n",
    "base.set_facecolor('none')      # Set the colour of the 'face' (i.e. interior of the polygon)\n",
    "base.set_linewidth(0.5)         # Set the width of the edge\n",
    "base.set_edgecolor('0.8')       # Set the colour of the edge (this is a light grey)\n",
    "```\n",
    "\n",
    "Now we'll create another map, this one containing the AirBnB points loaded from the listings shapefile that we created earlier:\n",
    "```python\n",
    "lyr = maps.map_point_shp(listings) # Create a point map from the AirBnB listings object\n",
    "lyr = maps.base_choropleth_classif(lyr, pricesq5.yb) # Turn this into a choropleth \n",
    "lyr.set_alpha(0.5)                 # Set the transparency (a.k.a. alpha) to 50%\n",
    "lyr.set_linewidth(0.)              # Set the line width around each point to 0\n",
    "lyr.set_sizes(np.repeat(2.5, len(listingsdata))) # Set the size of the points to 2.5 (points?)\n",
    "```\n",
    "\n",
    "The final step is to add both 'layers' to our figure (so, technically, we could have moved the `fig` line down to here, but it's easier to figure out what's going on if we leave it up above):\n",
    "```python\n",
    "ax = maps.setup_ax([base, lyr], [lsoas.bbox, listings.bbox]) # Add layer _and_ bounding box\n",
    "fig.add_axes(ax) # Specify the axes\n",
    "show() # Done!\n",
    "\n",
    "```\n",
    "Notice that you have to set up the x and y axes (`setup_ax`) using two parameters: the first array contains the layer objects created using `maps.map_<type>_shp(...)`, and the second contains the bounding boxes (`bbox`) derived as part of that process. **Note:** you *cannot* find the `bbox` on the object created using `maps...`, it is associated with the original object created using PySAL (`ps.open(...)`).\n",
    "\n",
    "Another example, this time using Fisher-Jenks. Notice how the results look very, very different and how the boundaries between bins have been radically altered but also capture a bit better our intuition as to how the data can be represented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "prices = np.array(listingsdata.by_col['price']) # Retrieve pricing data\n",
    "pricesjk = ps.esda.mapclassify.Fisher_Jenks(prices, k=5) # Classify into 5 quantiles\n",
    "print(pricesjk) # Show the classification result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make it easier to see what's happening in Central London why don't we try zooming in a little bit? Notice that in the previous example we passed in the _bounding box_ as part of the axes command. Have a look at how the `listings.bbox` is set up and see if you can't figure out how to zoom in to something like:\n",
    "\n",
    "Dimension | Max/Min | Value \n",
    "--------- | ------- | -----\n",
    "Northing | Max | 190000\n",
    "Northing | Min | 170000\n",
    "Easting | Max | 540000\n",
    "Easting | Min | 520000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(listings.bbox)\n",
    "my_bbox = ???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# To find out more about setting colours:\n",
    "# help(plt.cm)\n",
    "# help(maps.base_choropleth_classif)\n",
    "\n",
    "fig = figure(figsize=(7,7)) # Why do you think I changed the figure size this way?\n",
    "\n",
    "base = maps.map_poly_shp(lsoas)\n",
    "base.set_facecolor('none')\n",
    "base.set_linewidth(0.5)\n",
    "base.set_edgecolor('0.8')\n",
    "\n",
    "lyr = maps.map_point_shp(listings)\n",
    "lyr = maps.base_choropleth_classif(???, cmap='viridis')\n",
    "lyr.set_alpha(0.5)\n",
    "lyr.set_linewidth(0.)\n",
    "lyr.set_sizes(np.repeat(5, len(listingsdata))) # Sice of the dots\n",
    "\n",
    "plt.title(\"Fisher-Jenks Classification of AirBnB Sample\")\n",
    "ax = maps.setup_ax([base, lyr], [???, ???]) # You should get a 'crop' of London\n",
    "fig.add_axes(ax)\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at how the Fisher-Jenks classification relates to the underlying distribution by plotting the prices and bin boundaries derived by PySAL:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.kdeplot(prices, shade=True)\n",
    "for b in pricesjk.bins:\n",
    "    plt.vlines(b, 0, 0.011, color='red', linestyle='--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing Your Understanding\n",
    "\n",
    "Over to you! I'd like you to test [a number of the `mapclassify` options](http://pysal.readthedocs.io/en/latest/library/esda/mapclassify.html) and make up your own mind as to which is the 'best' for classifying the AirBnB data for London. To do this you will need to think about how you can automate the steps we've taken above, including:\n",
    "\n",
    "* Loading the data\n",
    "* Running the classifer\n",
    "* Outputting it to a map\n",
    "* Outputting the distribution with bin boundaries\n",
    "* Looking at 'goodness-of-fit' metrics (see the documentation)\n",
    "\n",
    "To make sure that you've done all of this correctly you should be able to click on 'Kernel' > 'Restart' and then run _only_ the code in the block below fully and successfully. So you'll need to factor in things like loading the libraries that you'll need, telling the notebook where to find the data, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [spats]",
   "language": "python",
   "name": "Python [spats]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}