We'd like to visualize some DEM data from Mt St Helens. See this blog post by Evan Bianco.
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
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.
after.shape
(1400, 979)
after.ptp()
41134.0
print np.amin(before), np.amax(before)
print np.amin(after), np.amax(after)
-32767.0 2951.0 -32767.0 8367.0
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)
nan nan nan nan
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.
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
(468, 327) (1400, 979)
The grids are different sizes. We need to resample.
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
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
difference = before_re - after_re
plt.imshow(difference, cmap='gist_earth')
plt.colorbar()
plt.show()