CoordTransform Cadastral & High-Precision Python

Parallelizing Coordinate Transforms with multiprocessing

Turning a single-core pyproj batch into a multi-core one is the concrete speed win that Concurrent pyproj Transformation Pipelines in Python frames in general terms, and it belongs to the wider automation practice in Batch Transformation & Automation for Cadastral Coordinate Pipelines; this page is the self-contained script that does it correctly, using concurrent.futures.ProcessPoolExecutor with an initializer that builds exactly one pyproj.Transformer per worker process.

Chunk, transform, and write-back by offset A single (N,2) array is divided into three contiguous chunks labelled by index. Each chunk goes to a worker process that has its own Transformer. Each worker returns its transformed chunk tagged with a start offset, and the collector writes it into the preallocated output array at that offset so output order matches input order exactly. (N, 2) array chunk 0 → w0 chunk 1 → w1 chunk 2 → w2 write back by start offset output

Figure — split by index, transform per-worker, write back by offset.

Why the initializer is the whole trick

A pyproj.Transformer is not a plain value you can hand to a subprocess. It wraps a PROJ context — a C structure holding file handles and a grid cache — that is bound to the process that created it. On Linux the default process-start method is fork(), and a forked child inherits a copy of the pointer to that context without inheriting a valid, independent context. Transform through it and the child races the parent’s memory: the visible outcome is either a crashed worker or coordinates that are quietly wrong. The remedy is to never let a transformer cross the boundary. ProcessPoolExecutor accepts an initializer callable that runs once inside each worker before it handles any task; building the transformer there means the context is native to the child. We keep it in a module-global and read it back from the task function. The parallel efficiency you can expect follows Amdahl’s law: if a fraction pp of the work is parallelizable across PP workers, the achievable speedup is

S(P)=1(1−p)+pP,S(P) = \frac{1}{(1 - p) + \dfrac{p}{P}},

so the serial remainder — pickling chunks, reassembly — caps the gain and is exactly why chunk sizing matters. The same context-isolation rule, applied to threads rather than processes, is what drives the design in thread-safe pyproj transformer caching.

Complete Runnable Implementation

The script is self-contained: it defines the per-worker initializer, a chunk transform task, a splitter, and a driver that preserves order and can time itself against a single-process baseline.

from __future__ import annotations

import time
from concurrent.futures import ProcessPoolExecutor
from dataclasses import dataclass

import numpy as np
from pyproj import Transformer

# One Transformer per worker, created by the initializer, read by the task.
_TF: Transformer | None = None


def _init(source_epsg: str, target_epsg: str) -> None:
    """Run once per worker: build the Transformer INSIDE the child so its PROJ
    context is owned locally and never inherited across fork()."""
    global _TF
    _TF = Transformer.from_crs(source_epsg, target_epsg, always_xy=True)


@dataclass(frozen=True)
class _Chunk:
    start: int          # first row of this chunk in the source array
    coords: np.ndarray  # (m, 2) float64 slice


def _work(chunk: _Chunk) -> tuple[int, np.ndarray]:
    """Transform one chunk with the worker-local Transformer. Vectorized:
    each column is handed to PROJ in a single call, not looped point by point."""
    if _TF is None:
        raise RuntimeError("worker Transformer not initialized")
    x, y = _TF.transform(chunk.coords[:, 0], chunk.coords[:, 1])
    return chunk.start, np.column_stack([x, y]).astype(np.float64, copy=False)


def _split(coords: np.ndarray, n_chunks: int) -> list[_Chunk]:
    """Split an (N,2) array into `n_chunks` contiguous, offset-tagged pieces."""
    if coords.ndim != 2 or coords.shape[1] != 2:
        raise ValueError(f"expected (N, 2), got {coords.shape}")
    bounds = np.array_split(np.arange(coords.shape[0]), max(1, n_chunks))
    return [_Chunk(int(b[0]), np.ascontiguousarray(coords[b[0]:b[-1] + 1]))
            for b in bounds if b.size]


def parallel_transform(coords: np.ndarray, source_epsg: str, target_epsg: str,
                       n_workers: int = 4) -> np.ndarray:
    """Transform an (N,2) array across `n_workers` processes, preserving order.
    Output row i is the transform of input row i, regardless of finish order."""
    coords = np.ascontiguousarray(coords, dtype=np.float64)
    out = np.empty_like(coords)
    chunks = _split(coords, n_workers * 4)          # a few chunks per worker
    with ProcessPoolExecutor(max_workers=n_workers, initializer=_init,
                             initargs=(source_epsg, target_epsg)) as pool:
        for start, block in pool.map(_work, chunks):  # map preserves submit order
            out[start:start + block.shape[0]] = block  # but we place by offset anyway
    return out


def timed_speedup(coords: np.ndarray, source_epsg: str, target_epsg: str,
                  n_workers: int = 4) -> dict[str, float]:
    """Measure wall-clock speedup of the parallel run over a single-process run."""
    t0 = time.perf_counter()
    _init(source_epsg, target_epsg)                 # build once for the baseline
    serial = np.column_stack(_TF.transform(coords[:, 0], coords[:, 1]))
    t_serial = time.perf_counter() - t0

    t1 = time.perf_counter()
    par = parallel_transform(coords, source_epsg, target_epsg, n_workers)
    t_par = time.perf_counter() - t1

    assert np.allclose(serial, par, atol=1e-6), "parallel result diverged from serial"
    return {"serial_s": round(t_serial, 4), "parallel_s": round(t_par, 4),
            "speedup": round(t_serial / t_par, 2) if t_par else float("inf")}

Parameter Reference

Name Type Units Valid range Notes
coords np.ndarray degrees / metres (N, 2) float64 row order is the join key; must be preserved
source_epsg str — e.g. "EPSG:4326" source CRS handed to every worker’s Transformer
target_epsg str — e.g. "EPSG:32633" target CRS; wrong UTM zone offsets easting grossly
n_workers int processes 1 … CPU count begin at physical core count
_Chunk.start int row offset 0 … N−1 the offset each transformed block is written back to
return np.ndarray metres (N, 2) float64 aligned row-for-row with coords

Minimal Worked Example

Transform four points from EPSG:4326 (WGS 84 longitude/latitude) to EPSG:32633 (WGS 84 / UTM zone 33N) and confirm the parallel path matches the serial one.

import numpy as np

pts = np.array([
    [13.4050, 52.5200],   # Berlin
    [14.4378, 50.0755],   # Prague
    [11.5820, 48.1351],   # Munich
    [16.3738, 48.2082],   # Vienna
], dtype=np.float64)      # (longitude, latitude)

out = parallel_transform(pts, "EPSG:4326", "EPSG:32633", n_workers=2)
print(np.round(out[0], 2))
# [ 392788.16 5819943.4 ]   Berlin easting/northing in metres

Validation Check

Gate the parallel output against the single-process result before trusting it — this catches a mis-built worker context and a reassembly bug at once.

import numpy as np

serial = np.column_stack(
    Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True)
    .transform(pts[:, 0], pts[:, 1])
)
assert np.allclose(serial, out, atol=1e-6), "parallel diverged from serial baseline"
assert out.shape == pts.shape, "row count or order changed"
print("ok — parallel matches serial within 1e-6 m")

Common Mistakes

Forking a Transformer built in the parent process
Building Transformer.from_crs(...) once at module top level and relying on fork() to share it gives each child an invalid PROJ context — a segfault or silently garbage coordinates. Build the transformer inside the initializer so every worker owns its own context, and never pass a transformer as a task argument.
Losing input↔output order
Appending results as they arrive scrambles the row-to-parcel join because workers finish out of order. Either use pool.map, which yields results in submit order, or write each block back by its start offset into a preallocated array. Do both here for defence in depth.
Chunks so small the overhead dominates
Splitting an array into thousands of tiny chunks spends more time pickling and dispatching than transforming, and the parallel run ends up slower than serial. Target a handful of chunks per worker so each task carries tens of thousands of points; per Amdahl's law the serial pickling remainder is what caps your speedup.