Masking Rasters with Polygons in Python

Cutting a raster down to a polygon looks like a one-line call, and it is — right up to the moment the output comes back the same size as the input, or entirely nodata, or written to disk with a transform that still points at the original corner. This guide is for anyone extracting a study area from a scene with rasterio.mask.mask: it sits under Zonal Statistics & Raster Sampling in Spatial Analysis & Advanced Query Techniques, and it assumes you can already open a dataset and read a band, as covered in Raster Data Handling with Rasterio.

Why This Approach / What Goes Wrong

rasterio.mask.mask fuses two operations that people habitually confuse. The first is masking: burning the geometry into a boolean array the same shape as the pixels, then hiding every pixel on the wrong side of it. The second is cropping: narrowing the read window to the geometry's bounding box so the array you get back is smaller than the file. They are controlled by different arguments, they fail in different ways, and only one of them changes the affine transform.

With the default crop=False the returned array has exactly the source shape, out_transform is byte-for-byte src.transform, and every pixel outside the polygon has been replaced by the nodata value. You have a clipped picture inside an unclipped frame — convenient when several masks must stay pixel-aligned with each other, wasteful when the polygon covers two percent of a 640-by-400 grid, because the whole band was decoded to produce it. With crop=True the window shrinks to the geometry's bounds rounded outward to whole pixels, so out_image has new height and width and out_transform carries a new origin. The pixel size terms are untouched; only the translation terms move. That distinction is the whole reason the function returns a transform at all: it is the piece of metadata you must carry into the output profile, and forgetting it is how a perfectly good array ends up georeferenced somewhere it never was.

The second surprise is that nothing in the call inspects a coordinate reference system. mask receives an open dataset and an iterable of GeoJSON-like mappings; geometries arriving through __geo_interface__ have already lost their CRS, so the function compares raw numbers against raw numbers. Hand it a floodplain in degrees and a scene in metres and the geometry's bounding box lands nowhere near the grid. What happens next depends on crop: with crop=True you get ValueError: Input shapes do not overlap raster., which at least stops the pipeline, and with the default crop=False you get a UserWarning about shapes being outside the bounds of the raster, a mask that is True everywhere, and an output array in which every pixel is the nodata fill. A warning in a log nobody reads, followed by a plausible-looking GeoTIFF full of -9999, is the single most common way this function fails in production.

Prerequisites

• rasterio>=1.3,<2 — provides rasterio.mask.mask, raster_geometry_mask, and the pad_width argument
• geopandas>=1.0 — reads the polygon layer and carries its CRS
• shapely>=2.0 — the geometry objects, each exposing __geo_interface__
• pyproj>=3.6 — the transformation engine behind to_crs()
• numpy>=1.26 — masked-array arithmetic on the result

python -m pip install "rasterio>=1.3,<2" "geopandas>=1.0" \
  "shapely>=2.0" "pyproj>=3.6" "numpy>=1.26"

On Windows, take the whole set from one channel — conda install -c conda-forge rasterio geopandas — rather than mixing a pip rasterio wheel with a conda gdal. Two GDAL builds on the DLL search path surface as ImportError: DLL load failed at the first import, never as a message about GDAL itself.

Step-by-Step Implementation

1. Open the raster, then move the vector onto its CRS.

The raster defines the working coordinate system for this operation; the polygons come to it. Reprojecting a few hundred vertices is exact and cheap, whereas reprojecting the grid would resample every pixel. Pin the check with an explicit comparison rather than trusting that both files "are in UTM" — the mechanics of these definitions are in Coordinate Systems with PyProj.

import geopandas as gpd
import rasterio

floodplain = gpd.read_file("floodplain_boundary.gpkg", layer="floodplain")

with rasterio.open("dem_utm32n.tif") as src:
    print(src.crs, src.count, src.shape, src.nodata)
    # EPSG:32632 1 (400, 640) -9999.0

    if floodplain.crs is None:
        raise ValueError("floodplain_boundary.gpkg carries no CRS — set_crs() first")
    if not floodplain.crs.equals(src.crs):
        floodplain = floodplain.to_crs(src.crs)          # move the vector, not the grid

    # shapely geometries implement __geo_interface__, which is all mask() wants
    shapes = list(floodplain.geometry)
    raster_bounds = src.bounds

Two habits pay for themselves here. Pass list(gdf.geometry), not the GeoDataFrame — the frame itself is not an iterable of geometries and raises before it does anything useful. And use crs.equals() rather than ==: two CRS objects can describe the same system through an EPSG code and a WKT string and compare unequal as objects while being identical in effect. Note also that a projected CRS is what you want for anything downstream that measures area; use the local UTM zone or an equal-area projection, never EPSG:3857, whose scale factor grows with latitude.

2. Make the call, and name every argument you rely on.

mask returns a two-tuple of (out_image, out_transform). The array is three-dimensional — (bands, rows, cols) — unless you pass a single integer to indexes, in which case it comes back two-dimensional.

import rasterio
from rasterio.mask import mask

with rasterio.open("dem_utm32n.tif") as src:
    fill = src.nodata if src.nodata is not None else -9999.0
    profile = src.profile.copy()          # capture while the dataset is open

    out_image, out_transform = mask(
        src,
        shapes,
        crop=True,          # window down to the shape bounds
        all_touched=False,  # pixel-centre rule
        invert=False,       # keep what is INSIDE the shapes
        filled=True,        # ndarray with nodata substituted
        nodata=fill,        # explicit — the default is 0 when the file declares none
    )

print(out_image.shape, out_image.dtype)   # (1, 96, 128) float32
print(out_transform)
# | 10.00, 0.00, 383250.00|
# | 0.00,-10.00, 5699560.00|
# | 0.00, 0.00, 1.00|

That nodata=fill line is not decoration. When a dataset declares no nodata value, mask falls back to 0, and on a digital elevation model zero is a perfectly ordinary elevation — the masked margin becomes indistinguishable from sea level, and every mean you compute afterwards is quietly dragged toward it.

3. Flip the three behavioural switches deliberately.

invert=True keeps what lies outside the shapes — punching a hole rather than cutting a patch, which is how you blank a cloud footprint or an exclusion zone. It cannot be combined with crop=True, and rasterio says so: ValueError: crop and invert cannot both be True. The reasoning is geometric, not arbitrary — the complement of a polygon has no meaningful bounding box to crop to.

all_touched=True includes every pixel the polygon boundary so much as clips, instead of only those whose centre falls inside. filled=False returns a numpy.ma.MaskedArray instead of substituting the fill value, which is what you want whenever the next step is arithmetic.

import numpy as np
import rasterio
from rasterio.mask import mask

with rasterio.open("dem_utm32n.tif") as src:
    # blank the exclusion zone, keep the rest of the scene at full extent
    with_hole, hole_transform = mask(src, shapes, crop=False, invert=True, nodata=fill)

    # a masked array: no sentinel to remember, no risk of averaging the fill
    band = mask(src, shapes, crop=True, filled=False)[0][0]

print(type(band))                       # <class 'numpy.ma.MaskedArray'>
print(band.count(), band.size)          # 7841 12288
print(f"mean elevation {band.mean():.1f} m")   # mean elevation 214.6 m

# the same statistic from the filled array needs the sentinel excluded by hand
filled_band = np.ma.masked_equal(out_image[0], fill)
print(f"mean elevation {filled_band.mean():.1f} m")   # mean elevation 214.6 m

Masked arrays are the safer default for analysis and the wrong choice for output: dst.write() wants a plain array, so filled=True — or band.filled(fill) — is what reaches disk.

4. Write the result with a profile that matches the array.

The output profile has to be updated in three places at once — height, width and transform — because crop=True changed all three. Copy the source profile so band count, dtype and compression settings survive, then override.

import rasterio
from rasterio.mask import mask

with rasterio.open("dem_utm32n.tif") as src:
    fill = src.nodata if src.nodata is not None else -9999.0
    out_image, out_transform = mask(src, shapes, crop=True, nodata=fill)
    profile = src.profile.copy()

profile.update(
    driver="GTiff",
    height=out_image.shape[1],     # rows — NOT src.height
    width=out_image.shape[2],      # cols — NOT src.width
    transform=out_transform,       # the new origin
    nodata=fill,
    compress="deflate",
    predictor=3,                   # floating-point predictor; use 2 for integer bands
)

with rasterio.open("dem_floodplain.tif", "w", **profile) as dst:
    dst.write(out_image)

Leaving height and width at their source values raises immediately, which is the friendly failure. Leaving the transform behind does not: the file writes cleanly and lands its 96-by-128 patch at the corner of the original scene, offset by thousands of metres. If the result is destined for a tile service or object storage, add overviews and internal tiling on the way out, following Resampling & Overviews When Writing COGs.

Verification

Three assertions catch every failure described above: the transform must describe the array, the crop must not have changed any value, and the result must not be empty.

import numpy as np
from rasterio.transform import array_bounds

# 1. the transform genuinely describes this array, and it covers the polygon
height, width = out_image.shape[1], out_image.shape[2]
left, bottom, right, top = array_bounds(height, width, out_transform)
minx, miny, maxx, maxy = floodplain.total_bounds
assert left <= minx and right >= maxx and bottom <= miny and top >= maxy, \
    "cropped window does not cover the geometry — wrong transform"

# 2. cropping trims the frame; it must not move a single value
kept = np.ma.masked_equal(out_image[0], fill)
assert np.isclose(kept.mean(), band.mean()), "crop changed the data — impossible"

# 3. an all-nodata result is the CRS-mismatch signature, not a small polygon
assert kept.count() > 0, "no pixels survived — check CRS and extent overlap"
print(f"{kept.count()} of {kept.size} pixels kept, mean {kept.mean():.1f} m")
# 7841 of 12288 pixels kept, mean 214.6 m

Promote that third assertion into any batch job. mask will happily produce an empty raster for every feature in a layer, and the only visible difference between "this polygon is genuinely off the scene" and "every polygon is in the wrong projection" is that the second case yields zero kept pixels for all of them.

Edge Cases & Debugging

• ValueError: crop and invert cannot both be True. — pick one. To blank an area and trim the result, call mask(..., invert=True, crop=False) first, then crop separately with a window derived from the geometry you actually want to keep.
• Output is entirely nodata and only a UserWarning appeared. — degrees against metres. Re-run with crop=True to convert the warning into ValueError: Input shapes do not overlap raster., then fix the CRS with to_crs(src.crs).
• nodata silently became 0. — the source file declares none, so mask defaults to zero. Pass nodata= explicitly, or use filled=False and keep a masked array until the moment you write.
• An integer band cannot hold NaN. — nodata=np.nan on a uint16 raster writes garbage. Either choose a sentinel inside the dtype's range or cast with out_image.astype("float32") before writing and set the profile dtype to match.
• TypeError when passing a GeoDataFrame. — mask needs an iterable of geometries: use list(gdf.geometry), or [mapping(g) for g in gdf.geometry] with shapely.geometry.mapping if you want plain dicts.
• Edge pixels vanish on a thin polygon. — the bounds-derived window rounds outward, but the centre rule still drops slivers. Use crop=True, pad=True, pad_width=0.5 to widen the window by half a pixel, and all_touched=True if the feature is narrower than a cell.
• A self-intersecting boundary burns wrongly. — the rasterizer does not repair topology. Clean the layer first, following Fixing Self-Intersecting Polygons Programmatically.

Frequently Asked Questions

Does crop=True change my pixel values? No. It changes only which pixels are returned, along with the height, width and origin needed to describe them. Any statistic computed over the unmasked pixels is identical either way — the assertion in the verification block proves it. What crop=True does change is memory: crop=False decodes the entire band before masking it, so on a large scene it is the difference between a few hundred kilobytes and several gigabytes. For scenes served over HTTP the same logic scales further, as described in Windowed Reads from Cloud-Optimized GeoTIFF.

When is all_touched=True the right choice? When the feature is narrow relative to the cell size and the centre rule would return almost nothing — a river centreline buffered by five metres over a thirty-metre grid, a road corridor, a field margin. Everywhere else it over-samples the boundary and pulls the result toward whatever surrounds the polygon. If the answer has to be defensible rather than merely non-empty, weight partial pixels by their actual overlap instead, which is the comparison drawn in Zonal Statistics with rasterstats vs exactextract.

Should I use rasterio.mask or rioxarray's clip? Use rasterio.mask when the data is a single scene you are reading from disk and writing back, and you want direct control over the fill value and the output profile. Use rio.clip when the data is already an xarray object — a multi-band or time-series cube — because it applies one mask across every coordinate and keeps the labelled dimensions intact; that path is worked through in Clipping Rasters by Vector Geometry with rioxarray. Both burn the geometry through the same rasterization code underneath.

How do I mask one raster per polygon instead of one for the whole layer? mask treats the shapes you pass as a single set and returns their union, so a layer of two hundred parcels yields one array. For per-feature output, loop and call it once per geometry with crop=True — each iteration reads only that feature's window, so peak memory tracks the largest polygon rather than the scene. If the features overlap and you only need one result per group, dissolve them first using the patterns in Dissolving & Aggregating Features by Attribute.