Chapter 8: Raster Data with NumPy and rasterio¶
Part II — Geographic Data Science
Learning Objectives¶
By the end of this chapter you will be able to:
- Read and write GeoTIFF files with rasterio
- Apply raster algebra and focal operations with NumPy
- Understand the NumPy array as a simulation substrate
# Standard imports — add chapter-specific imports below
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import rasterio as rio
from rasterio.plot import show
Maranhão's coastline, the running example since Chapter 6, is almost perfectly flat — a poor advertisement for a chapter about elevation. This chapter borrows terrain from a hillier part of Brazil instead: an SRTM (Shuttle Radar Topography Mission) elevation grid over a microregion in Minas Gerais, centered on the town of Conselheiro Lafaiete. Chapter 16 brings the coast back once raster and vector start working together.
Raster Data Model¶
A raster's defining structure is simple: a regular grid of rows and columns, one or more bands of values stacked on top of it (a single elevation band here; Chapter 9's climate data stacks dozens of time steps), and an affine transform — six numbers that map "row 340, column 12" to a real-world coordinate, so the grid can be placed on the Earth without storing a coordinate for every single cell.
Open a GeoTIFF with rasterio.open() and that structure is immediately inspectable, before reading a single pixel of actual data:
raster_path = "srtm_micro_lafaiete.tif"
src = rio.open(raster_path)
print("Bands:", src.indexes)
print("Data type:", src.dtypes)
print("Shape (rows, cols):", src.shape)
print("CRS:", src.crs)
print("Transform:", src.transform)
Bands: (1,)
Data type: ('int16',)
Shape (rows, cols): (1895, 3671)
CRS: EPSG:4326
Transform: | 0.00, 0.00,-44.42|
| 0.00,-0.00,-20.42|
| 0.00, 0.00, 1.00|
That transform is the raster analogue of the CRS work Chapter 6 did for vector data — it's what lets rasterio answer "where is pixel (340, 12)" the same way a GeoDataFrame's geometry column answers "where is this polygon."
Reading GeoTIFFs with rasterio¶
Metadata alone doesn't get you a NumPy array — that needs an explicit .read(), naming which band (bands are 1-indexed, not 0-indexed, one of rasterio's few surprises for anyone coming from NumPy):
elevation = src.read(1)
print(type(elevation), elevation.shape, elevation.dtype)
fig, ax = plt.subplots(figsize=(6, 6))
show(elevation, cmap="terrain", transform=src.transform, ax=ax)
ax.set_title("SRTM elevation — Conselheiro Lafaiete microregion")
plt.show()
<class 'numpy.ndarray'> (1895, 3671) int16
Real rasters almost always have a nodata value — a sentinel marking "no valid measurement here" (ocean, a sensor gap, a clipped edge), and it is not automatically NaN. SRTM tiles commonly use 0 for missing elevation, which is a genuinely bad default: 0 is also a perfectly plausible elevation at sea level, so a naive read can't tell "no data" apart from "this pixel is at sea level" without being told which is which:
print("nodata as read:", src.nodata)
# Read with the mask applied — returns a masked array instead of raw numbers
masked_elevation = src.read(1, masked=True)
print(type(masked_elevation))
print("Valid pixels:", masked_elevation.count(), "of", masked_elevation.size)
nodata as read: None <class 'numpy.ma.MaskedArray'> Valid pixels: 6956545 of 6956545
If src.nodata came back None even though you know certain pixels are placeholders, the fix is the same pattern Chapter 4 used on messy tabular data — replace the sentinel with NaN explicitly, this time with np.where instead of pd.to_numeric(errors="coerce"):
clean_elevation = np.where(elevation == 0, np.nan, elevation)
fig, ax = plt.subplots(figsize=(6, 6))
show(clean_elevation, cmap="terrain", transform=src.transform, ax=ax)
ax.set_title("After masking the 0-elevation sentinel")
plt.show()
Raster Algebra¶
Once a raster is a NumPy array, every vectorized trick from Chapter 2 and Chapter 3 applies directly — arithmetic on an entire grid at once, no pixel-by-pixel loop required. A hillshade-style relief exaggeration is one line:
exaggerated = clean_elevation * 2.5 # exaggerate relief for visualization
diff_from_mean = clean_elevation - np.nanmean(clean_elevation)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].imshow(exaggerated, cmap="terrain")
axes[0].set_title("Elevation × 2.5")
axes[1].imshow(diff_from_mean, cmap="RdBu_r")
axes[1].set_title("Deviation from mean elevation")
plt.show()
Boolean masks work exactly as they did on a pandas column in Chapter 4 — clean_elevation > 900 returns a same-shaped array of True/False, directly usable both as a mask and, cast to int, as its own raster:
highlands = clean_elevation > 900
print("Highland pixels:", np.nansum(highlands), "of", highlands.size)
plt.imshow(highlands, cmap="gray")
plt.title("Terrain above 900 m")
plt.show()
Highland pixels: 2865817 of 6956545
Focal Operations with NumPy¶
Raster algebra combines a cell only with itself, across layers. A focal (or moving-window) operation instead combines each cell with its neighbors within the same layer — a smoothing filter, a slope estimate, anything where the output at one pixel depends on a small neighborhood around it. scipy.ndimage implements the common ones efficiently, entirely in terms of arrays rasterio and NumPy already produced — with one catch: a moving-window mean naively applied to clean_elevation would let a single NaN neighbor (from the nodata masking two sections ago) contaminate every window that touches it, silently blanking out entire regions near any nodata edge. Splitting the mean into a sum-of-valid-values and a count-of-valid-values, each computed with a fast uniform_filter pass, avoids that: a window with 20 valid neighbors and 5 nodata ones still averages correctly over the 20, instead of returning NaN for the whole window. The example below uses a 15×15 window and zooms into the single most rugged patch of the raster — a small window, or a flat patch of terrain, both make a moving-window mean's effect genuinely hard to see by eye even when it's working correctly, since a shared, wide-range colormap like terrain compresses a modest local smoothing into a barely-visible shade difference. A window near a large nodata patch might still end up with only one or two valid neighbors, though — technically a defined average, but really just that one raw value wearing a smoothed-looking label. Discarding windows below a minimum valid fraction (min_valid_fraction) keeps only cells with enough real neighbors to make the average meaningful, at the cost of a slightly wider NaN border around every nodata patch.
from scipy.ndimage import uniform_filter
valid = ~np.isnan(clean_elevation)
filled = np.where(valid, clean_elevation, 0)
window = 15
min_valid_fraction = 0.5 # require at least half the window to be real data
sum_filter = uniform_filter(filled, size=window) * window**2
count_filter = uniform_filter(valid.astype(float), size=window) * window**2
with np.errstate(invalid="ignore", divide="ignore"):
smoothed = sum_filter / count_filter
smoothed[count_filter < min_valid_fraction * window**2] = np.nan
# Zoom into the most rugged patch, where a 15x15 mean has the most visible effect.
# Local variance via uniform_filter (fast, vectorized) instead of generic_filter
# (which calls a Python function per pixel and is far too slow on a full raster):
# Var(X) = E[X²] - E[X]²
rugged_window = 30
mean_ = uniform_filter(filled, size=rugged_window)
mean_sq = uniform_filter(filled**2, size=rugged_window)
local_std = np.sqrt(np.clip(mean_sq - mean_**2, 0, None))
cy, cx = np.unravel_index(np.argmax(local_std), local_std.shape)
half = 40
crop = (slice(max(cy - half, 0), cy + half), slice(max(cx - half, 0), cx + half))
vmin, vmax = np.nanmin(clean_elevation[crop]), np.nanmax(clean_elevation[crop])
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].imshow(clean_elevation[crop], cmap="terrain", vmin=vmin, vmax=vmax)
axes[0].set_title("Original (rugged crop)")
im = axes[1].imshow(smoothed[crop], cmap="terrain", vmin=vmin, vmax=vmax)
axes[1].set_title(f"{window}×{window} moving-window mean (rugged crop)")
fig.colorbar(im, ax=axes, shrink=0.7, label="Elevation (m)")
plt.show()
Cropping to the most rugged patch and sharing a color scale between both panels already makes the effect visible — but a map is still an indirect way to prove a filter is doing what it claims. A cross-section: a plain line plot of elevation along a single row of pixels, original against smoothed, removes the colormap from the question entirely.
row = cy # same rugged row used to find the crop above
cols = range(crop[1].start, crop[1].stop)
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(cols, clean_elevation[row, crop[1]], label="Original", alpha=0.7)
ax.plot(cols, smoothed[row, crop[1]], label=f"{window}×{window} smoothed", linewidth=2)
ax.set_xlabel("Column (pixel)")
ax.set_ylabel("Elevation (m)")
ax.legend()
ax.set_title(f"Elevation profile — row {row}")
plt.show()
The smoothed line should sit visibly rounder than the original: peaks shaved down, valleys filled in, exactly what averaging every pixel with its neighbors is supposed to do. That's a useful lesson on its own, independent of elevation data specifically: a moving-window operation generalizes directly to anything Part III's cellular automata need — a fire-spread model checking whether any neighbor is burning, a diffusion model averaging a concentration across neighbors each step — and getting the visualization right matters just as much as getting the array math right, since a silently-wrong plot and a silently-too-subtle-to-see-correct plot look identical at first glance. Chapter 18 reaches for this same neighborhood logic, just with rules instead of a mean.
Multi-band Imagery¶
Every array so far has had one band — elevation, a single number per pixel. Satellite imagery is usually stacked: a Landsat scene ships red, green, blue, near-infrared, and more as separate bands sharing the same grid and transform. rasterio reads a multi-band file into a 3D array, shaped (bands, rows, cols):
# with a multi-band file open as `src`:
# stack = src.read() # shape: (n_bands, rows, cols)
# red, nir = stack[2], stack[3] # band order depends on the sensor
# ndvi = (nir - red) / (nir + red) # a common per-pixel band algebra
That ndvi line is raster algebra again — Chapter 9 is where genuinely large, many-band stacks like this get a data structure built specifically for them (xarray, with named, labeled bands instead of a bare integer index).
NumPy Array as a Simulation Grid¶
Chapter 7 closed by treating a GeoDataFrame as a simulation substrate, one row per cell. A NumPy array is the same idea in raster form — and, for a regular grid, usually the faster one: no geometry objects to intersect, just array indexing.
DisSModel's RasterBackend, which Chapter 8's counterpart in Part III (Chapter 21) introduces properly, wraps exactly this: a 2D (or 3D, for multiple state variables) NumPy array as a model's state, with each cell's eight neighbors reachable by simple index arithmetic — grid[i-1:i+2, j-1:j+2] — rather than a spatial join. The focal operation above is a preview of the mechanism every cellular automaton in Part III depends on: rule() reads a cell's neighborhood, computes a new value, and the whole grid updates via one vectorized NumPy expression, not a Python loop over cells. Chapter 20 comes back to explain exactly why that vectorization step matters as a grid grows from Lafaiete's few hundred pixels to a continental one.
Exercises¶
- Nodata, correctly. Confirm that
src.nodataand the0-sentinel behave differently: count how many pixels equal exactly0inelevation, and compare that toelevation.size - masked_elevation.count(). Are they the same number? Should they be? - A different focal size. Rerun the Focal Operations smoothing with
window=5instead ofwindow=15, keeping the same rugged crop and shared color scale. Is the smoothing effect still visible? At what window size does it stop being visible in the cropped, shared-scale plot — and does that match whatnp.nanmax(np.abs(smoothed - clean_elevation))reports? - Threshold and count. Using
clean_elevation, find the elevation threshold above which exactly 10% of valid pixels fall (hint:np.nanpercentile). How does that threshold compare to the900m used in Raster Algebra? - Sketch the neighbor math. Without running any code, write out in words what
grid[i-1:i+2, j-1:j+2]selects for a cell at position(i, j), and how many cells that is in total. This is the exact slice Chapter 18's cellular automata use for a Moore neighborhood.
# Your code here
Summary¶
Key concepts introduced¶
- The raster data model: a regular grid, one or more bands, and an affine transform placing it in real-world coordinates
- Reading a GeoTIFF with
rasterio.open()and.read(), and why nodata sentinels (like SRTM's0) need explicit handling, not automatic detection - Raster algebra: NumPy arithmetic and boolean masks applied to a whole grid at once, carrying Chapter 2's vectorization habit forward
- Focal (moving-window) operations with
scipy.ndimage, as a preview of the neighborhood logic Part III's cellular automata run on every step - Multi-band imagery as a 3D array, and a first look at per-pixel band algebra (NDVI)
- The NumPy array as a simulation substrate — DisSModel's
RasterBackend, the raster counterpart to Chapter 7'sGeoDataFrame-as-substrate
Chapter 9 stays in raster territory but swaps this chapter's single 2D grid for a full multidimensional, labeled stack — the shape real climate and land-cover time series actually come in.
Further Reading¶
- rasterio documentation, Reading and writing data: https://rasterio.readthedocs.io/en/stable/topics/reading.html
- rasterio documentation, Nodata masking: https://rasterio.readthedocs.io/en/stable/topics/masking-by-shapefile.html
- NASA JPL, Shuttle Radar Topography Mission (SRTM): https://www2.jpl.nasa.gov/srtm/
- scipy documentation,
scipy.ndimage— the module behind this chapter's focal filter: https://docs.scipy.org/doc/scipy/reference/ndimage.html