Reprojecting CAD Coordinates with pyproj Transformer

To reproject CAD coordinates in Python, build a single pyproj.Transformer.from_crs(src, dst, always_xy=True), then call .transform(x, y) on numpy arrays of the whole point set at once. The always_xy=True flag is mandatory: it guarantees (easting/longitude, northing/latitude) order regardless of how the CRS authority defines its axes, which is the difference between correct output and silently swapped coordinates. This page is an implementation reference within the CRS Normalization Workflows topic. It assumes the CAD data already sits in a known projected CRS and in metres — if it does not, resolve that first, because reprojection cannot fix an unknown datum or a millimetre-scale drawing.

How pyproj Transformer Handles Reprojection

pyproj.Transformer wraps a PROJ coordinate operation: given a source and destination CRS, PROJ selects a concrete pipeline (projection, datum shift, and any grid corrections) and evaluates it. Two properties of that design govern how you should use it.

First, axis order follows the CRS authority, not intuition. EPSG:4326 is officially latitude-first. Since pyproj 2.0 the library respects this, so a naive Transformer.from_crs("EPSG:25832", "EPSG:4326") returns (lat, lon). Every CAD, GeoJSON, and shapefile toolchain expects (x, y) / (lon, lat). Passing always_xy=True normalises both input and output to that traditional order, eliminating an entire class of “my points landed in the wrong hemisphere” bugs.

Second, the Transformer is a reusable, thread-safe object. Constructing it resolves the CRS pair and selects an operation — work you do not want to repeat per point. Build it once, cache it, and call .transform() as many times as you like, including from multiple threads. The .transform() call also accepts numpy arrays directly, dispatching the whole batch into PROJ in a single C call rather than paying Python-loop overhead per coordinate.

pyproj does not read CAD files or infer their CRS — it only transforms coordinates you hand it. Establishing that the source CRS is correct, and that the geometry is in metres, is a precondition handled upstream. What pyproj guarantees is a correct, datum-aware mapping between two known CRSs.

Production-Ready Script

The script reprojects arrays of CAD points from a projected CRS (EPSG:25832, ETRS89 / UTM 32N) to EPSG:4326, in memory-bounded chunks, using one cached Transformer. It also includes a helper that auto-selects the appropriate WGS84 UTM CRS for a site.

# pyproj>=3.5.0, numpy>=1.24.0, Python 3.9+
from __future__ import annotations

from functools import lru_cache

import numpy as np
from pyproj import Transformer
from pyproj.aoi import AreaOfInterest
from pyproj.database import query_utm_crs_info


@lru_cache(maxsize=64)
def get_transformer(src: str, dst: str) -> Transformer:
    """Return a cached, thread-safe Transformer for a CRS pair.

    Cached so repeated calls with the same CRS pair reuse one instance
    rather than re-resolving the operation on every batch.
    """
    return Transformer.from_crs(src, dst, always_xy=True)


def reproject_points(
    xy: np.ndarray,
    src_crs: str,
    dst_crs: str = "EPSG:4326",
    chunk_size: int = 500_000,
) -> np.ndarray:
    """Reproject an (N, 2) array of (x, y) coordinates.

    Args:
        xy:       (N, 2) array of coordinates in src_crs (x/easting, y/northing).
        src_crs:  Source CRS, e.g. "EPSG:25832". Must be a known projected CRS.
        dst_crs:  Destination CRS. Defaults to WGS84 geographic.
        chunk_size: Rows processed per PROJ call to bound peak memory.

    Returns:
        (N, 2) array in dst_crs order (lon, lat for EPSG:4326). Out-of-domain
        inputs yield non-finite rows; the caller should filter with isfinite.
    """
    xy = np.atleast_2d(np.asarray(xy, dtype=np.float64))
    if xy.shape[1] != 2:
        raise ValueError(f"Expected (N, 2) array, got shape {xy.shape}.")

    transformer = get_transformer(src_crs, dst_crs)
    out = np.empty_like(xy)
    for start in range(0, xy.shape[0], chunk_size):
        stop = start + chunk_size
        x_out, y_out = transformer.transform(
            xy[start:stop, 0], xy[start:stop, 1]
        )
        out[start:stop, 0] = x_out
        out[start:stop, 1] = y_out
    return out


def finite_mask(xy: np.ndarray) -> np.ndarray:
    """Boolean mask of rows with all-finite coordinates (drops PROJ inf/nan)."""
    return np.isfinite(xy).all(axis=1)


def auto_utm_crs(lon: float, lat: float) -> str:
    """Return the WGS84 UTM EPSG code for a site's approximate lon/lat.

    Uses pyproj's CRS database for a robust result near zone boundaries;
    falls back to the standard zone formula if the query returns nothing.
    """
    utm_list = query_utm_crs_info(
        datum_name="WGS 84",
        area_of_interest=AreaOfInterest(
            west_lon_degree=lon, south_lat_degree=lat,
            east_lon_degree=lon, north_lat_degree=lat,
        ),
    )
    if utm_list:
        return f"EPSG:{utm_list[0].code}"
    zone = int((lon + 180.0) // 6.0) + 1
    return f"EPSG:{(32600 if lat >= 0 else 32700) + zone}"


if __name__ == "__main__":
    # CAD points already in ETRS89 / UTM 32N (metres). Replace with real data.
    cad_xy = np.array([
        [500000.0, 5570000.0],
        [500250.5, 5570180.2],
        [499800.0, 5569920.7],
    ])

    lonlat = reproject_points(cad_xy, "EPSG:25832", "EPSG:4326")

    mask = finite_mask(lonlat)
    if not mask.all():
        print(f"[WARN] Dropped {(~mask).sum()} out-of-domain point(s).")
    lonlat = lonlat[mask]

    # Validate: reprojected extent must sit inside the destination domain.
    lon_min, lat_min = lonlat.min(axis=0)
    lon_max, lat_max = lonlat.max(axis=0)
    assert -180.0 <= lon_min and lon_max <= 180.0, "Longitude out of range"
    assert -90.0 <= lat_min and lat_max <= 90.0, "Latitude out of range"

    print("Site centroid CRS:", auto_utm_crs(float(lonlat[:, 0].mean()),
                                             float(lonlat[:, 1].mean())))
    for lon, lat in lonlat:
        print(f"  lon={lon:.7f}  lat={lat:.7f}")

Key implementation notes:

• get_transformer is wrapped in lru_cache, so a pipeline that reprojects thousands of batches with the same CRS pair pays the construction cost once. The cached instance is safe to call from multiple worker threads.
• reproject_points passes numpy arrays straight into transformer.transform. This is the vectorized path — PROJ processes the whole chunk in one call. Never loop transform(x, y) per coordinate; that is orders of magnitude slower.
• chunk_size bounds peak memory for very large point clouds (millions of vertices from survey drawings) without changing results. Tune it to your worker’s memory budget.
• auto_utm_crs prefers query_utm_crs_info, which consults the PROJ CRS database and is correct at zone edges, and falls back to the closed-form zone formula only if the query yields nothing.
• The bounding-box asserts are a cheap regression guard: if a datum or axis mistake creeps in, reprojected coordinates leave the valid longitude/latitude range and the run fails loudly instead of writing garbage.

Compatibility Matrix

Fallback Strategies

1. Swapped coordinates from axis-order defaults

If reprojected points appear mirrored across the equator or land at impossible latitudes, always_xy=True was omitted and pyproj returned (lat, lon). Set the flag on every Transformer.from_crs call. There is no downside for (x, y) toolchains, and it makes input and output order explicit and identical.

2. Missing datum-shift grids degrade accuracy silently

When the source-to-destination operation needs an NTv2 or geoid grid that is not installed, PROJ falls back to a lower-accuracy operation without raising an error. Inspect the available operations and their accuracy with TransformerGroup, and pick deliberately:

# pyproj>=3.5.0
from pyproj.transformer import TransformerGroup

tg = TransformerGroup("EPSG:25832", "EPSG:4326", always_xy=True)
for t in tg.transformers:
    print(t.description, "| accuracy(m):", t.accuracy)
if tg.unavailable_operations:
    print("Missing grids:", len(tg.unavailable_operations))
    # tg.download_grids(verbose=True)  # fetch when network access is allowed

Confirm the active grid search path with import pyproj; print(pyproj.datadir.get_data_dir()), and pre-stage grids in air-gapped deployments.

3. Out-of-area points return inf or nan

A point far outside the source or destination projection domain (for example an outlier CAD vertex, or data from the wrong UTM zone) reprojects to inf. Writing that into PostGIS corrupts the layer’s bounding box and spatial index. Always filter with the finite_mask helper and log the dropped inputs for review rather than passing them downstream.

4. Per-point Transformer creation — the performance anti-pattern

Calling Transformer.from_crs(...) inside the loop that processes each coordinate re-resolves the CRS pair and re-selects an operation on every iteration, which can dominate runtime by two or three orders of magnitude. Build the Transformer once (cache it, as above) and pass arrays. This single change routinely turns a multi-minute reprojection of a large survey drawing into a sub-second one.

5. Reprojecting before units are metres

pyproj assumes projected coordinates are in the CRS’s linear unit — metres for UTM. A DXF drawing authored in millimetres passed straight into EPSG:25832 is off by a factor of 1000 and lands far outside the zone, returning inf. Normalise units first; the Unit Conversion Pipelines topic covers reading $INSUNITS and scaling geometry to metres before any reprojection.

For the datum-parameter side of alignment — when the source is not a registered CRS but a set of published transform parameters — see Applying a Helmert 7-Parameter Transform in Python. Once points are in WGS84, serialise them with Converting CAD Polylines to GeoJSON.

FAQ

Related Pages

• CRS Normalization Workflows — parent topic: full pipeline from CRS detection through reprojection and validation
• Applying a Helmert 7-Parameter Transform in Python — datum-shift parameters for cases where the source is not a registered CRS
• Converting CAD Local Coordinates to EPSG:4326 — the upstream step that registers an arbitrary CAD grid to a projected CRS from survey control
• Unit Conversion Pipelines — scaling CAD geometry to metres before any reprojection
• Converting CAD Polylines to GeoJSON — serialising reprojected WGS84 geometry for GIS consumption