Matplotlib Basemap tutorial 02 : Let’s add some earthquakes !

Now, let’s imagine we have a dataset containing latitude/longitudes of earthquakes, plus their depth and magnitude. Of course, you don’t always have this dataset available, so let’s build a fake one :

import numpy as np
lon = np.random.random_integers(11,79,1000)/10.
lat = np.random.random_integers(491,519,1000)/10.
depth = np.random.random_integers(0,300,1000)/10.
magnitude = np.random.random_integers(0,100,1000)/10.

We create random integer values, and we devide each of the the 1000-element array by 10.0.

Then, time is to convert from lat/lon world to matplotlib canvas world :

# with x,y=m(lon,lat) we calculate the x,y position of each quake in the map coordinates
x,y = m(lon,lat)

And finaly, to plot the earthquakes, with a color scale linked to the depth, and the size of the scatter point linked to the magnitude (in fact, 10 times the magnitude, in order to see something):

x,y = m(lon,lat)
m.scatter(x,y,s=10*magnitude,c=depth)

See the full code after the break ;

#

# BaseMap example by geophysique.be

# tutorial 02

from mpl_toolkits.basemap import Basemap

import matplotlib.pyplot as plt

import numpy as np

fig = plt.figure(figsize=(11.7,8.3))

#Custom adjust of the subplots
plt.subplots_adjust(left=0.05,right=0.95,top=0.90,bottom=0.05,wspace=0.15,hspace=0.05)
ax = plt.subplot(111)

#Let's create a basemap around Belgium
m = Basemap(resolution='i',projection='merc', llcrnrlat=49.0,urcrnrlat=52.0,llcrnrlon=1.,urcrnrlon=8.0,lat_ts=51.0)
m.drawcountries(linewidth=0.5)
m.drawcoastlines(linewidth=0.5)
m.drawparallels(np.arange(49.,53.,1.),labels=[1,0,0,0],color='black',dashes=[1,0],labelstyle='+/-',linewidth=0.2) # draw parallels
m.drawmeridians(np.arange(1.,9.,1.),labels=[0,0,0,1],color='black',dashes=[1,0],labelstyle='+/-',linewidth=0.2) # draw meridians

lon = np.random.random_integers(11,79,1000)/10.
lat = np.random.random_integers(491,519,1000)/10.
depth = np.random.random_integers(0,300,1000)/10.
magnitude = np.random.random_integers(0,100,1000)/10.
# with x,y=m(lon,lat) we calculate the x,y position of each quake in the map coordinates
x,y = m(lon,lat)
m.scatter(x,y,s=10*magnitude,c=depth) # we will scale the dots by 10 time the magnitude
c = plt.colorbar(orientation='horizontal')
c.set_label("Depth")
plt.show()

Leave a Reply

Your email address will not be published. Required fields are marked *


*