CoordTransform Cadastral & High-Precision Python

Benchmarking NumPy Vectorized vs Looped Transforms

Quantifying how much faster a single vectorized array call is than a per-point Python loop — honestly, without the measurement artefacts that inflate or deflate the number — is the empirical companion to Vectorizing Coordinate Transforms with NumPy and Dask, under Batch Transformation & Automation for Cadastral Coordinate Pipelines: the benchmark must warm the PROJ grid cache, exclude transformer construction from the timed region, hold both paths at float64, and prove the two results are identical to survey tolerance before any speedup figure is believed. A benchmark that reports a fast loop usually just forgot to warm the cache; one that reports a suspiciously huge win usually timed the object construction only once.

Benchmark harness for looped versus vectorized transforms One coordinate array feeds two timed branches. The upper branch is a per-point Python loop calling transform once per coordinate; the lower branch is a single vectorized transform call over the whole array. Each branch's elapsed time becomes a points-per-second throughput, and the two output arrays are checked for element-wise agreement within survey tolerance. Warm array float64, cache hot Looped: N calls perf_counter timed Vectorized: 1 call perf_counter timed points/sec + results equal?

Figure — one warm array feeds a looped and a vectorized path; both are timed, converted to throughput, and checked for identical output.

What a Fair Benchmark Measures

Throughput is the honest currency, not raw seconds, because it normalises across dataset sizes. For NN points transformed in Δt\Delta t seconds,

T=NΔt[points⋅s−1],S=TvecTloop=ΔtloopΔtvec,T = \frac{N}{\Delta t}\quad[\text{points·s}^{-1}], \qquad S = \frac{T_{\text{vec}}}{T_{\text{loop}}} = \frac{\Delta t_{\text{loop}}}{\Delta t_{\text{vec}}},

where SS is the speedup. Two systematic errors dominate naive benchmarks. The first is a cold PROJ cache: the very first transform after constructing a Transformer may resolve and memory-map grid files, so timing that first call charges one-time I/O to the transform. The fix is a warm-up call before the clock starts. The second is including Transformer.from_crs inside the timed region — object construction is amortised over the whole batch in production, so timing it per iteration slanders the loop and flatters nothing. The construction cost belongs to setup; the precision setup that makes both paths comparable is established in setting up high-precision coordinate reference systems, and the production concurrency numbers this feeds appear in concurrent pyproj transformation pipelines in Python.

Complete Runnable Implementation

The harness builds the transformer once, warms the cache, times each path with time.perf_counter over repeated trials taking the minimum (the least noise-contaminated estimate), and returns a structured result. Both paths consume the same float64 array, so any output difference is a real defect rather than a dtype artefact.

from __future__ import annotations

import time
from dataclasses import dataclass

import numpy as np
from pyproj import Transformer


@dataclass(frozen=True)
class BenchmarkResult:
    n_points: int
    loop_seconds: float
    vector_seconds: float
    loop_throughput: float     # points per second
    vector_throughput: float
    speedup: float
    max_abs_diff_m: float      # vectorized vs looped, metres
    results_match: bool


def _time_loop(tf: Transformer, x: np.ndarray, y: np.ndarray
               ) -> tuple[np.ndarray, np.ndarray, float]:
    """Per-point path: one transform call per coordinate, timed."""
    out_e = np.empty_like(x)
    out_n = np.empty_like(y)
    start = time.perf_counter()
    for i in range(x.size):
        out_e[i], out_n[i] = tf.transform(float(x[i]), float(y[i]))
    return out_e, out_n, time.perf_counter() - start


def _time_vector(tf: Transformer, x: np.ndarray, y: np.ndarray
                 ) -> tuple[np.ndarray, np.ndarray, float]:
    """Vectorized path: one transform call over the whole array, timed."""
    start = time.perf_counter()
    e, n = tf.transform(x, y)
    elapsed = time.perf_counter() - start
    return np.asarray(e, np.float64), np.asarray(n, np.float64), elapsed


def benchmark_transform(
    src_epsg: int,
    dst_epsg: int,
    x: np.ndarray,
    y: np.ndarray,
    repeats: int = 3,
    tolerance_m: float = 1e-4,   # 0.1 mm: the two paths must agree exactly
) -> BenchmarkResult:
    """Compare looped vs vectorized transform throughput on identical inputs.

    The Transformer is built once outside the timed region and the grid
    cache is warmed before any measurement, so neither construction nor
    first-call I/O contaminates the timings.
    """
    tf = Transformer.from_crs(src_epsg, dst_epsg, always_xy=True)  # setup, untimed
    x = np.ascontiguousarray(x, dtype=np.float64)
    y = np.ascontiguousarray(y, dtype=np.float64)
    _ = tf.transform(x[:1], y[:1])  # warm the PROJ grid cache

    loop_t = min(_time_loop(tf, x, y)[2] for _ in range(repeats))
    vec_e, vec_n, vec_t = min(
        (_time_vector(tf, x, y) for _ in range(repeats)), key=lambda r: r[2]
    )
    ref_e, ref_n, _ = _time_loop(tf, x, y)

    diff = np.hypot(vec_e - ref_e, vec_n - ref_n)
    max_diff = float(diff.max())
    return BenchmarkResult(
        n_points=int(x.size),
        loop_seconds=loop_t,
        vector_seconds=vec_t,
        loop_throughput=x.size / loop_t,
        vector_throughput=x.size / vec_t,
        speedup=loop_t / vec_t,
        max_abs_diff_m=max_diff,
        results_match=max_diff <= tolerance_m,
    )

Parameter Reference

Name Type Units Valid range Significance
src_epsg / dst_epsg int — valid EPSG codes the CRS pair benchmarked; both paths use the identical operation
x / y np.ndarray degrees (or CRS units) finite float64 shared inputs; feeding float32 here would make the paths differ spuriously
repeats int count ≥ 1 trials per path; the minimum is reported to suppress OS scheduling noise
tolerance_m float metres 1e-4 typical the two paths must agree within this; they should match to float noise
loop_throughput / vector_throughput float points·s⁻¹ > 0 comparable across dataset sizes, unlike raw seconds
speedup float — > 1 expected vectorized advantage; routinely one to two orders of magnitude
results_match bool — — True proves vectorizing did not change the maths

Minimal Worked Example

Benchmark a NAD83(2011) geographic to UTM zone 14N projection, EPSG:6318 to EPSG:6343, over one hundred thousand points — small enough that the looped path finishes quickly, large enough that the throughput gap is unmistakable.

rng = np.random.default_rng(seed=7)
lon = rng.uniform(-102.0, -96.0, size=100_000).astype(np.float64)
lat = rng.uniform(28.0, 34.0, size=100_000).astype(np.float64)

result = benchmark_transform(6318, 6343, lon, lat)
print(int(result.vector_throughput), int(result.loop_throughput),
      round(result.speedup, 1), result.results_match)
# e.g. 9500000 140000 68.0 True
#   vectorized ~ millions/sec, loop ~ hundred-thousands/sec, ~50-100x, exact

Absolute numbers vary with hardware, but the shape is invariant: vectorized throughput lands in the millions of points per second while the loop crawls in the hundred-thousands, a speedup of roughly one to two orders of magnitude, and results_match is True because both paths ran the same PROJ operation at float64.

Validation Check

The benchmark is only meaningful if the two paths are numerically identical; a speedup on non-matching results is worthless. Assert equality within tolerance and confirm the vectorized path actually won.

assert result.results_match, "vectorized and looped outputs disagree"
assert result.max_abs_diff_m < 1e-4, "unexpected drift between paths"
assert result.speedup > 1.0, "vectorized path was not faster; check the harness"

Common Mistakes

Timing the first call with a cold PROJ or grid cache
The first transform after constructing a Transformer can resolve and memory-map grid files, charging one-time I/O to the measurement. Always issue a warm-up transform on a slice of the data before starting the clock, as the harness does, and report the minimum of several repeats to suppress scheduler noise.
Including Transformer construction inside the timed loop
Calling Transformer.from_crs inside the timed region charges object construction to every iteration, which never happens in production where the transformer is built once and reused. Build it once outside the clock so the benchmark reflects the real per-point transform cost.
Comparing a float32 loop against a float64 vectorized call
If one path down-casts to float32, the two outputs diverge by up to decimetres and results_match fails for a reason that has nothing to do with speed. Coerce both inputs to float64 up front so the comparison isolates performance, not precision.