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

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')
after = np.loadtxt('data/st-helens_after.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()