Batch Transformation & Automation for Cadastral Coordinate Pipelines

Transforming a single parcel corner between datums is an exercise in geodetic arithmetic; transforming the millions of coordinates that make up a statewide cadastre is an exercise in determinism, reproducibility, and evidence. When a land-records agency reprojects an entire county from a legacy State Plane realization onto NAD83(2011) — or a national mapping authority migrates a dynamic reference frame forward by an epoch — the operation must produce byte-identical results on every re-run, preserve survey-grade tolerances across every point, and emit an audit trail that survives a boundary-retracement dispute years later. ISO 19111 governs the metadata every one of those transformed coordinates must carry: an unambiguous source and target coordinate reference system (CRS), the coordinate operation and its accuracy, and the coordinate epoch where the datum is time-dependent. This reference establishes the architecture for running those transformations at scale — for the surveyors, GIS engineers, and government technical teams who are answerable for the numbers. It covers how to cache and parallelize a pyproj.Transformer safely, how to vectorize the arithmetic across NumPy and Dask without leaking precision, and how to wrap the whole run in a compliance report an agency will accept. The foundational single-point mathematics this section builds on — CRS resolution, grid selection, and ISO 19111 metadata enforcement — is set out in the companion reference on core transformation fundamentals and standards.

The distinction that matters at scale is between a correct transformation and a defensible one. A research script that reprojects a GeoDataFrame in one line can be arithmetically correct and still be inadmissible: it may silently pick a lower-accuracy operation when a grid file is absent, reorder its output rows non-deterministically under parallelism, accumulate float32 truncation across a hundred million points, or leave no record of which PROJ operation actually ran. Each of those is a failure of process, not of geodesy, and each is what turns a bulk reprojection into a liability. The sections below treat the batch pipeline as an engineered, testable system with explicit gates rather than a convenience call.

The End-to-End Batch Pipeline at a Glance #

A production batch run is a fixed sequence of stages with two decision nodes: one that selects and caches the coordinate operation, and one that gates the transformed output against survey tolerance before anything is written to an agency deliverable. Coordinates enter, are partitioned into chunks that fit in memory, are transformed through a single cached operation per chunk, are checked against withheld control, and only then are hashed and packaged.

Figure — the end-to-end batch pipeline: ingest, cache one operation per CRS pair, partition into ordered chunks, transform, gate against control, then hash and package for submission.

Each stage in this diagram is treated in depth by a dedicated section of this reference. Concurrency and the safe reuse of the cached operation are covered in concurrent pyproj transformation pipelines in Python; the array mathematics that make the per-chunk transform fast without sacrificing precision are covered in vectorizing coordinate transforms with NumPy and Dask; the audit-hash and metadata stage that produces the deliverable is covered in compliance report generation for agency submission; and the decision of which transformation interface to build the pipeline on in the first place is covered in choosing a transformation API: PROJ, pyproj, or manual grids.

Why One Cached Transformer Per CRS Pair #

The single most consequential architectural decision in a batch pipeline is made before any coordinate is transformed: how the coordinate operation is constructed and reused. Building a fresh pyproj.Transformer for every point, or every chunk, is not merely slow — it is a correctness hazard. When PROJ resolves a source-to-target CRS pair it may find several candidate operations of differing accuracy, and the operation it selects depends on which grid files are present in the PROJ_DATA search path at that instant. Constructing the transformer once, up front, freezes that selection so that every chunk in the run shares an identical, reproducible operation. It also amortizes the non-trivial cost of parsing the CRS definitions and resolving the operation graph, which for a national CRS pair can dominate the per-point arithmetic by orders of magnitude.

Caching interacts directly with thread-safety, and here PROJ’s execution model must be respected rather than assumed. A PROJ context (PJ_CONTEXT) is not safe to share across threads: each thread that calls into PROJ needs its own context. A pyproj.Transformer object wraps such a context, so a single Transformer instance must not be called concurrently from multiple threads without synchronization. The pattern that scales is one cached Transformer per CRS pair per worker — a thread-local or process-local instance keyed on the ordered pair of CRS identifiers — so that the expensive operation-selection happens once per worker while no context is ever touched by two threads at once. The trade-offs between thread-local caching and full multiprocessing, and the exact locking discipline for sharing a cache, are worked through in thread-safe pyproj transformer caching and the multiprocessing variant in parallelizing coordinate transforms with multiprocessing.

Two further parameters must be pinned at construction. First, always_xy=True fixes the axis order to longitude-then-latitude (easting-then-northing) regardless of the authority’s declared axis order, so a latitude-first geographic CRS such as EPSG:4326 is never accidentally fed coordinates in the wrong order — a mistake that relocates a parcel by hundreds of kilometres and, crucially, does not raise an error. Second, the input arrays must be explicitly numpy.float64: PROJ operates in double precision internally, and feeding it a float32 array forces a promotion whose truncation is exactly the low-order millimetres that cadastral tolerance is graded against. The rationale for resolving each EPSG code to an explicit WKT2:2019 definition before caching — the deterministic anchoring that makes the cached operation auditable — is detailed in setting up high-precision coordinate reference systems.

Batch Architecture: Chunking, Ordering, and the Transformer Class #

The core of the pipeline is a class that owns the cached operation and streams the batch through it in bounded chunks. Chunking serves two independent purposes. It caps peak memory, so a batch of arbitrary size runs in fixed working set rather than materializing every intermediate array at once; and it provides the natural unit of parallelism, since each chunk is an independent transform. The non-negotiable constraint is deterministic ordering: the transformed output must appear in exactly the same row order as the input, on every run and regardless of how many workers processed it, because downstream a coordinate is joined back to its parcel record by position. The implementation below guarantees this by writing each chunk’s result back at its original index range rather than appending in completion order.

from __future__ import annotations

import logging
from dataclasses import dataclass
from typing import Iterator, Sequence

import numpy as np
from pyproj import CRS, Transformer

logger = logging.getLogger(__name__)


@dataclass(frozen=True)
class BatchResult:
    """Ordered transformation output for one batch run."""

    easting: np.ndarray
    northing: np.ndarray
    chunk_count: int
    point_count: int


class BatchTransformer:
    """Deterministic, chunked coordinate transformer for cadastral batches.

    A single pyproj.Transformer is built once per CRS pair and reused across
    every chunk, so the coordinate operation (ISO 19111 §C.5) is selected once
    and is identical for the whole run. always_xy=True pins the axis order and
    all arithmetic is carried in float64.
    """

    def __init__(
        self,
        source_epsg: int,
        target_epsg: int,
        chunk_size: int = 100_000,
    ) -> None:
        if chunk_size <= 0:
            raise ValueError("chunk_size must be a positive integer")
        self.source_crs: CRS = CRS.from_epsg(source_epsg)
        self.target_crs: CRS = CRS.from_epsg(target_epsg)
        self.chunk_size: int = chunk_size
        # Cache the operation once; the same PJ pipeline serves every chunk.
        self._transformer: Transformer = Transformer.from_crs(
            self.source_crs, self.target_crs, always_xy=True
        )

    @property
    def operation_wkt(self) -> str:
        """The selected coordinate operation as WKT, for the audit record."""
        return self._transformer.to_wkt()

    def _iter_chunks(self, n: int) -> Iterator[slice]:
        """Yield contiguous index slices covering [0, n) in input order."""
        for start in range(0, n, self.chunk_size):
            yield slice(start, min(start + self.chunk_size, n))

    def transform(
        self,
        x: Sequence[float],
        y: Sequence[float],
    ) -> BatchResult:
        """Transform an (x, y) batch chunk-by-chunk, preserving row order.

        Each chunk is written back at its ORIGINAL index range, so the output
        order is independent of chunk size and of any parallel scheduling.
        """
        xa = np.ascontiguousarray(x, dtype=np.float64)
        ya = np.ascontiguousarray(y, dtype=np.float64)
        if xa.ndim != 1 or xa.shape != ya.shape:
            raise ValueError("x and y must be 1-D arrays of equal length")

        n = int(xa.shape[0])
        out_x = np.empty(n, dtype=np.float64)
        out_y = np.empty(n, dtype=np.float64)

        chunk_count = 0
        for sl in self._iter_chunks(n):
            tx, ty = self._transformer.transform(xa[sl], ya[sl])
            out_x[sl] = np.asarray(tx, dtype=np.float64)
            out_y[sl] = np.asarray(ty, dtype=np.float64)
            chunk_count += 1
            logger.info(
                "chunk %d: transformed %d points at rows [%d:%d]",
                chunk_count, sl.stop - sl.start, sl.start, sl.stop,
            )
        return BatchResult(
            easting=out_x,
            northing=out_y,
            chunk_count=chunk_count,
            point_count=n,
        )


# Minimal worked example: NAD83(2011) geographic -> UTM zone 14N (EPSG:6343).
if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    bt = BatchTransformer(source_epsg=6318, target_epsg=6343, chunk_size=50_000)
    lon = [-99.1, -99.2, -99.3]
    lat = [31.0, 31.1, 31.2]
    result = bt.transform(lon, lat)
    assert result.point_count == 3
    assert result.easting.dtype == np.float64
    print(result.easting.round(3), result.northing.round(3))

The structured logging line is not decoration; it is the seam through which a batch run is observed and, later, reconstructed. Emitting the chunk index and the exact row range per chunk means an operator can confirm coverage — that every input row was transformed exactly once — and can locate the chunk responsible for any anomaly without re-running the batch. When the transform is distributed across a Dask cluster, the same slice-based write-back discipline is what keeps the assembled result deterministic even though partitions complete out of order; the partition-aware variant for very large datasets is developed in Dask-partitioned datum shift for large cadastral datasets.

A note on the throughput that justifies vectorizing at all: PROJ transforms whole arrays in a single call far faster than it transforms points one at a time, because the per-call overhead of crossing into the PROJ pipeline is paid once per chunk rather than once per point. Passing the full chunk array to transform — as the class above does — rather than looping in Python is the difference between a batch that finishes in seconds and one that finishes in hours. The measured speedups, and the array shapes at which they materialize, are quantified in benchmarking NumPy vectorized vs looped transforms.

Survey-Grade Precision & Tolerance Thresholds #

A batch does not pass or fail against a single number; it is graded against the tolerance budget of its survey class, and the pipeline’s gate enforces that budget per operation. The table below lists the thresholds the routing logic checks against and the residuals a correctly configured chain produces. Angular limits are converted to ground distance at mid-latitudes (1 arcsecond of latitude ≈ 30.9 m, so a 0.001″ limit ≈ 31 mm).

The gate in the pipeline diagram maps directly onto these rows. A batch that meets the geodetic-control row commits to the deliverable without review; a batch that meets the engineering row but not the cadastral row is routed to the quarantine set for manual adjudication rather than silently written. Because the criterion is an aggregate statistic over potentially millions of points, the pipeline computes it against a withheld set of control coordinates rather than the whole batch — the sampling and gating methodology is the same one calibrated in tuning transformation thresholds for survey grade.

The root-mean-square error the gate tests is defined over the per-point residual magnitudes $r_i$ across the control sample of size $n$:

where $(\hat{x}_i, \hat{y}_i)$ is the published control coordinate and $(x_i, y_i)$ is the pipeline’s transformed value. The maximum residual $\max_i r_i$ is tracked alongside the RMSE because a single gross outlier — one point transformed through a wrong operation or falling outside a grid’s extent — can pass an RMSE test while being individually inadmissible.

Compliance & Audit-Trail Generation #

A batch deliverable is only as defensible as the evidence attached to it. For legal admissibility the record must be tamper-evident and must bind together the three things a later reviewer needs to reproduce the run: the exact inputs, the exact operation parameters, and the exact outputs. The routine below folds all three into a single SHA-256 token. Hashing the raw float64 bytes of the input and output arrays — rather than a rounded string representation — captures full double precision, so any later edit to a single coordinate anywhere in the batch changes the digest.

from __future__ import annotations

import hashlib
import json
from datetime import datetime, timezone
from typing import Any

import numpy as np


def batch_audit_record(
    source_epsg: int,
    target_epsg: int,
    operation_wkt: str,
    inputs_xy: np.ndarray,
    outputs_xy: np.ndarray,
    epoch: float,
    rmse_m: float,
    max_residual_m: float,
) -> dict[str, Any]:
    """Bind inputs, parameters, and outputs into one reproducible token.

    ISO 19111 §C.4 — the CoordinateOperation and its evidence must be
    recoverable for every result. The SHA-256 is folded over the operation
    definition and the float64 byte views of the input and output arrays, so
    editing any single coordinate after the fact is detectable.
    """
    src = np.ascontiguousarray(inputs_xy, dtype=np.float64)
    dst = np.ascontiguousarray(outputs_xy, dtype=np.float64)
    if src.shape != dst.shape or src.ndim != 2 or src.shape[1] != 2:
        raise ValueError("inputs_xy and outputs_xy must be matching Nx2 arrays")

    hasher = hashlib.sha256()
    hasher.update(f"EPSG:{source_epsg}>EPSG:{target_epsg}".encode("utf-8"))
    hasher.update(operation_wkt.encode("utf-8"))
    # tobytes() preserves full precision; endianness is fixed to little so the
    # digest is reproducible across machines of the same byte order.
    hasher.update(src.astype("<f8").tobytes())
    hasher.update(dst.astype("<f8").tobytes())

    metadata: dict[str, Any] = {
        "source_crs": f"EPSG:{source_epsg}",
        "target_crs": f"EPSG:{target_epsg}",
        "coordinate_epoch": round(epoch, 3),
        "point_count": int(src.shape[0]),
        "rmse_m": round(rmse_m, 6),
        "max_residual_m": round(max_residual_m, 6),
        "generated_utc": datetime.now(timezone.utc).isoformat(timespec="seconds"),
        "audit_sha256": hasher.hexdigest(),
    }
    return metadata


if __name__ == "__main__":
    src = np.array([[-99.1, 31.0], [-99.2, 31.1]], dtype=np.float64)
    dst = np.array([[490000.12, 3430000.55], [480500.31, 3441000.90]])
    rec = batch_audit_record(
        source_epsg=6318, target_epsg=6343,
        operation_wkt="PROJCRS[...]",  # from BatchTransformer.operation_wkt
        inputs_xy=src, outputs_xy=dst,
        epoch=2010.00, rmse_m=0.0091, max_residual_m=0.0181,
    )
    assert len(rec["audit_sha256"]) == 64
    print(json.dumps(rec, indent=2))

The metadata dictionary this returns is the skeleton of the submission package. The fields an agency typically requires are the source and target CRS as authority codes, the coordinate epoch for a dynamic datum, the point count, the batch RMSE and maximum residual, a UTC generation timestamp, and the SHA-256 token that seals the run. Alongside these, the deliverable must version the artifacts the operation depended on — the WKT2:2019 CRS definitions and any grid binaries — because reproducing the run requires the exact grids that were on the PROJ_DATA path when it executed. The full submission-report schema, including the per-jurisdiction fields and the ISO 19111 metadata export, is developed in exporting ISO 19111 metadata for cadastral deliverables, and the hashing scheme is extended for incremental and re-run batches in generating audit hashes for transformation batches.

Validation & Control-Point Testing Across the Batch #

No batch is trustworthy until it has been checked against coordinates the pipeline never touched. The validation workflow withholds a set of independently surveyed control points that fall inside the batch extent, transforms them through the identical cached operation, and compares the result against their published values. The aggregate statistic is the RMSE of the residual magnitudes; the maximum residual is tracked separately to catch a single gross error that an average would mask. The function below computes both and returns the gate decision as a single structured result.

from __future__ import annotations

from dataclasses import dataclass

import numpy as np


@dataclass(frozen=True)
class ControlCheck:
    """Outcome of comparing a transformed batch to withheld control."""

    rmse_m: float
    max_residual_m: float
    n_control: int
    passed: bool


def validate_batch_against_control(
    transformed_xy: np.ndarray,
    control_xy: np.ndarray,
    tolerance_m: float = 0.02,
) -> ControlCheck:
    """Compare a transformed batch to published control coordinates.

    The pass rule follows the cadastral-boundary row: batch RMSE within
    tolerance and no single residual worse than twice that tolerance.
    """
    t = np.ascontiguousarray(transformed_xy, dtype=np.float64)
    c = np.ascontiguousarray(control_xy, dtype=np.float64)
    if t.shape != c.shape or t.ndim != 2 or t.shape[1] != 2:
        raise ValueError("both inputs must be matching Nx2 float arrays")

    residuals = np.hypot(t[:, 0] - c[:, 0], t[:, 1] - c[:, 1])
    rmse = float(np.sqrt(np.mean(residuals ** 2)))
    max_r = float(residuals.max())
    passed = rmse <= tolerance_m and max_r <= 2.0 * tolerance_m
    return ControlCheck(
        rmse_m=rmse,
        max_residual_m=max_r,
        n_control=int(t.shape[0]),
        passed=passed,
    )


# Worked example: State Plane metres, cadastral tolerance of 0.02 m.
transformed = np.array([[500123.412, 4649871.205],
                        [500987.004, 4650102.881],
                        [501450.660, 4650540.117]])
control = np.array([[500123.418, 4649871.199],
                    [500986.998, 4650102.889],
                    [501450.651, 4650540.126]])
check = validate_batch_against_control(transformed, control, tolerance_m=0.02)
assert check.passed, f"batch failed: RMSE={check.rmse_m:.4f} m"
print(f"RMSE {check.rmse_m:.4f} m over {check.n_control} control points")

A histogram of the residual array should be approximately symmetric about zero. A shifted mean signals a systematic error — a wrong operation, a missing coordinate epoch, or an unmodelled tectonic shift — rather than the random noise of measurement; the diagnostic separation between those causes is formalized in error distribution modeling in Python, and the estimation of the transformation parameters that a failed control check sends you back to is the subject of least squares adjustment for control networks. The RMSE-to-submission handoff — turning a passing ControlCheck into the numbers an agency form expects — is walked through end to end in the RMSE-to-agency-submission workflow in Python.

Failure Modes & Deterministic Mitigations #

The difference between a batch script and an audit-ready pipeline is how it behaves when something is wrong at scale. Each failure below is one that specifically emerges — or is amplified — when running in bulk and in parallel, and each has a single deterministic mitigation rather than a heuristic recovery.

• Transformer sharing across processes (not fork-safe). A pyproj.Transformer wraps a native PROJ context that must not be inherited across a fork() and used in the child. Constructing the transformer in the parent and relying on the worker to reuse it produces intermittent corruption or crashes. Build the cached transformer inside each worker after the fork — using an initializer that keys it on the CRS pair — so every process owns its context outright.
• CRS axis order. A latitude-first geographic CRS fed longitude-first coordinates transforms to a plausible-looking but grossly wrong location, and it does so silently. Pass always_xy=True at construction and assert the axis order from the resolved WKT2:2019 definition before the first chunk runs; do not rely on the input file’s column order matching the CRS.
• Silent lower-accuracy fallback when PROJ_DATA is unset. If the grid files a high-accuracy operation depends on are not on the PROJ_DATA search path, PROJ quietly selects a coarser operation and the batch completes without error at decimetre bias. Verify the required grids are present before the run and record the selected operation’s accuracy in the audit token; the routing discipline for a genuinely missing grid is set out in understanding NTv2 grid shift files in Python.
• Memory blow-up on huge arrays. Materializing every intermediate array for a hundred-million-point batch exhausts memory and triggers the OOM killer mid-run, leaving a partial, undated deliverable. Bound the working set with a fixed chunk_size and stream results to disk per chunk; for datasets that exceed a single machine, partition across a Dask cluster rather than growing the chunk.
• Non-deterministic output order under parallelism. Appending chunk results in completion order scrambles the row-to-parcel join and makes the run irreproducible. Write each chunk back at its original index range, as the BatchTransformer above does, so the output order is fixed regardless of scheduling.

Frequently Asked Questions #

Related #

• Concurrent pyproj transformation pipelines in Python — thread-safe caching and multiprocessing for the per-chunk transform.
• Vectorizing coordinate transforms with NumPy and Dask — array-level and partitioned transforms for datasets beyond one machine.
• Compliance report generation for agency submission — audit hashes, RMSE handoff, and ISO 19111 metadata export.
• Choosing a transformation API: PROJ, pyproj, or manual grids — selecting the interface the whole pipeline is built on.
• Core transformation fundamentals and standards — the single-point CRS, grid, and ISO 19111 foundations this section scales up.
• Algorithmic math and geodetic workflows — the estimation and error-propagation mathematics behind the tolerance gate.