CoordTransform Cadastral & High-Precision Python

Python Script for NADCON Datum Transformation

Transforming a NAD27 coordinate to NAD83 with NADCON grids is a deterministic latitude/longitude shift that must reproduce the published .las/.los node values to grid precision and agree with a known NGS control mark to within a survey-grade tolerance — typically a radial residual of ±0.01 m for Class I cadastral control under ISO 19111 coordinate-operation accuracy reporting. This page is the runnable NADCON implementation for the parent guide on choosing between NADCON and NTv2 datum-shift methods; where that comparison decides which grid family to use, this one builds the bilinear shift engine for the NADCON ASCII case and gates it the way an auditor would.

NADCON .las/.los transform decision flow An EPSG:4267 (NAD27) latitude/longitude pair enters a bounds check; coordinates outside the grid route to a deterministic fallback. In-bounds points locate their bilinear cell, interpolate the arc-second shifts from the .las and .los grids, apply the divide-by-3600 conversion, then pass a control-residual gate. Residuals over 0.01 metre route to the fallback; residuals within tolerance are rounded to 1e-8 degree and returned as an EPSG:4269 (NAD83) coordinate. Input (lat, lon) EPSG:4267 — latitude-first Inside grid bounds? out of grid in bounds Locate bounding cell row/col index fractions Bilinear interpolate Δφ from .las, Δλ from .los (arc-sec) Apply shift ÷ 3600 arc-seconds → degrees Control residual ≤ tolerance? > 0.01 m Deterministic fallback null-shift or abort within tolerance Round to 1e-8° deterministic Decimal quantize Return (lat, lon) EPSG:4269 — NAD83

Concept: What the NADCON Grid Shift Computes

NADCON stores the latitude shift Δφ and longitude shift Δλ as two co-registered grids of arc-second values — the .las file for latitude, the .los file for longitude — sampled on a regular geographic mesh. For an input coordinate that falls inside a cell, the shift is recovered by bilinear interpolation of the four surrounding nodes. With grid-index fractions trt_r (row) and tct_c (column) measured from the lower-left node, the interpolated shift is:

Δ(ϕ,λ)=Δ00(1−tr)(1−tc)+Δ10 tr(1−tc)+Δ01(1−tr)tc+Δ11 trtc \Delta(\phi,\lambda) = \Delta_{00}(1-t_r)(1-t_c) + \Delta_{10}\,t_r(1-t_c) + \Delta_{01}(1-t_r)t_c + \Delta_{11}\,t_r t_c

where tr=(ϕ−ϕ0)/δϕt_r = (\phi - \phi_0)/\delta\phi and tc=(λ−λ0)/δλt_c = (\lambda - \lambda_0)/\delta\lambda for grid origin (ϕ0,λ0)(\phi_0,\lambda_0) and node spacing (δϕ,δλ)(\delta\phi,\delta\lambda). The shifts are published in arc-seconds, so applying them to the source coordinate requires the divide-by-3600 conversion to degrees:

ϕ83=ϕ27+Δϕ′′3600,λ83=λ27+Δλ′′3600 \phi_{83} = \phi_{27} + \frac{\Delta\phi''}{3600}, \qquad \lambda_{83} = \lambda_{27} + \frac{\Delta\lambda''}{3600}

Two semantics matter for compliance. First, axis order is fixed: EPSG:4267 (NAD27) and EPSG:4269 (NAD83) are both latitude-first, longitude-second, and inverting them silently corrupts every record. Second, the shift magnitude (tens of metres for NAD27→NAD83) is not the residual — the residual is the leftover distance between the transformed point and an independent control coordinate in the target frame, and that is the quantity the ±0.01 m tolerance gates. Conflating the two makes every transform appear to fail. Epoch and vertical consistency are out of scope here and belong to validating datum alignment with control points; resolve those first so the residual reflects shift error rather than registration error.

Complete Runnable Implementation

The following self-contained pipeline parses the NADCON ASCII grids, interpolates the shift in deterministic Decimal arithmetic (eliminating the IEEE 754 drift that creeps into naive float rounding), applies the arc-second conversion, and — when a control coordinate is supplied — gates the result against an explicit survey-grade tolerance. Out-of-grid coordinates route to a deterministic fallback rather than silently extrapolating, the same defensive posture documented in fallback routing strategies for missing grid files. It runs as-is on Python 3.10+ with only the standard library.

import math
import logging
from dataclasses import dataclass
from decimal import Decimal, ROUND_HALF_UP, getcontext
from pathlib import Path
from typing import List, Optional, Union

# Deterministic high-precision context prevents IEEE 754 drift during rounding.
getcontext().prec = 28
getcontext().rounding = ROUND_HALF_UP

logger = logging.getLogger(__name__)

# 1 arc-second of latitude on the reference ellipsoid is ~30.8706 m.
_M_PER_ARCSEC_LAT = Decimal("30.8706")

Number = Union[float, str, Decimal]


@dataclass(frozen=True)
class NadconHeader:
    nrows: int
    ncols: int
    origin_lat: Decimal
    origin_lon: Decimal
    dlat: Decimal
    dlon: Decimal


@dataclass(frozen=True)
class TransformationResult:
    lat_out: Decimal
    lon_out: Decimal
    shift_lat_sec: Decimal
    shift_lon_sec: Decimal
    residual_m: Optional[Decimal]
    status: str  # "SUCCESS" | "FALLBACK_NULL_SHIFT"


class NadconPipeline:
    """ISO 19111-compliant NADCON (.las/.los) datum shift pipeline.

    Enforces EPSG:4267 (NAD27) -> EPSG:4269 (NAD83) axis ordering
    (latitude, longitude) and deterministic Decimal arithmetic. Shift grids
    are NADCON ASCII files expressed in arc-seconds.
    """

    def __init__(
        self,
        las_path: Path,
        los_path: Path,
        tolerance_m: Decimal = Decimal("0.01"),
        fallback_null_shift: bool = False,
    ) -> None:
        self.tolerance_m = tolerance_m
        self.fallback_null_shift = fallback_null_shift
        self.header = self._parse_grid_header(las_path)
        self.las_matrix = self._load_grid_matrix(las_path)
        self.los_matrix = self._load_grid_matrix(los_path)

    def _parse_grid_header(self, grid_path: Path) -> NadconHeader:
        """Parse the 3-line NADCON ASCII header: dims, origin, spacing."""
        with open(grid_path, "r", encoding="ascii") as fh:
            lines = [ln.strip() for ln in fh if ln.strip()]
        if len(lines) < 3:
            raise ValueError(f"Malformed NADCON header in {grid_path}")
        ncols, nrows = map(int, lines[0].split())
        origin_lon, origin_lat = map(Decimal, lines[1].split())
        dlon, dlat = map(Decimal, lines[2].split())
        return NadconHeader(nrows, ncols, origin_lat, origin_lon, dlat, dlon)

    def _load_grid_matrix(self, grid_path: Path) -> List[List[Decimal]]:
        """Load the arc-second shift values, skipping the 3-line header."""
        with open(grid_path, "r", encoding="ascii") as fh:
            body = [ln.strip() for ln in fh if ln.strip()][3:]
        matrix = [[Decimal(v) for v in ln.split()] for ln in body]
        if len(matrix) != self.header.nrows or any(
            len(row) != self.header.ncols for row in matrix
        ):
            raise ValueError(f"Grid dimension mismatch in {grid_path}")
        return matrix

    def _bilinear(self, lat: Decimal, lon: Decimal,
                  grid: List[List[Decimal]]) -> Decimal:
        """Bilinear interpolation of one shift grid (ISO 19111 sec. 9.3.2)."""
        row_idx = (lat - self.header.origin_lat) / self.header.dlat
        col_idx = (lon - self.header.origin_lon) / self.header.dlon
        # Clamp the lower node so the (r0, r0+1) cell stays inside the mesh.
        r0 = max(0, min(int(row_idx), self.header.nrows - 2))
        c0 = max(0, min(int(col_idx), self.header.ncols - 2))
        r1, c1 = r0 + 1, c0 + 1
        wr, wc = row_idx - r0, col_idx - c0  # interpolation weights
        one = Decimal(1)
        return (
            grid[r0][c0] * (one - wr) * (one - wc)
            + grid[r1][c0] * wr * (one - wc)
            + grid[r0][c1] * (one - wr) * wc
            + grid[r1][c1] * wr * wc
        )

    def _residual_to_control(self, lat_out: Decimal, lon_out: Decimal,
                             control_lat: Decimal,
                             control_lon: Decimal) -> Decimal:
        """Radial metric residual between the transformed point and a control mark."""
        lat_rad = float(lat_out) * math.pi / 180.0
        m_per_sec_lon = _M_PER_ARCSEC_LAT * Decimal(str(math.cos(lat_rad)))
        dy = (lat_out - control_lat) * Decimal(3600) * _M_PER_ARCSEC_LAT
        dx = (lon_out - control_lon) * Decimal(3600) * m_per_sec_lon
        return (dy * dy + dx * dx).sqrt()

    def _fallback(self, lat: Decimal, lon: Decimal,
                  reason: str) -> TransformationResult:
        """Deterministic fallback for out-of-bounds or tolerance violations."""
        logger.warning("NADCON fallback: %s at (%s, %s)", reason, lat, lon)
        if self.fallback_null_shift:
            return TransformationResult(
                lat, lon, Decimal("0"), Decimal("0"), None, "FALLBACK_NULL_SHIFT"
            )
        raise RuntimeError(f"Transformation aborted: {reason}. Fallback disabled.")

    def transform(self, lat_in: Number, lon_in: Number,
                  control_lat: Optional[Number] = None,
                  control_lon: Optional[Number] = None) -> TransformationResult:
        """Apply the NADCON shift. Axis order is (latitude, longitude).

        When a target-frame control coordinate is supplied, the residual is
        gated against ``tolerance_m``; otherwise only the bounds check applies.
        """
        lat, lon = Decimal(str(lat_in)), Decimal(str(lon_in))

        lat_max = self.header.origin_lat + (self.header.nrows - 1) * self.header.dlat
        lon_max = self.header.origin_lon + (self.header.ncols - 1) * self.header.dlon
        if not (self.header.origin_lat <= lat <= lat_max
                and self.header.origin_lon <= lon <= lon_max):
            return self._fallback(lat, lon, "COORDINATE_OUT_OF_GRID_BOUNDS")

        shift_lat = self._bilinear(lat, lon, self.las_matrix)
        shift_lon = self._bilinear(lat, lon, self.los_matrix)

        # Arc-seconds -> degrees, then deterministic round to 1e-8 deg.
        q = Decimal("0.00000001")
        lat_out = (lat + shift_lat / Decimal(3600)).quantize(q)
        lon_out = (lon + shift_lon / Decimal(3600)).quantize(q)

        residual: Optional[Decimal] = None
        if control_lat is not None and control_lon is not None:
            residual = self._residual_to_control(
                lat_out, lon_out, Decimal(str(control_lat)), Decimal(str(control_lon))
            )
            if residual > self.tolerance_m:
                return self._fallback(
                    lat, lon, f"RESIDUAL_EXCEEDED: {residual} m > {self.tolerance_m} m"
                )

        return TransformationResult(
            lat_out, lon_out, shift_lat, shift_lon, residual, "SUCCESS"
        )

Parameter and return reference

Name Type Units Valid range / meaning
las_path / los_path Path — NADCON ASCII latitude / longitude shift grids
tolerance_m Decimal metres Max radial residual gate; 0.01 Class I, 0.05 Class II
lat_in / lon_in float/str/Decimal degrees Source NAD27 coordinate, latitude-first (EPSG:4267)
control_lat / control_lon optional degrees Known NAD83 control mark for residual gating
shift_lat_sec / shift_lon_sec Decimal arc-seconds Interpolated Δφ, Δλ from the grids
lat_out / lon_out Decimal degrees NAD83 result rounded to 1e-8° (EPSG:4269)
residual_m Decimal / None metres Distance to control; None when no control supplied
status str — SUCCESS or FALLBACK_NULL_SHIFT

Minimal Worked Example

Synthetic 3×3 .las/.los grids over a Colorado test extent make the run reproducible without shipping the multi-megabyte NGS files. A linear node field means the bilinear interpolant is exact, so the expected shift can be predicted analytically and used to build the control point:

import tempfile
from decimal import Decimal
from pathlib import Path

ORIGIN_LON, ORIGIN_LAT, DLON, DLAT = -105.0, 39.0, 0.25, 0.25

def _lat_shift(r: float, c: float) -> float:   # arc-seconds, linear field
    return -0.10 - 0.01 * r - 0.005 * c

def _lon_shift(r: float, c: float) -> float:
    return -2.90 - 0.02 * r + 0.010 * c

def _write_grid(path: Path, fn) -> None:
    rows = [" ".join(f"{fn(r, c):.6f}" for c in range(3)) for r in range(3)]
    path.write_text(
        f"3 3\n{ORIGIN_LON} {ORIGIN_LAT}\n{DLON} {DLAT}\n" + "\n".join(rows) + "\n",
        encoding="ascii",
    )

tmp = Path(tempfile.mkdtemp())
las, los = tmp / "co.las", tmp / "co.los"
_write_grid(las, _lat_shift)
_write_grid(los, _lon_shift)

pipe = NadconPipeline(las, los, tolerance_m=Decimal("0.01"))

lat27, lon27 = 39.375, -104.625          # cell centre -> row_idx = col_idx = 1.5
q = Decimal("0.00000001")
ctrl_lat = (Decimal(str(lat27)) + Decimal(str(_lat_shift(1.5, 1.5))) / 3600).quantize(q)
ctrl_lon = (Decimal(str(lon27)) + Decimal(str(_lon_shift(1.5, 1.5))) / 3600).quantize(q)

res = pipe.transform(lat27, lon27, control_lat=ctrl_lat, control_lon=ctrl_lon)
print(res.lat_out, res.lon_out, res.shift_lat_sec, res.status)
# -> 39.37496597 -104.62580972 -0.12250000 SUCCESS

The interpolated latitude shift is exactly -0.1225″ (the mean of the four surrounding nodes, because the field is linear), giving a NAD83 latitude of 39.37496597°; the longitude shift -2.915″ moves the point by roughly 69 m west — far larger than the 0.01 m tolerance, which is precisely why the gate is applied to the control residual and not to the shift itself.

Validation Check

Gate the result with an explicit survey-grade assertion before any coordinate is written downstream:

assert res.status == "SUCCESS"
assert res.residual_m is not None and res.residual_m <= Decimal("0.01"), (
    f"Residual {res.residual_m} m exceeds survey-grade tolerance"
)
assert abs(res.shift_lat_sec - Decimal("-0.1225")) < Decimal("1e-9"), (
    "Bilinear interpolation did not reproduce the analytic node field"
)

For production cadastral runs, persist an audit record alongside the result carrying the grid filenames and version, the source/target EPSG codes (4267 → 4269), the interpolated Δφ/Δλ, the residual, and the pass/fail status — the reproducibility metadata that makes the transform legally defensible. The same control-residual discipline underpins validating coordinate precision to millimetre standards.

Common Mistakes

Inverting the axis order (longitude, latitude)
EPSG:4267 and EPSG:4269 are both latitude-first. Feeding a `(lon, lat)` tuple into a latitude-first pipeline indexes the wrong grid axis and applies Δφ where Δλ belongs — a systematic bias of tens of metres that still produces plausible-looking output. Keep the axis order explicit at every boundary and label coordinates by EPSG code, not by position in a tuple.
Forgetting the arc-second to degree conversion
NADCON `.las`/`.los` values are arc-seconds, not degrees. Adding a `-2.915` shift directly to a longitude in degrees moves the point by nearly three degrees instead of three arc-seconds. Always divide the interpolated shift by `3600` before applying it, and unit-test that a known node reproduces its published metre-scale displacement.
Gating the shift magnitude instead of the control residual
NAD27→NAD83 shifts are routinely tens of metres, so comparing the shift itself against a `0.01 m` tolerance makes every transform fail and trip the fallback. The tolerance belongs to the residual — the leftover distance between the transformed coordinate and an independent control mark in the target frame — which should collapse to a few millimetres on a correctly seated monument.