Fixing Rasterio CRS Mismatch Errors

TL;DR: A rasterio mask or overlay that returns an all-nodata array is almost always a CRS mismatch — reproject the vector to src.crs (or align the second raster with WarpedVRT) and place an assert src.crs == other.crs guard before every masking call so the failure raises instead of returning empty data.

Why a CRS Mismatch Produces Empty Output, Not an Error

rasterio.mask.mask and the windowed clip patterns underneath it work by intersecting a geometry with the raster’s pixel grid. When the geometry and the raster share a CRS, that intersection is the field. When they do not — a field boundary in geographic degrees (EPSG:4326) clipped against a raster in UTM metres (EPSG:32615) — the geometry’s coordinates (roughly −93, 42) fall nowhere near the raster’s (roughly 448000, 4651000). The intersection is empty, and here is the trap: rasterio does not treat an empty intersection as an error. It treats it as a perfectly valid clip that happens to contain zero valid pixels, and hands back an array filled with the nodata value.

Nothing raises. A downstream np.nanmean over that array returns nan, a zonal-statistics loop records 0.0 for the field, and the corrupt value flows into the season’s record indistinguishable from a real measurement. On a vegetation-index pipeline this is the difference between “this field is stressed” and “this field has no data” — and the two look identical until someone audits the raw arrays. The same mechanism strikes when two rasters are combined: subtracting a NIR raster in one CRS from a red raster in another produces a shape-compatible but spatially nonsensical result, or a broadcast error if the grids also differ in size. This guide is the raster-specific fix within the broader debugging CRS and projection errors reference.

Why the clip is empty instead of wrong On the left, a boundary in degrees near minus 93 by 41 and a raster in UTM metres near 500000 by 4540000 occupy completely different coordinate ranges, so their intersection is empty and the clip returns no pixels without raising an error. On the right, after reprojecting the boundary into the raster's CRS, the two overlap and the clip returns the field. Before — different coordinate ranges, empty intersection boundary EPSG:4326 x ≈ −93.4, y ≈ 41.9 raster EPSG:32615 x ≈ 500000, y ≈ 4540000 clip → 0 pixels no exception raised to_crs(src.crs) After — one CRS, real overlap raster · EPSG:32615 field clip → 61,842 pixels Reproject the vector, not the raster a few coordinates versus millions of resampled pixels

Prerequisites

Only one package differs in emphasis from the parent section — you need rasterio with its warp and vrt submodules, which ship in the standard wheel.

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 multi-band or single-band GeoTIFF with a populated .crs and .transform.
  • A field-boundary vector (GeoPackage, GeoJSON, or shapefile) with a populated .crs. If its CRS is None, assign it with set_crs first — see the parent guide, this fix assumes both CRS are known.

Step-by-Step

Four lines that turn an empty result into a message Four diagnostic steps: print the coordinate reference system of both the raster and the vector, compare them with an equality method rather than string comparison, reproject the vector into the raster's CRS, and assert that the extents overlap before attempting the clip. Print both CRSs src.crs and gdf.crs Compare properly CRS.equals, not str() Reproject the vector not the raster Assert overlap before clipping Comparing WKT strings reports a mismatch between two identical CRSs written by different tools, which sends you looking for a bug that is not there.

Step 1 — Detect the mismatch with a guard

Never mask before this assertion. It converts the silent empty-array failure into an immediate, informative exception.

PYTHON 
import rasterio
import geopandas as gpd
from pyproj import CRS

def assert_crs_match(raster_crs, vector_crs) -> None:
    a, b = CRS.from_user_input(raster_crs), CRS.from_user_input(vector_crs)
    assert a.to_epsg() == b.to_epsg(), (
        f"CRS mismatch: raster EPSG:{a.to_epsg()} != vector EPSG:{b.to_epsg()}. "
        "Reproject the geometry to the raster CRS before masking — otherwise the "
        "clip returns an all-nodata array with no error."
    )

Step 2 — Reproject the clip geometry to the raster CRS

The cheap, correct fix for raster/vector mismatch is to move the (small) vector into the (large) raster’s CRS, not the reverse. Reprojecting a boundary is a handful of coordinates; reprojecting a 10 GB orthomosaic is expensive and lossy.

PYTHON 
def load_geometry_in_raster_crs(vector_path: str, raster_crs) -> list:
    gdf = gpd.read_file(vector_path)
    assert gdf.crs is not None, "Vector CRS is None — set_crs before reprojecting"
    gdf = gdf.to_crs(raster_crs)          # match the raster exactly
    assert_crs_match(raster_crs, gdf.crs)
    return list(gdf.geometry)

Step 3 — Align a second raster on the fly with WarpedVRT

When two rasters disagree (a red band in EPSG:32615 and a NIR band delivered in EPSG:32616 from a neighbouring UTM zone), a WarpedVRT presents the second raster as if it were in the target CRS and grid, without writing a new file. Reads through the VRT are reprojected on demand.

PYTHON 
from rasterio.vrt import WarpedVRT
from rasterio.enums import Resampling

def open_aligned(secondary_path: str, reference: rasterio.DatasetReader) -> WarpedVRT:
    """Return the secondary raster warped to the reference CRS, grid, and size."""
    src = rasterio.open(secondary_path)
    return WarpedVRT(
        src,
        crs=reference.crs,
        transform=reference.transform,
        width=reference.width,
        height=reference.height,
        resampling=Resampling.bilinear,   # continuous data (reflectance/NDVI)
    )

Step 4 — Materialise a reprojected raster with rasterio.warp.reproject

When a raster must be permanently moved into the working CRS (so downstream tools open it repeatedly without re-warping), write a new GeoTIFF with rasterio.warp.reproject. calculate_default_transform derives the destination grid.

Step 5 — Complete runnable script

This combines the pattern end to end: guard, reproject the geometry, clip, and — as the reusable utility — reproject a raster to a target EPSG on disk.

PYTHON 
import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio.mask import mask
import geopandas as gpd
from pyproj import CRS


def assert_crs_match(raster_crs, vector_crs) -> None:
    a, b = CRS.from_user_input(raster_crs), CRS.from_user_input(vector_crs)
    assert a.to_epsg() == b.to_epsg(), (
        f"CRS mismatch: raster EPSG:{a.to_epsg()} != vector EPSG:{b.to_epsg()}"
    )


def reproject_raster(src_path: str, dst_path: str, dst_epsg: int,
                     resampling: Resampling = Resampling.bilinear) -> str:
    """Materialise src_path into EPSG:dst_epsg as a new GeoTIFF."""
    dst_crs = CRS.from_epsg(dst_epsg)
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds
        )
        profile = src.profile.copy()
        profile.update(crs=dst_crs, transform=transform,
                       width=width, height=height)
        with rasterio.open(dst_path, "w", **profile) as dst:
            for b in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, b),
                    destination=rasterio.band(dst, b),
                    src_transform=src.transform, src_crs=src.crs,
                    dst_transform=transform, dst_crs=dst_crs,
                    resampling=resampling,
                )
    return dst_path


def clip_raster_to_field(raster_path: str, vector_path: str) -> np.ndarray:
    """Clip a raster to a field boundary, guarding against the empty-clip trap."""
    with rasterio.open(raster_path) as src:
        gdf = gpd.read_file(vector_path)
        assert gdf.crs is not None, "Vector CRS is None — set_crs first"
        gdf = gdf.to_crs(src.crs)               # move vector to raster CRS
        assert_crs_match(src.crs, gdf.crs)      # guard BEFORE masking

        clipped, out_transform = mask(
            src, list(gdf.geometry), crop=True, filled=True, nodata=src.nodata
        )

    # Fail loudly if the clip is empty — the classic silent CRS-mismatch symptom
    valid = clipped != (src.nodata if src.nodata is not None else 0)
    assert valid.any(), (
        "Clip returned all-nodata — geometry does not overlap the raster even "
        "after reprojection; check for a datum relabel or wrong source CRS"
    )
    return clipped


if __name__ == "__main__":
    arr = clip_raster_to_field("ndvi_20240715.tif", "field_boundaries.gpkg")
    print(f"clipped shape: {arr.shape}, valid pixels: {(arr != 0).sum():,}")

Inline verification — confirm alignment and non-empty output separately from the clip:

PYTHON 
with rasterio.open("ndvi_20240715.tif") as src:
    fields = gpd.read_file("field_boundaries.gpkg").to_crs(src.crs)
    assert src.crs.to_epsg() == fields.crs.to_epsg()
    # bounding boxes must physically overlap
    b, fb = src.bounds, fields.total_bounds
    assert min(b.right, fb[2]) > max(b.left, fb[0]), "No X overlap — still mismatched"
    print("raster and vector overlap; safe to mask")

Gotchas & Edge Cases

  • Reprojecting the raster when you only needed to reproject the vector. Warping a large orthomosaic to match a boundary is slow and resamples every pixel; move the geometry to the raster CRS instead (Step 2). Only materialise a reprojected raster (Step 4) when downstream tools genuinely need the raster itself in a new CRS.
  • Wrong resampling on categorical rasters. WarpedVRT and reproject default matters: use Resampling.nearest for management-zone or mask-code rasters. bilinear/average interpolates class integers into meaningless fractional values, silently mislabeling zones fed to threshold mapping for crop health.
  • Same EPSG code, no overlap. If the guard passes but the clip is still empty, the raster or vector was relabelled with a CRS it is not actually in (a NAD83 layer tagged EPSG:4326). The bounding-box overlap assertion in the verification snippet catches this; the fix is a datum transform, covered in resolving pyproj datum shift warnings.
  • WarpedVRT left open. A WarpedVRT holds an open handle to the underlying dataset; use it inside a with block or call .close() on both the VRT and its source, or a batch loop will exhaust file descriptors on large flight sets.
Reproject the cheap side, and only once A table comparing what it costs to reproject the vector boundary, the raster, or the raster twice. Reprojecting the boundary is effectively free and leaves pixels untouched; reprojecting the raster resamples every pixel; doing it twice compounds the resampling error. What you reproject Cost What changes The field boundary a few dozen coordinates Microseconds Nothing — pixels untouched The raster millions of pixels Seconds to minutes Every pixel resampled once The raster, twice clip then export Twice over Resampling error compounds

Frequently Asked Questions

Why does rasterio mask return all nodata instead of raising an error?

The clip geometry and the raster are in different coordinate systems, so the geometry falls entirely outside the raster extent and the intersection is empty. Rasterio treats a non-overlapping geometry as a valid clip that happens to contain no data, filling the output with nodata rather than raising. Reproject the geometry to the raster CRS before masking, and assert the two CRS are equal first.

Should I use WarpedVRT or rasterio.warp.reproject to fix a mismatch?

Use WarpedVRT when you need a raster to appear in a target CRS on the fly without writing a new file, for example to align two rasters inside a windowed loop. Use rasterio.warp.reproject when you want to materialise and save a reprojected GeoTIFF that downstream tools will open repeatedly. WarpedVRT avoids disk writes but recomputes on every read.

Which resampling method should I use when reprojecting NDVI rasters?

Use bilinear resampling for continuous data such as NDVI or reflectance so reprojection does not introduce blocky artefacts. Use nearest-neighbour for categorical rasters such as management-zone or mask codes, because averaging class integers produces meaningless in-between values. Choosing average or cubic on a class raster silently corrupts the zone labels.

Parent Guide

This guide is part of Debugging CRS and Projection Errors in Python — see there for the full diagnosis flowchart and the assertion guards that cover missing CRS, datum shift, and axis-order failures alongside this one.

Parent section Debugging CRS and Projection Errors in Python Related Resolving pyproj Datum Shift Warnings