from IPython.display import Image
url = 'http://upload.wikimedia.org/wikipedia/commons/d/dc/MSH82_st_helens_plume_from_harrys_ridge_05-19-82.jpg'
Image(url=url, width=400)

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

before = np.loadtxt('data/st-helens_before.txt')
after = np.loadtxt('data/st-helens_after.txt')

after

after.shape

after.ptp()

print np.amin(before), np.amax(before)
print np.amin(after), np.amax(after)

plt.imshow(before, cmap="gist_earth")
plt.colorbar()
plt.show()

before[before==-32767.] = np.nan
after[after==-32767.] = np.nan

plt.imshow(before, cmap="gist_earth")
plt.colorbar()
plt.show()

print np.amin(before), np.amax(before)
print np.amin(after), np.amax(after)

print np.nanmin(before), np.nanmax(before)
print np.nanmin(after), np.nanmax(after)

from mpl_toolkits.mplot3d.axes3d import Axes3D

x = np.arange(before.shape[1])
y = np.arange(before.shape[0])
x, y = np.meshgrid(x, y)

fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(1, 1, 1, projection='3d')
p = ax.plot_surface(x, y, before, cmap="gist_earth", linewidth=0, vmin=676, vmax=2951)

after /= 3.28084 # convert to metres
x = np.arange(after.shape[1])
y = np.arange(after.shape[0])
x, y = np.meshgrid(x, y)

fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(1, 1, 1, projection='3d')
p = ax.plot_surface(x, y, after, cmap="gist_earth", linewidth=0, vmin=676, vmax=2951)

print before.shape, after.shape

import scipy.ndimage
before_re = scipy.ndimage.zoom(before, 3, order=0)
print before_re.shape

before_re = before_re[2:-2,1:-1]
print before_re.shape

after_re = after
print after_re.shape

difference = before_re - after_re

plt.imshow(difference, cmap='gist_earth')
plt.colorbar()
plt.show()

x = np.arange(difference.shape[1])
y = np.arange(difference.shape[0])
x, y = np.meshgrid(x, y)

fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.plot_surface(x, y, difference, cmap="gist_earth", linewidth=0, vmin=-50, vmax=1000 )
plt.show()

volume = 100 * np.nansum(difference)
print "Approximate expelled volume: {0} m³ or {1:.3} km³".format(int(round(volume,-5)), volume/1e9)