Clipping Rasters to Field Boundaries with rasterio.mask

TL;DR: Reproject the field-boundary polygon into the raster’s CRS first (a mismatch silently yields an all-nodata output), pass the geometries to rasterio.mask.mask with crop=True and an explicit nodata, write the returned out_transform into the profile, and save a Cloud-Optimized GeoTIFF — then assert the result actually contains valid pixels.

Why Clipping to a Boundary Arises

An orthomosaic or a derived index raster covers the whole flight block: the target field plus headlands, roads, treelines, and the neighbour’s crop. Every downstream statistic — mean NDVI, canopy cover, a yield-zone histogram — is contaminated the moment it includes pixels outside the field. Clipping to the boundary polygon is the step that turns a flight-block raster into a field-specific dataset, and doing it correctly is what makes a zonal mean trustworthy rather than diluted by whatever surrounds the field.

rasterio.mask.mask is the workhorse, but it has one failure mode that catches almost everyone: it does no CRS reconciliation. Hand it a polygon in EPSG:4326 (degrees) and a raster in EPSG:32615 (metres) and it does not raise — the polygon’s tiny decimal-degree coordinates fall nowhere near the metre-scale raster, the intersection is empty, and you get a raster full of nodata with no error. Aligning coordinate systems up front is therefore non-negotiable; if the concept is unfamiliar, read understanding CRS in precision agriculture before running this. The boundary polygon itself typically comes out of field boundary extraction with GeoPandas.

The diagram below shows the four gates every clip passes through.

Field-boundary raster clipping pipeline Four-stage flow: align the boundary CRS to the raster, run rasterio mask with crop and nodata, update the profile with the clipped transform, write a Cloud-Optimized GeoTIFF. Align CRS polygon → raster CRS mask(crop=True) nodata fill Update profile out_transform Write COG tiled + overviews 1. Reconcile 2. Clip 3. Reprofile 4. Export

Prerequisites

Only what differs from the Field Boundary Extraction with GeoPandas guide:

TEXT 
rasterio==1.3.10
geopandas==0.14.4
numpy==1.26.4

Install with:

BASH 
pip install rasterio==1.3.10 geopandas==0.14.4 numpy==1.26.4

Input requirements:

  • A raster (orthomosaic or single-band index such as an NDVI GeoTIFF) with a defined CRS and a set nodata value.
  • A field-boundary polygon as a GeoPackage, GeoParquet, or shapefile readable by GeoPandas — a Polygon or MultiPolygon, valid geometry, any CRS (the script reprojects it).
  • The two extents must genuinely overlap in the raster CRS; a boundary for the wrong field produces a legitimately empty clip that the final assert will reject.

Step-by-Step

One flag that changes what the field's area means Two panels comparing rasterio's all_touched setting. With it off, a pixel is included when its centre falls inside the polygon and the clipped area tracks the boundary's area, which suits statistics. With it on, any touched pixel is included, the area is inflated and mixed edge pixels enter the result. all_touched = False (the default) A pixel is in if its centre is inside The clipped area is close to the polygon's area. Edge pixels that are mostly outside are excluded. Right for statistics and for area-weighted rates. all_touched = True A pixel is in if the polygon touches it at all The clipped area exceeds the polygon's area. Mixed canopy-and-soil pixels enter the statistics. Right for masks, wrong for measurement.

Step 1 — Load both and align the CRS

Open the raster to read its CRS, then reproject the boundary GeoDataFrame into that exact CRS with to_crs. Never assume they match, even when both files “came from the same project” — a boundary digitised in QGIS often lands in EPSG:4326 while the mosaic is in a UTM zone. This single reprojection is what prevents the silent all-nodata output.

Step 2 — Extract geometries for the mask

rasterio.mask.mask wants an iterable of GeoJSON-like geometry mappings, not a GeoDataFrame. Dissolve multiple rows into one if the field is split across features, then pass gdf.geometry (rasterio accepts Shapely objects directly in 1.3). Guard against empty or invalid geometry here rather than debugging a blank raster later.

Step 3 — Run the mask with crop and nodata

Call mask(src, geoms, crop=True, nodata=<value>, all_touched=False, filled=True). crop=True trims the output to the polygon’s bounding box and returns a fresh out_transform. filled=True replaces outside pixels with the nodata value; filled=False instead returns a masked array if you want to keep computing before writing. Leave all_touched=False for area-accurate zonal statistics; flip it on only when you must retain every pixel the polygon grazes.

Step 4 — Update the profile with the clipped transform

The clipped array has new dimensions and a new origin, so the source profile is stale. Copy it and overwrite height, width, transform (from out_transform), and nodata. Skipping this is why some clipped rasters open shifted or squashed — the pixels are right but the georeferencing points at the old window.

Step 5 — Write a Cloud-Optimized GeoTIFF and verify

Add tiling and internal overviews so the output streams efficiently. The complete, directly runnable script:

PYTHON 
import numpy as np
import geopandas as gpd
import rasterio
from rasterio.mask import mask
from rasterio.enums import Resampling

RASTER_PATH = "ndvi_flight_block.tif"
BOUNDARY_PATH = "field_boundary.gpkg"
OUT_PATH = "ndvi_field_clip.tif"
NODATA = -9999.0
ALL_TOUCHED = False   # True keeps every pixel the polygon touches

with rasterio.open(RASTER_PATH) as src:
    raster_crs = src.crs
    assert raster_crs is not None, "Source raster has no CRS — cannot align the mask"

    # ── 1. Load boundary and align CRS to the raster (prevents empty output) ──
    gdf = gpd.read_file(BOUNDARY_PATH)
    assert gdf.crs is not None, "Boundary has no CRS — set it before reprojecting"
    if gdf.crs != raster_crs:
        gdf = gdf.to_crs(raster_crs)

    # Dissolve multi-feature boundaries into one geometry set; drop invalids.
    gdf = gdf[gdf.geometry.notna() & gdf.geometry.is_valid]
    assert not gdf.empty, "No valid boundary geometry after cleaning"

    # Fail fast if the extents do not overlap (wrong field / wrong CRS).
    rb = src.bounds
    minx, miny, maxx, maxy = gdf.total_bounds
    overlaps = (minx < rb.right and maxx > rb.left and
                miny < rb.top and maxy > rb.bottom)
    assert overlaps, "Boundary and raster extents do not intersect in the raster CRS"

    geoms = list(gdf.geometry)

    # ── 2 & 3. Clip: crop to bbox, fill outside with nodata ──────────────────
    out_image, out_transform = mask(
        src,
        geoms,
        crop=True,
        nodata=NODATA,
        all_touched=ALL_TOUCHED,
        filled=True,
    )

    # ── 4. Update the profile with the new window ────────────────────────────
    profile = src.profile.copy()

profile.update(
    height=out_image.shape[1],
    width=out_image.shape[2],
    transform=out_transform,
    nodata=NODATA,
    driver="GTiff",
    tiled=True,
    blockxsize=512,
    blockysize=512,
    compress="deflate",
)

# ── 5. Write a Cloud-Optimized GeoTIFF (tiled + internal overviews) ──────────
with rasterio.open(OUT_PATH, "w", **profile) as dst:
    dst.write(out_image)
    dst.build_overviews([2, 4, 8, 16], Resampling.average)
    dst.update_tags(ns="rio_overview", resampling="average")

# ── Verification: the clip must contain valid (non-nodata) pixels ───────────
with rasterio.open(OUT_PATH) as chk:
    band1 = chk.read(1, masked=True)
    valid = int(band1.count())            # count of unmasked (valid) pixels
    total = band1.size
    print(f"Clipped: {chk.width}×{chk.height} px, {valid}/{total} valid "
          f"({100 * valid / total:.1f}% inside boundary)")
    assert valid > 0, "Clip produced zero valid pixels — check CRS alignment and extent"
    assert chk.crs == raster_crs, "Output CRS drifted from the source raster"

print(f"OK: wrote {OUT_PATH}")

Inline verification: the assert valid > 0 above is the load-bearing check — an all-nodata result (the classic CRS-mismatch symptom) fails here instead of silently propagating into your zonal statistics. The reported percentage of valid pixels also flags a boundary that only partly overlaps the raster.

Gotchas and Edge Cases

  • Silent empty output from a CRS mismatch. If you skip the to_crs reproject in Step 1, mask returns a full-nodata array with no exception. The extent-overlap assert and the final valid > 0 assert both catch it, but the root fix is always aligning the polygon to the raster CRS first — never the other way around, since reprojecting the raster resamples pixels.
The resampling choice follows from what the pixels mean A decision diagram for resampling during a clip. Categorical rasters such as management zone identifiers must use nearest-neighbour so class labels stay whole; continuous surfaces such as reflectance or elevation should use bilinear or cubic interpolation. Raster to resample during the clip are the values categories? yes Nearest neighbour zone identifiers, soil classes and crop stage codes must stay whole numbers no Bilinear or cubic reflectance, elevation and index surfaces are continuous, so interpolation is meaningful and reduces staircase artefacts Averaging class 2 and class 3 into 2.5 produces a zone that does not exist, and nothing downstream will reject it.

  • crop=False leaves a giant sparse raster. Without crop=True, a 2-hectare field clipped from a 200-hectare mosaic keeps the full mosaic dimensions with nearly everything set to nodata — huge on disk and slow to read. Use crop=True and remember it hands back a new out_transform you must write into the profile.

  • all_touched=True inflates area statistics. It keeps every pixel the polygon so much as grazes, adding up to half a pixel of border all the way around. For a mean-NDVI or a total-canopy-area number that shifts the result; keep the default centre-in-polygon rule for statistics and reserve all_touched=True for preserving thin features you cannot afford to drop.

  • Forgetting nodata on the output profile. If the source had no nodata and you do not set one, the filled border pixels become a real data value (often 0), which then counts as valid NDVI of zero and drags your mean down. Always set an explicit nodata both in the mask call and in the written profile.

Frequently Asked Questions

Why does rasterio.mask return an all-nodata raster?

Almost always a CRS mismatch. If the boundary polygon and the raster are in different coordinate systems, their coordinates do not overlap in space, so the mask geometry falls entirely outside the raster and every pixel is filled with nodata. Reproject the polygon into the raster CRS before calling mask, and confirm the two extents actually intersect.

What does all_touched do when clipping to a field boundary?

By default a pixel is kept only if the polygon covers its centre, which shaves a thin rim of edge pixels. Setting all_touched to true keeps every pixel the polygon touches at all, growing the retained area by up to half a pixel around the perimeter. Keep the default for area-accurate statistics and enable all_touched when you must not drop any pixel along a narrow feature.

Do I need crop set to true when masking?

Not strictly, but you almost always want it. With crop set to false the output keeps the full raster dimensions and only replaces outside pixels with nodata, wasting space on a small field in a large mosaic. With crop set to true rasterio trims the array to the polygon bounding box and returns a new transform, which you must write into the output profile so the clipped raster stays georeferenced.

Parent Guide

This guide is part of Field Boundary Extraction with GeoPandas — see there for the full pipeline that produces and validates the boundary polygon this clip consumes.

Parent section Field Boundary Extraction with GeoPandas