Debugging CRS and Projection Errors in Python

A prescription map that is silently offset by two metres, a raster clip that returns an all-nodata array with no error, a yield layer that plots in the Gulf of Guinea instead of Iowa — these are the daily coordinate-reference-system failures of a Python ag-GIS pipeline. They rarely raise exceptions. Instead they produce output that looks plausible, passes a naive shape check, and only reveals itself when a planter drives the wrong rows or an agronomist notices a field boundary hanging in open water. This guide is a systematic diagnosis-to-fix reference for the CRS and projection bugs that corrupt geopandas and rasterio workflows, with copy-ready guard functions that convert silent misalignment into loud, early assertion failures.

The failure modes covered here are: a raster and a vector disagreeing on CRS (the empty-clip trap), a layer whose CRS is None, running metric operations on EPSG:4326 degrees, the ~1–2 m NAD83-versus-WGS84 datum offset that ruins RTK-grade prescriptions, axis-order (latitude/longitude versus longitude/latitude) swaps, UTM zone-boundary crossings, and pyproj’s “ballpark transformation” datum-grid warnings. This page is part of the Ag-GIS Data Fundamentals & Spatial Reference Systems section; start with understanding CRS in precision agriculture if you need the conceptual grounding before debugging a live failure.

Diagnosis flowchart: my two layers do not line up A decision tree starting from two misaligned layers. It first checks whether either CRS is None, then whether the two CRS codes match, then whether the offset is a small constant datum shift, then whether coordinates appear swapped, ending in the correct fix for each branch. Layers do not line up clip empty / offset / wrong hemisphere Is either CRS None? src.crs / gdf.crs yes Assign the known CRS set_crs — never reproject no Do EPSG codes match? a.to_epsg() == b.to_epsg() no Reproject to one CRS to_crs / warp.reproject yes Offset small + constant? ~1–2 m everywhere yes Datum shift NAD83↔WGS84 transform no Coords swapped? lat/lon vs lon/lat yes Axis order always_xy=True no Field spans UTM zones? crosses 6° meridian yes Single-zone or UTM pick one zone / use CONUS Albers no Alignment is correct — inspect transform, resampling, or nodata next Work top to bottom — the first "yes" is usually the whole bug.

Prerequisites

Python packages (exact versions tested):

  • rasterio==1.3.10
  • geopandas==0.14.4
  • pyproj==3.6.1
  • shapely==2.0.4
  • numpy==1.26.4

Install with:

BASH 
pip install rasterio==1.3.10 geopandas==0.14.4 pyproj==3.6.1 shapely==2.0.4 numpy==1.26.4

Input and environment requirements:

  • Every raster carries a populated .crs (a rasterio.crs.CRS) and an affine .transform; every vector layer carries a populated gdf.crs (a pyproj.CRS). Layers with crs is None are the first thing this guide diagnoses.
  • A designated working CRS for each field — a projected, metric system, almost always the local UTM zone (for example EPSG:32615 for UTM zone 15 N over central Iowa). All metric operations (distance, area, buffering, kriging, rasterio.mask) run in this CRS, never in EPSG:4326.
  • The PROJ data directory reachable by pyproj. Datum-shift transforms between NAD83 and WGS84 need PROJ transformation grids; the ballpark-transformation warning is resolved in the dedicated guide linked below.
  • Field extents small enough to sit inside one UTM zone (a 6°-wide band ≈ 400+ km). Multi-zone operations are flagged as an explicit failure mode.
  • RTK-sourced points (±1–3 cm) and yield-monitor points (±1–3 m) — knowing which accuracy tier a layer belongs to decides whether a 1–2 m datum offset is a real bug or below noise.

1. Concept & Why CRS Bugs Stay Silent

A coordinate reference system binds a pair of numbers to a location on the Earth. It has three parts that each fail differently: a datum (the reference ellipsoid and its realisation, e.g. NAD83 versus WGS84), a coordinate system (geographic degrees versus projected metres, and the axis order), and, for projected systems, a projection (e.g. Transverse Mercator for a UTM zone). A bug in any one of the three produces a different signature, and the reason these bugs are so costly in agriculture is that none of them raises a Python exception on its own.

Consider the four dominant patterns:

CRS mismatch (the empty-clip trap). When you clip a raster to a field boundary and the two layers declare different CRS, their coordinate ranges do not overlap. rasterio.mask.mask computes the intersection of the geometry with the raster window, finds nothing, and returns an array full of the nodata value — no error, just empty data. Downstream, an NDVI mean over that array is nan or zero and quietly enters the record for that field. This exact failure and its WarpedVRT fix are the subject of fixing rasterio CRS mismatch errors.

Missing CRS (None). A shapefile without a .prj sidecar or a CSV of GPS points loaded straight into geopandas has crs is None. The geometry numbers are fine, but no library knows what they mean. The dangerous fix is to reproject it (which no-ops or errors), when the correct fix is to assign the CRS you know the data is already in with set_crs. Confusing set_crs with to_crs is one of the most common mistakes and one of the easiest to guard.

Datum shift (the RTK killer). NAD83 and WGS84 were coincident in the 1980s but the North American plate has since carried NAD83 roughly 1–2 m away from WGS84. If a controller logs planting in WGS84 and your field boundary is NAD83, treating them as “close enough” offsets every row by that amount. At the ±1–3 cm accuracy an RTK planter delivers, a 1–2 m shift is catastrophic; at the ±1–3 m accuracy of a yield monitor it may be within noise. Whether the shift matters is an agronomic judgement, not a software default, which is exactly why resolving pyproj datum shift warnings treats it as a deliberate decision rather than a silent transform.

Axis order (lat/lon versus lon/lat). EPSG:4326’s authority definition orders axes latitude-then-longitude. GeoJSON, shapefiles, and virtually every raster store longitude-then-latitude. When pyproj honours the authority order, a point at (42.0 N, 93.6 W) round-trips to coordinates that place the field near the equator off the coast of West Africa. The always_xy=True flag forces pyproj to speak longitude-then-latitude and eliminates the whole class of swap.

The correct working posture is defensive: pick one metric working CRS per field, reproject everything into it at ingest, and assert equality at every boundary so that a mismatch fails loudly instead of producing a plausible wrong number. The rest of this page builds that posture into reusable functions.

2. Step-by-Step Diagnosis and Guards

The first question splits the problem in half A decision diagram for misaligned layers. If the extents do not overlap at all, the layers are in different coordinate reference systems and every operation returns empty without error. If they do overlap but are offset, the projection is not the fault and the datum, epoch and positional accuracy are what to check. Two layers will not line up do their extents overlap at all? yes Same CRS, real offset check the datum and epoch, then the positional error budget — the projection is not the fault no Different CRSs, or one is None the coordinate ranges are incompatible, so the intersection is empty and every clip returns nothing without raising

Step 1 — Inventory the declared CRS of every input

Before touching geometry, print what each layer claims to be. Most misalignments are visible in this inventory alone.

PYTHON 
import rasterio
import geopandas as gpd

def describe_crs(label: str, crs) -> None:
    """Print a compact CRS fingerprint for one layer."""
    if crs is None:
        print(f"{label:<22} CRS = None  <-- UNDEFINED, must set_crs before use")
        return
    from pyproj import CRS
    c = CRS.from_user_input(crs)
    print(
        f"{label:<22} EPSG:{c.to_epsg()}  "
        f"{'projected' if c.is_projected else 'GEOGRAPHIC(deg)':<15} "
        f"datum={c.datum.name if c.datum else '?'}"
    )

raster_path = "ndvi_20240715.tif"
field_path = "field_boundaries.gpkg"

with rasterio.open(raster_path) as src:
    describe_crs("raster (ndvi)", src.crs)

fields = gpd.read_file(field_path)
describe_crs("vector (fields)", fields.crs)

A typical smoking gun looks like raster EPSG:32615 projected next to vector EPSG:4326 GEOGRAPHIC(deg) — a guaranteed empty clip. Two identical EPSG codes but different datum names point at a datum-shift bug instead.

Step 2 — Classify the failure with a diagnosis function

Encode the flowchart as code so the diagnosis is repeatable rather than a manual eyeball each time.

PYTHON 
from pyproj import CRS

def diagnose_pair(crs_a, crs_b) -> str:
    """Return the most likely CRS failure class for two layers."""
    if crs_a is None or crs_b is None:
        return "MISSING_CRS: one layer has crs=None — assign with set_crs, do not reproject"

    a, b = CRS.from_user_input(crs_a), CRS.from_user_input(crs_b)

    if a.to_epsg() != b.to_epsg():
        # Same underlying datum + axes but different projection is still a mismatch
        return f"CRS_MISMATCH: EPSG:{a.to_epsg()} vs EPSG:{b.to_epsg()} — reproject to one"

    if a.datum and b.datum and a.datum.name != b.datum.name:
        return "DATUM_SHIFT: same code, different datum realisation — transform explicitly"

    if a.is_geographic and b.is_geographic:
        return "CHECK_AXIS_ORDER: geographic pair — confirm always_xy=True in transformers"

    return "CRS_OK: codes and datums agree — inspect transform, resampling, or nodata next"

print(diagnose_pair("EPSG:4326", "EPSG:32615"))
print(diagnose_pair("EPSG:32615", "EPSG:32615"))

Step 3 — Reproject to a single metric working CRS

Once the class is known, collapse everything into one projected working CRS with explicit EPSG codes. Never let a library pick a CRS for you.

PYTHON 
WORKING_EPSG = 32615  # UTM 15N — replace with your field's zone

def to_working_crs(gdf: gpd.GeoDataFrame, working_epsg: int = WORKING_EPSG) -> gpd.GeoDataFrame:
    """Reproject a vector layer into the metric working CRS, guarding missing CRS."""
    assert gdf.crs is not None, (
        "Vector CRS is None. Assign the CRS the data was recorded in with "
        "gdf.set_crs('EPSG:4326') BEFORE calling this function — set_crs relabels, "
        "to_crs reprojects. Using the wrong one silently corrupts coordinates."
    )
    out = gdf.to_crs(epsg=working_epsg)
    assert out.crs.is_projected, "Working CRS must be projected/metric"
    return out

fields_m = to_working_crs(fields)
print(f"fields now in EPSG:{fields_m.crs.to_epsg()}, projected={fields_m.crs.is_projected}")

The set_crs-versus-to_crs distinction is the single most consequential line here. set_crs attaches a label without moving any coordinates; use it when the CRS is missing or wrong-but-known. to_crs recomputes every coordinate; use it to move correctly-labelled data between systems. The mechanics of the WGS84→UTM move are covered in depth in how to convert WGS84 to UTM for farm mapping.

Step 4 — Guard raster/vector boundaries with an assertion

The empty-clip trap is defeated by one assertion placed immediately before any masking, sampling, or overlay call.

PYTHON 
from rasterio.crs import CRS as RioCRS

def assert_same_crs(raster_crs, vector_crs, *, tol_context: str = "") -> None:
    """Raise before an operation that requires two layers to share a CRS."""
    assert raster_crs is not None, f"Raster CRS is None {tol_context}"
    assert vector_crs is not None, f"Vector CRS is None {tol_context}"
    a = CRS.from_user_input(raster_crs)
    b = CRS.from_user_input(vector_crs)
    assert a.to_epsg() == b.to_epsg(), (
        f"CRS mismatch {tol_context}: raster EPSG:{a.to_epsg()} != "
        f"vector EPSG:{b.to_epsg()}. Reproject the geometry before masking — "
        "otherwise the clip returns an all-nodata array with no error."
    )

with rasterio.open(raster_path) as src:
    fields_r = fields.to_crs(src.crs)          # match the raster
    assert_same_crs(src.crs, fields_r.crs, tol_context="(pre-mask)")
    # ... rasterio.mask.mask(src, fields_r.geometry) is now safe

Step 5 — Verify against a known ground control coordinate

The final guard is empirical: transform a coordinate whose true position you know (a surveyed field corner, an RTK base station, a road intersection you can read off imagery) and assert the result lands where it should, within the accuracy tier of the data.

PYTHON 
from pyproj import Transformer

# A surveyed field corner: WGS84 lon/lat and its known UTM 15N easting/northing
KNOWN_LONLAT = (-93.6250, 42.0250)       # (lon, lat) — note lon first
KNOWN_UTM15N = (448_236.0, 4_651_890.0)  # (easting, northing), metres

tr = Transformer.from_crs("EPSG:4326", f"EPSG:{WORKING_EPSG}", always_xy=True)
e, n = tr.transform(*KNOWN_LONLAT)
print(f"transformed: E={e:.1f}  N={n:.1f}")

# RTK-grade tolerance is centimetres; here allow 2 m for a hand-read control point
assert abs(e - KNOWN_UTM15N[0]) < 2.0 and abs(n - KNOWN_UTM15N[1]) < 2.0, (
    "Control point missed target — check always_xy, datum, and EPSG code"
)
print("control-point check passed")

If this assertion fires with the coordinates swapped (easting and northing roughly transposed, or the point in the wrong hemisphere), the cause is axis order — the transformer was built without always_xy=True. If it fires with a consistent ~1–2 m miss, the cause is a datum shift. If it fires wildly, the EPSG code is wrong.

3. Key Parameters & Tuning

These are the parameters that most often decide whether a CRS operation is correct. Get them wrong and the pipeline still runs.

Parameter Type Default Agronomic Effect
always_xy bool False Passed to pyproj.Transformer.from_crs. Left False, EPSG:4326 uses lat/lon order and swaps coordinates, throwing points to the wrong hemisphere. Set True for every transformer touching farm data so pyproj speaks lon/lat like GeoJSON and shapefiles.
allow_ballpark bool True Passed to Transformer.from_crs. When True, pyproj falls back to a datum-free approximation (the “ballpark” warning) that can be off by 1–2 m. Set False on RTK-grade pipelines so a missing NAD83↔WGS84 grid fails loudly instead of silently degrading planting accuracy.
accuracy float (m) None Rejects any transform whose stated accuracy is worse than this many metres. Set to 0.05 for RTK guidance layers; leave None for ±1–3 m yield data where sub-metre precision is meaningless.
working CRS (WORKING_EPSG) int — The single metric EPSG all layers are reprojected into. Must be the field’s UTM zone (or a regional equal-area CRS). A geographic default (4326) breaks every distance, area, and buffer calculation downstream.
set_crs vs to_crs method — set_crs relabels without moving coordinates (use for missing/wrong CRS); to_crs recomputes coordinates (use to change a correct CRS). Swapping them either corrupts coordinates or no-ops the reprojection.
resampling enum nearest Passed to rasterio.warp.reproject. Use nearest for categorical rasters (management zones, mask codes) to preserve class values; bilinear for continuous NDVI to avoid blocky reprojection artefacts.
grid tolerance float grid-defined The NAD83↔WGS84 transform grid embeds its own accuracy (typically a few cm). Absent the grid, the transform degrades to ballpark ±1–2 m — the difference between an on-row and off-row prescription.

4. Edge Cases & Failure Modes

The all-nodata clip that never errors. A raster in EPSG:32615 clipped against a boundary in EPSG:4326 returns an array filled with nodata, not an exception. Any np.nanmean over it yields nan or 0.0 and silently poisons that field’s record. The fix is the Step 4 assertion plus reprojecting the geometry to the raster CRS first; the deeper reprojection-based recovery using WarpedVRT is detailed in fixing rasterio CRS mismatch errors.

Four failures that never raise an exception A table of four coordinate-system failures that produce no error: the deprecated init syntax bypassing datum grids, an assumed axis order placing data in the wrong hemisphere, GeoJSON written in a projected CRS, and nodata padding added during reprojection drifting index statistics. Silent failure What you see Guard init=epsg: syntax deprecated Datum grids never consulted CRS.from_epsg, pyproj 3+ Axis order assumed lat/lon vs lon/lat Data lands in the wrong hemisphere Let the library order it GeoJSON written in UTM RFC 7946 requires 4326 Web maps misplace it silently Reproject before writing Nodata padding after warp empty edge pixels Index means drift downward Mask after reprojection

Metric operations on EPSG:4326 degrees. Calling gdf.buffer(30) or .area on a layer still in EPSG:4326 buffers by 30 degrees and computes area in square degrees. There is no error — just a buffer that engulfs a continent. A one degree of longitude spans ~78 km at 45° latitude, so any distance or area result off a geographic CRS is wrong by four to five orders of magnitude. Guard every metric call with assert gdf.crs.is_projected. This is also why kriging and interpolation must run in UTM, as enforced in interpolating sparse yield monitor data with kriging.

NAD83 relabelled as WGS84. A field boundary digitised in NAD83 (EPSG:4269) but tagged EPSG:4326 in its metadata carries a hidden 1–2 m offset. Because both are geographic degrees, every CRS-equality check passes — the codes even look almost interchangeable to careless code. The only detection is a control-point check (Step 5) or noticing that RTK-logged operations consistently miss the boundary by a fixed metre-scale amount. Transform, don’t relabel, and read resolving pyproj datum shift warnings for the grid setup.

UTM zone-boundary crossing. A field or a farm straddling a 6°-wide UTM meridian (e.g. 90°W between zones 15N and 16N in the US Corn Belt) cannot be represented without distortion in a single zone — points far into the neighbouring zone accumulate scale error that grows to several metres near the far edge. Delivered as two tiles in two zones, they will not merge cleanly and produce a seam. Options: pick one zone and accept bounded distortion for a field that only slightly overhangs, or adopt a regional equal-area CRS such as CONUS Albers (EPSG:5070) for multi-zone operations. Validate the choice with validating coordinate systems for variable-rate maps.

The axis mapping strategy gotcha in mixed toolchains. rasterio and shapely are always lon/lat; raw pyproj.Transformer obeys authority axis order unless always_xy=True; GDAL command-line tools depend on the OAMS_TRADITIONAL_GIS_ORDER setting. A pipeline that hands coordinates between all three can swap axes at exactly one hop. Standardise on always_xy=True everywhere pyproj is constructed directly, and never hand-build a transformer without it.

Float32 coordinate truncation in projected CRS. UTM eastings and northings are six- and seven-digit numbers (e.g. 448236.0, 4651890.0). Stored as float32 they lose ~0.5 m of precision — invisible for yield data, fatal for RTK guidance. Keep projected coordinates in float64 end-to-end; only cast raster values (not coordinates) to float32.

5. Verification & Output Validation

Correctness is confirmed empirically, not by the absence of exceptions. Run this consolidated check after any reprojection or before any cross-layer operation.

PYTHON 
import numpy as np
import rasterio
import geopandas as gpd
from pyproj import CRS

def validate_alignment(raster_path: str, vector: gpd.GeoDataFrame) -> dict:
    """Confirm a raster and a vector actually occupy the same space."""
    with rasterio.open(raster_path) as src:
        assert src.crs is not None, "Raster CRS is None"
        assert vector.crs is not None, "Vector CRS is None"

        rc, vc = CRS.from_user_input(src.crs), CRS.from_user_input(vector.crs)
        assert rc.to_epsg() == vc.to_epsg(), (
            f"CRS mismatch: raster EPSG:{rc.to_epsg()} != vector EPSG:{vc.to_epsg()}"
        )
        assert rc.is_projected, "Work in a projected CRS for area/distance validity"

        # Bounding-box overlap is the empirical alignment test
        rb = src.bounds
        vb = vector.total_bounds  # (minx, miny, maxx, maxy)
        overlap_x = min(rb.right, vb[2]) - max(rb.left, vb[0])
        overlap_y = min(rb.top, vb[3]) - max(rb.bottom, vb[1])
        assert overlap_x > 0 and overlap_y > 0, (
            "Bounding boxes do not overlap — layers are in the same CRS code but "
            "different space; check for a datum-relabel or a bad reprojection"
        )
        return {
            "epsg": rc.to_epsg(),
            "projected": rc.is_projected,
            "overlap_m2": float(overlap_x * overlap_y),
        }

fields = gpd.read_file("field_boundaries.gpkg")
with rasterio.open("ndvi_20240715.tif") as src:
    fields = fields.to_crs(src.crs)
print(validate_alignment("ndvi_20240715.tif", fields))

A passing result guarantees three things simultaneously: neither CRS is missing, the codes agree, and the layers physically overlap. A layer pair that agrees on EPSG code but fails the overlap test is the fingerprint of a datum relabel or a reprojection that used the wrong source CRS — precisely the bugs that pass a naive equality check.

For a visual cross-check, load both layers into QGIS with the project CRS forced to the working EPSG, and confirm the vector boundary traces the raster’s field edges. A constant, direction-consistent offset of one to two metres between an RTK layer and its boundary is the visual signature of an unresolved datum shift.

6. Integration with the Pipeline

CRS debugging is not a stage; it is a guard clause you insert at every seam of the broader ag-GIS pipeline.

At ingest. The moment drone or satellite imagery lands, confirm and record its CRS — the ingesting multispectral drone imagery and orthomosaic stitching workflows both assume a known, projected CRS before any radiometric work begins. A missing or geographic CRS at ingest cascades into every downstream error catalogued here.

Before masking and index math. Every rasterio.mask and band-math operation sits behind the Step 4 assert_same_crs guard. The empty-clip trap is the most common way a clean-looking band math and raster algebra result turns out to be all-nodata; clipping rasters to field boundaries with rasterio.mask applies the same guard in context.

Before prescription export. A prescription written in the wrong CRS or with a residual datum shift steers real equipment off-target. Validate the CRS one final time before exporting prescription maps to John Deere GreenStar format, where a 1–2 m offset becomes a physically mis-applied input rate.

Two focused deep-dives. The two most involved failures each have a dedicated guide: fixing rasterio CRS mismatch errors for the raster/vector empty-clip and WarpedVRT reprojection path, and resolving pyproj datum shift warnings for the ballpark-transformation warning and PROJ grid installation.

Frequently Asked Questions

Why is my rasterio mask returning an empty array instead of raising an error?

The raster and the clip geometry are in different coordinate systems, so their coordinates do not overlap and the intersection is empty. Most masking functions treat a non-overlapping geometry as a valid zero-area clip rather than an error. Assert that the geometry CRS equals the raster CRS before masking, and reproject the geometry to match if they differ.

Does the difference between NAD83 and WGS84 matter for farm mapping?

For sub-metre RTK work it does. NAD83 and WGS84 have drifted apart by roughly one to two metres in the continental United States, which exceeds the one to three centimetre accuracy of RTK GPS used for guidance and planting. Treating the two datums as identical shifts every planted row and every prescription cell by that offset, so transform between them explicitly with pyproj instead of relabelling the CRS.

Why do my coordinates come out swapped as latitude and longitude?

EPSG:4326 formally defines its axis order as latitude then longitude, but most file formats and libraries store longitude then latitude. When pyproj follows the authority axis order it returns coordinates in the opposite order to what GeoJSON and shapefiles expect. Construct transformers with always_xy set to True so pyproj always uses longitude then latitude.

Parent section Ag-GIS Data Fundamentals & Spatial Reference Systems Related Field Boundary Extraction with GeoPandas Related Geospatial File Formats for Farm Data Pipelines Related Ingesting Multispectral Drone Imagery in Python Related Orthomosaic Stitching for Agricultural Drone Imagery Related RTK GPS Accuracy & Positional Error Budgets Related Understanding CRS in Precision Agriculture