Forestry GIS Python · Forestry · Ecology

How to Fix CRS Mismatches in geopandas with set_crs and to_crs

Coordinate Reference System (CRS) mismatches are the most common failure vector in spatial Python workflows, particularly when integrating field-collected plot inventories with agency-provided administrative boundaries or multispectral raster stacks. When projection metadata diverges, geometries fail to align, polygon area calculations return distorted values, spatial joins silently drop records, and raster extractions sample the wrong pixels. Fixing the mismatch comes down to one decision — does a layer need a label (set_crs) or a transform (to_crs) — followed by a bounds check. This page covers that narrow, high-frequency sub-task. It sits within the Coordinate Reference Systems for Forestry workflow, which is part of the broader Ecological GIS Data Foundations in Python framework for reproducible conservation analysis.

When to Use set_crs vs to_crs

The single decision that resolves almost every mismatch is whether the coordinate values in the GeoDataFrame already represent the system you want. If they do but the metadata is missing or wrong, you only need to attach a label. If the values themselves are in the wrong system, you must mathematically move them. Choosing wrong is the classic way to displace an entire dataset by hundreds of metres without raising a single error.

Situation .crs state Correct method What it does
Coordinates already correct, metadata absent None / empty set_crs(epsg) Attaches a definition; geometry values unchanged
Coordinates already correct, metadata wrong wrong EPSG set_crs(epsg, allow_override=True) Overwrites the label; no points move
Coordinates in the wrong projection known but wrong to_crs(target_epsg) Reprojects points through a pyproj pipeline
Untagged data that also needs reprojecting None set_crs then to_crs Label the source first, then transform

GeoPandas delegates projection handling to the pyproj library and stores the definition in the .crs attribute. A mismatch typically manifests when two GeoDataFrame objects report different EPSG codes, or when one returns None. Legacy forestry shapefiles exported from desktop GIS, or raw CSV exports of GPS waypoints, frequently lack embedded projection metadata entirely.

import geopandas as gpd

plots = gpd.read_file("field_plots.shp")
boundaries = gpd.read_file("provincial_forest_zones.gpkg")

# Inspect CRS metadata before doing anything spatial
print(plots.crs)                  # may be None
print(boundaries.crs.to_epsg())   # e.g. 3005

If plots.crs evaluates to None, the dataset is geographically unanchored: every overlay, distance, or area call against it is meaningless until you attach a definition.

Decision flow for resolving a CRS mismatch in GeoPandas A flowchart starting from a loaded GeoDataFrame. First decision: does the .crs attribute exist? If it is None or empty, use set_crs only when the coordinate values are already correct, otherwise label then transform. If a CRS is present but wrong, use to_crs to mathematically reproject. Both correct paths converge on validating total_bounds. A warning box notes that chaining set_crs and to_crs without checking the initial state double-transforms and corrupts geometry. Loaded GeoDataFrame inspect gdf.crs Is .crs present and correct? None / missing Coordinates already in target system? yes .set_crs(epsg, ...) label only — no math no present, wrong .to_crs(target_epsg) reproject coordinates via pyproj Validate alignment total_bounds match? Do NOT chain set_crs + to_crs on unchecked data — double-transform corrupts.
Choosing set_crs (attach a label without moving points) versus to_crs (mathematically reproject), then converging on a bounds check.

Minimal Reproducible Example

The snippet below covers the three situations end to end: label an untagged layer, reproject a layer that is in the wrong projection, and converge on a shared target CRS so the two datasets can be overlaid. For forestry work, settle on a regional equal-area projection (e.g. EPSG:3005 BC Albers, or EPSG:5070 for the contiguous US) so .area stays accurate for biomass and canopy metrics.

import geopandas as gpd

TARGET = "EPSG:3005"  # BC Albers — equal-area, metres

plots = gpd.read_file("field_plots.shp")
boundaries = gpd.read_file("provincial_forest_zones.gpkg")

# 1. Untagged plots whose coordinates are already WGS84 lon/lat:
#    attach the label only — do NOT transform.
if plots.crs is None:
    plots = plots.set_crs("EPSG:4326")

# 2. Now both layers have a known source CRS. Reproject each to TARGET.
plots = plots.to_crs(TARGET)
boundaries = boundaries.to_crs(TARGET)

# 3. They now share one projected CRS in metres and can be overlaid.
assert plots.crs == boundaries.crs
joined = gpd.sjoin(plots, boundaries, how="left", predicate="within")

The two methods do fundamentally different things. set_crs() only writes the .crs attribute — the (x, y) values are untouched — so it is correct only when the coordinates already match the system you are naming. to_crs() invokes a pyproj transformation pipeline that accounts for ellipsoid parameters, datum shifts, and grid interpolation, so the coordinate values change. Use allow_override=True with set_crs() only when you are deliberately replacing an incorrect existing label:

# Coordinates are already in WGS84 but the file carries a wrong EPSG tag
plots = plots.set_crs("EPSG:4326", allow_override=True)
What set_crs and to_crs do to coordinate values versus the CRS label Two panels compare the effect of each method on a single forest plot point. In the set_crs panel the point stays at the same map location while only its CRS label changes from None to EPSG:4326 — the (x, y) values are identical before and after. In the to_crs panel the CRS label changes from EPSG:4326 to EPSG:3005 and the point physically moves to a new location because pyproj recomputes the coordinates, with the values changing from geographic degrees to projected metres. .set_crs("EPSG:4326") relabel only — point does not move same (x, y) .crs: None → EPSG:4326 (values unchanged) .to_crs("EPSG:3005") reproject — pyproj recomputes the point before new (x, y) degrees → metres (values change)
set_crs rewrites only the .crs label and leaves every coordinate exactly where it was; to_crs sends each point through pyproj, so the values — and the on-map position — change.

This CRS discipline is the non-negotiable prerequisite for every downstream operation — it is the same alignment step required before merging shapefiles and rasters in Python and before any zonal extraction.

Parameter Reference

The arguments below are the ones that change the outcome of a CRS repair. Defaults are tuned for the common case; the failure cases in forestry data come from the edges.

Argument Method Type Default Recommended Rationale
crs set_crs / to_crs str / int / CRS — EPSG string e.g. "EPSG:3005" Pass an explicit EPSG; never rely on an assumed default
allow_override set_crs bool False True only to replace a known-wrong tag False raises if a CRS already exists, guarding against accidental relabel
inplace set_crs / to_crs bool False False Keep the source layer immutable so a mis-step is recoverable
always_xy Transformer.from_crs bool False True Forces lon/lat (x,y) ordering, avoiding silent axis swaps
accuracy Transformer.from_crs float (m) None 0.05 for datum shifts Raises instead of falling back to a coarse Helmert approximation

Expected Output and Verification

After transformation, confirm spatial congruence before trusting any join, buffer, or raster sample. The cheapest reliable check is the bounding box: two layers covering the same ground must report overlapping total_bounds, and projected coordinates in metres are numerically far larger than geographic degrees.

print("Plots bounds:    ", plots.total_bounds)
print("Boundaries bounds:", boundaries.total_bounds)

# Geographic degrees stay within [-180, 180] / [-90, 90];
# projected metres are much larger. A mix means a transform was skipped.
assert not plots.crs.is_geographic
assert plots.crs == boundaries.crs

# The two extents must actually overlap, not sit oceans apart.
px0, py0, px1, py1 = plots.total_bounds
bx0, by0, bx1, by1 = boundaries.total_bounds
assert px0 < bx1 and bx0 < px1 and py0 < by1 and by0 < py1

If bounds differ by orders of magnitude — [-180, -90, 180, 90] against [300000, 500000, 400000, 600000] — a transformation was skipped or applied to only one layer. For raster–vector work, also assert src.crs.to_epsg() == gdf.crs.to_epsg() before any sampling call.

Legacy PROJ strings and datum-grid accuracy

Older datasets may carry deprecated +init=epsg:XXXX strings, which modern pyproj warns about. Normalize them to a clean EPSG code:

from pyproj import CRS

legacy = CRS.from_string("+init=epsg:26910")   # NAD83 / UTM zone 10N
plots = plots.set_crs(legacy.to_epsg(), allow_override=True)

Transformations between NAD27, NAD83, and WGS84 need grid-shift files (.gsb). If pyproj cannot find them it silently falls back to a low-accuracy Helmert approximation, displacing geometry by a metre or more. Enable network access so pyproj can fetch official grids, or install the proj-data conda package to bundle them locally:

import pyproj
pyproj.network.set_network_enabled(active=True)

To require (rather than tolerate) the grid shift, build an explicit pyproj.Transformer with an accuracy budget. Setting accuracy=0.05 (metres) makes pyproj raise instead of approximating, so silent datum failures surface at development time:

from pyproj import Transformer
from shapely.ops import transform

tfm = Transformer.from_crs("EPSG:26910", "EPSG:3005",
                           always_xy=True, accuracy=0.05)
plots["geometry"] = plots.geometry.apply(lambda g: transform(tfm.transform, g))
plots = plots.set_crs("EPSG:3005", allow_override=True)

Detailed pipeline configuration lives in the official pyproj documentation.

Common Pitfalls

  • Using set_crs where to_crs is needed. Labelling UTM coordinates as EPSG:4326 does not convert them — it asserts the wrong truth and silently misplaces every feature by thousands of kilometres. Confirm the coordinate magnitudes match the system before labelling.
  • Double-transforming. Chaining set_crs() then to_crs() on data whose initial state you did not check reprojects already-correct coordinates a second time, irreversibly corrupting topology. Inspect gdf.crs first; only label untagged data.
  • Geographic CRS leaking into .area. Calling .area on an EPSG:4326 layer returns squared degrees, not square metres. Assert not gdf.crs.is_geographic before any area or distance metric.
  • Silent Helmert fallback. A NAD83↔WGS84 transform without datum grids quietly degrades to a coarse approximation. Use an explicit Transformer with an accuracy budget to force a failure you can see.

Frequently Asked Questions

What is the difference between set_crs and to_crs in geopandas?

set_crs() only writes the .crs metadata attribute and leaves the coordinate values untouched — use it when the points are already in the target system but the label is missing or wrong. to_crs() mathematically reprojects the coordinates through a pyproj pipeline, so the (x, y) values change. Picking the wrong one silently displaces your data instead of raising an error.

Why does my GeoDataFrame have a .crs of None?

The source file carried no projection metadata — common with shapefiles missing their .prj sidecar, GeoPackages built from raw coordinates, or CSV/GPS exports. The coordinate values still exist; they are just unlabelled. Determine the true system from the data provider, then attach it with set_crs() (label only) or follow it with to_crs() if the values also need reprojecting.

How do I know whether a transform actually happened?

Compare total_bounds before and after, and check crs.is_geographic. Geographic coordinates stay within roughly [-180, 180] longitude and [-90, 90] latitude; a projected CRS in metres produces values in the hundreds of thousands. If the magnitudes did not change, to_crs() was a no-op (the source and target were identical) or you accidentally called set_crs() instead.

Why are my reprojected coordinates off by about a metre?

pyproj could not find the datum-grid (.gsb) files for a NAD27/NAD83/WGS84 shift and fell back to a coarse Helmert approximation. Enable pyproj.network.set_network_enabled(active=True) or install the proj-data package, and build a Transformer with accuracy=0.05 so the missing grid raises an exception instead of degrading silently.