Vectorizing Coordinate Transforms with NumPy and Dask

Transforming a cadastral dataset one point at a time is the single most common performance defect in production coordinate pipelines, and eliminating it is the first optimisation covered under Batch Transformation & Automation for Cadastral Coordinate Pipelines: pyproj.Transformer.transform already accepts NumPy arrays and returns arrays, so a Python for loop that calls it per coordinate pays the C-extension crossing cost millions of times over for no accuracy benefit. This page narrows that parent topic to one concrete discipline — replacing per-point loops with a single vectorized array call, then scaling that same call beyond the memory of one machine by mapping it across partitioned Dask arrays and dataframes — while holding every value at float64 so results stay bit-for-bit reproducible against an agency-certified reference. The rule is simple: build the transformer once, hand it whole arrays, and let PROJ iterate in compiled code instead of the interpreter.

Figure — a coordinate array either transforms whole in one NumPy call or splits into Dask partitions that map, transform, and concatenate back to one result.

Why the Loop Is the Bottleneck #

Transformer.transform is a thin Python wrapper over PROJ’s compiled proj_trans_generic. When you pass two scalars it still crosses the Python/C boundary, allocates result objects, and returns them — a fixed overhead of a few microseconds that is invisible once but ruinous when repeated. For $N$ points the looped cost is

where $c_{\text{call}}$ is the per-call interpreter and boundary overhead and $c_{\text{trans}}$ is the actual projection arithmetic. The vectorized call collapses that to

paying the boundary cost exactly once. Because $c_{\text{call}}$ typically dominates $c_{\text{trans}}$ for a single point, the achievable speedup approaches $c_{\text{call}}/c_{\text{trans}}$, routinely one to two orders of magnitude. The full measurement discipline — timing both forms and reporting throughput in points per second — is the dedicated procedure in benchmarking NumPy vectorized versus looped transforms.

Vectorization alone stops at the memory ceiling: an array of $N$ float64 coordinate pairs needs $16N$ bytes for the input and as much again for the output, so a hundred-million-point cadastral extract wants several gigabytes just for the live arrays. When the working set exceeds RAM you keep the identical vectorized kernel but apply it to bounded partitions with Dask, which streams chunks through workers and never materialises the whole dataset at once. The out-of-core datum-shift case is worked end to end in Dask-partitioned datum shift for large cadastral datasets. Note that vectorizing is orthogonal to parallelism: a single vectorized call already saturates PROJ’s inner loop, and true multi-core distribution — including the transformer-reuse and thread-safety concerns — is the separate concern of concurrent pyproj transformation pipelines in Python.

Prerequisites and the always_xy Contract #

Two invariants make a vectorized transform reproducible. First, axis order: construct the transformer with always_xy=True so inputs and outputs are unconditionally (x=longitude/easting, y=latitude/northing), defeating the CRS-defined axis order that otherwise silently swaps latitude and longitude for EPSG:4326. Second, dtype: pass genuine float64 arrays. PROJ computes in double precision internally, so feeding float32 truncates your inputs to roughly seven significant digits — about a decimetre of longitude at mid-latitude — before the transform even runs. The mantissa budget is unforgiving: a coordinate near easting 500000 m carried in float32 cannot represent millimetres at all, because $2^{-23}\times 5\times10^{5}\approx 0.06\,\text{m}$. Everything below keeps arrays at float64 from ingestion to audit.

Step-by-Step Implementation #

Step 1 — Vectorize an (N,) array pair in one call #

The transformer is built once and reused; transform receives the full longitude and latitude arrays and returns easting and northing arrays of the same shape. No Python-level loop touches an individual coordinate.

from __future__ import annotations

import numpy as np
from pyproj import Transformer


def build_transformer(src_epsg: int, dst_epsg: int) -> Transformer:
    """Construct a reusable Transformer with deterministic axis order.

    always_xy=True fixes I/O to (x, y) = (lon/easting, lat/northing),
    overriding the CRS authority axis order (ISO 19111 coordinate system).
    """
    return Transformer.from_crs(src_epsg, dst_epsg, always_xy=True)


def transform_arrays(
    transformer: Transformer,
    lon: np.ndarray,
    lat: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
    """Transform paired coordinate arrays in a single vectorized call.

    Inputs are coerced to contiguous float64 so PROJ never sees a
    down-cast mantissa; the output dtype matches and is returned as
    (easting, northing).
    """
    x = np.ascontiguousarray(lon, dtype=np.float64)
    y = np.ascontiguousarray(lat, dtype=np.float64)
    if x.shape != y.shape:
        raise ValueError(f"lon/lat shape mismatch: {x.shape} vs {y.shape}")
    east, north = transformer.transform(x, y)  # PROJ iterates in C, not Python
    return np.asarray(east, dtype=np.float64), np.asarray(north, dtype=np.float64)

A subtle point: Transformer.transform returns -inf (not an exception) for points that fall outside the operation’s area of use, so the array call degrades gracefully rather than aborting the whole batch. You detect and quarantine those rows in verification rather than losing the run.

Step 2 — Scale the same kernel with Dask map_partitions #

When the array no longer fits in memory, wrap the coordinates in a Dask dataframe and apply the identical vectorized function to each partition. The transformer is rebuilt inside the mapped function so every worker owns its own PROJ context — a Transformer object holds a pointer into a native PROJ context that must not be pickled and shipped from the driver.

from __future__ import annotations

import dask.dataframe as dd
import pandas as pd


def transform_partition(
    part: pd.DataFrame,
    src_epsg: int,
    dst_epsg: int,
    lon_col: str = "lon",
    lat_col: str = "lat",
) -> pd.DataFrame:
    """Transform one partition; the Transformer is built here, on the worker.

    Building inside the mapped callable guarantees each worker gets a fresh
    native PROJ context instead of a serialised one from the driver.
    """
    transformer = build_transformer(src_epsg, dst_epsg)
    x = part[lon_col].to_numpy(dtype=np.float64)
    y = part[lat_col].to_numpy(dtype=np.float64)
    east, north = transformer.transform(x, y)
    out = part.copy()
    out["easting"] = np.asarray(east, dtype=np.float64)
    out["northing"] = np.asarray(north, dtype=np.float64)
    return out


def transform_dask(
    ddf: dd.DataFrame,
    src_epsg: int,
    dst_epsg: int,
) -> dd.DataFrame:
    """Map the vectorized transform across every partition lazily.

    `meta` declares the output schema so Dask never triggers an eager
    inference pass over real data, and float64 columns are pinned.
    """
    template = ddf._meta.copy()
    template["easting"] = np.float64()
    template["northing"] = np.float64()
    return ddf.map_partitions(
        transform_partition,
        src_epsg,
        dst_epsg,
        meta=template,
    )

The meta argument is the reproducibility hinge: without it Dask runs the function on a dummy frame to infer column dtypes and can guess float32 or object, corrupting precision before a single real row is processed. Declaring it makes the output schema explicit and the graph fully lazy until .compute() or .to_parquet() is called.

Parameter and Return-Value Reference #

Worked Example — Two Million Parcel Centroids #

Consider a statewide parcel layer whose centroids are stored as NAD83(2011) geographic coordinates, EPSG:6318, and must be delivered in NAD83(2011) / UTM zone 14N, EPSG:6343, the projected frame a Texas county records office ingests. We synthesise two million centroids across the zone, transform them whole in memory, then repeat through Dask to confirm the two paths agree exactly.

import numpy as np

rng = np.random.default_rng(seed=20260714)  # fixed seed → reproducible extract
n = 2_000_000

# Parcel centroids scattered across UTM zone 14N's longitude band.
lon = rng.uniform(-102.0, -96.0, size=n).astype(np.float64)
lat = rng.uniform(28.0, 34.0, size=n).astype(np.float64)

tf = build_transformer(6318, 6343)               # NAD83(2011) -> UTM 14N
east, north = transform_arrays(tf, lon, lat)

# A parcel exactly on the zone's central meridian (-99 deg) lands at
# easting 500000 m by definition of the transverse Mercator false easting.
e0, n0 = transform_arrays(tf, np.array([-99.0]), np.array([31.0]))
print(round(float(e0[0]), 3), round(float(n0[0]), 3))
# 500000.0 3429693.657  (northing approximate; easting exact on the CM)

Because the central meridian of UTM zone 14N is exactly $-99°$, any point there returns the false easting of 500000.0 m to the last representable digit — a free, exact correctness check baked into the projection definition. The northing depends on the latitude’s meridional arc and lands near 3.43 million metres at 31° N. Driving the same two million centroids through transform_dask and calling .compute() yields a frame whose easting/northing columns are element-wise identical to the in-memory arrays, because both branches run the same PROJ operation at float64.

Verification and Residual Analysis #

Vectorization must not change a single result bit versus the trusted per-point path, and out-of-extent points must be caught rather than shipped. The check below re-transforms a sample through a scalar loop, confirms the arrays agree within a sub-millimetre tolerance, counts any non-finite results, and emits a structured audit record. The residual against the reference is the planar separation

import json
import logging

logger = logging.getLogger("vectorize")


def verify_vectorized(
    transformer: Transformer,
    lon: np.ndarray,
    lat: np.ndarray,
    east: np.ndarray,
    north: np.ndarray,
    sample: int = 1000,
    tolerance_m: float = 1e-4,   # 0.1 mm: vectorized vs scalar must be exact
) -> dict[str, object]:
    """Compare a random sample of the vectorized result against per-point
    reference transforms and emit an audit record. Raises on drift."""
    idx = np.random.default_rng(0).choice(lon.size, size=min(sample, lon.size),
                                           replace=False)
    ref_e = np.empty(idx.size, dtype=np.float64)
    ref_n = np.empty(idx.size, dtype=np.float64)
    for k, i in enumerate(idx):  # scalar reference path, deliberately looped
        ref_e[k], ref_n[k] = transformer.transform(float(lon[i]), float(lat[i]))
    resid = np.hypot(east[idx] - ref_e, north[idx] - ref_n)
    non_finite = int(np.count_nonzero(~np.isfinite(east)))
    record = {
        "sample": int(idx.size),
        "max_residual_m": round(float(resid.max()), 6),
        "tolerance_m": tolerance_m,
        "non_finite_points": non_finite,
        "passed": bool(resid.max() <= tolerance_m),
    }
    logger.info("vectorize_verify %s", json.dumps(record))
    if not record["passed"]:
        raise ValueError(f"Vectorized drift {resid.max():.6f} m exceeds tolerance")
    return record


rec = verify_vectorized(tf, lon, lat, east, north)
assert rec["passed"] and rec["non_finite_points"] == 0

The max_residual_m should be 0.0 to floating-point noise, because the vectorized and scalar calls invoke the same operation on the same double-precision inputs — any non-zero value points to a dtype leak, not a maths difference. Pinning that residual against surveyed monuments rather than the solver’s own output belongs to the control-network discipline documented across Algorithmic Math & Geodetic Workflows, which supplies the residual statistics and error-ellipse machinery a submission ultimately cites.

Troubleshooting #

Frequently Asked Questions #

Related #

• Dask-partitioned datum shift for large cadastral datasets — the out-of-core datum shift worked end to end with partition sizing.
• Benchmarking NumPy vectorized vs looped transforms — measuring the speedup honestly and verifying identical results.
• Concurrent pyproj transformation pipelines in Python — distributing the same kernel across cores and reusing transformers safely.
• Algorithmic Math & Geodetic Workflows — the residual, error-ellipse, and adjustment machinery that certifies a batch.
• Batch Transformation & Automation for Cadastral Coordinate Pipelines — the parent reference on scaling, scheduling, and auditing transform runs.