%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)

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])


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

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')

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()

bmap.shadedrelief(scale=0.25)

background = '/Volumes/OptiHDD/data/natural_earth/medium_res/GRAY_50M_SR_OB/GRAY_50M_SR_OB.jpg'
bmap.warpimage(image=background, scale=0.25)

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")

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

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

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

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

huayna = SA[SA.volcano_name == 'Huaynaputina']
pprint_records(huayna)