Rasterio is to rasters as fiona is to shapefiles. Following this example.
In [1]:
import numpy as np
import rasterio
In [2]:
with rasterio.drivers(CPL_DEBUG=True):
with rasterio.open('data/Image2_b_TMIresidual.tif') as src:
r, g, b = src.read()
In [3]:
r.shape
Out[3]:
In [4]:
import matplotlib.pyplot as plt
%matplotlib inline
In [5]:
plt.imshow(r, cmap="Greys")
Out[5]:
Let's try plotting a section across here. Thanks Joe Kington for the tips!
In [6]:
x0, y0 = 500,500
x1, y1 = 1500,1500
length = int(np.hypot(x1-x0, y1-y0))
x, y = np.linspace(x0, x1, length), np.linspace(y0, y1, length)
zi = r[x.astype(np.int), y.astype(np.int)]
fig, axes = plt.subplots(ncols=2)
axes[0].imshow(r)
axes[0].plot([x0, x1], [y0, y1], 'ro-')
axes[0].axis('image')
axes[1].plot(zi)
plt.show()
Make a fake elevation map.
In [7]:
u, v = np.mgrid[-1000:893, -1000:1318]
elev = np.sin((u/100.) + (v/100.))
In [8]:
elevi = elev[x.astype(np.int), y.astype(np.int)]
fig, axes = plt.subplots(ncols=2)
# Plot the elevation map.
axes[0].imshow(elev)
axes[0].plot([x0, x1], [y0, y1], 'ro-')
axes[0].axis('image')
# Plot the extracted profile.
axes[1].plot(elevi)
plt.show()
Make a 'transect x', the x-coordinate along the transect:
In [9]:
transx = np.linspace(0, length, elevi.size)
In [10]:
plt.scatter(transx, elevi, c=zi, edgecolor='', s=2)
Out[10]:
In [10]: