Fixing EPSG Mismatches in rasterio.open

Direct answer: call src.crs.to_epsg() immediately after opening the file, compare it against your expected EPSG code, then either override the metadata tag via MemoryFile (zero resampling, milliseconds) or reproject the pixel grid with rasterio.warp.reproject (necessary only when coordinates truly diverge).

import rasterio
from rasterio.crs import CRS
from rasterio.io import MemoryFile

with rasterio.open("scene.tif") as src:
    if src.crs is None or src.crs.to_epsg() != 32633:
        target_crs = CRS.from_epsg(32633)
        memfile = MemoryFile()
        profile = src.profile | {"crs": target_crs}
        with memfile.open(**profile) as dst:
            dst.write(src.read())
        # keep `memfile` alive — freeing it invalidates `ds`
        ds = memfile.open()

This approach is part of the Mastering CRS Transformations in Rasterio workflow. For the broader context of raster ingestion architecture, see Core Raster Fundamentals & STAC Mapping.

Why EPSG Mismatches Occur During Ingestion

Remote sensing pipelines routinely ingest heterogeneous data: legacy GeoTIFFs, cloud-optimized archives, and STAC catalog assets. Mismatches concentrate around three root causes.

Header vs. data divergence. The TIFF tags declare one EPSG (e.g., EPSG:4326), but the pixel grid coordinates and affine transform clearly follow a projected system like EPSG:32633. This pattern is common when data is manually edited in desktop GIS software without updating projection metadata. You can catch it quickly: if src.crs claims EPSG:4326 (degree-valued) but src.transform.c is a seven-digit number, the label is wrong.

Deprecated or ambiguous WKT strings. Older GDAL versions wrote custom Well-Known Text strings that map to multiple EPSG identifiers. Modern PROJ rejects ambiguous definitions unless explicitly resolved, so src.crs.to_epsg() returns None even when a valid CRS is embedded.

Catalog metadata overrides. proj:epsg fields in STAC items frequently conflict with actual raster headers, especially after mosaicking or manual curation. Relying solely on catalog JSON without validating the binary payload introduces silent spatial drift — a particularly dangerous failure mode when feeding data into spectral index calculation pipelines or spatial joins.

The diagram below shows where a mismatch can enter the pipeline and the two correction branches available:

EPSG mismatch correction decision flow A flowchart showing rasterio.open leading to CRS validation, then branching into a fast metadata-only override path (MemoryFile, no resampling) and a slow reproject path (calculate_default_transform + warp). rasterio.open(path) read header src.crs.to_epsg() cross-check bounds EPSG correct? Use as-is no action needed coords diverge? MemoryFile override tag fix, no resample warp.reproject pixel grid transform yes no tag only diverge

Environment & Setup

Package Minimum version Why required
rasterio 1.3 rasterio.open, MemoryFile, warp.reproject
pyproj 3.4 CRS authority lookups, Transformer
numpy 1.24 Array operations on band data 
pip install "rasterio>=1.3" "pyproj>=3.4" "numpy>=1.24"

Complete Working Example

The function below handles both the fast path (metadata tag wrong, coordinates correct) and the slow path (pixel grid genuinely needs transformation). It returns a (dataset, memfile) tuple so the caller can control the lifetime of the in-memory buffer.

Relabel or reproject? Compare the dataset's bounds against the declared CRS. If the numbers are impossible for that CRS — degrees outside plus or minus 180, or metre coordinates in a geographic system — the label is wrong and the fix is to overwrite the CRS without touching pixels. If the bounds are plausible but disagree with another layer, the data really is in a different CRS and needs a warp. Do the bounds make sense for the CRS the file claims? print(src.crs, src.bounds) the first command, every time impossible numbers plausible numbers the label is wrong e.g. EPSG:4326 with bounds of 499980 … 609780 the data really differs correct UTM bounds, but a different zone from your grid rewrite the CRS tag only — no warp calculate_default_transform + reproject Choosing the wrong one moves data 
import rasterio
from rasterio.crs import CRS
from rasterio.io import MemoryFile
from rasterio.warp import reproject, Resampling, calculate_default_transform


def open_with_verified_crs(
    filepath: str,
    target_epsg: int,
    force_reproject: bool = False,
) -> tuple[rasterio.DatasetReader, MemoryFile | None]:
    """
    Open a raster and guarantee the returned dataset carries the correct CRS.

    Parameters
    ----------
    filepath : str
        Path to a local GeoTIFF or any rasterio-supported format.
    target_epsg : int
        The EPSG code the caller expects (e.g. 32633 for UTM Zone 33N).
    force_reproject : bool
        If True, warp the pixel grid even when only the metadata tag differs.
        Use this when bounds validation confirms true coordinate divergence.

    Returns
    -------
    (dataset, memfile)
        Keep `memfile` alive in the caller's scope. If it is garbage-collected,
        the in-memory buffer is freed and `dataset` becomes invalid — which
        causes a segmentation fault, not a Python exception.
    """
    target_crs = CRS.from_epsg(target_epsg)

    with rasterio.open(filepath) as src:
        if src.crs is None:
            raise ValueError(
                f"No CRS metadata found in '{filepath}'. "
                "Assign a CRS explicitly before ingestion."
            )

        src_epsg = src.crs.to_epsg()  # returns None if ambiguous

        # --- Fast path: EPSG already matches; return a direct file handle ---
        if src_epsg == target_epsg and not force_reproject:
            return rasterio.open(filepath), None

        memfile = MemoryFile()

        if not force_reproject:
            # Metadata-only override: copy raw pixel data without any
            # resampling. Safe only when coordinates are already in the
            # target system but the TIFF tag is wrong (e.g. labelled 4326
            # while the affine values are clearly in UTM metres).
            profile = src.profile | {"crs": target_crs}
            with memfile.open(**profile) as dst:
                dst.write(src.read())  # single bulk read, no per-band loop
            return memfile.open(), memfile

        # --- Slow path: true coordinate transformation ---
        # calculate_default_transform derives the optimal output extent and
        # resolution so no source pixels are clipped or over-sampled.
        transform, width, height = calculate_default_transform(
            src.crs,
            target_crs,
            src.width,
            src.height,
            *src.bounds,  # left, bottom, right, top in src CRS units
        )
        profile = src.profile | {
            "crs": target_crs,
            "transform": transform,
            "width": width,
            "height": height,
        }
        with memfile.open(**profile) as dst:
            for band_idx in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, band_idx),
                    destination=rasterio.band(dst, band_idx),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=target_crs,
                    resampling=Resampling.bilinear,
                    # bilinear is appropriate for continuous data (reflectance,
                    # elevation); use Resampling.nearest for categorical data
                    # (land cover, cloud masks).
                )
        return memfile.open(), memfile

Usage:

ds, memfile = open_with_verified_crs("sentinel2_B04.tif", target_epsg=32633)
print(ds.crs.to_epsg())   # 32633
print(ds.bounds)
# keep ds and memfile in scope until you are done reading
ds.close()
if memfile:
    memfile.close()

Variant Patterns

Variant 1 — Sentinel-2 STAC asset with proj:epsg conflict

Reading the CRS family from the numbers alone Geographic coordinates stay within plus or minus 180 and 90. UTM eastings run from 166021 to 833978 with northings up to ten million. Web Mercator reaches twenty million in both axes. Recognising the magnitude tells you what the file really is before you look at any metadata tag. Bounds magnitude is a fingerprint looks like typical bounds so the CRS is small decimals 36.1 … 37.0, −1.4 … −0.6 geographic — EPSG:4326 six and seven digits 499980 … 609780, 9890200 UTM — EPSG:326xx / 327xx seven and eight digits, both axes −20037508 … 20037508 Web Mercator — EPSG:3857 degrees, but src.crs says UTM 36.1 … 37.0 tagged EPSG:32636 mislabelled — fix the tag, not the pixels A northing above 10,000,000 in a northern-hemisphere UTM zone is the other classic tell: it belongs to a southern zone.

When loading assets via a STAC catalog query, the proj:epsg field in the item’s properties may disagree with the actual GeoTIFF header. This variant reads the STAC-declared EPSG, compares it against the header, and applies the override only when necessary:

import rasterio
from rasterio.crs import CRS
from rasterio.io import MemoryFile
import pystac


def open_stac_asset_verified(
    item: pystac.Item,
    asset_key: str = "B04",
) -> tuple[rasterio.DatasetReader, MemoryFile | None]:
    """
    Open a STAC asset, using the raster header as ground truth for the CRS.
    The STAC proj:epsg field is treated as a hint only.
    """
    asset_href = item.assets[asset_key].href
    stac_epsg: int | None = item.properties.get("proj:epsg")

    with rasterio.open(asset_href) as src:
        header_epsg = src.crs.to_epsg() if src.crs else None

    if stac_epsg and header_epsg and stac_epsg != header_epsg:
        # STAC metadata disagrees with the header; trust the header.
        # Log the discrepancy for audit trails.
        print(
            f"CRS conflict for {item.id}/{asset_key}: "
            f"STAC says EPSG:{stac_epsg}, header says EPSG:{header_epsg}. "
            f"Using header value."
        )

    # Validate bounds against the header EPSG regardless of STAC claim
    expected_epsg = header_epsg or stac_epsg
    if expected_epsg is None:
        raise ValueError(f"Cannot determine CRS for {item.id}/{asset_key}")

    return open_with_verified_crs(asset_href, target_epsg=expected_epsg)

Variant 2 — Batch pipeline with pyproj bounds validation

For batch ingestion of heterogeneous archives, cross-referencing src.bounds against the official EPSG geographic extent catches mislabelled files before they corrupt downstream outputs. The handling pixel resolution and scaling cluster covers how affine transform parameters interact with the CRS choice.

from pyproj import CRS as ProjCRS
from pyproj.aoi import AreaOfInterest
from pyproj.database import query_utm_crs_info
import rasterio


def validate_crs_by_bounds(filepath: str) -> dict:
    """
    Return a validation report: whether the embedded CRS is plausible
    given the raster's geographic extent.
    """
    with rasterio.open(filepath) as src:
        if src.crs is None:
            return {"valid": False, "reason": "no_crs"}

        epsg = src.crs.to_epsg()
        bounds = src.bounds

        # Convert bounds to WGS 84 degrees for the plausibility check
        proj_crs = ProjCRS.from_user_input(src.crs)
        wgs84_bounds = proj_crs.area_of_use

        if wgs84_bounds is None:
            return {"valid": None, "reason": "no_aou_data", "epsg": epsg}

        # A crude sanity check: the raster's geographic footprint must
        # overlap the CRS's declared area of use.
        in_range = (
            wgs84_bounds.west <= bounds.right
            and wgs84_bounds.east >= bounds.left
            and wgs84_bounds.south <= bounds.top
            and wgs84_bounds.north >= bounds.bottom
        )
        return {
            "valid": in_range,
            "epsg": epsg,
            "reason": "bounds_ok" if in_range else "bounds_outside_aou",
        }


# Example batch run
import pathlib

for tif in pathlib.Path("./scenes").glob("*.tif"):
    report = validate_crs_by_bounds(str(tif))
    if not report["valid"]:
        print(f"WARNING {tif.name}: {report}")

Variant 3 — Dask-friendly variant for large mosaics

When working with large mosaics processed via rasterio.vrt.WarpedVRT, the MemoryFile pattern does not scale. Instead, pass the correct CRS directly to WarpedVRT, which performs lazy on-demand warping without materialising the entire array:

import rasterio
from rasterio.crs import CRS
from rasterio.vrt import WarpedVRT

def open_warped_vrt(filepath: str, target_epsg: int) -> WarpedVRT:
    """
    Return a WarpedVRT that exposes the dataset in `target_epsg` without
    loading all pixels into memory. Compatible with rasterio's window-read
    interface used in Dask pipelines.
    """
    src = rasterio.open(filepath)
    # WarpedVRT wraps the source and handles on-demand pixel transformation.
    # No data is read until a window or array is requested.
    return WarpedVRT(src, crs=CRS.from_epsg(target_epsg))

Common Errors

AttributeError: 'NoneType' object has no attribute 'to_epsg' src.crs is None — the file has no CRS metadata at all. Assign one via the MemoryFile pattern above or reject the file early with a clear error message rather than propagating None downstream.

rasterio.errors.CRSError: Input is not a CRS CRS.from_epsg() received an invalid code (e.g. None or a deprecated authority string). Verify the integer before calling; EPSG codes are positive integers in the range 1024–32767 for common projected and geographic systems.

Segmentation fault after MemoryFile.close() The MemoryFile object was garbage-collected while a dataset handle opened from it was still in use. Always return (or store) both (dataset, memfile) and call memfile.close() only after dataset.close().

Frequently Asked Questions

Q: Why does rasterio.open return None for src.crs? The embedded WKT or GeoTIFF tag is either missing, malformed, or too ambiguous for PROJ to resolve to a single authority code. Use CRS.from_wkt() or assign a known EPSG explicitly.

Q: When should I use a metadata override versus a full reproject? Use the metadata override (copy pixels, change only the CRS tag) when bounds validation confirms the coordinates are already in the target system — e.g. the file is labelled EPSG:4326 but the affine values are clearly in UTM metres. Use a full reproject only when the pixel coordinates genuinely need transformation.

Q: Why does my MemoryFile dataset become invalid after the function returns? The MemoryFile buffer is freed when its reference is garbage-collected. Return the MemoryFile object alongside the dataset handle and keep both alive in the calling scope.