How to Read, Transform, and Write Rasters with Rasterio

Rasterio is a Python library that simplifies geospatial raster data processing by wrapping the Geospatial Data Abstraction Library (GDAL) in idiomatic Python interfaces. This article explains how Rasterio leverages GDAL's C/C++ engine to read spatial raster bands into NumPy arrays, apply geometric and array-level transformations, and write the modified georeferenced datasets to disk while preserving critical spatial metadata.

The GDAL and Rasterio Connection

GDAL is the industry standard for translating and manipulating geospatial raster formats such as GeoTIFF, HDF, and NetCDF. While GDAL includes native Python bindings, they closely mirror its C++ implementation, leading to verbose and non-Pythonic code. Rasterio solves this by using Cython to interface directly with GDAL’s C API. It handles GDAL data structures behind the scenes, mapping GDAL raster bands directly to standard NumPy N-dimensional arrays and organizing spatial metadata—such as Coordinate Reference Systems (CRS) and affine transform matrices—into intuitive Python dictionaries.

Reading Geospatial Raster Bands

Rasterio opens datasets using rasterio.open(). When a file is opened, GDAL reads the file header without immediately loading the entire pixel dataset into memory. This lazy-loading model allows you to inspect metadata before processing.

import rasterio

# Open the georeferenced file
with rasterio.open("input.tif") as src:
    # Access metadata
    print("Dimensions:", src.width, src.height)
    print("Band count:", src.count)
    print("CRS:", src.crs)
    print("Affine Transform:", src.transform)
    
    # Read the first band into a NumPy array (1-indexed)
    band1 = src.read(1)
    
    # Read all bands as a 3D NumPy array (bands, rows, columns)
    all_bands = src.read()

When calling read(), GDAL translates the raw raster values into a NumPy array. Rasterio also supports windowed reading via rasterio.windows.Window, which permits loading specific bounding boxes or pixel sub-regions directly, preventing memory exhaustion when handling massive datasets.

Transforming Raster Data

Transformations in Rasterio fall into two primary categories: array-level (radiometric) manipulation and coordinate-level (spatial) transformations.

Array-Level Manipulations

Because raster bands are loaded as NumPy arrays, mathematical operations on pixel values are straightforward:

import numpy as np

# Example: Calculating Normalized Difference Vegetation Index (NDVI)
with rasterio.open("multispectral.tif") as src:
    red = src.read(1).astype(float)
    nir = src.read(2).astype(float)
    
    # Avoid division by zero
    ndvi = np.where((nir + red) == 0, 0, (nir - red) / (nir + red))

Spatial Transformations and Reprojection

Spatial transformations alter the raster's coordinate system, pixel resolution, or geographic alignment. Rasterio handles this using rasterio.warp.reproject and GDAL's internal warping algorithms:

from rasterio.warp import calculate_default_transform, reproject, Resampling

dst_crs = "EPSG:4326"

with rasterio.open("input.tif") as src:
    # Calculate target transform and dimensions for the new CRS
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds
    )
    
    # Allocate destination array
    destination = np.zeros((src.count, height, width), dtype=src.dtypes[0])
    
    # Reproject using bilinear interpolation
    reproject(
        source=rasterio.band(src, list(range(1, src.count + 1))),
        destination=destination,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=transform,
        dst_crs=dst_crs,
        resampling=Resampling.bilinear
    )

GDAL handles the complex coordinate transformations, while Rasterio ensures the resulting array aligns with the newly calculated affine transform.

Writing Geospatial Raster Bands

Writing raster data back to disk requires both the pixel arrays and an accurate geospatial profile. The profile specifies the driver (e.g., 'GTiff'), data dimensions, data type, band count, CRS, and affine transform matrix.

with rasterio.open("input.tif") as src:
    # Copy and update the source file's metadata profile
    profile = src.profile
    profile.update(
        dtype=rasterio.float32,
        count=1,
        compress="lzw"
    )

# Write the processed array to a new GeoTIFF
with rasterio.open("ndvi_output.tif", "w", **profile) as dst:
    dst.write(ndvi.astype(rasterio.float32), 1)

When calling rasterio.open() in 'w' mode, Rasterio passes the configuration options to GDAL's file creation drivers. GDAL allocates the dataset container, compresses the data (if specified, such as using LZW or Deflate), builds the header containing the geographic coordinate reference, and writes the array values into their corresponding raster bands. Opening files within Python context managers (with blocks) ensures that GDAL flushes all data buffers to disk and properly closes the file handles upon completion.