This guide is part of Raster Alignment & Resampling Techniques, within the broader Automated Vector & Raster Cleaning Workflows reference.

Masking Rasters with Vector Geometries Using rasterio

Clipping a raster to a boundary is a two-line operation with four decisions hidden inside it, and three of them are silent when they go wrong.

Why Masking Produces Surprising Output

• A CRS mismatch returns an empty result, not an error. The geometry simply intersects nothing, and the output is all nodata.
• Pixel-centre semantics under-count small features. A three-pixel-wide watercourse can lose most of its pixels to the default rule.
• The transform must be updated. Masking without writing the returned transform produces a correctly-clipped array with the wrong georeferencing.
• Nodata has to exist. Masking sets excluded pixels to a fill value; if none is declared, zero is used and becomes indistinguishable from real data.

Version and Environment Compatibility

pip install "rasterio>=1.3" "geopandas>=1.0" "shapely>=2.0"

What all_touched Decides

The mask_to_geometry Recipe

from __future__ import annotations

import logging

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.mask import mask as rio_mask

logger = logging.getLogger(__name__)


def mask_to_geometry(
    raster_path: str,
    boundaries: gpd.GeoDataFrame,
    output_path: str,
    all_touched: bool = False,
    crop: bool = True,
    nodata: float | int | None = None,
) -> dict:
    """Clip a raster to vector geometry, writing a correctly-referenced output."""
    with rasterio.open(raster_path) as src:
        if boundaries.crs is None:
            raise ValueError("mask geometry has no CRS — cannot align it to the raster")

        # Reproject the geometry, never the raster: vertices are cheap, pixels are not.
        aligned = boundaries.to_crs(src.crs)
        geometry = [aligned.geometry.union_all()]

        minx, miny, maxx, maxy = aligned.total_bounds
        if (minx > src.bounds.right or maxx < src.bounds.left
                or miny > src.bounds.top or maxy < src.bounds.bottom):
            raise ValueError("mask geometry does not overlap the raster extent")

        fill = nodata if nodata is not None else src.nodata
        if fill is None:
            raise ValueError(f"{raster_path} declares no nodata; pass one explicitly")

        data, transform = rio_mask(
            src, geometry, crop=crop, all_touched=all_touched,
            nodata=fill, filled=True,
        )

        profile = src.profile.copy()
        profile.update(
            height=data.shape[1], width=data.shape[2],
            transform=transform, nodata=fill,
            compress="deflate", tiled=True, blockxsize=512, blockysize=512,
        )

    with rasterio.open(output_path, "w", **profile) as dst:
        dst.write(data)

    valid = int(np.count_nonzero(data[0] != fill))
    total = int(data[0].size)
    logger.info("%s: %d of %d pixels retained (%.1f%%)", output_path, valid, total,
                100 * valid / total if total else 0.0)
    return {"valid_pixels": valid, "total_pixels": total, "transform": transform}

Key Implementation Notes

• The geometry is reprojected, not the raster. Warping a raster to match a polygon resamples every pixel and changes the values being measured; moving the polygon changes nothing.
• union_all collapses a multi-feature mask. Passing many geometries works but rasterises each separately; the union is one pass and produces identical output.
• The overlap check runs before masking. An empty result from a non-overlapping mask is otherwise indistinguishable from a raster that is genuinely all nodata.
• A missing nodata raises rather than defaulting to zero. Zero is a legitimate elevation, reflectance and class code, and using it as a fill value makes masked and real pixels identical.
• The profile is updated with the returned transform. This is the step most often omitted, and it produces a file that opens, displays and is georeferenced to the wrong place.
• The output is tiled and compressed. A masked product is usually consumed remotely, so writing it in a cloud-friendly layout costs nothing here — the settings in writing cloud-optimized GeoTIFFs with rio-cogeo go further.

Masking a Raster Larger Than Memory

rasterio.mask.mask reads the masked window into memory, which is fine for an administrative area and not for a national mosaic clipped to a coastline. For those, mask by window.

Compute the geometry’s window with rasterio.features.geometry_window, iterate the raster’s internal blocks intersecting it, and mask each block against the geometry before writing it to the output. Peak memory then tracks the block size rather than the clipped extent, and the operation parallelises across blocks with no coordination.

The one subtlety is that each block must be masked against the same geometry rather than against a clipped copy of it, because clipping the polygon per block introduces edges at block boundaries that then appear in the output as seams — the same class of problem as the chunk-boundary artefacts described in handling precision and coordinate rounding.

Troubleshooting Raster Masks

Integration Note

Masking belongs after alignment and before any statistic is computed, so that every measurement covers exactly the area of interest — the stage ordering in raster alignment and resampling techniques. Record the mask geometry’s source and version alongside the output, because a boundary that changes between runs shifts every derived figure, and that is otherwise indistinguishable from a change in the data itself.

Masking Versus Zonal Statistics

Clipping a raster to a boundary and computing a statistic over that boundary look like the same operation and are not, and choosing the wrong one costs either accuracy or a great deal of storage.

Mask when the clipped raster is the product. A tile handed to an analyst, an input to a model that expects a rectangular array, an image published for a specific area — these need pixels written to a file, and the extent matters.

Use zonal statistics when only the number is wanted. Mean elevation per catchment, total built area per district, land-cover proportions per parcel: these read the pixels within each zone and emit one row per zone. Nothing is written per zone, which for ten thousand small polygons is the difference between ten thousand files and one table.

The mistake worth avoiding is masking per polygon in a loop to compute a statistic. It produces the right answers, writes or holds one array per feature, and runs one to two orders of magnitude slower than a single pass that rasterises the zones once and aggregates by zone identifier. Where a library such as rasterstats or an xarray grouping is available, prefer it, and reserve explicit masking for the cases where the raster itself is the deliverable.

Related

• Raster Alignment & Resampling Techniques
• Reprojecting Rasters to a Common Grid with rasterio.warp
• Writing Cloud-Optimized GeoTIFFs with rio-cogeo
• Choosing Resampling Methods for Categorical and Continuous Rasters