Masking Rasters to Stand Boundaries with rasterio.mask

Clipping a raster to a stand boundary is the most-run operation in forest GIS and the one with the most silent variants. Whether a pixel on the boundary belongs to the stand, whether the result keeps the original grid or gets a new one, and what happens to pixels outside the polygon are three independent decisions, and different defaults across tools give measurably different hectare figures. This guide settles them, as a focused operation inside Raster-Vector Overlay Techniques, part of Ecological GIS Data Foundations in Python.

When to use this approach

If the end product is a table of means and areas per stand, do not clip at all — zonal statistics is a single pass and avoids writing a raster per stand. Clipping earns its place when the masked pixels themselves are the product, or when a subsequent operation needs a raster.

The decision that changes your numbers

all_touched decides what happens to a pixel the boundary crosses. With the default False, a pixel is included when its centre falls inside the polygon; with True, when any part of it is covered.

The difference scales with the perimeter-to-area ratio, which is why it matters far more in forestry than in most raster work. A compact 100 ha stand at 20 m resolution has a boundary about 4 km long, touching roughly 200 pixels out of 2,500 — so the two settings can differ by several percent. A narrow riparian buffer 30 m wide is nearly all boundary, and the two settings can differ by half.

The rule that keeps results defensible: use all_touched=False when you are measuring (area, mean, biomass), because centre-based selection is unbiased — boundary pixels are as likely to fall outside as inside. Use all_touched=True when you are extracting for display or for a subsequent operation that must not have holes, such as a mask that will be eroded anyway.

Minimal reproducible example

import geopandas as gpd
import rasterio
from rasterio.mask import mask


def clip_to_stands(raster_path: str, stands_path: str, out_path: str,
                   all_touched: bool = False, nodata: float = -9999.0):
    """Clip a raster to a stand layer, writing a cropped GeoTIFF."""
    stands = gpd.read_file(stands_path)
    with rasterio.open(raster_path) as src:
        shapes = stands.to_crs(src.crs).geometry.values          # CRS match is mandatory
        data, transform = mask(
            src, shapes,
            crop=True,               # shrink the output to the geometry's extent
            all_touched=all_touched,
            filled=True,
            nodata=nodata,
        )
        profile = src.profile | {
            "height": data.shape[1],
            "width": data.shape[2],
            "transform": transform,
            "nodata": nodata,
            "compress": "deflate",
            "predictor": 2,
            "tiled": True,
        }
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(data)
    return out_path

filled=True with an explicit nodata is deliberate. The alternative, filled=False, returns a masked array, which is more expressive in memory but loses the mask the moment anyone calls .data or writes it without the mask — a class of bug that produces stands whose “outside” pixels are quietly zero.

To exclude rather than keep — cutting lakes, roads, or a retention area out of a stand raster — invert:

def exclude_features(raster_path: str, exclusions_path: str, out_path: str,
                     buffer_m: float = 0.0, nodata: float = -9999.0):
    """Set every pixel inside the exclusion layer to NoData, keeping the rest."""
    excl = gpd.read_file(exclusions_path)
    with rasterio.open(raster_path) as src:
        shapes = excl.to_crs(src.crs).geometry
        if buffer_m:
            shapes = shapes.buffer(buffer_m)
        data, transform = mask(src, shapes.values, invert=True, crop=False,
                               filled=True, nodata=nodata)
        profile = src.profile | {"nodata": nodata, "transform": transform}
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(data)
    return out_path

Note crop=False here. When excluding, the output should keep the original extent and grid; cropping to the exclusion layer’s bounds would be exactly backwards.

Parameter reference

Expected output and verification

The masked raster’s pixel count should reconstruct the polygon area to within a boundary pixel or so. That comparison is the check worth automating, because it catches CRS mismatches, all_touched surprises, and NoData confusion in one number.

import geopandas as gpd
import numpy as np
import rasterio


def verify_mask(clipped_path: str, stands_path: str, equal_area_epsg: int = 3035) -> dict:
    """Compare masked pixel area against the polygon area from an equal-area CRS."""
    stands = gpd.read_file(stands_path)
    poly_ha = stands.to_crs(epsg=equal_area_epsg).area.sum() / 10_000.0

    with rasterio.open(clipped_path) as src:
        arr = src.read(1, masked=True)
        cell_m2 = abs(src.transform.a) * abs(src.transform.e)
    valid_px = int((~arr.mask).sum()) if np.ma.isMaskedArray(arr) else int(arr.size)
    raster_ha = valid_px * cell_m2 / 10_000.0

    return {
        "polygon_ha": round(poly_ha, 2),
        "raster_ha": round(raster_ha, 2),
        "diff_pct": round(100.0 * (raster_ha - poly_ha) / poly_ha, 2),
        "valid_px": valid_px,
    }


report = verify_mask("chm_stands.tif", "stands.gpkg")
assert abs(report["diff_pct"]) < 5.0, f"masked area disagrees with the polygons: {report}"

A difference of a percent or two is the boundary discretisation and is expected. A difference of 30% or more means all_touched is not what you thought, or the raster is in a different CRS than the polygons were reprojected to. A difference near 100% means the mask returned nothing and the geometries do not overlap the raster at all.

Common pitfalls

• Passing geometries in the wrong CRS. mask does not reproject. Mismatched CRSs usually raise “Input shapes do not overlap raster”, but when the two happen to overlap numerically you get a plausible clip of the wrong ground.
• Leaving nodata at 0. Zero is a real canopy height, a real reflectance, and a real NDVI. Use a sentinel outside the data range and set it in the profile.
• crop=True when excluding. Inverting to remove lakes and cropping at the same time returns a raster the size of the lakes.
• Assuming the clip preserves the grid. With crop=True the origin moves to the geometry’s bounding window. That is fine, unless you are stacking several clips that must align — in which case crop them all to one common window instead.
• Masking band by band in a loop. mask handles all bands at once and reads the window once; looping multiplies the I/O.

Frequently Asked Questions

Related

• Raster-Vector Overlay Techniques — the wider treatment of combining the two data models
• Zonal Statistics for Forest Inventory with rasterstats — the alternative when you only need summaries
• Choosing an Equal-Area CRS for Forest Area Statistics — the frame the verification step measures in
• Reading Cloud Optimized GeoTIFFs with rasterio Windows — clipping a remote raster without downloading it

Up: Raster-Vector Overlay Techniques · Ecological GIS Data Foundations in Python