Concurrent pyproj Transformation Pipelines in Python

Transforming a few thousand coordinates is instantaneous; transforming a national cadastral extract of tens of millions is the throughput problem that Batch Transformation & Automation for Cadastral Coordinate Pipelines exists to solve, and concurrency is the lever that turns an overnight job into a coffee-break one. This page narrows that framework to a single, well-scoped operation: running pyproj.Transformer across many CPU cores at once without corrupting a single coordinate. The central difficulty is not writing a parallel loop — it is that a Transformer wraps a live PROJ C context that cannot be shared across a fork() boundary, so a naive process pool that inherits a parent-built transformer will segfault or, far worse, return plausible-looking garbage that passes straight into a survey deliverable. The correct pattern builds one transformer inside each worker, fans a large coordinate array out as ordered chunks, fans the transformed chunks back in while preserving input order exactly, and gates the reassembled result against a control-point residual before anything is written.

Figure — one input array fans out as indexed chunks, each worker builds its own Transformer, and results fan back in order-preserved before the residual gate.

Why concurrency here is a correctness problem, not just a speed one #

CPython runs one bytecode stream at a time under the Global Interpreter Lock, so the choice between threads and processes is not cosmetic. A pyproj.Transformer.transform call over a NumPy array does release the GIL while the underlying PROJ C routine runs, which means threads can overlap the C work — but the surrounding array marshalling, the Python-level chunk bookkeeping, and any pure-Python pre/post-processing remain serialized. For a transform-dominated batch you get real parallelism from processes because each process has its own interpreter, its own GIL, and — critically — its own PROJ context. The decision table below is the one that governs every pipeline you build.

The failure that trips up almost everyone comes from the interaction of os.fork() (the default process-start method on Linux) with the PROJ context. A Transformer is not just Python state: it holds a pointer into a PJ_CONTEXT, a C structure that owns file handles, an in-memory grid cache, and thread-affinity state. When you build a transformer in the parent and let a forked child inherit it, the child receives a copy of the pointer but not a valid, independent context. Using it races the parent, double-frees on teardown, or reads stale grid memory. The symptom is either an outright segfault that kills the worker, or — the dangerous case — silently wrong coordinates that are off by the exact magnitude of a datum shift and still look like coordinates. Transformer objects are therefore neither fork-safe nor picklable in a way you should rely on. The rule is absolute: never build a transformer in the parent and hand it to workers. Build one inside each worker, once, in an initializer. The same context-isolation discipline underpins the thread-scoped variant in thread-safe pyproj transformer caching.

Prerequisites and the ordering guarantee #

Concurrency reorders completion, not data. Workers finish in whatever order the scheduler and chunk sizes dictate, so if you append results as they arrive, chunk 7 may land before chunk 3 and every downstream row is silently misaligned against its parcel ID. The pipeline must therefore carry an explicit index with each chunk and reassemble by that index, not by arrival. Formally, if the input array $X$ of $N$ rows is partitioned into contiguous chunks $C_0, C_1, \dots, C_k$ where chunk $C_j$ spans rows $[s_j, e_j)$, then the output array $Y$ must satisfy

regardless of the order in which $T(C_j)$ is computed. Preserving this invariant is what lets a transformed coordinate keep its join to the surveyed attribute it belongs to. Choosing the chunk size $m$ is a trade-off: with $P$ workers and $N$ points, too-large chunks ($m \approx N/P$) starve late-arriving workers of work and unbalance the finish, while too-small chunks drown the run in pickling and dispatch overhead. A robust default targets a few chunks per worker, e.g. $m = \lceil N / (8P) \rceil$, bounded to a sane floor so tiny inputs never spawn thousands of micro-tasks.

The worked pipeline uses EPSG:4326 (WGS 84 geographic, latitude/longitude) as the source and EPSG:32633 (WGS 84 / UTM zone 33N, easting/northing in metres) as the target — a projection that a Central-European cadastral agency uses for planar parcel geometry. Because 4326 is a latitude-first CRS and 32633 is easting-first, we pin always_xy=True so every worker agrees that the array columns are (longitude, latitude) on input and (easting, northing) on output; axis-order drift across workers is a classic source of transposed deliverables, and the broader CRS-setup discipline is covered in setting up high-precision coordinate reference systems.

Step-by-Step: building a concurrent transform runner #

Each step below is a runnable, type-hinted unit. Together they form a ConcurrentTransformer that you can drop into a batch job. Build them in order.

Step 1 — Define the chunk task and result carriers #

The carrier objects are deliberately small and picklable: a chunk task ships an index, a start offset, and the coordinate slice; a result ships the same index and offset with the transformed slice. Nothing PROJ-related crosses the process boundary.

from __future__ import annotations

import logging
from dataclasses import dataclass

import numpy as np

logger = logging.getLogger("concurrent_transform")


@dataclass(frozen=True)
class ChunkTask:
    """One unit of work sent to a worker. `index` fixes reassembly order."""
    index: int          # 0-based chunk ordinal, monotonic with input order
    start: int          # first row of this chunk in the source array
    coords: np.ndarray  # (m, 2) float64 slice: columns (lon, lat) for EPSG:4326


@dataclass(frozen=True)
class ChunkResult:
    """A transformed chunk, tagged so the collector can restore order."""
    index: int          # echoes ChunkTask.index
    start: int          # echoes ChunkTask.start
    coords: np.ndarray  # (m, 2) float64: columns (easting, northing) for EPSG:32633

Step 2 — Build one Transformer per worker in an initializer #

ProcessPoolExecutor calls its initializer exactly once per worker process, before that worker handles any task. This is the only correct place to construct the Transformer: it runs inside the child, so the PROJ context is native to that process and never crosses a fork boundary. We stash it in a module-global keyed to the worker; the task function reads it back.

from pyproj import Transformer

# Process-global, populated once per worker by the initializer below.
_WORKER_TF: Transformer | None = None


def init_worker(source_epsg: str, target_epsg: str) -> None:
    """Runs once per worker process. Builds the Transformer INSIDE the child so
    the PROJ context (PJ_CONTEXT) is owned by this process and never inherited
    across a fork() — the invariant that keeps results deterministic."""
    global _WORKER_TF
    # always_xy=True forces (x, y) = (lon, lat) in / (easting, northing) out,
    # so every worker agrees on column order regardless of CRS axis metadata.
    _WORKER_TF = Transformer.from_crs(source_epsg, target_epsg, always_xy=True)
    logger.debug("worker Transformer built: %s -> %s", source_epsg, target_epsg)

Step 3 — Transform a chunk with the worker-local Transformer #

The task function is pure with respect to its input: it reads the process-global transformer, transforms the two columns as vectors (PROJ releases the GIL for this call), and returns a tagged result. It never mutates shared state.

def transform_chunk(task: ChunkTask) -> ChunkResult:
    """Transform one chunk using this worker's own Transformer. Vectorized:
    the whole (m,) column is passed to PROJ in a single call, not looped."""
    if _WORKER_TF is None:  # defensive: initializer must have run first
        raise RuntimeError("worker Transformer is not initialized")
    lon = np.ascontiguousarray(task.coords[:, 0], dtype=np.float64)
    lat = np.ascontiguousarray(task.coords[:, 1], dtype=np.float64)
    easting, northing = _WORKER_TF.transform(lon, lat)
    out = np.column_stack([easting, northing]).astype(np.float64, copy=False)
    return ChunkResult(index=task.index, start=task.start, coords=out)

Step 4 — Split the array into indexed chunks #

The chunker computes a chunk size from the worker count and yields ChunkTask objects carrying their own start offset. Contiguous, index-labelled chunks are what make order-preserving reassembly a one-liner later.

from collections.abc import Iterator


def chunk_array(coords: np.ndarray, n_workers: int,
                chunks_per_worker: int = 8, min_chunk: int = 1024) -> Iterator[ChunkTask]:
    """Split an (N,2) array into contiguous, index-tagged chunks. Target a few
    chunks per worker for load balance without drowning in dispatch overhead."""
    if coords.ndim != 2 or coords.shape[1] != 2:
        raise ValueError(f"expected (N, 2) array, got {coords.shape}")
    n = coords.shape[0]
    target = max(min_chunk, -(-n // max(1, n_workers * chunks_per_worker)))  # ceil div
    for index, start in enumerate(range(0, n, target)):
        stop = min(start + target, n)
        yield ChunkTask(index=index, start=start,
                        coords=np.ascontiguousarray(coords[start:stop]))

Step 5 — Fan out, fan in, and reassemble in order #

The runner ties it together. It preallocates the output array, dispatches every chunk to the pool with the initializer wired in, and writes each returned chunk back into its original slice by start offset. Because placement is by offset and not by arrival, the output order is deterministic no matter how the scheduler interleaves the workers.

from concurrent.futures import ProcessPoolExecutor, as_completed


def run_concurrent_transform(coords: np.ndarray, source_epsg: str,
                             target_epsg: str, n_workers: int = 4) -> np.ndarray:
    """Transform an (N,2) array concurrently, preserving input↔output order.
    Returns a new (N,2) float64 array aligned row-for-row with `coords`."""
    coords = np.ascontiguousarray(coords, dtype=np.float64)
    out = np.empty_like(coords)                       # preallocated, order-preserving
    tasks = list(chunk_array(coords, n_workers))
    with ProcessPoolExecutor(
        max_workers=n_workers,
        initializer=init_worker,                      # one Transformer per worker
        initargs=(source_epsg, target_epsg),
    ) as pool:
        futures = [pool.submit(transform_chunk, t) for t in tasks]
        for fut in as_completed(futures):             # completion order is arbitrary
            res = fut.result()
            out[res.start:res.start + res.coords.shape[0]] = res.coords  # placed by offset
    logger.info("transformed %d points across %d workers in %d chunks",
                coords.shape[0], n_workers, len(tasks))
    return out

Parameter and Return-Value Reference #

Worked Example: WGS 84 geographic to UTM zone 33N #

Take three points in central Europe and transform them from EPSG:4326 to EPSG:32633. The array is (lon, lat) because always_xy=True is set; the output is (easting, northing) in metres.

import numpy as np

pts = np.array([
    [13.4050, 52.5200],   # Berlin
    [12.4964, 41.9028],   # Rome (still within zone 33N usable span)
    [14.4378, 50.0755],   # Prague
], dtype=np.float64)      # columns: (longitude, latitude)

result = run_concurrent_transform(pts, "EPSG:4326", "EPSG:32633", n_workers=2)
print(np.round(result, 3))
# [[ 392788.163 5819943.399]
#  [ 292013.148 4642166.567]
#  [ 458337.727 5548439.019]]

The Berlin easting near 392 788 m and northing near 5 819 943 m is the reproducible expected value for zone 33N; a small discrepancy would signal an axis-order or datum-pipeline mismatch between workers. Selecting the transform operation deterministically — so every worker resolves the same pipeline rather than PROJ silently choosing a lower-accuracy path when a grid is absent — is the concern addressed in fallback routing strategies for missing grid files; pinning the operation is what makes a concurrent run bit-for-bit repeatable.

Verification and Residual Analysis #

A concurrent pipeline earns trust only when its aggregate output is pinned against an independent reference. The cheapest robust check is a round-trip: transform forward, transform back, and confirm every point returns to within a tight tolerance of where it started. Combined with a control-point comparison, this catches both a mis-built worker transformer and a reassembly bug in one pass. The per-point residual is the planar separation, and the batch statistic is the root-mean-square:

import json
import numpy as np


def verify_roundtrip(coords: np.ndarray, source_epsg: str, target_epsg: str,
                     n_workers: int = 4, tolerance_m: float = 0.001) -> dict[str, object]:
    """Forward then inverse transform; gate the max residual against tolerance.
    A concurrent bug (bad worker context, misordered reassembly) shows up here
    as either a huge residual or a row-shifted pattern."""
    fwd = run_concurrent_transform(coords, source_epsg, target_epsg, n_workers)
    back = run_concurrent_transform(fwd, target_epsg, source_epsg, n_workers)
    # Compare in a metric-ish sense on the geographic degrees round-trip:
    residual = np.abs(back - coords)
    max_deg = float(residual.max())
    record = {
        "points": int(coords.shape[0]),
        "workers": n_workers,
        "max_roundtrip_deg": round(max_deg, 12),
        "passed": bool(max_deg < 1e-9),   # sub-nanodegree ≈ sub-mm at these latitudes
    }
    logger.info("concurrent_verify %s", json.dumps(record))
    if not record["passed"]:
        raise ValueError(f"round-trip residual {max_deg:.3e} deg exceeds gate")
    return record


rec = verify_roundtrip(pts, "EPSG:4326", "EPSG:32633", n_workers=2)
assert rec["passed"]

The structured record — point count, worker count, residual, pass flag — is the minimal audit payload a batch stage should emit per run, and it composes directly into the certification trail built in the compliance workflows for agency submission. Pinning the same residual against surveyed monuments rather than a round-trip is the stricter check described in validating datum alignment with control points.

Troubleshooting and Gotchas #

Frequently Asked Questions #

Related References #

• Parallelizing Coordinate Transforms with multiprocessing — a self-contained ProcessPoolExecutor script with speedup measurement.
• Thread-Safe pyproj Transformer Caching — the thread-local variant for I/O-bound reuse.
• Fallback routing strategies for missing grid files — pin one operation so every worker resolves the same pipeline.
• Setting up high-precision coordinate reference systems — axis order and unit setup that keeps workers in agreement.
• Batch Transformation & Automation for Cadastral Coordinate Pipelines — the parent reference on scaling, scheduling, and audit for coordinate batches.