We'd like to visualize some DEM data from Mt St Helens. See this blog post by Evan Bianco.
In [1]:
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)
Out[1]:
In [2]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [19]:
before = np.loadtxt('data/st-helens_before.txt')
after = np.loadtxt('data/st-helens_after.txt')
In [20]:
after
Out[20]:
Those weird values are NaNs. They will cause trouble, but we will deal.
In [21]:
after.shape
Out[21]:
In [22]:
after.ptp()
Out[22]:
In [23]:
print np.amin(before), np.amax(before)
print np.amin(after), np.amax(after)
In [24]:
plt.imshow(before, cmap="gist_earth")
plt.colorbar()
plt.show()
In [25]:
before[before==-32767.] = np.nan
after[after==-32767.] = np.nan
plt.imshow(before, cmap="gist_earth")
plt.colorbar()
plt.show()
In [26]:
print np.amin(before), np.amax(before)
print np.amin(after), np.amax(after)
In [27]:
print np.nanmin(before), np.nanmax(before)
print np.nanmin(after), np.nanmax(after)
Seems like before is in metres, while after is in feet. Of course.
In [12]:
from mpl_toolkits.mplot3d.axes3d import Axes3D
In [42]:
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)
In [43]:
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)
In [44]:
print before.shape, after.shape
The grids are different sizes. We need to resample.
In [45]:
import scipy.ndimage
before_re = scipy.ndimage.zoom(before, 3, order=0)
print before_re.shape
Remove 2 elements from top and bottom of first dimension, and 1 element from second
In [46]:
before_re = before_re[2:-2,1:-1]
print before_re.shape
after_re = after
print after_re.shape
The "after" data remains unchanged
Now this matches the shape of the "after" array, so we can do math on them
In [49]:
difference = before_re - after_re
In [50]:
plt.imshow(difference, cmap='gist_earth')
plt.colorbar()
plt.show()
In [51]:
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()
In [77]:
volume = 100 * np.nansum(difference)
print "Approximate expelled volume: {0} m³ or {1:.3} km³".format(int(round(volume,-5)), volume/1e9)
Mayavi doesn't perform well in the notebook so it is omitted. Try it on your machine... Installation instructions.
Note that VTK is a prerequisite, and installing it can be a bit involved.