How to use matplotlib's basemap and Python 3.x to plot earthquake data on a world map

matplotlib-basemap-python3x-quakes.md

Mapping earthquakes in Python 3.x using matplotlib and matplotlib's basemap

(with a big assist from Anaconda)

I currently use Python for nearly all of my data science and wrangling work these days but usually find myself switching to R to visualize data using ggplot2. This is due in part to ggplot2's general excellence, but also because I had a lot of trouble learning Python's most popular viz library, matplotlib on my own...its homepage is decent enough...but its variety of plotting APIs (--pylab? OOP? %matplotlib???) has led to widely differing examples and best practices among the many online matplotlib guides (not dissimilar to the general problem of trying to practice either Python 2.x or 3.x).

That changed yesterday when I stumbled across Nicolas P. Rougier's beautifully designed and comprehensive -- but plainly-titled, Matplotlib tutorial, which includes easy-to-follow and interesting geovisualization examples, among many other elegant demonstrations. That was enough to get me to try matplotlib again.

Another big problem for me has been the lack of well-documented best practices for GIS work for Python 3.x (there are plenty of 2.x examples). However, I recently found this mapping tutorial from the Intro to Python guide. The best part of it for me was the installation instructions, which alluded to the possibility that the Anaconda Python 3.x framework includes a recipe for Basemap, a plugin for matplotlib for plotting data on map projections.

Lo and behold, it was as easy as this on OS X:

   $ conda install basemap

-- which added this to my Anaconda package listing:

   basemap                   1.0.7               np110py35_0    defaults

(the value and convenience that Anaconda brings to casual Python programmers like me cannot be overstated...and I highly recommend it if you're new to Python)

And now, map plotting just works, in the way that I have long since done it with ggplot2: render a shapefile and then plot the geocoordinates like any other scatterplot. Here's my result:

You can see Rougier's tutorial on mapping the earthquakes here: http://www.labri.fr/perso/nrougier/teaching/matplotlib/matplotlib.html#earthquakes

His example shows how to animate the earthquakes, similar to an animated torque map at CartoDB.

I've opted to show less matplotlib-specific code, including pretending that numpy and pandas don't exist (I want to teach matplotlib initially by showing that lists of values can serve as the plotting arguments, and that numpy/pandas just expedite the process).

The code snippets below obviously should go into one file, but I've broken them up by function:

1. Downloading and parsing the data from USGS
2. Rendering the map and data

Download the USGS quakes data

The USGS earthquakes feed can be found at its Earthquake Hazards Program:

http://earthquake.usgs.gov/earthquakes/feed/v1.0/csv.php

import requests
from csv import DictReader
DATA_URL = 'http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/4.5_month.csv'
print("Downloading", DATA_URL)
resp = requests.get(DATA_URL)
quakes = list(DictReader(resp.text.splitlines()))
# ...avoiding numpy/pandas Array() for now, and can't care enough to do this less awkwardly...
lngs = [float(q['longitude']) for q in quakes]
lats = [float(q['latitude']) for q in quakes]
mags = [2 ** float(q['mag']) for q in quakes]

Create the map

# %matplotlib  # in iPython
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
plt.figure(figsize=(14, 8))
earth = Basemap()
earth.bluemarble(alpha=0.42)
earth.drawcoastlines(color='#555566', linewidth=1)
plt.scatter(lngs, lats, mags, 
            c='red',alpha=0.5, zorder=10)
plt.xlabel("M4.5 earthquakes in the past 30 days from March 18, 2016 (USGS)")
plt.savefig('usgs-4.5quakes-bluemarble.png', dpi=350)

Note: I'm deliberately not using matplotlib's OOP style in favor of the global plt object, just to show the relative bare minimum of code to produce a visualization in matplotlib...however, I am avoiding the use of pylab, which seems to be a legacy practice on the verge of being deprecated.

The USGS earthquakes data file is provided below in this gist: usgs-4.5M-month-quakes.csv

Saving as SVG

I always thought you needed Cairo -- which is currently not trivial to install on 3.x and OS X, as far as I can tell -- to export as SVG. But it turns out it's as easy as just telling matplotlib to save as '.svg'...To produce a SVG file (which I've included below), just modify this single line:

plt.savefig('usgs-4.5quakes.svg')

I also took out the Blue Marble raster since I don't know how that would render as SVG:

# earth.bluemarble(alpha=0.42)

Using plt.subplots()

Call Basemap using the ax argument:

fig, ax = plt.subplots()
earth = Basemap(ax=ax)
earth.drawcoastlines(color='#555566', linewidth=1)
ax.scatter(lngs, lats, mags, c='red',alpha=0.5, zorder=10)
ax.set_xlabel("Hello there")
fig.savefig('hello-ax.png')

Using pandas

import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap


DATA_URL = 'http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/4.5_month.csv'
print("Downloading", DATA_URL)
df = pd.read_csv(DATA_URL)

fig, ax = plt.subplots()
earth = Basemap(ax=ax)
earth.drawcoastlines(color='#556655', linewidth=0.5)
ax.scatter(df['longitude'], df['latitude'], df['mag'] ** 2, 
           c='red', alpha=0.5, zorder=10)
ax.set_xlabel("This month's 4.5M+ earthquakes")
fig.savefig('usgs-monthly-4.5M.png')

usgs-4.5M-month-quakes.csv

usgs-4.5M-month-quakes.svg

what projection did you use ? how did you get the coordinates to match?

@johnmisquita the coordinates of the quakes are latlong, so the matching happens within basemap at the same time it's fitting the background to the projection, afaik.