Coordinate Systems with PyProj: Precision Transformations & Pipeline Integration

Spatial accuracy begins with rigorous coordinate reference system (CRS) management, and PyProj is the component of Mastering Core Geospatial Python Libraries that owns it — the thin, fast wrapper over the PROJ engine that every other tool defers to for datum-safe conversions. This guide covers how PyProj models a CRS, how the Transformer executes conversions, and how to wire both into reproducible pipelines that feed GeoPandas DataFrames and Shapely geometry operations without silent datum shifts. The patterns below target pyproj>=3.4 on Python>=3.10, eliminate axis-flip bugs, and enforce strict spatial validation across vector and raster workflows.

PyProj transformation pipeline A source CRS in EPSG 4326 passes through a reusable Transformer built with always_xy true, producing coordinates in a projected target CRS such as UTM zone 33N, with a note that always_xy enforces longitude-latitude order. EPSG:4326 lon/lat degrees Transformer always_xy=True build once, reuse EPSG:32633 UTM 33N, metres always_xy=True enforces (lon, lat) order — the fix for the most common axis-flip bug
Build the Transformer once with always_xy=True and reuse it — re-parsing CRS objects inside loops is the classic performance trap.

Architecture & Data Structures

PyProj exposes two objects that carry almost every workflow: CRS, an immutable description of a coordinate reference system, and Transformer, a compiled operation that maps coordinates from one CRS to another. PyProj 3+ sits on PROJ 8+ and speaks WKT2:2019 and EPSG authority codes natively, so the legacy +proj= PROJ-4 strings and the old Proj/transform free functions are deprecated in favour of these two classes. A CRS never transforms anything on its own — it is metadata (datum, ellipsoid, axis order, units, area of use). The Transformer is where the arithmetic lives, and it is the object you cache.

from pyproj import CRS, Transformer

# CRS = immutable metadata; Transformer = the compiled operation
wgs84 = CRS.from_epsg(4326)          # geographic, lon/lat in degrees
utm33n = CRS.from_epsg(32633)        # projected, easting/northing in metres

print(wgs84.is_geographic, utm33n.is_projected)   # True True
print(utm33n.axis_info[0].abbrev)                  # 'E' (easting first)

to_utm = Transformer.from_crs(wgs84, utm33n, always_xy=True)
easting, northing = to_utm.transform(12.4924, 41.8902)   # (lon, lat) in, metres out

Underneath the flat CRS façade is a small object graph, and knowing its shape turns most "why did this transform behave like that" questions into a two-line inspection. A CRS owns a datum (which owns an ellipsoid and a prime meridian), a coordinate system (which owns the ordered axes and their units), and an area of use. A projected CRS additionally exposes the conversion that produced it, and a compound CRS exposes its components through sub_crs_list. Every one of those is reachable without leaving Python.

from pyproj import CRS

bng = CRS.from_epsg(27700)          # OSGB36 / British National Grid

print(bng.datum.name)                # OSGB36
print(bng.ellipsoid.name, bng.ellipsoid.semi_major_metre)   # Airy 1830 6377563.396
print(bng.prime_meridian.name)       # Greenwich
print([ax.unit_name for ax in bng.axis_info])               # ['metre', 'metre']
print(bng.area_of_use.bounds)        # (-9.01, 49.75, 2.01, 61.01) — lon/lat degrees
print(bng.coordinate_operation.method_name)                 # Transverse Mercator

# Compound (2D + vertical) CRSs decompose into their parts
osgb_with_height = CRS.from_epsg(7405)                      # OSGB36 / BNG + ODN height
print(osgb_with_height.is_compound, len(osgb_with_height.sub_crs_list))   # True 2

Two equality rules follow from that graph and both catch people out. crs_a == crs_b compares the definitions, so a CRS built from an EPSG code and the same CRS parsed from a shapefile .prj can compare unequal because one carries an authority and the other does not. crs_a.equals(crs_b, ignore_axis_order=True) is the comparison you usually want, and is_exact_same is the strict one for reproducibility tests. GeoPandas layers that "should" already match but produce empty joins almost always fail on this distinction rather than on the geometry — the serialization side of it is covered in EPSG code vs PROJ string vs WKT. CRS.from_user_input() is the permissive constructor that accepts any of those forms — an int, an EPSG: string, WKT of any dialect, PROJJSON, or another CRS — and is the right entry point when you are ingesting CRS metadata you did not write.

The critical property to internalize is that a CRS knows its own authority-defined axis order, and for most geographic systems (including EPSG:4326) that order is latitude-then-longitude. Passing coordinates in the intuitive (lon, lat) order without accounting for this is the single most common PyProj bug. always_xy=True normalizes the Transformer so it always accepts and returns coordinates in (x, y) = (lon, lat) / (easting, northing) order, regardless of what the underlying authorities declare. When a transform still fails or returns inf, the diagnosis flow lives in fixing PyProj CRS transformation errors.

Authority axis order versus always_xy normalization The identical coordinate pair 12.49, 41.89 is passed to a Transformer twice. Under the default authority axis order, EPSG:4326 reads it as latitude then longitude and the point lands in the Gulf of Aden — wrong. With always_xy set to True the Transformer reads it as longitude then latitude and the point lands in Rome — correct. Default — authority axis order CRS.from_epsg(4326) declares (lat, lon) transform(12.49, 41.89) 12.49 → latitude 41.89 → longitude Lands in the Gulf of Aden ✗ wrong location vs always_xy=True — (x, y) order Transformer normalizes to (lon, lat) transform(12.49, 41.89) 12.49 → longitude 41.89 → latitude Lands in Rome, Italy ✓ correct location
The bytes are identical — only the declared axis order differs. always_xy=True is the one-line guarantee that (lon, lat) means what you think it means.

Environment Configuration & Dependency Resolution

PyProj bundles its own copy of the PROJ C library and the core proj.db authority database inside the wheel, so a pip install pyproj on a modern platform is self-contained — you do not need a system PROJ. Version alignment still matters: PyProj carries a hard minimum PROJ version, and a stale conda proj package on the path can shadow the wheel's bundled one, producing subtle discrepancies. Pin explicitly and verify at build time.

# Reproducible install — pin both the wrapper and let it carry its PROJ
python -m pip install "pyproj>=3.6,<4" numpy
import pyproj

# Run this in CI; it prints PyProj, PROJ, and the data-directory paths
pyproj.show_versions()
# Confirm the bundled data dir is the one in use (not a stray system PROJ_LIB)
print(pyproj.datadir.get_data_dir())

High-accuracy transforms — NADCON/NTv2 datum grids, geoid models, ITRF realizations — are not shipped in the wheel because they are large. PyProj can fetch them on demand from the PROJ CDN when network access is enabled:

import pyproj

# Allow on-demand grid downloads (sub-metre datum-shift accuracy)
pyproj.network.set_network_enabled(True)   # or set PROJ_NETWORK=ON in the environment

In containers, prefer the environment variable PROJ_NETWORK=ON so the setting is declarative, and consider baking the grids you actually need into the image to avoid a runtime dependency on the CDN. Never set PROJ_LIB by hand on a modern PyProj install — it exists only for legacy system-PROJ setups and is the usual cause of "proj.db not found" errors covered under platform debugging below.

One version detail is worth committing to memory because it silently changes which override wins: PROJ 9.1 renamed the data-directory variable from PROJ_LIB to PROJ_DATA. Both are honoured for now, PROJ_DATA taking precedence, but a Dockerfile that exports the old name against a new PROJ, or a CI runner that inherits PROJ_LIB from a base image while your code sets PROJ_DATA, produces two different answers on two machines from identical source. Assert the resolved path rather than the variable.

The pairing between the wrapper and the engine is equally worth pinning. PyProj carries a hard minimum PROJ version and bundles a matching build in its wheel, so the practical rule is that you choose a pyproj version and inherit a PROJ version rather than selecting them independently — pyproj 3.4 ships PROJ 9.x, 3.6 moved to PROJ 9.3, and 3.7 raised the floor to Python 3.10. Because the operation catalogue in proj.db is versioned too, an unpinned upgrade can re-rank the operations chosen for a datum pair and change your output coordinates without a line of your code changing.

import pyproj

print(pyproj.__version__)          # the Python wrapper
print(pyproj.__proj_version__)     # the PROJ engine actually loaded
print(pyproj.datadir.get_data_dir())        # where proj.db was found
print(pyproj.datadir.get_user_data_dir())   # where downloaded grids land

# Fail the build if the engine is not the one the results were validated against
assert pyproj.__proj_version__.startswith("9."), pyproj.__proj_version__

Two installation routes work and they should not be mixed. A pure pip environment is self-contained because every wheel bundles its own PROJ and GEOS. A conda-forge environment shares one system PROJ across pyproj, gdal, and rasterio, which is what you want when several libraries must agree on the same authority database. Layering a pip pyproj on top of a conda gdal gives you two PROJ builds and a coin-flip about which proj.db gets read — the mechanics of that failure on Windows are laid out in how to install and configure GeoPandas on Windows.

Vectorized Operations & Core Workflow

The Transformer.transform() method is array-aware: hand it NumPy arrays and it dispatches to a single vectorized PROJ call rather than looping in Python. This is the difference between transforming ten million points in a fraction of a second versus minutes. The rule is simple — build the Transformer once, outside any loop, and feed it whole arrays.

import numpy as np
from pyproj import Transformer

# Build ONCE — reuse across every batch. Re-creating this per call is the trap.
to_utm = Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True)

# A batch of sensor stations near Rome, in lon/lat degrees
station_lon = np.array([12.45, 12.46, 12.47, 12.48])
station_lat = np.array([41.90, 41.91, 41.92, 41.93])

# One vectorized C call — not a Python loop
easting, northing = to_utm.transform(station_lon, station_lat)

# Round-trip check: transform back and compare within tolerance
back = Transformer.from_crs("EPSG:32633", "EPSG:4326", always_xy=True)
lon_rt, lat_rt = back.transform(easting, northing)
assert np.allclose(station_lon, lon_rt, atol=1e-7)

For arrays that exceed comfortable memory (tens of millions of points), split the work with numpy.array_split and transform chunk-by-chunk with the same cached Transformer, keeping peak memory bounded. Use radians=False (the default) for degree input, and pass a fourth zz array when transforming 3D coordinates through a compound or geoid-aware CRS. Avoid calling transform() element-by-element inside a Python for loop — that pays the dispatch cost per point and forfeits the entire performance advantage.

# Chunked transform for very large point clouds — bounded peak memory
def transform_chunked(transformer, xs, ys, chunk=5_000_000):
    out_x, out_y = [], []
    for xi, yi in zip(np.array_split(xs, max(1, len(xs) // chunk)),
                      np.array_split(ys, max(1, len(ys) // chunk))):
        tx, ty = transformer.transform(xi, yi)
        out_x.append(tx)
        out_y.append(ty)
    return np.concatenate(out_x), np.concatenate(out_y)

Three further switches on transform() change its memory and throughput profile, and each solves a different problem. inplace=True writes the results back into the input arrays instead of allocating new ones, which halves peak memory on a point cloud that is already at the edge of RAM — at the cost of destroying the source coordinates, so only use it when the input is a scratch copy. direction=TransformDirection.INVERSE runs the same compiled pipeline backwards, meaning one cached transformer covers both directions and a round-trip check costs one PROJ lookup rather than two. And itransform() consumes an iterable of coordinate tuples and yields transformed ones, which is the right shape when coordinates arrive from a database cursor or a streaming reader and never exist as a single array.

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

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

# Scratch arrays — inplace=True overwrites them and allocates nothing
lidar_x = np.array([12.45, 12.46, 12.47], dtype="float64")
lidar_y = np.array([41.90, 41.91, 41.92], dtype="float64")
to_utm.transform(lidar_x, lidar_y, inplace=True)      # now metres, in the same buffers

# ONE transformer, both directions — no second from_crs() call
lon, lat = to_utm.transform(lidar_x, lidar_y, direction=TransformDirection.INVERSE)
assert np.allclose(lon, [12.45, 12.46, 12.47], atol=1e-7)

# Streaming shape: tuples in, tuples out, nothing materialised
station_stream = ((12.45, 41.90), (12.46, 41.91), (12.47, 41.92))
for easting, northing in to_utm.itransform(station_stream):
    pass

Know where the approach stops scaling. Vectorized transform() is a single C call over a contiguous float64 array, so throughput is bounded by memory bandwidth and typically lands in the millions of points per second for a simple projection; a grid-based datum shift is slower because it interpolates a raster per point, and a time-dependent transformation slower again. What actually breaks first is not speed but allocation: transform() returns new arrays, so a naive call on a billion-point cloud needs roughly 32 GB of headroom for the inputs and outputs together. Chunk before you reach that point, and reach for inplace=True only after profiling confirms the copy is the constraint. Above roughly a hundred million points the honest answer is to stop transforming coordinates in Python at all and reproject at the storage layer instead — see reprojecting large datasets without memory errors.

Geometry / Data Processing Details

Coordinate transformation is only correct when it happens at the right stage of a geometry workflow. Every metric operation — buffering, area, length, distance, and most overlays — assumes a planar, equal-scale coordinate space. Running them in a geographic CRS silently interprets degrees as if they were a flat Cartesian unit, which is wrong by a factor that varies with latitude. The discipline is: project to a suitable metric CRS, run the Shapely geometry operation, then transform the result back only for display or export.

PyProj integrates with Shapely 2.x through shapely.ops.transform, which walks a geometry's coordinates and applies any callable — including a Transformer.transform — while preserving structure and ring order.

from shapely.geometry import Point
from shapely.ops import transform as shp_transform
from shapely import make_valid
from pyproj import Transformer

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

# A monitoring station; buffer it by 1 km — a METRIC operation
station = Point(12.4924, 41.8902)                 # lon, lat
station_utm = shp_transform(to_utm.transform, station)

catchment_utm = station_utm.buffer(1000)          # 1000 metres, correct in UTM
catchment_wgs = shp_transform(from_utm.transform, catchment_utm)

# Projection warping can introduce self-intersections — repair before use
catchment_wgs = make_valid(catchment_wgs)

Because always_xy=True is set on both transformers, the coordinate order stays consistent through the round trip and no axis flip creeps in. When a projection stretches a polygon across a zone boundary it can produce slivers or invalid rings; validating with make_valid (Shapely 2.x) after any transform is cheap insurance. Never compute .area or .distance on a geometry still in EPSG:4326 — the numbers come back in square-degrees and degrees, which are meaningless as physical measurements.

Measuring on the ellipsoid instead of projecting

Project-then-measure is the right default, but it has a hard limit: a projected CRS is only accurate inside its area of use, so a feature that spans a continent, crosses several UTM zones, or reaches into polar latitudes cannot be measured honestly in any single projection. PyProj's third major class, Geod, solves that case by computing directly on the ellipsoid using Karney's geodesic algorithms — no projection, no zone choice, and accuracy at the sub-millimetre level anywhere on Earth including antipodal pairs, where the older Vincenty formulation fails to converge.

from pyproj import Geod
from shapely.geometry import LineString, Polygon

geod = Geod(ellps="WGS84")

# 1. True ground distance between two airports, no projection involved
_, _, distance_m = geod.inv(-0.4543, 51.4700, 103.9915, 1.3644)   # LHR -> SIN
print(f"{distance_m / 1000:,.1f} km")
# 10,878.7 km

# 2. Move a known bearing and distance from a start point (the forward problem)
end_lon, end_lat, back_azimuth = geod.fwd(12.4924, 41.8902, 45.0, 25_000)

# 3. Ellipsoidal area and perimeter of a polygon still in lon/lat degrees
survey_block = Polygon([(12.45, 41.89), (12.50, 41.89), (12.50, 41.93), (12.45, 41.93)])
area_m2, perimeter_m = geod.geometry_area_perimeter(survey_block)
print(f"{abs(area_m2):,.0f} m2  ·  {perimeter_m:,.0f} m")

# 4. Length of a flight track that crosses many UTM zones
track = LineString([(-0.4543, 51.4700), (55.3644, 25.2528), (103.9915, 1.3644)])
print(f"{geod.geometry_length(track) / 1000:,.1f} km")

geometry_area_perimeter returns a signed area — negative for a clockwise ring — which is why the example wraps it in abs(); the sign is useful when you want to detect winding order but it will quietly turn a total into a subtraction if you sum a mixture of rings without normalizing. Note also that these methods take degrees, so the input geometry must still be geographic; feeding a projected geometry to Geod produces a large, confident, meaningless number.

The trade-off against project-then-measure is throughput. Each geodesic call solves an iterative problem per coordinate pair, so Geod runs one to two orders of magnitude slower than a planar .length on a projected GeoSeries. That makes it the wrong tool for measuring a million parcels and the right tool for the cases a projection cannot serve: long-haul distances, features crossing the antimeridian, validating a projected result against ground truth, and any single number that ends up in a report. A practical division of labour is to measure in bulk in the local metric CRS and spot-check a random sample with Geod.inv, asserting the two agree within the tolerance the projection promises.

CRS Alignment & Projection Pipeline

At the DataFrame layer, GeoPandas wraps PyProj so alignment is declarative: .to_crs() builds the right Transformer internally and applies it column-wide. The prerequisite is that the source frame actually carries a CRS — an untagged frame (gdf.crs is None) means GeoPandas has no idea what the numbers mean and cannot transform them. Assign the known source CRS with .set_crs() first; only then .to_crs() to reproject.

import geopandas as gpd

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

# Guard: a missing CRS is a data error, not something to guess around
if parcels.crs is None:
    parcels = parcels.set_crs("EPSG:4326")     # declare the KNOWN source, do not reproject

# Reproject the whole frame to a metric CRS for area/length work
parcels_utm = parcels.to_crs("EPSG:32633")
parcels_utm["area_m2"] = parcels_utm.geometry.area

For a national or continental dataset that spans several UTM zones, a single fixed zone distorts the edges. The robust pattern is to pick the zone per feature or per bounding box — the logic is laid out in choosing the right UTM zone automatically in Python, which uses the dataset extent to derive the EPSG code before projecting. When you only need the extent of a reprojection (for tiling or window math) rather than every vertex, Transformer.transform_bounds reprojects a bounding box while densifying the edges so the result is not clipped by nonlinear warping:

from pyproj import Transformer

to_utm = Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True)
# Densified bbox reprojection — safe for curved projection edges
minx, miny, maxx, maxy = to_utm.transform_bounds(12.3, 41.8, 12.6, 42.0)
Choosing a metric CRS by analysis scope Geographic input in EPSG:4326 branches by scope. A local or regional analysis projects to a UTM zone for metre-accurate area and length. A continental or national analysis uses an equal-area CRS such as Lambert Azimuthal Equal Area or Albers to preserve area. Web Mercator EPSG:3857 is reserved for web-tile display and must never be used for metric analysis. Geographic input EPSG:4326 — lon/lat degrees Match the CRS to analysis scope project before any metric operation Local / regional UTM zone — EPSG:326xx metre grid, <1 m in-zone ✓ area, length, buffer Continental / national Equal-area — LAEA / Albers preserves area over wide extent ✓ zonal statistics Web tiles / display only Web Mercator — EPSG:3857 scale grows with latitude ✗ never for metric work
Scope decides the projection: a UTM zone for local metric work, an equal-area CRS for wide-extent statistics, and EPSG:3857 strictly for the tile-render boundary — never for computing area or distance.

Reserve Web Mercator (EPSG:3857) strictly for the final web-tile render stage; its scale factor grows away from the equator, so any area or distance computed in it is wrong. For continental statistics choose an equal-area CRS (for example a Lambert Azimuthal Equal Area or Albers definition for your region) rather than a single UTM zone.

Vertical coordinates travel on a separate track and are the part of the pipeline most often dropped without anyone noticing. A 2D Transformer accepts a third zz argument and passes it through unchanged, so elevations survive the call while being silently wrong the moment the two CRSs use different height systems — ellipsoidal height above the reference ellipsoid is not the same quantity as orthometric height above a geoid, and the difference reaches tens of metres in many regions. Getting a real vertical transformation requires naming a CRS that includes the vertical component, at which point PROJ inserts a geoid model into the pipeline.

from pyproj import Transformer

# 2D transformer: z is carried through untouched, NOT converted
flat = Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True)
print(flat.transform(12.4924, 41.8902, 120.0)[2])       # 120.0 — unchanged

# Compound source: 3D geographic + EGM2008 orthometric height -> 3D geographic
# (requires the geoid model on disk or PROJ_NETWORK=ON)
vertical = Transformer.from_crs("EPSG:9518", "EPSG:4979", always_xy=True)
lon, lat, ellipsoidal_h = vertical.transform(12.4924, 41.8902, 120.0)

If the third value comes back byte-identical to what you passed in, no vertical transformation happened — that is the one-line diagnostic. Where a geoid model is genuinely needed and missing, the grid-resolution mechanics are the same as for horizontal datum shifts, covered in datum shifts and transformation grids in PyProj.

Production Export & Integration

In a long-running service, transformer construction — not the transform itself — dominates cost if done carelessly. Cache transformers by their CRS pair and warm them at startup. Since pyproj 3.7 the library manages one PROJ context per thread automatically, so the older pyproj.set_use_global_context() call is deprecated and unnecessary; a module-level lru_cache factory is enough.

from functools import lru_cache
import pyproj

@lru_cache(maxsize=128)
def get_transformer(src_epsg: int, dst_epsg: int) -> pyproj.Transformer:
    """Cached factory — avoids redundant PROJ pipeline compilation."""
    return pyproj.Transformer.from_crs(
        f"EPSG:{src_epsg}", f"EPSG:{dst_epsg}", always_xy=True
    )

def project_points(src_epsg, dst_epsg, coords):
    """CI-ready wrapper: coords is a list of (lon, lat) tuples."""
    transformer = get_transformer(src_epsg, dst_epsg)
    xs, ys = zip(*coords)
    return transformer.transform(list(xs), list(ys))

# Reproject a station list from WGS84 to UTM 33N for metric analysis
easting, northing = project_points(4326, 32633, [(12.4924, 41.8902)])
assert len(easting) == 1

For integration, the EPSG code is the contract across systems. When writing to PostGIS, the geometry column's SRID must match the CRS you reprojected to — gdf.crs.to_epsg() gives the integer to hand to the srid argument of GeoDataFrame.to_postgis(). Cloud-native formats behave the same way: GeoParquet stores the CRS in its metadata, and a .to_parquet() export from a correctly-tagged frame round-trips the CRS without a separate .prj sidecar. Rasterio-based pixel pipelines described in raster data handling with Rasterio consume the same CRS objects, so a single canonical EPSG code can govern both the vector and raster halves of a pipeline.

Deployment checklist:

Windows / Platform Edge Cases & Debugging

The dominant class of platform failure is a mismatched or shadowed PROJ data directory. A leftover PROJ_LIB environment variable — often set by an old GDAL, QGIS, or conda install — points PyProj at the wrong proj.db and surfaces as pyproj.exceptions.DataDirError or a "proj.db not found / version mismatch" message. On a modern pip-installed PyProj the fix is to unset it and let PyProj use its bundled directory.

import os, pyproj

# Diagnose: which proj.db is actually loaded?
print(pyproj.datadir.get_data_dir())
print("PROJ_LIB =", os.environ.get("PROJ_LIB"))   # ideally None on a wheel install

# If a stray PROJ_LIB is shadowing the wheel, point PyProj back at its own data
pyproj.datadir.set_data_dir(pyproj.datadir.get_data_dir())

Common platform gotchas and one-line fixes:

Frequently Asked Questions

When should I use Geod instead of projecting to a metric CRS? Use a projection whenever the data fits inside one projected CRS's area of use and you are measuring many features — it is one to two orders of magnitude faster and integrates with every vectorized GeoPandas and Shapely operation. Reach for Geod when no single projection is honest: transcontinental distances, features crossing the antimeridian or a polar region, and any measurement spanning more than a few UTM zones. The third case is validation — measure in bulk on the projection and spot-check a sample with Geod.inv, which gives you a ground-truth reference that does not depend on the projection being right.

Do I need PyProj at all if I only use GeoPandas? You are already using it — .to_crs(), .estimate_utm_crs(), and gdf.crs are all PyProj objects behind a DataFrame façade. Working with PyProj directly becomes necessary in three situations: when you need to inspect or assert which transformation operation PROJ selected, which .to_crs() gives you no access to; when coordinates arrive as raw arrays rather than geometries, where the array API is dramatically faster than building throwaway geometry; and when you need geodesic measurement, which has no GeoPandas equivalent. For everything else, staying at the GeoPandas level is the right altitude.

How much does building a Transformer actually cost? Enough to matter and not enough to panic about. Each from_crs() call is a proj.db query plus pipeline compilation — on the order of a millisecond, which is invisible once and ruinous a million times. The rule that follows is simple: the cost is per construction, not per coordinate, so caching on the CRS pair collapses it entirely. The pathological version is building the transformer inside a per-row apply(), which turns a sub-second array operation into minutes and is the single most common PyProj performance bug. A module-level lru_cache factory fixes it in three lines.

Is EPSG:4326 the same thing as WGS 84? Not exactly, and the gap is about two metres. EPSG:4326 names the WGS 84 ensemble, a family of realizations (G730 through G2296) that have drifted apart as the reference frame has been refined, and PROJ treats the ensemble as accurate to roughly 2 m internally. That is irrelevant for a web map and decisive for survey or GNSS work, where you should name the specific realization by its own EPSG code instead. It is also why a transform between EPSG:4326 and a modern national frame may report a stated accuracy of 2 m even with every grid present — the uncertainty is in the source definition, not in the operation.

Why does the same EPSG code give different results in two libraries? Because the code identifies the CRS, not the operation used to get there, and operation selection depends on the PROJ version, the available grids, and whether an area of interest was supplied. Two tools that both "use EPSG:25833" can pick different candidate operations and land a metre apart, entirely legitimately. Pin pyproj in your lockfile, provision the same grids everywhere, and record transformer.description alongside the CRS codes for anything reproducible — the description is what actually determines the numbers.

Can I define a CRS that is not in the EPSG registry? Yes, and there are two reasonable ways. CRS.from_wkt() with a hand-written WKT2 definition is the durable option for a local grid or an engineering coordinate system that must be stored and re-read. Transformer.from_pipeline() with an explicit +proj=pipeline definition is the option for a fixed, reproducible operation you want to lock down against any future re-ranking by PROJ. Both are legitimate; what is not legitimate is a bare +proj= string standing in for a CRS in stored metadata, because it carries no datum identity — the full argument is in EPSG code vs PROJ string vs WKT.

Related