Earthquake in Basemap Demo

In [1]:

    

import datetime
from dateutil.tz import tzutc, tzlocal

#External modules.
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()


    

In [2]:

    

#Load earthquake data (either live or from file).
LIVE_DATA = True
import json
if LIVE_DATA:
    import urllib
    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" % (
                            starttime, endtime))
    o = json.load(f)
    f.close()
    title = "Earthquake Magnitudes %s - %s" % (start_dt.strftime("%d %b %Y %H:%M %Z"), end_dt.strftime("%d %b %Y %H:%M %Z"))
else:
    title = "Earthquake Magnitudes - Sep. 25 2013 - Sep. 26 2013"
    f = open("./earthquake.json", "r")
    o = json.load(f)
    f.close()


    

In [3]:

    

#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[0], coords[1]
    props = feature['properties']
    mag = props['mag']
    data.append((lon, lat, mag))
del o


    

In [5]:

    

#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))
#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)
    lons.append(xpt)
    lats.append(ypt)
    mags.append(mag)
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)
plt.colorbar(sc, orientation='vertical')

#Title
r = plt.title(title)

#Day/night shading
class DateWrapper(object):
    """
    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)

    def utcoffset(self):
        return None
date = DateWrapper(datetime.datetime.today().replace(tzinfo=tz_local).astimezone(tz_utc))
CS=m.nightshade(date)

#Tight layout
fig.tight_layout()


    

In [4]: