Coordinate Reference System Transformations in Modern Python Workflows

Accurate spatial analysis depends entirely on correct Coordinate Reference System (CRS) transformations. Misaligned projections introduce metric distortion, break spatial relationships, and invalidate downstream analytics. This stage of Geospatial Data Ingestion & Processing Workflows is where every incoming layer is normalized to one authoritative projection before it reaches Spatial Joins & Merging or Topology Validation & Repair. It covers production CRS management with pyproj, GeoPandas, and Rasterio, from EPSG initialization through raster warping and batch coordinate conversion. For layers too large to reproject in RAM, the deep dive is Reprojecting Large Datasets Without Memory Errors.

Architecture & Data Structures

A CRS is not a single number; it is a stack of definitions. At the base sits a datum (an ellipsoid plus a reference frame such as WGS84 or ETRS89), on top of which a projection maps the curved surface to a plane, and around which a coordinate system fixes axis order and units. pyproj.CRS is the object that carries all three, wrapping the underlying PROJ database. GeoPandas and Rasterio both delegate to it: a GeoDataFrame.crs attribute and a raster dataset's .crs are pyproj.CRS instances, which is why a single mental model covers vector and raster alike.

The one rule that prevents most production incidents is to identify every CRS by its EPSG authority code rather than a hand-written PROJ string. Authority codes are unambiguous, versioned, and carry the correct axis metadata; free-form strings drift.

from pyproj import CRS

# Authoritative initialization — an EPSG code resolves the full definition
target_crs = CRS.from_epsg(25832)          # ETRS89 / UTM zone 32N, metres
assert target_crs.is_projected, "Metric analysis needs a projected CRS"
assert target_crs.axis_info[0].unit_name == "metre"

print(target_crs.name)                      # ETRS89 / UTM zone 32N
print(target_crs.to_authority())            # ('EPSG', '25832')

The is_projected, is_geographic, and axis_info attributes are the fields you interrogate before trusting any measurement. A geographic CRS (degrees) will happily return an "area" from .area, but the number is meaningless — an assertion at ingestion is cheaper than a silent error downstream.

A second structural fact catches people out: a pair of CRS objects does not imply one answer. When PROJ is asked to go from one to the other it searches proj.db for every registered coordinate operation linking the two datums, filters them by area of use, ranks the survivors by declared accuracy, and hands back one. The Transformer is therefore a third object with its own identity, its own accuracy figure and its own validity envelope — not a function mechanically derived from its endpoints. Which one you get depends on what grid files the machine has and, if you supply it, on where your data actually sits.

from pyproj import Transformer
from pyproj.aoi import AreaOfInterest

# Constrain the search to the extent the layer actually covers
HAMBURG = AreaOfInterest(
    west_lon_degree=9.6, south_lat_degree=53.3,
    east_lon_degree=10.4, north_lat_degree=53.8,
)

to_utm = Transformer.from_crs(
    "EPSG:4326", "EPSG:25832", always_xy=True, area_of_interest=HAMBURG
)
print(to_utm.description)          # the operation PROJ actually selected
print(to_utm.accuracy)             # declared metres of error, or -1 when unknown
print(to_utm.is_network_enabled)   # False here means grids must already be local

accuracy is the field worth logging next to every batch you ship: a pipeline that quietly falls back from a one-centimetre grid-based operation to a one-metre Helmert approximation produces output that still passes every assertion in this article. How PROJ ranks candidates, and how to force it to fail rather than degrade, is the subject of Datum Shifts and Transformation Grids in PyProj.

The third fact concerns height. Almost every EPSG code you will type — 4326, 25832, 27700 — is two-dimensional. Hand a third z array to a 2D transformer and the heights come back byte-identical, with no error and no warning, which is exactly wrong when the source carries ellipsoidal heights and the consumer expects orthometric ones. If elevation means anything in your pipeline, build a compound CRS so the vertical axis is transformed too.

from pyproj import CRS, Transformer
from pyproj.crs import CompoundCRS

# 2D target: the z value simply passes through untouched
flat = Transformer.from_crs("EPSG:4979", "EPSG:25832", always_xy=True)
print(flat.transform(9.9937, 53.5511, 42.0)[2])   # 42.0 — unchanged, silently

# 3D target: ellipsoidal height -> EGM2008 orthometric height (uses a geoid grid)
utm32_egm2008 = CompoundCRS(
    name="ETRS89 / UTM zone 32N + EGM2008 height",
    components=["EPSG:25832", "EPSG:3855"],
)
tall = Transformer.from_crs("EPSG:4979", utm32_egm2008, always_xy=True)
easting, northing, ortho_h = tall.transform(9.9937, 53.5511, 42.0)
print(round(ortho_h, 2))   # differs from 42.0 by the local geoid separation

CRS.to_3d() promotes a 2D authority code to its three-dimensional sibling where one exists, which is the cheaper route when both sides reference the same ellipsoid. What it cannot do is invent a vertical datum: a z column that survived a 2D GeoDataFrame.to_crs() is still in whatever vertical frame it arrived in, whatever the horizontal label now says.

Environment Configuration & Dependency Resolution

CRS transformations lean on PROJ, the C library that pyproj binds. Version skew between pyproj, the bundled PROJ, and the GDAL that GeoPandas and Rasterio compile against is the single most common cause of "works on my machine" reprojection bugs. Pin the stack from conda-forge, which builds all three against one PROJ.

conda install -c conda-forge \
  "pyproj>=3.6" "geopandas>=0.14" "rasterio>=1.3" "shapely>=2.0"

The legacy +init=epsg: PROJ string syntax was removed in PROJ 6, and any tutorial that still uses it will emit deprecation errors or silently swap axis order — replace it with CRS.from_epsg(). In containers, enable the PROJ network so datum-shift grids download on demand rather than failing mid-pipeline.

import pyproj

# CI / container sanity check — log the resolved environment
print("pyproj:", pyproj.__version__)
print("PROJ:", pyproj.proj_version_str)
print("grids dir:", pyproj.datadir.get_data_dir())

# Enable on-demand transformation-grid downloads (set PROJ_NETWORK=ON in prod)
pyproj.network.set_network_enabled(True)

Four environment variables decide where PROJ looks and whether it may write. PROJ_DATA replaced PROJ_LIB in PROJ 9.1 — both are honoured for now, but a container that sets only the old name against a new build will resolve proj.db by luck rather than by configuration. PROJ_NETWORK gates CDN downloads. PROJ_USER_WRITABLE_DIRECTORY is the one people discover the hard way: with the network enabled PROJ caches downloaded grids under the user's data directory, and on a read-only or ephemeral filesystem — a serverless function, a hardened container, a Spark executor — that write fails and every request re-downloads or silently degrades.

# Deterministic CRS behaviour inside a container image
export PROJ_DATA=/opt/conda/share/proj          # PROJ >= 9.1; was PROJ_LIB
export GDAL_DATA=/opt/conda/share/gdal
export PROJ_NETWORK=ON
export PROJ_USER_WRITABLE_DIRECTORY=/tmp/proj   # writable even on a read-only FS

Air-gapped and cost-sensitive deployments should invert that arrangement: leave the network off and bake the grids your area of operation needs into the image, so the operation PROJ selects is a property of the build rather than of the runtime's connectivity. Whichever way you go, assert it at startup rather than discovering it from a coordinate that looks a metre off.

import pyproj
from pyproj import Transformer

REQUIRED_ACCURACY_M = 0.1

def assert_transform_quality(src: str, dst: str) -> None:
    """Fail the boot, not the batch, when the environment cannot meet the spec."""
    tf = Transformer.from_crs(src, dst, always_xy=True)
    known = tf.accuracy is not None and tf.accuracy >= 0    # -1 means "not stated"
    if known and tf.accuracy > REQUIRED_ACCURACY_M:
        raise RuntimeError(
            f"{src} -> {dst} resolved to a {tf.accuracy} m operation "
            f"({tf.description}); required <= {REQUIRED_ACCURACY_M} m. "
            "A transformation grid is probably missing."
        )

assert_transform_quality("EPSG:4326", "EPSG:25832")

Concurrency has its own rules. A Transformer wraps a PROJ context that is not fork-safe, so build it inside each worker after the fork rather than constructing it in the parent and inheriting it through multiprocessing; the symptom of getting this wrong is a segfault or a stream of identical wrong coordinates rather than an exception. PyProj memoizes Transformer.from_crs on its argument tuple, so repeating the identical call is nearly free — but passing a freshly-generated WKT string each time produces a cache miss every call and rebuilds the pipeline from proj.db.

Vectorized Operations & Core Workflow

The core operation is GeoDataFrame.to_crs(), which reprojects an entire geometry column in one vectorized pass. The canonical ingestion workflow is: read, backfill a CRS only if the source lacks one, then reproject to the project's single target. Never assign a CRS with set_crs to fix misalignment — that relabels coordinates without moving them; only to_crs actually transforms them.

import geopandas as gpd

TARGET_EPSG = 25832  # project-wide metric CRS (ETRS89 / UTM 32N)

parcels = gpd.read_file("parcels.shp")

# A missing .prj gives crs=None — backfill the KNOWN source CRS, don't guess
if parcels.crs is None:
    parcels = parcels.set_crs("EPSG:4326")   # relabels only, no coordinates move

# Reproject the whole geometry column to the project target
parcels_metric = parcels.to_crs(epsg=TARGET_EPSG)

assert parcels_metric.crs.equals(CRS.from_epsg(TARGET_EPSG))
parcels_metric["area_m2"] = parcels_metric.area   # now metric and meaningful

For raw coordinate arrays or point clouds that live outside a GeoDataFrame, pyproj.Transformer is faster and avoids per-row Python overhead. Build the transformer once and reuse it — pyproj caches the compiled transformation pipeline, so re-instantiating inside a loop throws that cache away.

from pyproj import Transformer
import numpy as np

# always_xy=True forces (lon, lat) / (easting, northing) input-output order
to_utm = Transformer.from_crs("EPSG:4326", "EPSG:25832", always_xy=True)

lon = np.array([9.9937, 10.0014])   # Hamburg-area sensor longitudes
lat = np.array([53.5511, 53.5600])
easting, northing = to_utm.transform(lon, lat)   # vectorized over the arrays

A transform that runs off the edge of its projection does not raise — it returns inf. Transverse Mercator diverges as you move away from the central meridian, so a stray point from another continent that slipped into a national extract comes back as infinity, propagates into total_bounds, and turns every subsequent bounding-box filter into an empty result. Screen for it explicitly, or pass errcheck=True to make PROJ raise instead.

import numpy as np
from pyproj import Transformer
from pyproj.enums import TransformDirection

to_utm = Transformer.from_crs("EPSG:4326", "EPSG:25832", always_xy=True)

lon = np.array([9.9937, 10.0014, -152.4])   # third row: a mis-keyed sensor record
lat = np.array([53.5511, 53.5600, 61.2])

easting, northing = to_utm.transform(lon, lat)
off_grid = ~np.isfinite(easting)
print(off_grid.sum(), "coordinates fell outside the projection")   # 1

# The same object runs backwards — never build a second transformer for the inverse
lon_back, lat_back = to_utm.transform(
    easting[~off_grid], northing[~off_grid], direction=TransformDirection.INVERSE
)
assert np.allclose(lon_back, lon[~off_grid], atol=1e-9)

Performance is dominated by two very different costs. Building the transformer means a proj.db lookup and pipeline compilation — order of milliseconds, which is nothing once but ruinous inside a per-row loop. Executing it is a tight C loop over a coordinate array: on commodity hardware expect roughly a million coordinates per second for a projection-only conversion, and several times slower when a grid-based datum shift sits in the pipeline and every point costs an interpolation into a raster of offsets. GeoDataFrame.to_crs() adds the cost of unpacking each geometry's coordinate sequence and rebuilding it, so a layer of vertex-dense polygons reprojects far slower per feature than a point layer of the same row count — vertices, not rows, are the unit of work.

For datasets that exceed available RAM, partition the reprojection with dask-geopandas or stream it in record batches — the dedicated method is covered in Reprojecting Large Datasets Without Memory Errors.

Geometry & Data Processing Details

Two subtleties bite even experienced practitioners: axis order and datum shifts.

Axis order. EPSG defines geographic CRS such as EPSG:4326 with latitude first. PROJ honours that authority definition, so without always_xy=True a Transformer expects (lat, lon) and returns (lat, lon). Most Python code, GeoJSON (RFC 7946), and web maps assume (lon, lat). Setting always_xy=True on every Transformer normalizes to longitude-first and removes the ambiguity.

from pyproj import CRS, Transformer

wgs84 = CRS.from_epsg(4326)
print(wgs84.axis_info[0].direction)   # 'north'  → latitude is axis 0

# WITHOUT always_xy this silently swaps the coordinates:
naive = Transformer.from_crs(4326, 25832)
print(naive.transform(53.55, 9.99))   # interprets 53.55 as LONGITUDE — wrong

# WITH always_xy the (lon, lat) you pass is the (lon, lat) it uses:
safe = Transformer.from_crs(4326, 25832, always_xy=True)
print(safe.transform(9.99, 53.55))    # correct easting, northing

Datum shifts. Moving between datums (for example legacy NAD27 to WGS84, or a national grid to ETRS89) is not a formula but a grid-based correction. PROJ downloads the correct shift grid when the network is enabled; without it you get metre-level errors that pass every syntactic check. Define such transforms explicitly and confirm a residual on a known control point.

from pyproj import Transformer

# Explicit legacy → modern datum transform (grid-based under the hood)
nad27_to_wgs84 = Transformer.from_crs(
    "EPSG:4267",   # NAD27
    "EPSG:4326",   # WGS84
    always_xy=True,
)
lon, lat = nad27_to_wgs84.transform(-77.0369, 38.9072)   # Washington, DC
print(f"{lon:.6f}, {lat:.6f}")   # shifted by ~tens of metres vs. a naive copy

Edges are not transformed, only vertices are. Every reprojection maps the points of a geometry and then reconnects them with straight segments in the destination plane. That is fine for a building footprint, where the vertices are metres apart, and quietly wrong for anything long: a boundary drawn as a straight line between two monuments 300 km apart in UTM is a curve in geographic coordinates, and a two-vertex representation of it cuts the corner. The fix is to densify before transforming, not after — shapely.segmentize inserts intermediate vertices at a chosen spacing in the source CRS, and each of those then lands on the true curve.

import geopandas as gpd
import shapely
from shapely.geometry import LineString
from pyproj import Geod

# A survey baseline drawn straight in ETRS89 / UTM 32N
baseline = gpd.GeoSeries(
    [LineString([(400000, 5600000), (700000, 5900000)])], crs="EPSG:25832"
)

naive = baseline.to_crs("EPSG:4326")
dense = gpd.GeoSeries(
    shapely.segmentize(baseline.values, max_segment_length=1000),  # 1 km spacing
    crs="EPSG:25832",
).to_crs("EPSG:4326")

geod = Geod(ellps="WGS84")
a = naive.iloc[0].interpolate(0.5, normalized=True)
b = dense.iloc[0].interpolate(0.5, normalized=True)
_, _, gap_m = geod.inv(a.x, a.y, b.x, b.y)
print(f"midpoint divergence: {gap_m:.1f} m")   # a three-figure number of metres

Choose max_segment_length from the accuracy you owe: a spacing that keeps the sagitta of each segment under your tolerance. For continental polygons and flight or shipping tracks that means kilometres; for cadastral parcels densification is wasted work, because the vertices are already closer together than the curvature matters over. Note the ordering — densifying after the transform adds vertices to the wrong curve and fixes nothing.

The antimeridian and the poles break rings, not coordinates. Reproject a Pacific-spanning layer into EPSG:4326 and the individual longitudes are all correct while the rings that connect them are not: a polygon whose vertices run from 179.5° to −179.5° is interpreted as spanning the other 359 degrees of the globe, so it renders as a band across the entire map and its total_bounds reports the whole world. RFC 7946 requires such geometry to be cut at ±180° into two parts before serialization, so split it at the export boundary rather than hoping the renderer copes. Polar data has the mirror problem: a polygon that encloses a pole has no valid geographic ring at all, and needs a polar stereographic CRS such as EPSG:3413 or EPSG:3031 for any analysis.

Missing geometry passes through untouched. to_crs() transforms the geometries it can and leaves None entries as None, and empty geometries come back empty with the new CRS label attached. Neither raises. If your validation gate asserts on feature counts alone it will not notice, so check notna() and is_empty alongside the CRS assertion — the pruning pass in Automating Shapefile Cleanup with Python belongs before the reprojection, not after it.

CRS Alignment & Projection Pipeline

The pipeline that keeps a whole project consistent has three gates. First, an ingestion gate that detects or backfills the source CRS and logs anything ambiguous. Second, a reprojection step to the one project-wide target. Third, a validation gate that confirms units and axis order before any metric operation runs. Encapsulating this makes every incoming layer trustworthy by the time it reaches a spatial join.

import logging
import geopandas as gpd
from pyproj import CRS

logger = logging.getLogger("crs")
TARGET = CRS.from_epsg(25832)   # project-wide metric target

def normalize_crs(gdf: gpd.GeoDataFrame, fallback_epsg: int | None = None):
    """Detect, backfill and reproject a layer to the project target CRS."""
    if gdf.crs is None:
        if fallback_epsg is None:
            raise ValueError("Layer has no CRS and no fallback provided")
        logger.warning("No CRS on source — backfilling EPSG:%s", fallback_epsg)
        gdf = gdf.set_crs(epsg=fallback_epsg)

    if not gdf.crs.equals(TARGET):
        logger.info("Reprojecting %s → %s", gdf.crs.to_authority(),
                    TARGET.to_authority())
        gdf = gdf.to_crs(TARGET)

    # Validation gate: metric analysis MUST be projected, in metres
    assert gdf.crs.is_projected, "Target CRS is not projected"
    assert gdf.crs.axis_info[0].unit_name == "metre"
    return gdf

For raster stacks the alignment step is rasterio.warp.reproject. Unlike vectors, rasters need a resampling algorithm, because reprojection resamples the pixel grid: use Resampling.nearest for categorical layers (land cover, masks) and Resampling.bilinear or cubic for continuous surfaces (elevation, reflectance).

import rasterio
from rasterio.warp import reproject, calculate_default_transform, Resampling

DST_CRS = "EPSG:25832"

with rasterio.open("ortho.tif") 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("ortho_utm.tif", "w", **profile) as dst:
        for band in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, band),
                destination=rasterio.band(dst, band),
                src_transform=src.transform, src_crs=src.crs,
                dst_transform=transform, dst_crs=DST_CRS,
                resampling=Resampling.bilinear,   # continuous imagery
            )

calculate_default_transform derives a pixel size from the source extent, which yields awkward numbers like 0.4372 m and makes two scenes warped separately land on two different grids. Pass resolution= to snap the destination to a grid you chose, and set the nodata values explicitly so the areas outside the rotated footprint are flagged rather than filled with zeros that later average into your statistics.

from rasterio.warp import calculate_default_transform

with rasterio.open("ortho.tif") as src:
    transform, width, height = calculate_default_transform(
        src.crs, DST_CRS, src.width, src.height, *src.bounds,
        resolution=(0.5, 0.5),        # a grid you picked, not one PROJ derived
    )

Four reproject() arguments are worth setting deliberately on every warp: resolution (or an explicit dst_transform) to fix the grid, src_nodata/dst_nodata so the collar is masked instead of zero-filled, num_threads to use more than one core on large scenes, and warp_mem_limit to raise GDAL's default working buffer when the warp thrashes. When several rasters must line up pixel-for-pixel, do not warp them independently and hope — reproject each onto the first one's exact transform, or use the reproject_match pattern described in Reprojecting Raster Cubes with reproject_match.

The cheapest alignment is frequently no warp at all. Reprojecting a vector layer costs a few million coordinate transforms and is lossless; warping a raster resamples every pixel and permanently degrades it, and resampling twice degrades it twice. So when the job is sampling a raster under polygons, move the vectors into the raster's CRS — the direction taken in Zonal Statistics & Raster Sampling. Warp the raster only when it genuinely has to be co-registered with other rasters or served as tiles.

Production Export & Integration

Once a layer is projected and validated, the export format should preserve the CRS as metadata rather than force a lossy conversion. GeoParquet embeds the full CRS in its metadata and supports predicate pushdown, making it the preferred intermediate for downstream Cloud-Native Geospatial Formats and for analytics engines such as DuckDB spatial. When loading into a spatial database via PostGIS integration, the SRID must match the reprojected CRS or the server rejects the geometry.

Web mapping is the one place a geographic or Web Mercator CRS is correct: tile pipelines expect EPSG:4326 for geographic coordinates and EPSG:3857 for rendered tiles. Reproject to those standards only at the serialization boundary — never run metric analysis in them.

# Metric intermediate for analytics / PostGIS — CRS preserved in metadata
parcels_metric.to_parquet("parcels_utm.parquet")

# Serialization boundary: reproject to WGS84 only for the web / GeoJSON output
parcels_wgs84 = parcels_metric.to_crs(epsg=4326)
parcels_wgs84.to_file("parcels.geojson", driver="GeoJSON")   # RFC 7946, lon/lat

How much of the CRS survives export varies sharply by format, and that determines what the next reader can do. GeoParquet stores the full definition as PROJJSON in the file's metadata, so an authority code, its area of use and its datum all round-trip. GeoPackage stores WKT in a gpkg_spatial_ref_sys table and round-trips nearly as well. A Shapefile .prj holds WKT1, which keeps the datum name — enough for PROJ to still choose the right operation — but drops the usage and area-of-use blocks, and cannot express a compound or 3D CRS at all; the trade-offs between the dialects are unpacked in EPSG vs PROJ String vs WKT CRS Formats.

One naming detail causes recurring confusion at the web boundary. GeoJSON coordinates are longitude-latitude, which is OGC:CRS84, not EPSG:4326 — the two describe the same datum with opposite declared axis order. Writers that label GeoJSON output EPSG:4326 are relying on every reader applying the same longitude-first convention the spec mandates, which is true in practice and false in principle; the practical consequence is that a strict validator or a database import may transpose your coordinates. Where a CRS must be named alongside lon/lat data, OGC:CRS84 is the honest label.

from pyproj import CRS

crs84 = CRS.from_user_input("OGC:CRS84")
print(crs84.axis_info[0].direction)      # 'east'  → longitude first
print(CRS.from_epsg(4326).axis_info[0].direction)   # 'north' → latitude first
print(crs84.datum == CRS.from_epsg(4326).datum)     # True — same datum, other order

Loading into PostGIS raises the same question from the other side: reproject in Python before the insert, or store the source SRID and call ST_Transform in the query. Reproject on the way in when the table has a single analytical purpose — the geometry is then indexable in the units you measure in, and the GiST index actually helps. Keep the source SRID and transform per query only when several consumers need different targets and storage is cheaper than recomputation. What you must not do is mix SRIDs within one column: PostGIS enforces the declared SRID on insert and will reject the row rather than silently reproject it, which is the one place in this stack where a CRS mismatch is a loud error instead of a quiet one.

Production checklist

• Set PROJ_NETWORK=ON in every container so datum-shift grids download automatically.
• Reuse Transformer objects across calls — never instantiate them inside a loop.
• Keep intermediates as GeoParquet, not Shapefile; Shapefile does not round-trip a full CRS definition or long field names.
• Run a round-trip transform (WGS84 → UTM → WGS84) on control points and assert sub-millimetre residuals in CI to catch a misconfigured grid directory.
• Log pyproj.__version__ and pyproj.proj_version_str with each run for reproducibility.

Platform Edge Cases & Debugging

Most CRS failures are environmental, not logical, and fall into a few recurring causes:

• CRSError: Invalid projection or missing PROJ grids on Windows. The PROJ_DATA (older builds: PROJ_LIB) environment variable points at a stale directory, often when GDAL, Fiona, and pyproj were installed from mixed channels. Confirm the resolved path with pyproj.datadir.get_data_dir() and reinstall the stack from a single conda-forge environment.
• Coordinates land in the ocean off West Africa (0, 0-ish) or are transposed. Axis-order confusion — a Transformer built without always_xy=True, or GeoJSON read as (lat, lon). Rebuild transformers with always_xy=True.
• .prj-less Shapefile reads back as crs=None. GeoPandas cannot infer the CRS; set_crs the known authority code before anything else, as shown in Shapefile & GeoJSON Parsing.
• Areas or distances are off by a consistent factor. The geometry is still geographic (degrees) or in Web Mercator, where scale distorts with latitude. Reproject to a local UTM zone or national grid and assert crs.axis_info[0].unit_name == "metre" first.
• A sjoin silently returns no matches or wrong matches. The two operands carry different CRSs. Normalize both to the identical target before the join — see Spatial Joins & Merging.
• Results differ between your laptop and CI by a few decimetres. Different PROJ builds, so a different transformation grid was available on each. Log pyproj.proj_version_str and Transformer.description on both, pin the version in the lockfile, and ship the grids in the image rather than fetching them at runtime.
• Coordinates are inf after a transform, and every later bounding-box filter is empty. A point sat outside the projection's valid range — usually one mis-keyed record in an otherwise regional extract. Filter with np.isfinite immediately after the transform, or construct the transformer with errcheck=True so PROJ raises on the offending point.
• A segfault or frozen worker when reprojecting under multiprocessing. A Transformer was built in the parent and inherited across the fork. Construct it inside the worker, or switch the start method to spawn.
• The first transform in a container takes 30 seconds, then it is fast. PROJ_NETWORK=ON with an empty cache: PROJ is downloading a shift grid over HTTP on the first call. Bake the grids into the image, or point PROJ_USER_WRITABLE_DIRECTORY at a persistent volume so the download survives a restart.
• Elevations are unchanged after reprojecting a 3D layer. The target was a 2D CRS, so the z values passed straight through. Build a CompoundCRS with an explicit vertical component, and confirm the geoid grid is present.

Frequently Asked Questions

When do I use set_crs versus to_crs? set_crs only labels the geometry with a CRS — it moves no coordinates and is for backfilling a genuinely missing definition (a Shapefile with no .prj). to_crs actually reprojects the coordinates. If your layers are misaligned, you almost always want to_crs; reaching for set_crs to "fix" alignment relabels wrong data as right.

Why do I need always_xy=True everywhere? EPSG:4326 and many national CRS are defined latitude-first by their authority, and PROJ honours that. Python code, GeoJSON, and web maps assume longitude-first. always_xy=True normalizes every Transformer to (x, y) = (lon, lat), eliminating the most common silent coordinate swap.

Which projected CRS should I reproject to for analysis? A locally appropriate metric CRS: the correct UTM zone for your area, or a national grid such as ETRS89 / UTM (EPSG:25832) in Europe. Never do metric work in EPSG:4326 (degrees) or EPSG:3857 (Web Mercator distorts scale with latitude). Confirm the target reads metre before measuring.

How do I reproject a dataset that is larger than RAM? Stream it in record batches or partition it with dask-geopandas, building the transformer once and reusing it per chunk, so peak memory stays flat. The full recipe is Reprojecting Large Datasets Without Memory Errors.

My transform runs but the result is tens of metres off — what is wrong? That is a datum shift with no grid available. Enable the PROJ network (PROJ_NETWORK=ON) so the correct transformation grid downloads, and verify with a round-trip on a known control point that residuals are sub-millimetre.

Should I reproject in Python or let PostGIS do it with ST_Transform? Reproject in Python when the table serves one analytical purpose: the stored geometry is then already in the units you measure in, and the GiST index is built over the coordinates your queries actually use. Leave the source SRID in place and call ST_Transform per query only when several consumers need different targets. Reprojecting inside a WHERE clause is the case to avoid — a function applied to the indexed column defeats the index and forces a sequential scan of the whole table.

Does reprojecting a long line or a large polygon distort its shape? Yes, because only the vertices are transformed and the segments between them are redrawn straight in the destination plane. Over a few hundred metres the error is negligible; over hundreds of kilometres it can reach hundreds of metres at the midpoint of a segment. Densify with shapely.segmentize in the source CRS before calling to_crs(), choosing a spacing that keeps the deviation inside your tolerance.

How do I reproject elevation as well as position? A plain EPSG code is two-dimensional and passes z through unchanged. Compose a CompoundCRS from the horizontal code plus a vertical one such as EPSG:3855 (EGM2008 height), confirm the geoid grid is installed, and check that the returned height actually differs from the input — an unchanged value means the vertical component was ignored.

Is EPSG:4326 accurate enough as a target for survey-grade work? Not on its own. EPSG:4326 names the WGS 84 datum ensemble, whose member realizations differ by up to about two metres, so a transform into it can be internally consistent and still two metres from a specific realization. Where that matters, target a named realization or a plate-fixed frame such as ETRS89, and record which one in the dataset metadata.