Maps of earthquakes in Japan, Korea and China, 1970-2013

by David Taylor, www.prooffreader.com
Uses Python's pandas, matplotlib and basemap

Steps before this script was run:

Visited http://www.iris.edu/ieb/ -- the IRIS Earthquake Browser
Downloaded .xls files and opened them in Excel
kept columns named magnitude, lon and lat
parsed date to make always 10 characters, e.g. 1/1/1970 becomes 01/01/1970
parsed date to create month, day and year columns (integers)
created days_of_year column counting number of days since the beginning of that year
All of this could be done in Python; IRIS has a REST API from which the same data can be downloaded. This repo contains a script called earthquakes-rest-api.py to do so.



In [1]:

    
import pandas as pd
quakes = pd.read_csv('jp-earthquakes-5plus.csv')

quakes[['year', 'month', 'day']] = quakes[['year', 'month', 'day']].astype(int)
quakes.sort('magnitude', ascending=True, inplace=True)

quakes.to_pickle('jpquakes.pickle')

print quakes.tail()









    



      magnitude        date    lat     lon  month  day  year  days_of_year  \
6616        8.1  09/25/2003  41.75  143.87      9   25  2003           267   
7341        8.1  01/13/2007  46.23  154.50      1   13  2007            12   
5014        8.2  10/04/1994  43.66  147.38     10    4  1994           276   
7222        8.3  11/15/2006  46.68  153.21     11   15  2006           318   
8368        9.0  03/11/2011  38.30  142.37      3   11  2011            69   

      month_half  
6616           1  
7341           0  
5014           0  
7222           1  
8368           0



In [2]:

    
prefix = 'map_'  # for saved maps

import pandas as pd
quakes = pd.read_pickle('jpquakes.pickle')

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

alphalist = [0.95, 0.85, 0.75, 0.65, 0.55, 0.45, 0.35]

def get_alpha(magnitude):
    if magnitude < 5.5:
        return alphalist[0]
    elif magnitude < 6:
        return alphalist[1]
    elif magnitude < 6.5:
        return alphalist[2]
    elif magnitude < 7.0:
        return alphalist[3]
    elif magnitude < 7.5:
        return alphalist[4]
    elif magnitude < 8.0:
        return alphalist[5]
    else:
        return alphalist[6]

for this_year in range(2011,2012):
     for this_month in range(3,4): 
        
            plt.figure(figsize=(15,10))
            map = Basemap(projection='merc', lat_0 = 57, lon_0 = 120,
                          resolution = 'c', area_thresh = 0.1,
                          llcrnrlon=83, llcrnrlat=23,
                          urcrnrlon=156, urcrnrlat=52)
            
            map.warpimage(image='darkrelief.jpg')
            
            tempdf = quakes[(quakes.year == this_year) & (quakes.month == this_month)]
            latlist = list(tempdf.lat)
            lonlist = list(tempdf.lon)
            maglist = list(tempdf.magnitude)
            
            label=str(this_year) + "/"
            
            if this_month < 10: label += '0'
                
            label += str(this_month)
            
            plt.annotate(label, xy=(0.01, 0.97), ha = 'left', va='top', xycoords='axes fraction', color='#ffff33', 
                         fontsize=35, fontname = 'Consolas')#, bbox=dict(facecolor='#ffffaa', alpha=0.8))
            
            min_marker_size = 20
            for lon, lat, mag in zip(lonlist, latlist, maglist):
                x,y = map(lon, lat)
                msize = (mag - 4.5)**1.5 * min_marker_size
                thisalpha = get_alpha(mag)
                map.plot(x, y, color='#ffff00', marker='o', markersize=msize, alpha=thisalpha)
                            
            plt.savefig(prefix+label.replace('/', '')+'.png')
            plt.show()
            plt.close()
            print this_year, this_month



In [3]:

    
# to draw circles for legend

circles = pd.read_csv('circles.csv')

# circles has just five data points set in months 5-9 of the year 1900

circles[['year', 'month', 'day']] = circles[['year', 'month', 'day']].astype(int)
circles.sort('magnitude', ascending=True, inplace=True)

prefix = 'circles'


alphalist = [0.95, 0.85, 0.75, 0.65, 0.55, 0.45, 0.35]

def get_alpha(magnitude):
    if magnitude < 5.5:
        return alphalist[0]
    elif magnitude < 6:
        return alphalist[1]
    elif magnitude < 6.5:
        return alphalist[2]
    elif magnitude < 7.0:
        return alphalist[3]
    elif magnitude < 7.5:
        return alphalist[4]
    elif magnitude < 8.0:
        return alphalist[5]
    else:
        return alphalist[6]

for this_year in range(1900,1901):
     for this_month in range(5,10): 
        
            plt.figure(figsize=(15,10))
            map = Basemap(projection='merc', lat_0 = 57, lon_0 = 120,
                          resolution = 'c', area_thresh = 0.1,
                          llcrnrlon=83, llcrnrlat=23,
                          urcrnrlon=156, urcrnrlat=52)
            
            map.warpimage(image='darkrelief.jpg')
            
            tempdf = circles[(circles.year == this_year) & (circles.month == this_month)]
            latlist = list(tempdf.lat)
            lonlist = list(tempdf.lon)
            maglist = list(tempdf.magnitude)
            
            label='Magnitude '
                
            label += str(this_month)
            
            plt.annotate(label, xy=(0.01, 0.97), ha = 'left', va='top', xycoords='axes fraction', color='#ffff33', 
                         fontsize=35, fontname = 'Consolas')#, bbox=dict(facecolor='#ffffaa', alpha=0.8))
            
            min_marker_size = 20
            for lon, lat, mag in zip(lonlist, latlist, maglist):
                x,y = map(lon, lat)
                msize = (mag - 4.5)**1.5 * min_marker_size
                thisalpha = get_alpha(mag)
                map.plot(x, y, color='#ffff00', marker='o', markersize=msize, alpha=thisalpha)
                            
            
            
            plt.savefig(prefix+label.replace('/', '')+'.png')
            #plt.show()
            plt.close()



In [4]:

    
# blank map

prefix = 'good_dt_relief__blank'  # for saved maps

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np



plt.figure(figsize=(15,10))
map = Basemap(projection='merc', lat_0 = 57, lon_0 = 120,
              resolution = 'c', area_thresh = 0.1,
              llcrnrlon=83, llcrnrlat=23,
              urcrnrlon=156, urcrnrlat=52)

map.warpimage(image='darkrelief.jpg')
                
plt.savefig(prefix+'.png')
plt.show()
plt.close()



In [5]:

    
# all circles superimposed

prefix = 'all_'  

import pandas as pd
quakes = pd.read_pickle('jpquakes.pickle')

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

alphalist = [0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1]

def get_alpha(magnitude):
    if magnitude < 5.5:
        return alphalist[0]
    elif magnitude < 6:
        return alphalist[1]
    elif magnitude < 6.5:
        return alphalist[2]
    elif magnitude < 7.0:
        return alphalist[3]
    elif magnitude < 7.5:
        return alphalist[4]
    elif magnitude < 8.0:
        return alphalist[5]
    else:
        return alphalist[6]

plt.figure(figsize=(15,10))
map = Basemap(projection='merc', lat_0 = 57, lon_0 = 120,
              resolution = 'c', area_thresh = 0.1,
              llcrnrlon=83, llcrnrlat=23,
              urcrnrlon=156, urcrnrlat=52)

map.warpimage(image='darkrelief.jpg')
    
for this_year in range(1970,2014):
     for this_month in range(1,13): 
            
            tempdf = quakes[(quakes.year == this_year) & (quakes.month == this_month)]
            latlist = list(tempdf.lat)
            lonlist = list(tempdf.lon)
            maglist = list(tempdf.magnitude)
                    
            min_marker_size = 20
            for lon, lat, mag in zip(lonlist, latlist, maglist):
                x,y = map(lon, lat)
                msize = (mag - 4.5)**1.5 * min_marker_size
                thisalpha = get_alpha(mag)
                map.plot(x, y, color='#ffff00', marker='o', markersize=msize, alpha=thisalpha)
                            
label='1970-2013'

plt.annotate(label, xy=(0.01, 0.97), ha = 'left', va='top', xycoords='axes fraction', color='#000000', 
             fontsize=30, fontname = 'Arial')
            
plt.savefig(prefix+label.replace('/', '')+'.png')
plt.show()
plt.close()



In [ ]: