## Introduction¶

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
Image(url=url, width=400)

Out[1]:

## Import libraries and data¶

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

In [19]:
before = np.loadtxt('data/st-helens_before.txt')

In [20]:
after

Out[20]:
array([[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
...,
[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
[-32767., -32767., -32767., ..., -32767., -32767., -32767.]])

Those weird values are NaNs. They will cause trouble, but we will deal.

In [21]:
after.shape

Out[21]:
(1400, 979)
In [22]:
after.ptp()

Out[22]:
41134.0
In [23]:
print np.amin(before), np.amax(before)
print np.amin(after), np.amax(after)

-32767.0 2951.0
-32767.0 8367.0


## Plotting and cleaning data¶

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)

nan nan
nan nan

In [27]:
print np.nanmin(before), np.nanmax(before)
print np.nanmin(after), np.nanmax(after)

676.0 2951.0
2231.0 8367.0


Seems like before is in metres, while after is in feet. Of course.

## 3D plotting in matplotlib¶

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

(468, 327) (1400, 979)


The grids are different sizes. We need to resample.

## Resampling¶

In [45]:
import scipy.ndimage
before_re = scipy.ndimage.zoom(before, 3, order=0)
print before_re.shape

(1404, 981)


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

(1400, 979)
(1400, 979)


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)

Approximate expelled volume: 2004700000 m³ or 2.0 km³


## 3D views in Mayavi¶

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.