Rasters in Python with rasterio

Every GeoDataFrame in the last few lessons stored discrete things: a point, a line, a polygon, each with its own boundary. A satellite image has no boundaries to draw. It is a grid of numbers, one value per cell, covering a rectangle of the Earth at some fixed resolution. GeoPandas has no idea what to do with that shape of data, because a grid is not a table of geometries. The library built for grids is called rasterio, and it treats a raster file almost exactly the way NumPy treats an array, because under the surface, that is exactly what it is.

Opening a file without loading it

A raster file, usually a GeoTIFF with a .tif extension, can be a few kilobytes or several gigabytes. rasterio does not read the whole thing into memory just to answer "how big is this and what CRS is it in." Opening a dataset with rasterio.open() inside a with block gives a handle to the file's metadata immediately, and pixel data only gets pulled in when something explicitly asks to read it.

import rasterio

with rasterio.open("elevation.tif") as dataset:
    print(dataset.count, dataset.width, dataset.height)
    print(dataset.crs)

dataset.count is the number of bands stored in the file. A single elevation raster usually has one band, a plain grayscale value per cell. A satellite scene might carry four or more: red, green, blue, and near-infrared, each its own band stacked in the same file. dataset.width and dataset.height are the grid's size in pixels, columns and rows, not meters or degrees. dataset.crs reports the coordinate reference system the same way GeoPandas does, and it matters just as much here. A raster in EPSG:4326 has cells that are a fixed size in degrees of longitude and latitude, which are not the same physical distance everywhere on Earth. A raster reprojected into a local UTM zone has cells that are a fixed size in meters, which is usually what an analyst wants before measuring area or distance on the grid.

Reading a band pulls the pixel values out as a NumPy array, and the whole point of rasterio is that this array behaves like any other NumPy array from there on.

import rasterio

with rasterio.open("elevation.tif") as dataset:
    band1 = dataset.read(1)
    print(type(band1), band1.shape, band1.dtype)
    print(band1.min(), band1.max())

Band numbering starts at 1, not 0, because that is the convention GDAL uses underneath rasterio, and it is one of the few places rasterio breaks from Python's usual zero-based counting. The returned array has shape (height, width), rows first, which matches how NumPy indexes any 2D array: band1[row, col]. The min() and max() calls are plain NumPy methods, nothing raster-specific about them. This is the payoff of the design. Learning NumPy once means an analyst already knows how to slice a raster, mask it, average it, or reshape it, because a band is just an array like any other.

Two more properties round out what a dataset reports before any pixels move. dataset.transform is an affine transform, a small set of six numbers that maps a pixel's row and column to real-world coordinates. Without it, an array of elevation numbers has no idea where on Earth it sits, only that it is 3601 by 3601 cells. The transform is what lets rasterio answer "what coordinate is pixel (100, 250)" or the reverse, "which pixel covers this coordinate." dataset.nodata is the value, often something like -9999 or a masked value, that a sensor or a processing pipeline uses to mark a cell with no real reading: a cloud-covered pixel, a gap in a scan, an area outside a sensor's swath. An analyst has to know which value means "nothing here" before running any math on the grid, because a stray -9999 silently averaged into real elevation numbers wrecks a result quietly, with no error and no warning.

Reading a piece instead of the whole file

A satellite scene covering a whole state can run into the tens of gigabytes. Pulling that entire array into memory just to look at one county is wasteful, and on a laptop it may simply fail. rasterio can read a rectangular chunk of a file, called a window, without ever touching the bytes outside it. A window is defined in pixel coordinates: a column offset, a row offset, a width, and a height.

import rasterio
from rasterio.windows import Window

with rasterio.open("scene.tif") as dataset:
    chunk = dataset.read(1, window=Window(col_off=1000, row_off=500, width=256, height=256))
    print(chunk.shape)

This is not a convenience feature bolted on top of a full read. GeoTIFF files store their pixels in internal tiles, small blocks the format is built around, and a windowed read asks the underlying GDAL library for only the tiles that overlap the requested rectangle. A file built as a Cloud Optimized GeoTIFF, usually shortened to COG, takes this further by organizing those tiles with an index at the front of the file and by storing lower-resolution overviews alongside the full-resolution data. A windowed read against a COG sitting on a remote server, an S3 bucket for instance, pulls only the relevant bytes over the network instead of downloading the entire multi-gigabyte file first. The same window= argument works whether the file is local or remote, so nothing about the code changes, only how much data crosses the wire.

Sometimes the question is not "give me this rectangle" but "what is the elevation at this exact spot." dataset.sample() takes a list of (x, y) coordinates in the dataset's CRS and, for each one, returns the pixel values at whichever cell contains that point. It is a generator, so wrapping it in list() pulls out the actual results.

import rasterio

with rasterio.open("elevation.tif") as dataset:
    points = [(-105.2705, 40.0150), (-105.30, 40.02)]
    values = list(dataset.sample(points))
    print(values)

Each result is a small array, one entry per band, because sample() was built to work on multi-band files too. For a single-band elevation raster that array holds one number: the elevation at that coordinate. This is the pattern behind a huge amount of real work, checking a rainfall grid at a set of weather stations, pulling an elevation at every point along a hiking trail, matching a satellite reading to a list of field survey locations. The coordinates supplied to sample() must be in the same CRS as the raster, the same rule that governs every other spatial operation covered so far, and a mismatch there returns numbers that look fine and point at the wrong ground entirely.

Doing math on the array

Because a band is just a NumPy array, ordinary array arithmetic is the whole toolkit for combining bands. The most common example in remote sensing is the Normalized Difference Vegetation Index, or NDVI, a simple formula that highlights healthy plant growth using two bands most satellites already carry: red light and near-infrared light. Healthy vegetation reflects a lot of near-infrared and absorbs most red, so the difference between the two, scaled by their sum, comes out high wherever plants are thriving and low or negative over bare ground, water, or pavement.

import numpy as np
import rasterio

with rasterio.open("scene.tif") as dataset:
    red = dataset.read(3).astype("float32")
    nir = dataset.read(4).astype("float32")

with np.errstate(divide="ignore", invalid="ignore"):
    ndvi = (nir - red) / (nir + red)

print(ndvi.shape, ndvi.dtype)
print(np.nanmin(ndvi), np.nanmax(ndvi))

Two things in that code are doing real work beyond the formula itself. Both bands get cast to float32 before the subtraction, because the raw pixel values usually arrive as integers, and integer division in a formula like this would truncate results to whole numbers, throwing away the entire point of the index. The np.errstate block quiets a warning that fires whenever nir + red lands on exactly zero, which happens at nodata pixels or masked edges of a scene, producing a division by zero. Silencing the warning does not fix the bad pixel, it just stops it from cluttering the output. A careful version of this calculation would mask out dataset.nodata values on both bands before doing the division, so the zero-division cells never make it into a plotted map. That is one line of extra code, np.where(red == dataset.nodata, np.nan, red) on each band, and it is the difference between a clean NDVI map and one with speckled garbage values along its edges.

This is the shape of nearly all raster analysis. Read the needed bands as arrays, apply array math, and the result is another array in the same grid, ready to be written back out to a new GeoTIFF with rasterio.open(path, "w", **profile) using the same transform and CRS as the source. The geometry model changed completely between the last few lessons and this one, points and polygons on one side, a fixed grid on the other, but the tool for exploring it turned out to be the same array library used for anything else in Python.