F0C-13-Maps.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FLIP (00): Data Science \n",
    "**(Module 02: Data Visualization)**\n",
    "\n",
    "---\n",
    "- Materials in this module include resources collected from various open-source online repositories.\n",
    "- You are free to use,but NOT allowed to change and distribute this package.\n",
    "\n",
    "Prepared by and for \n",
    "**Student Members** |\n",
    "2006-2018 [TULIP Lab](http://www.tulip.org.au), Australia\n",
    "\n",
    "---\n",
    "\n",
    "# Session N: `Basemap`\n",
    "\n",
    "The matplotlib library provides hooks for creating extensions called toolkits. To date, several have been created to compliment matplotlib, and over the course of the next few lessons, we'll be using one in particular, called Basemap, for mapping data.\n",
    "\n",
    "The Basemap matplotlib toolkit is a library for plotting 2-dimensional data onto maps. As an extension to matplotlib, Basemap does not do any of the plotting itself, but instead it provides several map projections for transforming geographic coordinates into points on a 2-dimensional cartesian plane, in other words, these projections turn longitude and latitude into x and y coordinates. Basemap was created by scientists for scientists, and as such, it can be a bit cumbersome to work with at times for the layman, but it's extremely powerful and featureful, and is likely to be the only tool you'll ever need for mapping data with matplotlib.\n",
    "\n",
    "In this lesson, we'll see how to use the Basemap library to draw a simple map of the world complete with political boundaries and lines of longitude and latitude. \n",
    "\n",
    "Before we can draw our map, however, we first need to setup our environment and choose a projection method.\n",
    "\n",
    "<!---\n",
    "In this lesson, we'll take our first steps with Basemap. We'll be using the toolkit to draw a simple map of the world and style it to make it aesthetically pleasing. In this lesson\n",
    "\n",
    "In this lesson, we'll take our first steps with Basemap. We'll be using the toolkit to draw a simple map of the world, and along the way, we'll see how to do things like draw the meridians and parallels on a map, otherwise known as the longitude and latitude lines on a map. We'll see how to draw the coastlines and political borders of a map. And, how to style a map to make it aesthetically pleasing as well as informative. But, before we can do any of those things, our very first task at hand is choosing a map projection.\n",
    "--->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# The current version of NumPy available from conda is issuing a warning \n",
    "# message that some behavior will change in the future when used with the \n",
    "# current version of matplotlib available from conda. This cell just keeps\n",
    "# that warning from being displayed.\n",
    "import warnings\n",
    "warnings.simplefilter(action=\"ignore\", category=FutureWarning)\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "from IPython.display import set_matplotlib_formats\n",
    "\n",
    "# Turn on retina display mode\n",
    "set_matplotlib_formats('retina')\n",
    "\n",
    "# Turn off interactive mode\n",
    "plt.ioff()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choosing a Projection\n",
    "\n",
    "As we learned in the previous lesson, the act of transforming points on a 3-dimensional object to coordinates on a 2-dimensional surface is not something that can be done without some error. As a result, every map you've ever used is in some way incorrect, but even with inaccuracies, each has some advantage in particular use cases. The example we gave in the last section was that of the mercator projection which has advantages for nautical navigation at the expense of greatly distorting the size of some land masses. \n",
    "\n",
    "Picking the right map projection is a question of trade-offs. You need to think about your data, and the message you're trying to convey, and choose the right projection for both. However, since this is a course on matplotlib and not mapping, getting into the details of how to choose the right map projection is well outside of the scope of this course, but if you're interested in learning more, I've included a few links to some excellent articles below.\n",
    "\n",
    "> Some good articles on map projections, the problems associated with them, and choosing the right one.\n",
    "- [Choosing the Right Map Projection][1]\n",
    "- [Why Google Maps gets Africa wrong][2]\n",
    "\n",
    "Now, to get a list of the projections available in the Basemap library, just import the `basemap` module and print the `basemap.supported_projections` variable. This will print out a table with the abbreviated names that you'll pass into the `Basemap` constructor, and a brief description of the type of each map projection. The [Basemap documentation][3] also has a list of each projection with a link to a short explanation of each one and an example of how to use it. \n",
    "\n",
    "[1]: https://source.opennews.org/en-US/learning/choosing-right-map-projection/\n",
    "[2]: http://www.theguardian.com/world/2014/apr/02/google-maps-gets-africa-wrong\n",
    "[3]: http://matplotlib.org/basemap/users/mapsetup.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits import basemap\n",
    "print basemap.supported_projections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see from the list above, there are plenty of projections to choose from, but for this example, we'll be using the Robinson projection. \n",
    "\n",
    "### The Robinson Projection\n",
    "\n",
    "The [Robinson projection][2] was created by [Arthur H. Robinson][3] in 1963 for the Rand McNally company. We'll be using it here for a few reasons. First, it's a projection that's meant to show the entire world at once, and since we're planning on showing the entire world, this projection is exactly what we're looking for. \n",
    "\n",
    "Second, the Robinson projection is, in my opinion at least, just a really good looking map. The approach that Dr. Robinson took when designing his projection, according to a 1988 [New York Times article][4], was backwards from most mapmakers. Rather than starting with the math, he approached the problem more artistically, visualizing the best-looking shapes and sizes first, and then figuring out the mathematical formula later to produce the effect he wanted. The result was a map projection that was so good it became the general purpose map for the National Geographic Society for nearly all of the nineties. \n",
    "\n",
    "Finally, every projection has a set of required parameters that you must specify when creating the projection, for the Robinson projection that set is minimal. In fact, only one parameter is needed to create the Robinson projection and that's the longitude of the center of the map. We'll be using the [prime meridian][5] as our central point, so the only parameter we'll need to pass in is a longitude of $0^{\\circ}$.\n",
    "\n",
    "Creating our projection is easy. Simply call the `Basemap` constructor and pass in the name of the projection that you want to use---in our case that will be `'robin'`. Then we just specify the longitude of the center of the map by passing it in as the value of the `lon_0` parameter. \n",
    "\n",
    "<!---\n",
    "Finally, we're going to specify the resolution of the map. This setting determines how detailed the boundaries of the map (e.g., the coastlines) will be. The options available for the resolution are `'c'` (crude), `'l'` (low), `'i'` (intermediate), `'h'` (high), `'f'` (full) or `None`. The higher the resolution the longer it takes to draw the map. The default value for this is 'c', which we will go ahead and specify in our example. If None is chosen for this parameter, no boundary data will be read in and any calls to methods such as `drawcoastlines` will raise an error.\n",
    "--->\n",
    "\n",
    "[1]: https://en.wikipedia.org/wiki/Robinson_projection\n",
    "[2]: http://matplotlib.org/basemap/users/robin.html\n",
    "[3]: https://en.wikipedia.org/wiki/Arthur_H._Robinson\n",
    "[4]: http://www.nytimes.com/1988/10/25/science/the-impossible-quest-for-the-perfect-map.html\n",
    "[5]: https://en.wikipedia.org/wiki/Prime_meridian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create a Robinson projection that's centered around \n",
    "# longitude 0 (the prime meridian)\n",
    "m = Basemap(projection='robin', lon_0=0.)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a map projection object, let's go ahead and set up our map. The lines below will create a new figure, set the background color, and remove the outline from the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(20,10))\n",
    "\n",
    "# Fill the background of the map with a greyish color\n",
    "fig.patch.set_facecolor('#e6e8ec')\n",
    "\n",
    "# Draw the white map background and remove the boundary line\n",
    "m.drawmapboundary(color='none', fill_color='white');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to draw the map itself. The lines below will draw the coastlines, fill all of the continents, and draw the political boundaries between countries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Draw the continents and the countries\n",
    "\n",
    "# This line isn't necessarily needed, but it does add some extra definition \n",
    "# to some geographic elements (the Amazon river, for example).\n",
    "m.drawcoastlines(color='white')\n",
    "m.fillcontinents(color='black', lake_color='white')\n",
    "m.drawcountries(linewidth=1, color='white');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we'll draw the graticule---the network of lines representing the meridians and parallels, or longitude and latitude lines respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Draw the meridians (longitude) and parallels (latitude)\n",
    "\n",
    "# The space between each line in the graph\n",
    "graticule_width = 20\n",
    "# Create lists of meridians and parallels to draw on map\n",
    "parallels = np.arange(-90, 91, graticule_width)\n",
    "meridians = np.arange(-180., 181., graticule_width)\n",
    "# The pattern for the dashes in the graticule lines. The default is\n",
    "# [1,1], but we're going with 3 pixels on, 5 pixels off.\n",
    "dashes = [3, 5] \n",
    "# The color of the graticule lines\n",
    "graticule_color = 'grey'\n",
    "\n",
    "m.drawparallels(parallels, \n",
    "                dashes=dashes, \n",
    "                color=graticule_color)\n",
    "\n",
    "m.drawmeridians(meridians, \n",
    "                dashes=dashes, \n",
    "                color=graticule_color)\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "And, with the map above complete, we've come to the end of this part. In this lesson, we created our first map using the Basemap library. We learned how to create a map projection and how to use it to draw a map and style it slightly. \n",
    "\n",
    "In the coming lessons we'll build on these skills to create a couple of interesting data visualizations with maps. \n",
    "\n",
    "## References\n",
    "\n",
    "This notebook shows how to work with maps in Python using matplotlib and basemap. It is based on two articles: [\"Working With Maps in Python\"][1] by Max Berggren which is itself based on the second article [\"Working with shapefiles, projections and world maps in ggplot\"][2] by Kristoffer Magnusson. The latter uses R and ggplot and the former is a translation of the latter into Python using matplotlib and the [basemap][3] library.\n",
    "\n",
    "[1]: http://maxberggren.github.io/2015/08/04/basemap/\n",
    "[2]: http://rpsychologist.com/working-with-shapefiles-projections-and-world-maps-in-ggplot\n",
    "[3]: http://matplotlib.org/basemap/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choropleth\n",
    "\n",
    "In the previous part, we learned about the [Basemap][1] mapping toolkit for the matplotlib library. We also played around with it a bit to create a simple, aesthetically pleasing map of the world, but aside from being able to see the general location of each country in the world, our map wasn't necessarily packed with information. \n",
    "\n",
    "In this lesson, we'll use the basemap toolkit to create a chorlopleth map. \n",
    "\n",
    "A [choropleth][2] is a [thematic map][3] in which areas are shaded, patterned, or colored to represent some underlying value associated with a particular phenomenon. The specific phenomenon we'll be studying is unemployment. Specifically, we'll be looking at the unemployment rate per county in the US before, during, and after the [Great Recession][4].\n",
    "\n",
    "Let's not waste anymore time talking about the task at hand, let's dive right in and get set up.\n",
    "\n",
    "[1]: http://matplotlib.org/basemap/\n",
    "[2]: https://en.wikipedia.org/wiki/Choropleth_map\n",
    "[3]: https://en.wikipedia.org/wiki/Thematic_map\n",
    "[4]: https://en.wikipedia.org/wiki/Great_Recession_in_the_United_States"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from IPython.display import set_matplotlib_formats\n",
    "set_matplotlib_formats('retina')\n",
    "\n",
    "# The current version of NumPy available from conda is issuing a warning \n",
    "# message that some behavior will change in the future when used with the \n",
    "# current version of matplotlib available from conda. This cell just keeps\n",
    "# that warning from being displayed.\n",
    "import warnings\n",
    "warnings.simplefilter(action=\"ignore\", category=FutureWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the Unemployment Data\n",
    "\n",
    "Just as we did in our earlier lessons, in this lesson, we'll be making use of `namedtuple`s to wrap each unemployment record in the data making it easier to work with throughout the rest of the lesson. \n",
    "\n",
    "The code below creates the `UnemploymentRecord` data structure with the `namedtuple` function. Of the fields passed into the `namedtuple` function, the two that we are the most interested in, are the FIPS code and the unemployment rate. A [FIPS code][1] is a five-digit number that uniquely identifies each county in the US, and we'll be using it to match up the unemployment rates to the county objects in the shapefile that we'll read in later. This will allow us to color each county in the map according to its unemployment rate.\n",
    "\n",
    "[1]: https://en.wikipedia.org/wiki/FIPS_county_code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "\n",
    "UnemploymentRecord = namedtuple('UnemploymentRecord', ('laus_code', \n",
    "                                                       'fips_code', \n",
    "                                                       'name', \n",
    "                                                       'year',\n",
    "                                                       'labor_force',\n",
    "                                                       'employed',\n",
    "                                                       'unemployment_level',\n",
    "                                                       'unemployment_rate'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a data structure for each unemployment record, we're ready to read in the data. Since we're going to be loading data from multiple years, it makes sense to create a function for doing this task, so we don't end up copying and pasting the code into multiple places. \n",
    "\n",
    "The code below creates a function called `load_unemployment_rates` that takes one parameter, a year as an integer value, downloads the data for the given year, and returns a list of `UnemploymentRecord` objects. \n",
    "\n",
    "The data that we'll be using comes from the [Bureau of Labor Statistics][1]. Specifically, we'll be using the [county unemployment data][2] in the Excel format. The data's also available in a text format, but it isn't in a machine-friendly format such as a comma, or tab delimited format, but thankfully, the fantastic [xlrd][3] library is available as a conda install, and it makes it very convenient to read in the data from an Excel file. Given this requirement, before you can run the code in the `load_unemployment_rates` function below, you'll need to do a quick `conda install xlrd` to install the required package.\n",
    "\n",
    "[1]: http://www.bls.gov/\n",
    "[2]: http://www.bls.gov/lau/#cntyaa\n",
    "[3]: http://www.python-excel.org/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import urllib2\n",
    "import xlrd\n",
    "\n",
    "# The URL for the Labor Force Data by County, or Local Area \n",
    "# Unemployment Statistics (LAUS), from the Bureau of Labor Statistics\n",
    "#base_url = \"http://www.bls.gov/lau/laucnty%02d.xlsx\"\n",
    "\n",
    "# The following code pulls the local data files in instead of downloading them on the fly\n",
    "import os\n",
    "base_url = \"http://www.bls.gov/lau/laucnty%02d.xlsx\" # % os.getcwd()\n",
    "\n",
    "def load_unemployment_rates(year):\n",
    "    \"\"\"Loads the unemployment data for the given year from the BLS.\n",
    "    \n",
    "    Arguments:\n",
    "    year -- The year (as an int) to load unemployment data for\n",
    "    \"\"\"\n",
    "    url = base_url % (year - 2000)\n",
    "    response = urllib2.urlopen(url)\n",
    "    workbook = xlrd.open_workbook(file_contents=response.read())\n",
    "    sheet = workbook.sheet_by_index(0)\n",
    "    \n",
    "    data = []\n",
    "    for i in range(6, sheet.nrows-4):\n",
    "        row = sheet.row_values(i)\n",
    "        data.append(UnemploymentRecord(\n",
    "            laus_code           = row[0],\n",
    "            fips_code           = row[1] + row[2],\n",
    "            name                = row[3],\n",
    "            year                = int(row[4]),\n",
    "            # For some reason, they have an extra, empty column here\n",
    "            labor_force         = int(row[6]) if row[6] != 'N.A.' else 0,\n",
    "            employed            = int(row[7]) if row[7] != 'N.A.' else 0, \n",
    "            unemployment_level  = int(row[8]) if row[8] != 'N.A.' else 0,\n",
    "            unemployment_rate   = float(row[9]) if row[9] != 'N.A.' else 0.))\n",
    "        \n",
    "    return data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code for the `load_unemployment_rates` function is really quite straightforward. We simply download the data for the given year, parse it with the xlrd library, and wrap each row's data in an `UnemploymentRecord` tuple. The only bit of preprocessing that we need to do on the data can be seen on lines 24 and 26-31. On line 24, we concatenate the state and county FIPS codes into a single FIPS code string, and on lines 26-31 we cast some of the data, e.g., the unemployment rate, into `int`s and `float`s. We have to be a bit careful when doing this though since, in some records, we can get an `'N.A.'` for the value. In those cases, we simply set the value to 0. \n",
    "\n",
    "## Create a Color Mapper\n",
    "\n",
    "Alright, so now we have a function for loading in the unemployment data, next we'll need a way to map each unemployment rate to a color. Personally, I prefer as simple an interface as possible, e.g., I'd like to just call a function with a value and get back a color. Unfortunately, we just can't grab a color map from matplotlib and pass in an unemployment rate. The problem with doing so is that a color map expects a float between 0 and 1 and our data ranges from roughly 2-30, thus every value in our data falls outside of the range of the color map, and so we end up with a map where every county is \"painted\" black. \n",
    "\n",
    "A color map also supports passing in an integer value, and it maps that value to one of 256 colors in the map. Given this option, we could simply call `int` on each value before passing it into the color map, but unfortunately, this creates another problem. Since our data consists of values between 2 and 30, all of our values are on the low end of the spectrum when compared to the color map's range of 0-255. With this method, we would end up with an extremely \"watered down\" map where all of the counties look a nearly identical shade of pale yellow.\n",
    "\n",
    "The solution to both of these problems is to normalize each unemployment rate before passing it into our color map. Fortunately, matplotlib provides the `pyplot.Normalize` class to make doing so extremely easy. But, as I mentioned before, we want a simple interface that allows us to pass in a value and get a color, and we don't want to normalize each value before doing so. One way to accomplish this is through a callable object that encapsulates the creation of the normalization function for our data, as well as its application on each individual value. Again matplotlib provides what we need to make this really easy. The `pyplot.cm.ScalarMappable` class is a mixin class for adding color map functionality to custom classes. In the code snippet below, we use this class to create our `HeatMapper` helper class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class HeatMapper(plt.cm.ScalarMappable):\n",
    "    \"\"\"A callable that maps cold colors to low values, and hot to high.\n",
    "    \"\"\"\n",
    "    def __init__(self, data=None):\n",
    "        norm = plt.Normalize(vmin=min(data), vmax=max(data))\n",
    "        cmap = plt.cm.hot_r\n",
    "        super(HeatMapper, self).__init__(norm, cmap)\n",
    "        \n",
    "    def __call__(self, value):\n",
    "        return self.to_rgba(value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our `HeatMapper` class is very simple since it inherits from the `plt.cm.ScalarMappable` class that takes on the brunt of the hard work for us. All we have to do is create the normalize function from the data that's passed into the constructor. Then, we implement the callable interface by overriding the `__call__` method to return the result from calling `ScalarMappable`'s  `.to_rgba` function.\n",
    "\n",
    "## Drawing the Choropleth\n",
    "\n",
    "In this section we'll create a function that will draw the choropleth for a given set of unemployment rates. One of the things we'll need to create our `draw_choropleth` function, is a list of FIPS codes for all 48 contiguous states since we don't want to draw the counties for any of the states outside of the contiguous United States."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# List of the FIPS codes for the 48 contiguous states \n",
    "# source: http://www.bls.gov/cew/cewedr10.htm)\n",
    "contiguous_states_fips = [\n",
    "    '01',   # Alabama \n",
    "    '04',   # Arizona \n",
    "    '05',   # Arkansas    \n",
    "    '06',   # California  \n",
    "    '08',   # Colorado    \n",
    "    '09',   # Connecticut \n",
    "    '10',   # Delaware    \n",
    "    '11',   # District of Columbia    \n",
    "    '12',   # Florida \n",
    "    '13',   # Georgia \n",
    "    '16',   # Idaho   \n",
    "    '17',   # Illinois    \n",
    "    '18',   # Indiana \n",
    "    '19',   # Iowa    \n",
    "    '20',   # Kansas  \n",
    "    '21',   # Kentucky    \n",
    "    '22',   # Louisiana   \n",
    "    '23',   # Maine   \n",
    "    '24',   # Maryland    \n",
    "    '25',   # Massachusetts   \n",
    "    '26',   # Michigan    \n",
    "    '27',   # Minnesota   \n",
    "    '28',   # Mississippi \n",
    "    '29',   # Missouri    \n",
    "    '30',   # Montana \n",
    "    '31',   # Nebraska    \n",
    "    '32',   # Nevada  \n",
    "    '33',   # New Hampshire   \n",
    "    '34',   # New Jersey  \n",
    "    '35',   # New Mexico  \n",
    "    '36',   # New York    \n",
    "    '37',   # North Carolina  \n",
    "    '38',   # North Dakota    \n",
    "    '39',   # Ohio    \n",
    "    '40',   # Oklahoma    \n",
    "    '41',   # Oregon  \n",
    "    '42',   # Pennsylvania    \n",
    "    '44',   # Rhode Island    \n",
    "    '45',   # South Carolina  \n",
    "    '46',   # South Dakota    \n",
    "    '47',   # Tennessee   \n",
    "    '48',   # Texas   \n",
    "    '49',   # Utah    \n",
    "    '50',   # Vermont \n",
    "    '51',   # Virginia    \n",
    "    '53',   # Washington  \n",
    "    '54',   # West Virginia   \n",
    "    '55',   # Wisconsin   \n",
    "    '56',   # Wyoming \n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our `draw_choropleth` function takes 3 parameters: a `Basemap` class instance, a list of `UnemploymentRecord` tuples, and a callable that maps unemployment rates to colors.\n",
    "\n",
    "The code for our drawing function itself is relatively simple since it assumes that the Basemap object is completely setup at the time it's passed in. By \"completely setup\", I mean that the shapefile for the counties has already been loaded into the `Basemap` object, so we can simply loop through each of the counties in the map and draw each one with its associated color. To draw a county on our map, we create an instance of the `Polygon` class, with a county's coordinates, and pass it into the `Axes.add_patch` function, which we can see happening on line 27 below. \n",
    "\n",
    "So, the major portion of the code below, specifically lines 17-27, is just a big loop, where for each county in the map, we get a color for its unemployment rate, draw a polygon with the county's coordinates, fill it with the corresponding color, and add it to our map. This will give the map that familiar patchwork look of a choropleth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.basemap import Polygon\n",
    "\n",
    "def draw_choropleth(m, data, color):\n",
    "    \"\"\"Draws a choropleth.\n",
    "    \n",
    "    Arguments:\n",
    "    m -- The map projection\n",
    "    data -- The unemployment data\n",
    "    color -- A function that takes an unemployment rate and returns \n",
    "             a hex value representing a color.\n",
    "    \"\"\"\n",
    "    # Convert the unemployment data into a dict mapping FIPS codes to unemployment_rates\n",
    "    unemployment_rates = {r.fips_code: r.unemployment_rate for r in data}\n",
    "    \n",
    "    # Draw all of the counties to the current axes\n",
    "    ax = plt.gca()\n",
    "    for info, shape in zip(m.counties_info, m.counties):\n",
    "        # We're only showing unemployment rates for the 48 contiguous states, \n",
    "        # so skip anything that doesn't show up in the list.\n",
    "        if info['STATEFP'] not in contiguous_states_fips:\n",
    "            continue\n",
    "        try:\n",
    "            fips = info['STATEFP'] + info['COUNTYFP']\n",
    "            unemployment_rate = unemployment_rates[fips]\n",
    "        except KeyError:\n",
    "            continue\n",
    "        ax.add_patch(Polygon(np.array(shape), facecolor=color(unemployment_rate), \n",
    "                             edgecolor='#221e1f', linewidth=0.5))    \n",
    "\n",
    "    # Calling drawmapboundary will ensure that the x- and y-axis \n",
    "    # limits are set correctly for the map and draws either simple \n",
    "    # border (or no border as in our example below) in place of a \n",
    "    # typical x- and y-axis with tick marks.\n",
    "    m.drawmapboundary(color='none');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create the Map Object\n",
    "\n",
    "The next bit of code creates our `Basemap` instance. In the last lesson, we used the Robinson projection, which is great for showing a map of the world, but it distorts both the shape and size of the land masses closer to the poles. For this lesson, we'll use the [Lambert Conformal Conic][1] projection, which should show us a more accurate portrayal of the United States since it does a better job of preserving the shape of the continents.\n",
    "\n",
    "To use this projection, we have to pass in the lat/lon coordinates for the region of the map that we want to draw. We can get that information from Google Maps by simply ctrl-clicking (right clicking) on a point on the map, and choosing the \"What's here?\" option. \n",
    "\n",
    "Next, we need either 1 or 2 lines of latitude for the standard parallels. These are used by the projection to adjust where the largest amounts of distortion in the map arise, and the closer these are to the area that you want to display, the smaller the amount of distortion. We're going to use the coordinates for the [geographic center of the United States][2] to limit the distortion as much as possible. \n",
    "\n",
    "Finally, we need to choose the line of longitude for the center of the map, and again we'll use the geographic center of the US here as well.\n",
    "\n",
    "[1]: https://en.wikipedia.org/wiki/Lambert_conformal_conic_projection\n",
    "[2]: https://en.wikipedia.org/wiki/Geographic_center_of_the_United_States"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.basemap import Basemap\n",
    "\n",
    "# Lambert Conformal map of lower 48 states.\n",
    "m = Basemap(\n",
    "    # These four attributes define the bounding box, \n",
    "    # i.e., the area of the world that we want to show \n",
    "    # in our map.\n",
    "    llcrnrlon=-119,\n",
    "    llcrnrlat=22,\n",
    "    urcrnrlon=-64,\n",
    "    urcrnrlat=49,\n",
    "    # We've chosen the Lambert Conformal Conic projection. \n",
    "    projection='lcc',\n",
    "    # These two params define the 1st (and 2nd) standard parallel\n",
    "    # and the central meridian used by the projection.\n",
    "    lat_1=39,\n",
    "    lon_0=-98\n",
    "   )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the Counties Shapefile\n",
    "\n",
    "While `Basemap` instances have methods for drawing things like lakes and rivers, coastlines, and political boundaries, they don't have anything for drawing the boundaries of the counties of the US built in. `Basemap` instances do, however, have the ability to read in a [shapefile][1] to make it possible to draw the county boundaries ourselves. \n",
    "\n",
    "The next line of code reads in the [counties shapefile][2] from the United States Census website, which I've already downloaded and included in the repository to make things a tad bit easier.\n",
    "\n",
    "[1]: https://en.wikipedia.org/wiki/Shapefile\n",
    "[2]: https://www.census.gov/geo/maps-data/data/cbf/cbf_counties.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import wget\n",
    "\n",
    "link_to_data = 'http://www2.census.gov/geo/tiger/GENZ2016/shp/cb_2016_us_county_500k.zip'\n",
    "DataSet = wget.download(link_to_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!ls *.zip\n",
    "!ls ./"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set drawbounds to False to make sure the map isn't drawn in this step\n",
    "m.readshapefile('./cb_2016_us_county_500k', 'counties', \n",
    "                drawbounds=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Great Recession\n",
    "\n",
    "Well, helper functions written and shapefiles read in, we're finally ready to explore the data. The next little bit of code loads in three years of unemployment data and draws a choropleth map for each one. The years that we show are 2006 (the year right before the Great Recession), 2009 (the recession at its peak), and 2014 (the latest year for which we have census data). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# The Great Recession began in 2007 and hit its peak in 2009. \n",
    "# We'll look at the unemployment rate, before it began, at its\n",
    "# peak, and finally the most recent year that we have data for.\n",
    "years = [2006, 2009, 2014]\n",
    "\n",
    "plt.figure(figsize=(20, 10))\n",
    "for i, year in enumerate(years):\n",
    "\n",
    "    # Get the data from the file\n",
    "    data = load_unemployment_rates(year)\n",
    "    color = HeatMapper([r.unemployment_rate for r in data])\n",
    "\n",
    "    # Create a new subplot and add a title\n",
    "    plt.subplot(1, len(years), i+1)\n",
    "    plt.title('Unemployment for %s' % year)\n",
    "    \n",
    "    # Draw the map\n",
    "    draw_choropleth(m, data, color)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the maps above, it's easy to see the damage that the recession did to the country's employment rate. 2006 shows unemployment rates to be fairly low across nearly all of the states, with slightly higher unemployment rates on the west cost and the southeast. In the 2009 map, on the other hand, we see the west coast and the southeast really heat up as unemployment rates reach their peak. It also looks like Michigan's employment rate took a bit hit during the recession as well, which makes sense since saw such a steep decline in manufacturing jobs across the US. Finally, in the 2014 map, we can see that the country as a whole has made a complete recovery, and in some cases, the midwest in particular, we even see a decrease in unemployment rates compared to the pre-recession map of 2006.\n",
    "\n",
    "## The Peak of the Recession\n",
    "\n",
    "Lastly, just to get a better view of the peak of the recession, and to show off just how good our choropleth's look, we blow up the 2009 map in the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,10))\n",
    "plt.title(\"Unemployment for 2009\")\n",
    "data = load_unemployment_rates(2009)\n",
    "color = HeatMapper([r.unemployment_rate for r in data])\n",
    "draw_choropleth(m, data, color)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "In prior lessons, we first learned about map projections in general, and then we saw how to use the basemap library to create our own maps. In this lesson, we took our knowledge of the basemap library a bit further by using it to create a choropleth map that visualized the unemployment data in the US. Mapping data like this has a very distinct advantage over other visualization techniques, since it allows us to quickly see where, and to what extent, a particular phenomenon has on a particular geographic area of the population.\n",
    "\n",
    "In the next, and final basemap lesson, we'll see yet another example of just how useful the basemap library can be by mapping local crime data using maps from an ArcGIS online service."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}