I had previously written about the Python library maplotlib and the cool 2D and 3D plots you could create with it. Recently, I noticed that a lot of folks have been building tool kits on top of matplotlib for various scientific domains. One toolkit that caught my eye was Basemap, a library used to plot data on maps.
What is really cool about Basemap is that it offers 30 different map projections of the Earth, right out of the box. And since the underlying data plotting mechanisms are all based on matplotlib, it didn’t take me long to figure out how to put markers on my location.
I stumbled upon a really nice tutorial for Basemap. The earlier lessons were a great supplement to the Basemap documentation. The tutorial I linked to shows how the author plotted earthquake data on a map using basemap. This article was 8 months old, though, and the government web service seemed to have changed. Also, I had recently been playing with matplotlib’s image processing capabilities, and I wanted to put a color bar on my plots as a key to show what the colors meant data-wise. So I shamelessly “borrowed” some of the author’s code and updated the data extraction and parsing, and added some new features to the plot. The result is the map I generated below:
Earthquake data plotted with matplotlib and Basemap
I created this interactively using an IPython notebook. If you are not familiar with IPython, it is a really cool interactive shell for Python. And the notebook feature lets you run the shell as a back end process that runs a mini web server so you can have an interactive Python session in your web browser. Best of all, plots from matplotlib can be rendered right in the browser. These notebooks can be saved and shared.
The code from my notebook is below. If you want to run it from a script, you will need to add a statement to show or save the plot at the end of the script (e.g. plt.show()). I have a variable, LIVE_DATA, defined and set to True. This causes the the data to be loaded from the USGS government web service. If you want to experiment and don’t want to keep downloading the data, you can download the data once by using curl or even just placing the HTTP GET request in your browser and saving the result to a file. Set LIVE_DATA to False and make sure the code is attempting to load from the correct file (I was loading from ‘earthquake.json’ in the folder I launched the notebook from). You may also want to experiment with various map types. I chose ‘moll’ for this demonstration, but I also played with ‘robin’ and ‘ortho’. Some types of maps require different information, so you may need to add or change parameters.
# Plot Earthquake data.
from dateutil.tz import tzutc, tzlocal
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
from matplotlib.cm import get_cmap, ScalarMappable
from mpl_toolkits.basemap import Basemap
import numpy as np
#Time zone info.
tz_utc = tzutc()
tz_local = tzlocal()
#Set figure size.
fig = figure(figsize=(14,10))
#Choose map type.
# make sure the value of resolution is a lowercase L,
# for 'low', not a numeral 1
#m = Basemap(projection='ortho', lat_0=50, lon_0=-100, resolution='l', area_thresh=1000.0)
m = Basemap(projection='moll', lat_0=0, lon_0=-100, resolution='l', area_thresh=1000.0)
#Fill in map features.
r = m.drawcoastlines()
r = m.fillcontinents(color='coral',lake_color='aqua')
r = m.drawmapboundary(fill_color='aqua')
r = m.drawmeridians(np.arange(0, 360, 30))
r = m.drawparallels(np.arange(-90, 90, 30))
#Load earthquake data (either live or from file).
LIVE_DATA = True
end_dt = datetime.datetime.today().replace(tzinfo=tz_local)
start_dt = (end_dt - datetime.timedelta(days=1))
starttime = start_dt.strftime("%Y%m%dT%H:%M:%S%z")
endtime = end_dt.strftime("%Y%m%dT%H:%M:%S%z")
f = urllib.urlopen("http://comcat.cr.usgs.gov/fdsnws/event/1/query?starttime=%s&endtime=%s&format=geojson&eventtype=earthquake" % (
o = json.load(f)
title = "Earthquake Magnitudes %s - %s" % (start_dt.strftime("%d %b %Y %H:%M %Z"), end_dt.strftime("%d %b %Y %H:%M %Z"))
title = "Earthquake Magnitudes - Sep. 25 2013 - Sep. 26 2013"
f = open("./earthquake.json", "r")
o = json.load(f)
#Parse data and extract the values in which we are interested.
features = o['features']
data = 
for feature in features:
geo = feature['geometry']
coords = geo['coordinates']
lon, lat = coords, coords
props = feature['properties']
mag = props['mag']
data.append((lon, lat, mag))
#Set up a normalizer and color map for the scatter plot.
norm = Normalize()
cmap = get_cmap('jet')
#Create series for longitude, latitude, and magnitude.
lons = 
lats = 
mags = 
for lon, lat, mag in data:
xpt,ypt = m(lon,lat)
x = np.array(lons)
y = np.array(lats)
#Create normalized magnitudes (range [0, 1.0]) for color map.
mag_norms = norm(np.array(mags))
#Create marker sizes as a function of magnitude.
z = (mag_norms * 10.0)**2.0
#Plot the data and create a colorbar.
sc = m.scatter(x,y, s=z, zorder=4, cmap=cmap, c=mags)
r = plt.title(title)
This class is a hack-- basemap is expecting a UTC time to have `None` as the utcoffset attrib.
def __init__(self, dt):
self._dt = dt
def __getattr__(self, x):
return getattr(self._dt, x)
date = DateWrapper(datetime.datetime.today().replace(tzinfo=tz_local).astimezone(tz_utc))