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

Where the Transformer is built in a process pool The parent process builds nothing but the CRS pair. Each worker process runs the pool initialiser exactly once and builds its own Transformer there, storing it in a module global. Every task the worker then receives is a chunk of coordinates, not a Transformer, so nothing PROJ-related crosses the process boundary. The result is that the expensive operation lookup happens once per worker rather than once per chunk, and no unpicklable object is ever sent. Parent process holds only the CRS pair Pool initialiser — once per worker builds the Transformer Worker 1 own Transformer Worker 2 own Transformer Worker n own Transformer Tasks carry coordinates only never a Transformer

Figure — build the Transformer once per worker, not once per task.

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.

Total wall-clock against chunk size for one million points Wall-clock seconds against chunk size for one million points across eight workers. At one thousand points per chunk the job takes 4.1 seconds, dominated by dispatch overhead across a thousand tasks. The minimum, about 0.7 seconds, sits between fifty and one hundred thousand points per chunk. At five hundred thousand it rises again to 1.6 seconds because only two chunks exist and six workers sit idle. 0 1 2 3 4 1 50 100 250 500 chunk size (thousands of points) s wall-clock (s)

Figure — chunk size against total time: too small pays dispatch, too large loses balance.

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.

Frequently Asked Questions

Why can’t I just pass a Transformer to each worker?

Because it wraps a native PROJ context that does not survive being pickled and sent across a process boundary. Depending on the platform and start method you will get a pickling error, or — worse — an object that appears to arrive and then misbehaves. Building it inside each worker sidesteps the question entirely and is also faster, because the construction cost is paid once per worker rather than once per task.

How large should each chunk be?

Large enough that dispatch overhead disappears and small enough that every worker gets several chunks. For coordinate transforms, tens of thousands of points per chunk is usually the flat part of the curve. Two chunks across eight workers is the failure mode at the other end: six workers idle while two do all the work.

Does multiprocessing change the numerical result?

It must not, and if it does you have a bug. The arithmetic is per point and independent, so worker count cannot affect any individual coordinate. What does change results is rounding or aggregation performed per chunk — round after the transform, on the assembled output, and the result becomes independent of how the work was divided.