Clipping Rasters by Vector Geometry with rioxarray

Cutting a raster down to a study-area boundary looks like a one-liner, and it is — right up until the polygons arrive in a different projection, the result comes back empty, or the clipped GeoTIFF renders with a black rectangle where the outside should be transparent. This guide walks the whole rio.clip path for anyone working with xarray & rioxarray raster cubes inside Mastering Core Geospatial Python Libraries: CRS alignment first, then the four arguments that decide what the output actually contains, then getting it back out as a file other tools will read correctly.

Why This Approach / What Goes Wrong

.rio.clip() does two separate things in one call, and almost every surprise comes from confusing them. First it rasterizes your geometries onto the array's own grid, producing a boolean mask the same width and height as the raster. Second it applies that mask — setting everything outside to nodata — and optionally crops the array down to the mask's bounding window. Masking changes values; cropping changes the extent, the shape and the affine transform. drop is the switch between "same array, some cells blanked" and "smaller array".

The rasterization step is where the CRS bites. rioxarray builds the mask by feeding your geometry coordinates through the raster's affine transform, so the numbers in the geometry have to be in the raster's coordinate system. A GeoDataFrame in EPSG:4326 carries longitudes near 11.4 and latitudes near 50.9; a Sentinel-2 tile in a UTM zone has eastings near 412000. Hand the first to the second and the mask lands nowhere near the grid, and rioxarray raises rioxarray.exceptions.NoDataInBounds. That exception is the good outcome. The bad one is two projected CRSs that overlap numerically — a neighbouring UTM zone, or metres versus US survey feet — where the mask lands somewhere on the grid and you get a plausible-looking raster of the wrong place. clip accepts a crs= argument naming the CRS your geometries are in, and transforms them for you; nothing anywhere infers it.

The CRS gate in front of every rioxarray clip A parcels GeoDataFrame in EPSG 4326 and an NDVI cube in EPSG 32633 both feed a CRS gate that asks whether the geometry CRS equals the raster CRS. On the failure branch, where crs is not passed, degree coordinates are rasterized as if they were metres, the mask lands off the grid, zero pixels are selected and rioxarray raises NoDataInBounds. On the success branch, after calling to_crs or passing crs explicitly, the flow continues along the bottom of the figure: the geometries are rasterized onto the cube's grid with geometry_mask, the mask is applied with where so that outside cells become NaN, and drop equals True then crops the array to the mask window. One gate decides everything: are the geometries in the raster's CRS? parcels GeoDataFrame · 812 polygons EPSG:4326 ndvi_cube DataArray (time, y, x) EPSG:32633 CRS gate geoms.crs == cube.rio.crs rioxarray never guesses no yes no — and crs= was not passed degrees get rasterized as if they were metres the mask lands off the grid → zero pixels kept raise NoDataInBounds yes — after .to_crs(cube.rio.crs) cube.rio.clip(geoms, parcels.crs) or pass crs= and let rioxarray transform them rasterize onto the cube grid features.geometry_mask() apply the mask .where(mask) → NaN outside then crop, or don't drop=True | drop=False Masking changes values; cropping changes the extent, the shape and the transform.
The clip is really three stages — transform, rasterize, apply — and only the first one can go silently wrong.

Prerequisites

python -m pip install \
  "rioxarray>=0.15" "xarray>=2024.7" "rasterio>=1.3.9" \
  "geopandas>=1.0" "shapely>=2.0" "pyproj>=3.6"

Keep rioxarray and rasterio from one package source. The accessor calls rasterio's masking internals directly, so a pip rasterio wheel layered over a conda GDAL produces an import-time DLL failure rather than a clip error — the same rule laid out in Raster Data Handling with Rasterio.

Step-by-Step Implementation

1. Open the raster as a masked, chunked array and load the boundary.

masked=True is not cosmetic here: it promotes the array to float and turns the file's nodata into NaN, which gives the clip a value it can legally write outside the polygon. Skip it on integer data and you inherit the dtype problem covered in step 6.

import geopandas as gpd
import rioxarray  # registers the .rio accessor

ndvi_cube = rioxarray.open_rasterio(
    "ndvi_2023_summer.tif", masked=True, chunks={"x": 1024, "y": 1024}
).squeeze("band", drop=True)

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

print(ndvi_cube.rio.crs)   # EPSG:32633
print(floodplain.crs)      # EPSG:4326  <- different, and that is the whole problem

2. Align the CRS before anything else.

Two options, both correct. Reproject the vector layer with to_crs() and pass the geometries as-is, or leave the layer alone and tell clip which CRS the geometries are in via crs=. Reprojecting the vector side is almost always cheaper — a few hundred polygons versus tens of millions of pixels — and it leaves you with a geometry you can plot and sanity-check. Whichever you pick, never reproject the raster just to make a clip work: warping resamples every pixel, and the mechanics of doing it deliberately belong in reprojecting raster cubes with reproject_match. If you build a transformer by hand to check a corner coordinate, pass always_xy=True, because EPSG:4326 is formally latitude-first and PyProj honours that.

# Reproject the vector side — cheap, inspectable, and the mask is built from these numbers
floodplain = floodplain.to_crs(ndvi_cube.rio.crs)
assert floodplain.crs == ndvi_cube.rio.crs

# One dissolved geometry rasterizes faster than 800 separate rings
boundary = floodplain.union_all()          # GeoPandas 1.0+; .unary_union on 0.x

3. Clip, and choose drop deliberately.

clip takes an iterable of geometries — Shapely objects or GeoJSON-like mappings — so [boundary], floodplain.geometry.values and floodplain.geometry all work. drop=True, the default, crops the array to the window that still holds data: the shape shrinks, the x/y coordinates shrink with it, and the affine transform moves. drop=False keeps the original grid intact and only blanks the outside, which is what you want whenever the result has to stay pixel-aligned with another layer in the same stack.

drop=True versus drop=False on the same six-by-four raster The same polygon clips the same four-row by six-column raster twice. With drop equals True, the default, only the two-by-three window of cells the polygon touches survives; the original extent is shown as a dashed outline and the result has shape two by three with a new affine transform. With drop equals False the full four-by-six grid is returned, the same six cells hold their values, and the remaining eighteen cells are set to NaN, so the shape and the transform are unchanged. Same mask, same values — two different extents drop=True · the default original extent (gone) 0.62 0.58 0.61 0.60 0.64 0.59 shape (2, 3) · new transform cells outside the window no longer exist drop=False extent preserved nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan 0.62 0.58 0.61 0.60 0.64 0.59 shape (4, 6) · transform unchanged still stackable with its neighbours on the grid
Use drop=True when the clip is the final product and drop=False when the array has to keep lining up with the rest of the cube.
clipped = ndvi_cube.rio.clip(
    [boundary],
    crs=floodplain.crs,   # state it even after to_crs — explicit beats implicit
    drop=True,
)

print(ndvi_cube.shape)   # (10980, 10980)
print(clipped.shape)     # (1842, 2416)  <- cropped to the mask window

4. Reach for all_touched=True only when the geometry is thin.

By default a pixel is kept when its centre falls inside the geometry. A 15 m riparian corridor on a 30 m grid can cross a column of pixels without covering a single centre, so the default rule selects nothing and clip raises NoDataInBounds on a geometry that plainly overlaps the raster. all_touched=True keeps every pixel the boundary intersects at all, which fixes slivers and over-samples everything else — the same trade-off that governs masking rasters with polygons in Python.

A narrow stream buffer under the two all_touched settings A 15 metre stream buffer, drawn as a narrow slanted band, crosses a grid of 30 metre pixels whose centres are marked with dots. On the left, with all_touched left at its default of False, the band passes between two columns of pixel centres and covers none of them, so zero pixels are selected and rioxarray raises NoDataInBounds. On the right, with all_touched set to True, the eight pixels the band touches at all are highlighted and selected, at the cost of including cells that are mostly outside the buffer. A 15 m buffer on a 30 m grid misses every pixel centre all_touched=False · centroid rule 0 pixels selected NoDataInBounds: no data found in bounds all_touched=True 8 pixels selected every cell the band touches, mostly-outside ones included Dots mark pixel centres — the default rule tests those, not the pixel area.
The centroid rule is the right default for study areas hundreds of pixels wide; all_touched=True exists for the geometries that are narrower than a cell.
from rioxarray.exceptions import NoDataInBounds

stream_buffer = gpd.read_file("stream_15m_buffer.gpkg").to_crs(ndvi_cube.rio.crs)

try:
    corridor = ndvi_cube.rio.clip(stream_buffer.geometry, stream_buffer.crs)
except NoDataInBounds:
    # Narrower than a pixel: fall back to touched-cell selection
    corridor = ndvi_cube.rio.clip(
        stream_buffer.geometry, stream_buffer.crs, all_touched=True
    )

5. Invert the mask to punch holes instead of cutting shapes.

invert=True keeps what is outside the geometries — removing lakes from a terrain model, or blanking a restricted zone before publishing. Note that drop is still evaluated against the pixels that survive, and those surround the hole, so invert=True with drop=True normally crops nothing at all. Pass drop=False to make the intent obvious.

water = gpd.read_file("water_bodies.gpkg").to_crs(ndvi_cube.rio.crs)

# Everything except the lakes, on the original grid
land_only = ndvi_cube.rio.clip(
    [water.union_all()], water.crs, invert=True, drop=False
)

6. Crop with clip_box first, and keep nodata honest.

clip rasterizes the mask at the raster's full width and height, so clipping a small parcel out of a 10980 × 10980 scene builds a 120-megapixel boolean array to keep a few thousand cells. clip_box is the cheap rectangular pre-crop: it slices by coordinate window with no rasterization at all. Run it first, then clip the much smaller result. Because the box is expressed in the raster's CRS units, take the bounds from the already-reprojected GeoDataFrame — and pad by a pixel so a boundary sitting exactly on a cell edge is not lost to rounding.

minx, miny, maxx, maxy = floodplain.total_bounds     # already in the cube's CRS
pad = abs(ndvi_cube.rio.resolution()[0])             # one pixel of slack

subset = ndvi_cube.rio.clip_box(
    minx - pad, miny - pad, maxx + pad, maxy + pad,
    auto_expand=True,       # widen the box if it degenerates to a single row/column
)
clipped = subset.rio.clip([boundary], floodplain.crs, drop=True)

Now the nodata question. clip blanks the outside with .where(), which produces NaN, and then fills those cells with .rio.nodata if one is set. On a float array opened with masked=True there is nothing to do — NaN is already the fill and it round-trips into a COG cleanly. On an integer array such as a land-cover class raster there is no integer that means "absent", so declare one before clipping rather than after.

How the source nodata setting propagates through a clip into the output file Three rows trace nodata from the source array to the written COG. A float32 array opened with masked equals True carries NaN as its nodata, the clip writes NaN outside the polygon, and the COG is tagged nodata equals nan so viewers render it transparent. An int16 land-cover array with write_nodata minus 9999 called first gets minus 9999 outside, and the COG is tagged so statistics skip those cells. An int16 array with no nodata set has no value that means absent, so the outside cells are untagged and downstream tools count them as real class codes. Whatever .rio.nodata says is what lands outside the polygon source array what clip writes outside what the COG carries open_rasterio(masked=True) float32 · nodata is NaN NaN float array, nothing to convert nodata=nan renders transparent .rio.write_nodata(-9999) int16 land cover · declared first -9999 a value the dtype can hold nodata=-9999 statistics skip those cells .rio.nodata is None int16 land cover · nothing declared no value means “absent” the fill is undefined for the dtype no nodata tag outside counted as real classes Declare nodata before the clip, not after — the fill happens inside the call.
Every downstream statistic depends on this one attribute; an untagged clip is the reason zonal means come back suspiciously low.

7. Write the clip back out as a COG.

Clipping a Dask-backed cube stays lazy — the mask is built from the geometry, not from pixel values — so nothing is read until the write. Re-attach the CRS and nodata immediately before to_raster, because intermediate xarray operations drop non-dimension coordinates without warning.

import numpy as np

out = (
    clipped
    .rio.write_crs(ndvi_cube.rio.crs)
    .rio.write_nodata(np.nan, encoded=True)
)
out.rio.to_raster(
    "ndvi_floodplain_2023.tif",
    driver="COG",
    compress="DEFLATE",
    predictor=2,
)

The COG driver builds the tiling and overview pyramid itself; the overview resampling choice for categorical versus continuous data is worked through in resampling and overviews when writing COGs, and the wider storage picture in Cloud-Native Geospatial Formats.

Verification

A clip is correct when the CRSs matched going in, the extent came out inside the geometry's bounds, and some data actually survived. Check all three — a clip that returns an array of pure NaN is not an error, it is a silent misalignment.

import numpy as np

# 1. The mask was built in the raster's own coordinates
assert floodplain.crs == ndvi_cube.rio.crs, "Reproject the vector layer before clipping"

# 2. The clipped extent sits inside the geometry bounds, within one pixel
res = abs(ndvi_cube.rio.resolution()[0])
cminx, cminy, cmaxx, cmaxy = clipped.rio.bounds(recalc=True)
gminx, gminy, gmaxx, gmaxy = floodplain.total_bounds
assert cminx >= gminx - res and cmaxx <= gmaxx + res

# 3. Something survived
valid = int(np.isfinite(clipped.values).sum())
assert valid > 0, "All-NaN result — geometry and raster do not really overlap"

print(f"{clipped.shape} kept, {valid} valid cells, CRS {clipped.rio.crs.to_epsg()}")
# (1842, 2416) kept, 2214883 valid cells, CRS 32633
print(clipped.rio.nodata)      # nan

Round-trip the file too: reopen it and confirm the nodata tag and transform survived the write.

check = rioxarray.open_rasterio("ndvi_floodplain_2023.tif", masked=True)
assert check.rio.crs == clipped.rio.crs
assert np.isnan(check.rio.nodata)
print(check.rio.transform(recalc=True))
# | 10.00, 0.00, 418340.00|
# | 0.00,-10.00, 5637820.00|

Edge Cases & Debugging

Frequently Asked Questions

Does clip reproject my geometries automatically? Only when you tell it which CRS they are in. Passing crs=parcels.crs lets rioxarray transform the geometries onto the raster's coordinate system before rasterizing. Omit it and the geometries are assumed to already be in the raster's CRS — no check, no warning. Reprojecting the vector layer with to_crs() and passing crs= anyway costs nothing and makes both halves of the assumption visible.

Should I use clip or clip_box? clip_box for a rectangle, clip for a shape. The box version is a pure coordinate-window slice with no rasterization, so it is dramatically cheaper and is the right first step on a large scene even when a polygon clip follows. Use clip alone only when the raster is already small relative to the geometry's bounding box.

Why is my clipped raster black instead of transparent? The nodata tag did not make it into the file. to_raster writes what .rio.nodata holds at write time, and several xarray operations drop the attribute along the way. Call .rio.write_nodata(np.nan, encoded=True) on the final object immediately before writing, then reopen the file and assert on check.rio.nodata as in the verification block.

Does clipping a Dask-backed cube load it into memory? No. The mask is derived from the geometry and the affine transform, never from pixel values, so clip extends the task graph like any other lazy operation and the time dimension is untouched. Data is read at .compute() or at to_raster. The exception is from_disk=True, which deliberately goes back to the file with rasterio to read a smaller window.

Related