⚙️ How-To Raster & rasterio Geometry & topology Raster Analysis · step 5 of 14 Updated

How to Clip a Raster to a Polygon in Python

Problem statement

You have a national elevation raster and you need the part inside one catchment boundary. The obvious approach reads 4 GB to keep 40 MB:

import rasterio, geopandas as gpd

with rasterio.open("elevation_gb.tif") as src:
    data = src.read(1)          # the whole country, in memory
MemoryError

And when it does fit, the result is often wrong in a way that only shows up later:

clipped = data[rows, cols]      # sliced by index

That array has no transform, so it has no location. Written out, it lands at the origin of the coordinate system. Plotted against the catchment boundary, it is somewhere off the coast of Ghana.

Clipping a raster properly means keeping three things in step: the array, the transform that says where the array is, and the mask that says which cells are inside the polygon.

Quick answer 

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

boundary = gpd.read_file("catchment.gpkg")

with rasterio.open("elevation_gb.tif") as src:
    boundary = boundary.to_crs(src.crs)              # ← always, before anything else
    out_image, out_transform = mask(
        src, boundary.geometry, crop=True, filled=True, nodata=src.nodata
    )
    profile = src.profile.copy()

profile.update(
    height=out_image.shape[1],
    width=out_image.shape[2],
    transform=out_transform,
)

with rasterio.open("elevation_catchment.tif", "w", **profile) as dst:
    dst.write(out_image)
Three panels showing the full raster, the bounding-box crop, and the polygon-masked result with corners set to NoData.
`crop=True` shrinks the grid to the bounding box. The mask sets everything outside the polygon to NoData.
Argument Effect if you leave it out
to_crs(src.crs) the polygon misses the raster entirely — empty or ValueError
crop=True output keeps the full national extent, mostly NoData
nodata=src.nodata outside cells become 0, indistinguishable from sea level
profile.update(...) the output claims the wrong size and position

Step-by-step solution

1. Put the polygon into the raster's CRS, not the other way round 

boundary = gpd.read_file("catchment.gpkg")
print(boundary.crs)                      # EPSG:4326

with rasterio.open("elevation_gb.tif") as src:
    print(src.crs)                       # EPSG:27700
    boundary = boundary.to_crs(src.crs)

Reprojecting the polygon is cheap — a few hundred coordinates. Reprojecting the raster rebuilds millions of cells and resamples every value. Always move the vector.

If src.crs is None, stop. There is nothing to align to, and clipping will use whatever numbers happen to be in the transform. See raster and vector do not line up.

2. Check that the two actually overlap

rasterio.mask raises a bare ValueError: Input shapes do not overlap raster when they do not, which tells you nothing about why.

from shapely.geometry import box

with rasterio.open("elevation_gb.tif") as src:
    raster_box = box(*src.bounds)
    boundary = boundary.to_crs(src.crs)

if not raster_box.intersects(boundary.union_all()):
    print("raster bounds   ", tuple(round(v) for v in src.bounds))
    print("boundary bounds ", tuple(round(v) for v in boundary.total_bounds))
    raise SystemExit("no overlap — check the CRS of both inputs")

Printing both bounding boxes usually diagnoses it instantly. Numbers like (-2.7, 53.3, -1.9, 53.9) next to (320000, 405000, 420000, 480000) mean one input is in degrees and the other in metres.

union_all() is the current GeoPandas spelling; older versions use unary_union.

3. Understand crop, filled and invert

Grid comparing crop, filled and invert arguments against what each produces.
Three independent switches. The defaults are rarely the ones you want.
out, transform = mask(src, shapes, crop=True, filled=True, invert=False, nodata=-9999)
  • crop=True shrinks the output grid to the polygon's bounding box. Without it you keep the full raster extent with everything outside the polygon set to NoData — a 4 GB file to hold a 40 MB catchment.
  • filled=True (the default) returns a plain array with outside cells set to nodata. filled=False returns a masked array instead, which is better if you are going straight to statistics.
  • invert=True keeps what is outside the polygon. Useful for erasing an area — a nodata region, a military exclusion zone — rather than extracting one.
  • all_touched=True includes every cell the polygon boundary touches, not only those whose centre falls inside. It grows the result by roughly one cell all round.

4. Update the profile from the actual output 

profile = src.profile.copy()
profile.update(
    height=out_image.shape[1],
    width=out_image.shape[2],
    transform=out_transform,
    nodata=src.nodata,
    compress="deflate",
    tiled=True,
)

Take the shape from out_image, never from your own arithmetic. mask(crop=True) snaps the crop to whole cells of the source grid, so the output bounding box is not exactly the polygon's bounding box — it is the smallest set of whole source cells containing it. Computing width and height from the polygon bounds gives a number that is off by one often enough to matter.

out_image is always 3-D, (bands, rows, cols), even for a single-band raster. That is why the indices are [1] and [2].

5. Verify 

import numpy as np

with rasterio.open("elevation_catchment.tif") as clipped:
    arr = clipped.read(1, masked=True)
    print("shape     ", arr.shape)
    print("valid     ", f"{100 * (1 - arr.mask.mean()):.1f}%")
    print("range     ", arr.min(), "…", arr.max())
    print("bounds in?", box(*clipped.bounds).intersects(boundary.union_all().buffer(1)))
shape      (1264, 986)
valid      63.2%
range      12 … 634
bounds in? True

valid 63.2% is the ratio of the catchment's area to its bounding box — a plausible number for an irregular shape. valid 100% means the mask did not apply, and valid 0% means the polygon and the raster overlap in bounds but not in reality.

Code examples

Example 1: clipping without loading the whole raster

rasterio.mask reads only the windows it needs when the source is tiled, but you can be explicit about it — useful when the raster is remote or enormous:

import rasterio
from rasterio.windows import from_bounds
from rasterio.mask import mask
from rasterio.features import geometry_mask
import numpy as np

def clip_by_window(raster_path, geom, dst_path):
    with rasterio.open(raster_path) as src:
        window = from_bounds(*geom.bounds, transform=src.transform)
        window = window.round_lengths().round_offsets()
        data = src.read(window=window, masked=True)          # only this window
        transform = src.window_transform(window)

        outside = geometry_mask(
            [geom], out_shape=data.shape[1:], transform=transform, invert=False
        )
        data.mask = data.mask | outside

        profile = src.profile.copy()
        profile.update(height=data.shape[1], width=data.shape[2], transform=transform)

    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(data.filled(profile["nodata"]))
    return dst_path

geometry_mask(..., invert=False) returns True for cells outside the geometry, which is exactly the polarity a NumPy mask wants. The naming catches everyone once: invert=True gives True inside.

This pattern is worth knowing because it separates the two things mask does — the window read and the polygon mask — so you can insert work between them, such as scaling values or applying a second mask.

Example 2: one output file per polygon

The common real task is not one clip but hundreds — a raster per ward, per catchment, per plot.

from pathlib import Path
import geopandas as gpd, rasterio
from rasterio.mask import mask
import re

def clip_per_feature(raster_path, gdf, out_dir, name_col, min_valid=0.01):
    out_dir = Path(out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    results = []

    with rasterio.open(raster_path) as src:
        gdf = gdf.to_crs(src.crs)
        for row in gdf.itertuples():
            name = re.sub(r"[^\w\-]+", "_", str(getattr(row, name_col))).strip("_")
            dst = out_dir / f"{name}.tif"
            try:
                arr, transform = mask(src, [row.geometry], crop=True,
                                      filled=False, nodata=src.nodata)
            except ValueError:
                results.append({"name": name, "status": "no overlap"})
                continue

            valid = 1 - arr.mask.mean()
            if valid < min_valid:
                results.append({"name": name, "status": "empty", "valid": valid})
                continue

            profile = src.profile.copy()
            profile.update(height=arr.shape[1], width=arr.shape[2],
                           transform=transform, compress="deflate")
            tmp = dst.with_suffix(".tmp.tif")
            with rasterio.open(tmp, "w", **profile) as out:
                out.write(arr.filled(src.nodata))
            tmp.replace(dst)
            results.append({"name": name, "status": "ok", "valid": valid,
                            "shape": arr.shape[1:]})

    for r in results:
        mark = {"ok": "✓", "empty": "·", "no overlap": "✗"}[r["status"]]
        print(f"  {mark} {r['name']:<24} {r['status']}")
    return results

wards = gpd.read_file("wards.gpkg")
clip_per_feature("elevation_gb.tif", wards, "out/wards", "ward_name")
  ✓ Ancoats_and_Beswick     ok
  ✓ Ardwick                 ok
  · Piccadilly              empty
  ✗ Isle_of_Man             no overlap

The three outcomes are deliberately distinct. "Empty" (the polygon overlaps the raster's bounds but every cell is NoData) and "no overlap" (the polygon is outside the raster entirely) have different causes and different fixes, and collapsing them into one "failed" bucket loses the diagnosis. This is the same reporting discipline as logging and summarising errors in a batch job.

Opening the source once outside the loop matters more than it looks: on a compressed, tiled GeoTIFF, reopening per feature discards the block cache and can cost more than the clipping itself.

Example 3: clipping to a buffered boundary, then trimming

Analyses near an edge need data from beyond it. Clip wide, analyse, then trim:

import geopandas as gpd, rasterio, numpy as np
from rasterio.mask import mask
from rasterio.features import geometry_mask
from scipy import ndimage

boundary = gpd.read_file("catchment.gpkg")

with rasterio.open("elevation_gb.tif") as src:
    boundary = boundary.to_crs(src.crs)
    core = boundary.union_all()
    wide = core.buffer(500)                       # 500 m of context

    arr, transform = mask(src, [wide], crop=True, filled=False, nodata=src.nodata)
    profile = src.profile.copy()

# an operation that looks at neighbours — smoothing needs data past the edge
smoothed = ndimage.uniform_filter(arr[0].filled(np.nan), size=9)

# now trim to the real boundary
outside_core = geometry_mask([core], out_shape=smoothed.shape,
                             transform=transform, invert=False)
smoothed[outside_core] = np.nan

profile.update(height=smoothed.shape[0], width=smoothed.shape[1],
               transform=transform, dtype="float32", nodata=np.nan)
with rasterio.open("elevation_smoothed.tif", "w", **profile) as dst:
    dst.write(smoothed.astype("float32"), 1)

Without the buffer, the smoothing window at the boundary averages real values with NoData and produces a rim of wrong numbers exactly where people look first. The read-wide, write-narrow pattern is the same one used when splitting a large layer into tiles, and for the same reason: the operation has a reach, and the clip must be larger than that reach.

Explanation

Scene showing a polygon boundary crossing cells, with centre-based inclusion compared against all_touched inclusion.
A cell is in or out; there is no half. Which rule you pick biases the area one way or the other.

Clipping a raster is two operations wearing one name, and most confusion comes from not separating them.

The first is cropping — reducing the grid to a smaller rectangle. This is cheap and lossless. The array is a subset of the original, the transform's origin moves to the new top-left corner, and no value changes. crop=True does this.

The second is masking — marking cells outside an arbitrary shape as NoData. The grid stays rectangular, because a raster grid must be; the polygon is expressed by which cells are valid. This changes values (to the sentinel) but not positions.

rasterio.mask.mask does both, which is convenient and slightly misleading. Calling it with crop=False gives you masking alone: a full-extent raster where only the catchment has data. Calling rasterio.windows alone gives cropping without masking: a rectangle of data around the catchment. Most of the time you want both, and the failure mode of forgetting crop=True is a file 100 times larger than it needs to be that still looks correct.

Cell inclusion is a threshold decision. By default a cell is inside if its centre is inside the polygon, so a boundary cutting through the middle of a cell either takes all of it or none. This means the clipped area is not exactly the polygon's area — it is a pixelated approximation, and the error is proportional to the perimeter times half the cell size. On a 25 m raster and a 10 km perimeter, that is around 12.5 hectares of ambiguity. For small polygons relative to cell size this dominates: clipping a 30 m raster to a 40 m plot can return one cell, four cells, or none, depending on where the plot falls. all_touched=True biases the other way — it over-includes rather than under-includes — and neither is more correct in general. When the difference matters, zonal statistics with area weighting is the honest tool.

The transform is the part that must not be lost. Slicing a NumPy array is easy; the hard part of clipping is that the result's position changed, and NumPy has no idea. mask returns out_transform for exactly this reason, and every bug where a clipped raster ends up in the Gulf of Guinea comes from writing the array with the original profile. If you take one habit from this article, take this one: the array and the transform travel together, always.

Edge cases or notes

  • ValueError: Input shapes do not overlap raster is almost always a CRS mismatch, not a genuine lack of overlap. Print both bounding boxes.
  • Pass a list of geometries, not a GeoDataFrame. mask(src, gdf.geometry, ...) works because a GeoSeries is iterable; mask(src, gdf, ...) does not.
  • Several polygons are unioned, not kept separate. Passing five shapes gives one raster covering all five. For one file each, loop — see Example 2.
  • crop=True snaps to source cells. The output bounds are the smallest whole-cell rectangle containing the polygon, never the polygon's exact bounds.
  • out_image is always 3-D. out_image[0] for a single band.
  • nodata=0 on unsigned data is a trap — zero is usually a real value. Use a sentinel outside the data range, or a float dtype with NaN.
  • filled=False returns a masked array, which is what you want if the next step is mean() or sum().
  • Invalid polygons make mask fail or produce nonsense. Run the geometry validity fix first.
  • gdalwarp -cutline boundary.gpkg -crop_to_cutline in.tif out.tif is the command-line equivalent and streams, so it handles rasters too large for memory.

FAQ

Why do I get "Input shapes do not overlap raster"?

The polygon and the raster are in different coordinate systems nine times out of ten. Reproject the polygon with gdf.to_crs(src.crs) and print both bounding boxes to confirm.

Should I reproject the raster or the polygon?

The polygon. It has a few hundred coordinates; the raster has millions of cells, and reprojecting it resamples every value.

What does crop=True actually do?

It shrinks the output grid to the polygon's bounding box, snapped to whole source cells. Without it the output keeps the full source extent with everything outside the polygon set to NoData.

How do I keep the area outside the polygon instead?

mask(src, shapes, invert=True). Leave crop=False, since cropping to the bounding box of the area you are removing makes no sense.

Why is my clipped raster in the wrong place?

The output was written with the source profile instead of out_transform. The array and its transform must be updated together.

Should I use all_touched=True?

Use it when under-inclusion is worse than over-inclusion — small polygons, thin features, coastal cells. It includes every cell the boundary touches instead of only those whose centre is inside.

How do I clip a raster that will not fit in memory?

rasterio.mask already reads only the needed windows on a tiled source. For very large jobs, gdalwarp -cutline -crop_to_cutline streams and never loads the whole file.