How to Clip a Raster to a Polygon in Python

On this page

Problem statement

You have an elevation raster covering a large area and you need the part inside one district boundary. Here the raster is 30 m Copernicus elevation for the southern Pennines and the district is Calderdale (the setup is under The data used on this page). The obvious approach reads everything to keep a small part:

import geopandas as gpd
import rasterio

with rasterio.open("elevation_gb.tif") as src:
    data = src.read(1)          # the whole raster, in memory
print(f"{data.shape[1]:,} × {data.shape[0]:,} cells, {data.nbytes / 1e6:,.0f} MB read")

4,474 × 3,725 cells, 67 MB read

On a national raster at the same resolution the same line wants gigabytes, and is a MemoryError on an ordinary machine. And when it does fit, the result is often wrong in a way that only shows up later:

clipped = data[2000:2800, 1500:2400]      # sliced by index
with rasterio.open("sliced.tif", "w", driver="GTiff", height=clipped.shape[0], width=clipped.shape[1],
                   count=1, dtype=clipped.dtype) as dst:
    dst.write(clipped, 1)
with rasterio.open("sliced.tif") as check:
    print(check.crs, check.bounds)

None BoundingBox(left=0.0, bottom=800.0, right=900.0, top=0.0)

That array has no transform and no CRS, so it has no location. Written out, it sits at pixel coordinates from the origin, unconnected to anywhere on Earth.

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)

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

What you get — from 1. Put the polygon into the raster's CRS, not the other way round

boundary: EPSG:4326
raster:   EPSG:27700

Polygon into the raster's CRS, check overlap, then mask

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

boundary = gpd.read_file("catchment.gpkg")
print("boundary:", boundary.crs)

with rasterio.open("elevation_gb.tif") as src:
    print("raster:  ", src.crs)
    boundary = boundary.to_crs(src.crs)

boundary: EPSG:4326
raster:   EPSG:27700

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)

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

in_degrees = gpd.read_file("catchment.gpkg")                     # what a missing to_crs() leaves you with
print("unprojected     ", tuple(round(float(v), 2) for v in in_degrees.total_bounds),
      "overlaps:", raster_box.intersects(in_degrees.union_all()))

raster bounds    (332972, 344725, 467192, 456475)
boundary bounds  (388663, 413197, 418114, 436596)
unprojected      (-2.17, 53.62, -1.73, 53.83) overlaps: False

Printing both bounding boxes usually diagnoses it instantly. Numbers like (-2.2, 53.6, -1.7, 53.8) next to six-digit metres 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`

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      (781, 983)
valid      52.7%
range      42.113815 … 480.30133
bounds in? True

The valid share is the ratio of the district'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.

Clipping without loading the whole raster

The data used on this page

The Copernicus 30 m digital elevation model (GLO-30), read straight from its public bucket on AWS: the two 1° tiles covering 53–54°N, 3–1°W, merged and reprojected to British National Grid at 30 m and saved as elevation_gb.tif — standing in for any large national raster. The boundary is the Calderdale district, from the ONS in WGS84, saved as catchment.gpkg; wards.gpkg holds Calderdale's wards plus one ward from Kingston upon Hull, well east of the raster, so that example 2 has something to miss.

import io

import geopandas as gpd
import numpy as np
import rasterio
import requests
from rasterio.merge import merge
from rasterio.warp import Resampling, calculate_default_transform, reproject

TILE = "https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_{0}_DEM/Copernicus_DSM_COG_10_{0}_DEM.tif"
tiles = [rasterio.open(TILE.format(t)) for t in ("N53_00_W003_00", "N53_00_W002_00")]
mosaic, mosaic_transform = merge(tiles)
src_crs = tiles[0].crs

transform, width, height = calculate_default_transform(src_crs, "EPSG:27700", mosaic.shape[2], mosaic.shape[1],
                                                       *rasterio.transform.array_bounds(mosaic.shape[1], mosaic.shape[2],
                                                                                        mosaic_transform),
                                                       resolution=30)
elevation = np.full((height, width), -9999, dtype="float32")
reproject(mosaic[0], elevation, src_transform=mosaic_transform, src_crs=src_crs, dst_transform=transform,
          dst_crs="EPSG:27700", resampling=Resampling.bilinear, dst_nodata=-9999)
with rasterio.open("elevation_gb.tif", "w", driver="GTiff", height=height, width=width, count=1, dtype="float32",
                   crs="EPSG:27700", transform=transform, nodata=-9999, tiled=True, compress="deflate") as dst:
    dst.write(elevation, 1)

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"


def ons(layer, where, fields):
    r = requests.get(f"{ONS}/{layer}/FeatureServer/0/query", timeout=120,
                     params={"where": where, "outFields": fields, "outSR": 4326, "f": "geojson"})
    return gpd.read_file(io.BytesIO(r.content))


ons("Local_Authority_Districts_DEC_2025_Boundaries_UK_BFC", "LAD25NM = 'Calderdale'", "LAD25NM").to_file("catchment.gpkg")
wards = ons("Wards_December_2024_Boundaries_UK_BFC", "LAD24NM = 'Calderdale'", "WD24NM")
hull = ons("Wards_December_2024_Boundaries_UK_BFC", "LAD24NM = 'Kingston upon Hull, City of'", "WD24NM").head(1)
(gpd.pd.concat([wards, hull], ignore_index=True).rename(columns={"WD24NM": "ward_name"})
   .to_file("wards.gpkg"))
print(f"elevation_gb.tif: {width:,} × {height:,} cells at 30 m; {len(wards)} Calderdale wards + {hull['WD24NM'].iloc[0]}")

elevation_gb.tif: 4,474 × 3,725 cells at 30 m; 17 Calderdale wards + Avenue

The cost of a clip is decided by this line, not by the geometry.

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 numpy as np
import rasterio
from rasterio.features import geometry_mask, geometry_window


def clip_by_window(raster_path, geom, dst_path):
    with rasterio.open(raster_path) as src:
        window = geometry_window(src, [geom])                # whole cells, rounded outward
        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


out = clip_by_window("elevation_gb.tif", boundary.union_all(), "calderdale_window.tif")
with rasterio.open(out) as a, rasterio.open("elevation_catchment.tif") as b:
    same = a.shape == b.shape and a.transform == b.transform and np.array_equal(a.read(1), b.read(1))
    print(f"window clip {a.shape}, mask() clip {b.shape}; identical: {same}")

window clip (781, 983), mask() clip (781, 983); identical: True

geometry_window is the step that is easy to get subtly wrong. The first version of this function built the window with from_bounds(...) and then round_offsets() and round_lengths(), which round to the nearest cell; on Calderdale that came out one row and one column smaller than mask()'s result, cutting off cells that are inside the district. geometry_window rounds outward to whole cells, exactly as mask(crop=True) does, and the check above confirms the two clips are now identical.

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")
results = clip_per_feature("elevation_gb.tif", wards, "out/wards", "ward_name")

  ✓ Brighouse                ok
  ✓ Calder                   ok
  ✓ Elland                   ok
  ✓ Greetland_and_Stainland  ok
  ✓ Hipperholme_and_Lightcliffe ok
  ✓ Illingworth_and_Mixenden ok
  ✓ Luddendenfoot            ok
  ✓ Northowram_and_Shelf     ok
  ✓ Ovenden                  ok
  ✓ Park                     ok
  ✓ Rastrick                 ok
  ✓ Ryburn                   ok
  ✓ Skircoat                 ok
  ✓ Sowerby_Bridge           ok
  ✓ Todmorden                ok
  ✓ Town                     ok
  ✓ Warley                   ok
  ✗ Avenue                   no overlap

The last row is the ward from Hull, added to show the third outcome: it lies outside this raster altogether. 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 — and compare with clipping tight:

from scipy import ndimage


def mean_filter(a, size):
    """A moving-window mean that ignores NaN, instead of spreading it."""
    valid = np.isfinite(a)
    total = ndimage.uniform_filter(np.where(valid, a, 0.0), size)
    count = ndimage.uniform_filter(valid.astype("float64"), size)
    out = total / count
    out[~valid] = np.nan
    return out


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, wide_t = mask(src, [core.buffer(500)], crop=True, filled=False, nodata=src.nodata)
    tight, tight_t = mask(src, [core], crop=True, filled=False, nodata=src.nodata)

plain = ndimage.uniform_filter(tight[0].filled(np.nan), size=9)
print(f"scipy uniform_filter on the tight clip: {np.isnan(plain).mean():.0%} of cells NaN")

smooth_wide = mean_filter(wide[0].filled(np.nan).astype("float64"), 9)
smooth_tight = mean_filter(tight[0].filled(np.nan).astype("float64"), 9)

# line the tight grid up inside the wide one (same source cells, different crop)
row0, col0 = map(int, rasterio.transform.rowcol(wide_t, tight_t.c + 1, tight_t.f - 1))
same_cells = smooth_wide[row0:row0 + tight.shape[1], col0:col0 + tight.shape[2]]
inside = ~tight[0].mask
diff = np.abs(smooth_tight - same_cells)[inside]
print(f"cells inside Calderdale: {inside.sum():,}")
print(f"tight clip vs buffered clip: {(diff > 1).sum():,} cells differ by more than 1 m, "
      f"worst {diff.max():.1f} m")

scipy uniform_filter on the tight clip: 100% of cells NaN
cells inside Calderdale: 404,404
tight clip vs buffered clip: 6,680 cells differ by more than 1 m, worst 23.0 m

Two lessons, one found by accident. First, scipy.ndimage.uniform_filter is built on running sums, so a single NaN does not stay local — it poisons every row and column it touches, and the whole district comes back NaN. A moving-window mean over data with holes needs the normalised form above: average the valid values and divide by how many there were.

Second, the reason for the buffer. With a tight clip, a cell near the boundary averages only the cells inside the district, so its smoothed value is pulled towards whatever lies inland — a rim of wrong numbers exactly where people look first. Reading 500 m wider gives those cells their real neighbours, and trimming to the boundary afterwards costs nothing. 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.

Cropping and masking are two operations

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 many times larger than it needs to be that still looks correct — for Calderdale in this raster, 4,474 × 3,725 cells instead of 983 × 781, twenty-two times the size.

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.

No overlap, geometry lists and nodata

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.

Internal links

The raster data model explained — why the transform must travel with the array
How to reproject a raster in Python with Rasterio — when the polygon's CRS is not negotiable
How to calculate zonal statistics in Python — summarising instead of extracting
Raster and vector do not line up in Python — diagnosing "shapes do not overlap"
How to clip spatial data in Python with GeoPandas — the vector equivalent
How to split a large layer into tiles — read-wide, write-narrow, on the vector side
How to fix invalid geometries in Python — before passing shapes to mask
How to merge and mosaic rasters in Python — reassembling clipped pieces

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.