Geophysics Deskwork

Scott's Non-Field Notebook

Mapping the largest Holocene Volcanic Eruptions

This iPython Notebook demonstrates how to generate a map of major Holocene volcanic eruptions. I think it's a good introduction to basic analysis and visualization of georeferenced data with Python tools.

Data is from the Smithsonian Global Volcanism Program. Recently the website was updated with improved search tools and the ability to download database search results to .csv format. I did a global search for eruptions with a volcano explosivity index (VEI) greater than 4. This index is a measure of the magnitude of an explosive volcanic eruption, and can be converted to erupted volume on a logarithmic scale:

VEI Description Plume Height Volume Classification How often Example
0 non-explosive < 100 m 1000s m3 Hawaiian daily Kilauea
1 gentle 100-1000 m 10,000s m3 Haw/Strombolian daily Stromboli
2 explosive 1-5 km 1,000,000s m3 Strom/Vulcanian weekly Galeras, 1992
3 severe 3-15 km 10,000,000s m3 Vulcanian yearly Ruiz, 1985
4 cataclysmic 10-25 km 100,000,000s m3 Vulc/Plinian 10's of years Galunggung, 1982
5 paroxysmal >25 km 1 km3 Plinian 100's of years St. Helens, 1980
6 colossal >25 km 10s km3 Plin/Ultra-Plinian 100's of years Krakatau, 1883
7 super-colossal >25 km 100s km3 Ultra-Plinian 1000's of years Tambora, 1815
8 mega-colossal >25 km 1,000s km3 Ultra-Plinian 10,000's of years Yellowstone, 2 Ma

table copied from Oregon State's Volcano World website

We'll use Numpy and Matplotlib for analyzing data and plotting. Load the libraries and make some changes to default plot appearance for starters:

In [20]: 
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['figure.figsize'] = (11.0, 8.5)
fontdict = dict(family='sans-serif',weight='bold',size=16) #NOTE: doesn't effect axes labels
plt.rc('font', **fontdict)

Reading .CSV files with python

An extremely useful utility comes with matplotlib to read csv data into a numpy record array. This automatically converts columns of data to strings, integers, and floats!

In [21]: 
from matplotlib.mlab import csv2rec, rec2txt #ported matlab functions
gvp = csv2rec('/Users/scott/Desktop/andes_seminar/GVP_VEI4_Global.csv',skiprows=1)

#print gvp.dtype.fields # uncomment to see all column headers
print 'There are {0} recorded eruptions in this database'.format(gvp.size)
volcanos = np.unique(gvp.volcano_name)
print 'The eruptions have occured at {0} volcanos'.format(volcanos.size)
ind = np.argmax(gvp.start_year)
print 'The most recent eruption was a VEI {0} at {1} in {2}'.format(gvp.vei[ind], gvp.volcano_name[ind], gvp.start_year[ind])

There are 713 recorded eruptions in this database
The eruptions have occured at 245 volcanos
The most recent eruption was a VEI 4 at Eyjafjallaj�kull in 2010

Note that many volcanos have names with strange characters, you can fix them manually this way:

In [22]: 
gvp.volcano_name[ind] = 'Eyjafjallajökull' #manual fix
print gvp.volcano_name[ind]

Eyjafjallajökull

Let's look at a histogram of this data set, binned by eruption size. The colors are set in a list to match the icons we'll use on the map:

In [23]: 
fig, ax = plt.subplots()

freq, bins, patches = plt.hist(gvp.vei, bins=np.arange(9),
                                lw=2,
                                alpha=0.5,
                                align='left')
colors = ['k','k','white', 'gray', 'yellow', 'orange', 'red', 'purple']
for p,c in zip(patches, colors): # change bar colors to match map markers
    print p,c
    p.set_facecolor(c)
ax.set_xlabel('VEI')
ax.set_xlim(0,8)
ax.grid()
ax.set_title('Histogram of Large Holocene Eruptions')
plt.savefig('histogram.png', bbox_inches='tight')

Rectangle(-0.5,0;1x0) k
Rectangle(0.5,0;1x0) k
Rectangle(1.5,0;1x0) white
Rectangle(2.5,0;1x0) gray
Rectangle(3.5,0;1x486) yellow
Rectangle(4.5,0;1x172) orange
Rectangle(5.5,0;1x49) red
Rectangle(6.5,0;1x6) purple

Plotting Data on a Map

Matplotlib is a fantastic plotting library, and the Basemap Toolkit is an excellent way to make publication quality maps.

In [24]: 
import mpl_toolkits.basemap as bm
import matplotlib.patheffects as PathEffects

# Initialize Global Map
bm.latlon_default = True
fig = plt.figure(figsize=(11, 8.5)) 
bmap = bm.Basemap(projection='robin',lon_0=0,resolution='c') #Robinson Projection
bmap.drawcoastlines()
Out[24]:
<matplotlib.collections.LineCollection at 0x10d57bdd0>

That's easy... but boring. Let's add a background. Basemap comes with some built-in background options. For example, shaded relief:

In [25]: 
bmap.shadedrelief(scale=0.25)
Out[25]:
<matplotlib.image.AxesImage at 0x10d749590>

The default shaded relief background image comes from the Natural Earth website, which has a variety of professional raster images. The "Gray Earth" backdrop is particularly good if you want data to stand out. As long as you have an image with global coverage it's easy to set the background with the warpimage() Basemap method.

In [26]: 
background = '/Volumes/OptiHDD/data/natural_earth/medium_res/GRAY_50M_SR_OB/GRAY_50M_SR_OB.jpg'
bmap.warpimage(image=background, scale=0.25)
Out[26]:
<matplotlib.image.AxesImage at 0x10d771850>

The scale keyword is used to downsample the resolution by 4x, since we are zoomed out so far. Now we can plot the locations of the eruptions on this background:

In [27]: 
def plot_volcanos(volcanos, ms=12, color='r',alpha=0.5):
    ''' draw volcanos on basemap '''
    for i in xrange(volcanos.size):
        v = volcanos[i]
        mark, = bmap.plot(v.longitude, v.latitude, marker='^', ms=ms, mfc=color, mec='white', mew=1, ls='none', alpha=alpha)
    return mark

l = plot_volcanos(gvp)
bmap.warpimage(image=background, scale=0.25)
t = plt.title("Holocene Eruptions VEI > 4")

We can convey more information by coloring and sizing the triangles by the size of the eruption:

In [28]: 
def print_info(volcanos):
    print 'Name, Eruption Year, Tephra Volume'
    for i in xrange(volcanos.size):
        v = volcanos[i]
        textstr = v.volcano_name + ', ' + str(v.start_year) + ', ' + str(v.tephra_volume)
        print textstr
In [29]: 
vei4 = (gvp.vei == 4) # Pelee 1902, Eyjafjallaiokull 2010
vei5 = (gvp.vei == 5) # St Helens 1980m Vesuvius 79
vei6 = (gvp.vei == 6) # Pinatubo 1991, Krakatoa 1883, Huaynaputina 1600
vei7 = (gvp.vei == 7) # Tambora 1815, Mazama 5600B.C.
vei8 = (gvp.vei == 8) # Yellowstone 640,000 BC, Toba 74,000 BC, Taupo 24,500 BC

plt.title('Holocene Eruptions')
bmap.warpimage(image=background, scale=0.25)
v4 = plot_volcanos(gvp[vei4],10,'yellow')
v5 = plot_volcanos(gvp[vei5],16,'orange')
v6 = plot_volcanos(gvp[vei6],22,'red')
v7 = plot_volcanos(gvp[vei7],30,'purple')

plt.legend( [v4,v5,v6,v7], ['4','5','6','7'], 
            numpoints=1, 
            loc='lower center', 
            title='VEI', 
            framealpha=0.7,
            shadow=True)
plt.savefig('map.png', bbox_inches='tight')


print_info(gvp[vei7]) # print info on the largest eruptions

Name, Eruption Year, Tephra Volume
Changbaishan, 1000, 9.6x10^10
Crater Lake, -5677, 1.5x10^11
Kikai, -4350, 1.5x10^11
Kurile Lake, -6440, 1.6x10^11
Santorini, -1610, 9.9x10^10
Tambora, 1812, 1.6x10^11

Notice that eruptions concentrate on active plate boundaries. Despite having many of the world's active volcanos, the Andes Cordillera on the western margin of South America has no documented Holocene eruption of VEI>6. What was the largest eruption in the Andes?

In [30]: 
SA = gvp[gvp.region == 'South America']
biggest = (SA.vei == 6)
print_info(SA[biggest])

Name, Eruption Year, Tephra Volume
Huaynaputina, 1600, 3.0x10^10
Hudson, Cerro, -4750, >1.8x10^10
Hudson, Cerro, -1890, >1.0x10^10
Quilotoa, 1280, 2.1x10^10

Turns out it was Huaynaputina Volcano in Peru in 1600. You can check out Huaynaputina's page on the Smithsonian website here Let's see what else we can find out about Huaynaputina from this database:

In [31]: 
def pprint_records(recarray):
    for key in recarray.dtype.names:
        #print '{0}:\t\t {0}'.format(key, recarray[key])
        print '{0:25}{1}'.format(key, recarray[key]) #fix first entry to 25 spaces
In [32]: 
huayna = SA[SA.volcano_name == 'Huaynaputina']
pprint_records(huayna)

volcano_name             ['Huaynaputina']
volcano_number           [354030]
activity_type            ['Confirmed Eruption']
activity_id              [11795]
country                  ['Peru']
region                   ['South America']
subregion                ['Per\x9c']
latitude                 [-16.608]
longitude                [-70.85]
elevation                [4850]
start_year               [1600]
start_uncertainty_years  [-1]
start_month              [2]
start_day                [17]
start_uncertainty_days   [1]
end_year                 [1600]
end_uncertainty_years    [-1]
end_month                [3]
end_day                  [6]
end_uncertainty_days     [-1]
evidence_method_dating   ['Historical']
event_1                  ['Eruption (Explosive)']
event_2                  ['Ash']
event_3                  ['Pyroclastic flow']
event_4                  ['Property damage']
event_5                  ['Fatalities']
event_6                  ['Lahar or mudflow']
event_7                  ['Evacuations']
event_8                  ['Earthquakes (undefined)']
event_9                  ['Earthquakes (undefined)']
event_10                 ['Lightning']
vent_description         ['Central vent, Flank vent']
area_of_activity         ['Summit and south flank']
vei                      [6]
total_volume             ['\xca']
tephra_volume            ['3.0x10^10']
lava_volume              ['\xca']
fatalities_count         ['1,500?']
fatalities_remarks       [ 'About 1,500 persons were killed, excluding fatalities from lahars down the Rio Tambo or earthquakes in Arequipa.']
evacuations_count        ['\xca']
evacuations_remarks      ['\xca']

That's it for now! The iPython notebook is a great way to interactively explore data sets and simultaneously keep a copy of your work. If you wanted to share the notebook as a regular python script, or pdf (requires LaTeX) you can do that too pretty easily with the following ipython terminal commands:

ipython nbconvert --to python mapping_volcanic_eruptions.ipynb
ipython nbconvert --to latex --template article --post PDF mapping_volcanic_eruptions.ipynb

This post was written in an IPython notebook, which can be downloaded here, or viewed statically here.

Comments