CoordTransform Cadastral & High-Precision Python

Dask-Partitioned Datum Shift for Large Cadastral Datasets

Applying a datum shift to tens of millions of cadastral coordinates that will not fit in RAM is the out-of-core case of Vectorizing Coordinate Transforms with NumPy and Dask, itself part of Batch Transformation & Automation for Cadastral Coordinate Pipelines: the same single vectorized Transformer.transform call is applied block by block with map_partitions, the transformer is built once per worker, and every column is pinned to float64 so the shift is never quietly degraded to metre-level noise. This page is a self-contained recipe for streaming a NAD27→NAD83(2011) datum shift across a dataset larger than memory without ever materialising the whole thing.

Out-of-core partitioned datum shift A Parquet dataset on disk is read lazily into three partitions. Each partition is handed to a worker that builds a pyproj Transformer once and applies the NAD27 to NAD83(2011) datum shift to its rows at float64. The shifted partitions stream straight back to a Parquet output on disk; the full dataset is never collected into the driver's memory. Parquet in > RAM Partition 0 · worker build transformer · shift Partition 1 · worker build transformer · shift Partition n · worker build transformer · shift Parquet out streamed

Figure — partitions read lazily, shift on their own worker, and stream back to disk without collecting in the driver.

What a Datum Shift Costs in Memory

A geographic datum shift such as NAD27 (EPSG:4267) to NAD83(2011) (EPSG:6318) is realised by PROJ as a grid-interpolated operation: each coordinate’s correction is read from a shift grid and added. Numerically the shift is small — typically tens of metres between NAD27 and NAD83 — but it is spatially non-uniform, which is exactly why it cannot be a single translation and why every point must be looked up. The compute is trivial; the constraint is memory. Storing NN coordinate pairs plus their shifted outputs at float64 needs roughly

M≈4×8×N  bytes=32N bytes,M \approx 4 \times 8 \times N \;\text{bytes} = 32N\,\text{bytes},

so 50 million parcels want about 1.6 GB just for the four live columns, before pandas overhead, index, and attribute payload. Dask’s answer is to never hold all of NN at once: it reads the source in partitions, applies the transform to each, and writes each result out, so peak memory scales with one partition, not the dataset. The identical concern of reusing a transformer safely when those partitions run truly concurrently is covered in concurrent pyproj transformation pipelines in Python.

Complete Runnable Implementation

The function below builds the transformer inside the mapped callable — once per partition, on the worker — applies the vectorized shift at float64, and declares an explicit meta so Dask never infers a float32 schema. It reads and writes Parquet so nothing but one partition is resident at a time.

from __future__ import annotations

import numpy as np
import pandas as pd
import dask.dataframe as dd
from pyproj import Transformer


def shift_partition(
    part: pd.DataFrame,
    src_epsg: int,
    dst_epsg: int,
    lon_col: str = "lon",
    lat_col: str = "lat",
) -> pd.DataFrame:
    """Apply a datum shift to one partition of coordinates.

    The Transformer is constructed here so each worker owns a native PROJ
    context (a PJ handle is not picklable and must not be shipped from the
    driver). Inputs are pinned to float64 to preserve the shift at the
    millimetre level demanded by cadastral tolerances.
    """
    transformer = Transformer.from_crs(src_epsg, dst_epsg, always_xy=True)
    x = part[lon_col].to_numpy(dtype=np.float64)
    y = part[lat_col].to_numpy(dtype=np.float64)
    shifted_lon, shifted_lat = transformer.transform(x, y)  # vectorized, in C
    out = part.copy()
    out["lon_shifted"] = np.asarray(shifted_lon, dtype=np.float64)
    out["lat_shifted"] = np.asarray(shifted_lat, dtype=np.float64)
    return out


def datum_shift_dataset(
    source_parquet: str,
    dest_parquet: str,
    src_epsg: int,
    dst_epsg: int,
    partition_rows: int = 1_000_000,
) -> None:
    """Stream a larger-than-RAM datum shift from one Parquet store to another.

    partition_rows sets the rows per partition; each partition's live arrays
    stay well under a worker's memory budget while the task count stays low
    enough to keep scheduler overhead negligible.
    """
    ddf = dd.read_parquet(source_parquet)
    ddf = ddf.repartition(partition_size=f"{partition_rows} rows")

    meta = ddf._meta.copy()
    meta["lon_shifted"] = np.float64()   # pin output dtype; never let Dask guess
    meta["lat_shifted"] = np.float64()

    shifted = ddf.map_partitions(
        shift_partition, src_epsg, dst_epsg, meta=meta,
    )
    # write_index=False keeps the on-disk schema clean for agency ingestion
    shifted.to_parquet(dest_parquet, write_index=False)

Nothing in datum_shift_dataset reads the whole dataset: read_parquet is lazy, map_partitions only records a graph node, and to_parquet drives the computation partition-by-partition, writing each result before the next is loaded.

Parameter Reference

Name Type Units Valid range Cadastral significance
source_parquet / dest_parquet str — path / URI on-disk stores; the dataset never lives fully in memory
src_epsg / dst_epsg int — valid EPSG codes the declared datum pair; wrong codes invalidate the whole deliverable
partition_rows int rows ≈ 0.5–5 M too large risks OOM per worker; too small floods the scheduler
lon_col / lat_col str — column names source coordinate columns; read as float64
lon_shifted / lat_shifted float64 degrees finite or -inf shifted coordinates; -inf marks a point outside the grid extent
meta pd.DataFrame — dtype template forces float64 outputs; omission risks a float32 down-cast

Minimal Worked Example

Synthesise a small stand-in for the real store, shift NAD27 to NAD83(2011), and inspect the magnitude of the correction. In production the same call runs against a multi-gigabyte Parquet directory unchanged.

rng = np.random.default_rng(seed=42)
frame = pd.DataFrame({
    "parcel_id": np.arange(8, dtype=np.int64),
    "lon": rng.uniform(-104.0, -95.0, size=8).astype(np.float64),
    "lat": rng.uniform(29.0, 36.0, size=8).astype(np.float64),
})
ddf = dd.from_pandas(frame, npartitions=2)      # two partitions for the demo

meta = ddf._meta.copy()
meta["lon_shifted"] = np.float64()
meta["lat_shifted"] = np.float64()
result = ddf.map_partitions(shift_partition, 4267, 6318, meta=meta).compute()

dlon_m = (result["lon_shifted"] - result["lon"]) * 111_320.0  # deg -> ~metres
print(round(float(dlon_m.abs().max()), 1))
# ~ tens of metres: the characteristic NAD27 -> NAD83 horizontal shift

The longitude difference scaled to metres is on the order of tens of metres, the well-known NAD27-to-NAD83 continental shift — the sanity signal that the datum operation actually ran rather than passing coordinates through unchanged.

Validation Check

Gate the streamed output before it is certified: the shift must be non-trivial (a zero shift means the operation silently no-op’d), every result must be finite, and the output dtype must still be float64.

assert result["lon_shifted"].dtype == np.float64, "output dtype down-cast"
assert np.isfinite(result[["lon_shifted", "lat_shifted"]].to_numpy()).all(), \
    "non-finite coordinate: point outside grid extent"
assert float(dlon_m.abs().max()) > 1.0, "shift is implausibly small; datum no-op?"

Common Mistakes

Building the Transformer in the driver and passing it into map_partitions
A Transformer wraps a native PROJ PJ handle that cannot be pickled; referencing a driver-built one inside the mapped function raises a serialisation error or, worse, corrupts state. Build it inside shift_partition so every worker constructs its own context exactly once per partition.
Partitions too large (OOM) or too small (scheduler overhead)
One partition's four float64 columns must fit comfortably in a worker's memory alongside the attribute payload, so a partition of hundreds of megabytes will crash a modest worker. Conversely, hundred-row partitions create millions of tasks whose scheduling dwarfs the transform cost. Target roughly half a million to a few million rows and tune from observed worker memory.
Losing precision to float32 through a missing meta
Without an explicit meta, Dask infers the output schema by running the function on a dummy frame and may assign float32 to the new columns, capping a metre-scale datum shift at roughly decimetre precision. Always declare meta with np.float64() columns and assert the dtype after compute.