pyproj vs GDAL for batch reprojection

Both libraries call PROJ, so a correctly configured pyproj transform and a correctly configured ogr2ogr transform produce the same coordinates to the last representable digit. The comparison is therefore not about accuracy. It is about throughput on a county-sized layer, whether the reprojection happens inside your process or as a subprocess, how each handles a missing grid-shift file, and which one lets you assert on the result. This guide measures those on a 214,000-parcel layer moving from a state plane CRS to EPSG:4326, extending CRS alignment strategies.

Diagnosis: which question are you actually asking?

Two distinct problems get labelled “batch reprojection” and they have different answers.

Reprojecting a file. A shapefile arrives in EPSG:2926 and you want a GeoPackage in EPSG:4326. Nothing else happens to the data. This is a file-format operation and GDAL’s command-line tools do it in one line, faster than anything you can write.

Wall time on a 214k-parcel layer, by where the work happens Grouped bar chart of wall-clock seconds to reproject a 214,000-parcel layer, comparing four approaches across two scenarios. For a whole-file conversion, ogr2ogr takes 31 seconds and pyproj with its own file I/O takes 74. Inside a pipeline where geometries are already in memory, vectorised pyproj takes 9 seconds and shelling out to ogr2ogr per batch takes 186. Each tool is roughly four times faster in its own scenario. Wall time on a 214k-parcel layer, by where the work happens 0 50 100 150 200 250 31 74 whole-file conversion 186 9 inside a pipeline Wall-clock seconds ogr2ogr subprocess pyproj, vectorised
Neither library is faster in general — each wins its own scenario by about 4×. The 186-second bar is the common mistake: a per-batch subprocess pays process startup and two format conversions per batch, and it looks like slow reprojection rather than slow plumbing.

Reprojecting inside a transformation. Parcels are being parsed, validated, normalised and upserted, and reprojection is one step among several. The geometries are already Python objects; writing them to a file so GDAL can read them back is pure overhead. This is pyproj’s case.

Choosing wrong costs about 4× either way, and the failure is quiet — a pipeline that shells out to ogr2ogr per batch spends most of its time on process startup and I/O.

Head-to-head

pyproj GDAL (ogr2ogr / Python bindings)
PROJ version whatever the wheel bundles whatever the system GDAL links
Reprojects coordinate arrays, in process files and layers
Best throughput shape one Transformer, vectorised arrays one process, whole layer
Grid shift missing raises or falls back per only_best warns, may silently use a lower-accuracy path
Format conversion none — bring your own I/O 80+ drivers
Geometry repair while reprojecting no -makevalid, -wrapdateline
Install weight a wheel with bundled PROJ a system library, or a large wheel
Assertable in tests directly via subprocess output or bindings

Step-by-step implementation

1. Build the transformer once, and transform arrays

The single largest pyproj performance mistake is constructing a Transformer per geometry. Construction parses the CRS definitions and selects an operation; reuse is what makes it fast.

import numpy as np
from pyproj import CRS, Transformer

SOURCE = CRS.from_epsg(2926)      # NAD83(HARN) / Washington North (ftUS)
TARGET = CRS.from_epsg(4326)

# always_xy keeps the argument order as (x, y) / (lon, lat) regardless of what the
# authority declares. Omitting it is the classic silent axis-order bug.
TF = Transformer.from_crs(SOURCE, TARGET, always_xy=True)


def reproject_rings(rings: list[np.ndarray]) -> list[np.ndarray]:
    """Transform every ring in one vectorised call.

    Concatenating first matters: one call over 400k points is roughly 30x faster
    than 8k calls over 50 points, because the per-call overhead dominates at that
    size."""
    if not rings:
        return []
    lengths = [len(r) for r in rings]
    flat = np.concatenate(rings)
    lon, lat = TF.transform(flat[:, 0], flat[:, 1])
    out = np.column_stack([lon, lat])
    bounds = np.cumsum([0] + lengths)
    return [out[bounds[i]:bounds[i + 1]] for i in range(len(rings))]
Where the time goes when a transform is called per geometry Lollipop chart of wall-clock seconds to reproject the same 400,000 coordinates four ways. Constructing a Transformer per geometry takes about 41 seconds. Reusing one Transformer but calling it per geometry takes about 2.8 seconds. Calling it per ring takes about 0.6 seconds. One concatenated vectorised call takes about 0.09 seconds. A dashed rule marks one second. Where the time goes when a transform is called per geometry 0 10 20 30 40 50 one second — above this the overhead dominates one vectorised call 0.09 s one call per ring 0.6 s one call per geometry 2.8 s new Transformer per geometry 41 s Seconds for 400k coordinates
Two independent mistakes, each worth roughly an order of magnitude: constructing the transformer inside the loop, and calling it per geometry instead of once over a concatenated array. Fixing both is about 450× on the same coordinates.

2. Make the grid-shift requirement explicit

A datum shift that needs a grid file will silently fall back to a lower-accuracy path if the grid is absent — a metre-scale error that never raises. State the requirement instead of hoping.

from pyproj.transformer import TransformerGroup


def assert_best_transform_available(source, target):
    """Fail at startup if the highest-accuracy operation is unavailable.

    An off-by-a-metre reprojection is worse than a crash: it passes every geometry
    check and puts parcel boundaries a metre from where the county has them."""
    group = TransformerGroup(source, target, always_xy=True)
    if not group.best_available:
        missing = [g.name for g in group.unavailable_operations]
        raise RuntimeError(
            f"best transform {source.to_epsg()}->{target.to_epsg()} unavailable; "
            f"missing grids for: {missing}. Run `projsync --source-id us_noaa` "
            f"or set PROJ_NETWORK=ON."
        )
    return group.transformers[0]
What a missing grid file does, and how each tool reports it A four-row by three-column matrix of grid-shift situations. Rows are the best grid present, the grid absent with network lookup off, the grid absent with PROJ network on, and the wrong source CRS declared. Columns are the resulting positional accuracy, what pyproj reports and what GDAL reports. The absent-grid case degrades to metre-scale accuracy while both tools report only a warning. What a missing grid file does, and how each tool reports it positional accuracy pyproj reports GDAL reports Best grid present centimetre best_available true nothing Grid absent, network off metre-scale a warning a warning Grid absent, network on centimetre, after a fetch best_available true nothing Wrong source CRS declared hundreds of metres nothing nothing accurate warns only silently wrong
The second row is the dangerous one: a metre-scale positional error that passes every geometry check, every bounds check and every test, reported as a warning nobody reads. That is why the startup assertion exists — a crash is preferable to a quiet metre.

3. Use GDAL when the unit of work is a file 

# One process, whole layer. -t_srs picks the operation; -makevalid repairs rings that
# only become invalid after the transform, which happens on long thin parcels.
ogr2ogr \
  -f GPKG parcels_4326.gpkg parcels_2926.shp \
  -s_srs EPSG:2926 -t_srs EPSG:4326 \
  -makevalid \
  -nlt PROMOTE_TO_MULTI \
  -lco SPATIAL_INDEX=YES \
  --config OGR_ENABLE_PARTIAL_REPROJECTION YES

-s_srs is explicit on purpose. Shapefile .prj files from municipal portals are frequently wrong or absent, and letting GDAL guess is how a layer ends up reprojected from the wrong source CRS — which produces coordinates in the right numeric range and the wrong place.

4. Decide once, at the pipeline boundary 

def reproject(source_path, target_path, *, in_pipeline: bool):
    """Files at the edges, arrays in the middle.

    The costly mistake is mixing them — a per-batch ogr2ogr subprocess inside an
    ingestion loop pays process startup and two format conversions per batch."""
    if in_pipeline:
        return [reproject_rings(g) for g in read_geometries(source_path)]
    return subprocess.run(
        ["ogr2ogr", "-f", "GPKG", target_path, source_path,
         "-s_srs", declared_crs(source_path), "-t_srs", "EPSG:4326", "-makevalid"],
        check=True, capture_output=True, text=True,
    )

Verification & testing

The test that matters is not “does it reproject” but “do both paths agree, and does a missing grid fail loudly.”

def test_both_paths_agree_within_a_millimetre(tmp_path, county_fixture):
    """The libraries share PROJ, so disagreement means a configuration difference —
    usually axis order or an implicit source CRS — not a numerical one."""
    via_pyproj = reproject_rings(read_geometries(county_fixture))
    out = tmp_path / "gdal.gpkg"
    subprocess.run(["ogr2ogr", "-f", "GPKG", str(out), str(county_fixture),
                    "-s_srs", "EPSG:2926", "-t_srs", "EPSG:4326"], check=True)
    via_gdal = read_geometries(out)

    for a, b in zip(via_pyproj, via_gdal):
        # ~1e-8 degrees is about a millimetre at this latitude.
        np.testing.assert_allclose(a, b, atol=1e-8)


def test_axis_order_is_pinned():
    """Without always_xy, EPSG:4326 yields (lat, lon) and every downstream bbox
    filter silently matches nothing."""
    tf = Transformer.from_crs(2926, 4326, always_xy=True)
    lon, lat = tf.transform(1_270_000, 230_000)
    assert -125 < lon < -116 and 45 < lat < 49


def test_missing_grid_raises_rather_than_degrading(monkeypatch):
    monkeypatch.setenv("PROJ_NETWORK", "OFF")
    monkeypatch.setenv("PROJ_DATA", str(EMPTY_GRID_DIR))
    with pytest.raises(RuntimeError, match="missing grids"):
        assert_best_transform_available(CRS.from_epsg(4267), CRS.from_epsg(4326))


def test_round_trip_is_stable_over_repeated_transforms(county_fixture):
    """A reprojection applied and reversed should return the original to within
    grid accuracy. Drift here means an operation without a defined inverse."""
    fwd = Transformer.from_crs(2926, 4326, always_xy=True)
    rev = Transformer.from_crs(4326, 2926, always_xy=True)
    for ring in read_geometries(county_fixture)[:200]:
        lon, lat = fwd.transform(ring[:, 0], ring[:, 1])
        x, y = rev.transform(lon, lat)
        np.testing.assert_allclose(np.column_stack([x, y]), ring, atol=1e-3)  # 1 mm in ftUS

Failure recovery

A layer reprojected from the wrong source CRS. The coordinates land in a plausible range, so nothing fails. Detect it by checking the reprojected envelope against the jurisdiction’s known bounds, and recover by reprojecting from the original file with an explicit -s_srs — never by reprojecting the already-wrong output, which compounds the error.

Rings that became invalid after the transform. Long thin parcels near a projection’s edge can self-intersect at the new precision. -makevalid in GDAL, shapely.make_valid after reproject_rings in Python; either way assert validity after the transform rather than before.

Two environments producing different coordinates. Almost always different PROJ data versions rather than different code. Pin the PROJ version and the grid set in the container, and log pyproj.proj_version_str and pyproj.datadir.get_data_dir() at startup so the difference is visible in a run’s own output.

Frequently asked questions

Which is faster on a county-sized layer?

For a whole-file conversion, GDAL — it does the reprojection and the format write in one pass, and there is nothing to beat about that. For geometries already in memory as arrays, pyproj with a reused Transformer and vectorised calls, because GDAL would require writing to a file and reading it back. The two numbers are not comparable; they answer different questions.

Do I need both installed?

Usually yes, and not for reprojection. GDAL is already in the stack for reading shapefiles, GeoPackages and GeoJSON, and pyproj arrives with pretty much any Python geospatial dependency. Having both is normal; the risk is that they link different PROJ versions, which is worth checking at startup rather than discovering from a coordinate discrepancy.

Is always_xy=True always right?

For this kind of work, yes — set it everywhere and stay in (x, y) / (lon, lat) order throughout. The authority-declared order for EPSG:4326 is latitude-first, and honouring it means every call site has to remember which convention applies. Pinning the order in one place and asserting it in a test removes an entire class of silent bug.

What about GeoPandas to_crs?

It is pyproj underneath, applied per geometry column, and it is the right choice when the data is already a GeoDataFrame — the vectorisation is handled for you. The caveats are the same: build the transform once (GeoPandas does), pin the axis order (it does), and check that the best operation is available (it does not — that check is still yours).