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.

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 $t_r$ (row) and $t_c$ (column) measured from the lower-left node, the interpolated shift is:

Figure — a NADCON transformation reads two files, and both must agree on their header.

where $t_r = (\phi - \phi_0)/\delta\phi$ and $t_c = (\lambda - \lambda_0)/\delta\lambda$ for grid origin $(\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:

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 #

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:

Figure — a NAD27 to NAD83 shift profile across one degree of longitude (typical CONUS magnitudes).

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 #

Frequently Asked Questions #

What does a NADCON shift value of exactly -999 mean? #

It is the no-data sentinel: that node lies outside the modelled area, and the correct response is to reject the query rather than interpolate through it. A reader that treats -999 as a shift in arc-seconds will move the point roughly 30 kilometres, which is obvious; a reader that averages it with three valid neighbours moves the point about 7.5 kilometres and produces a coordinate that still looks like a coordinate. Check every one of the four surrounding nodes before interpolating, not just the nearest.

Why is my longitude shift positive when I expected it negative? #

NADCON follows the geodetic convention that longitude is positive west, so a shift that moves a point east appears as a negative value and vice versa, and the sign flips again if your code stores longitude positive east. This is the single most common NADCON bug and it produces an error of roughly twice the shift — three or four metres in the conterminous United States, which is large enough to matter and small enough to be mistaken for a datum problem.

Can NADCON transform heights? #

No. NADCON models the horizontal difference between NAD27 and NAD83 only. Vertical datum transformations use a separate model — a geoid or vertical shift grid — and are a different operation with a different accuracy statement. Combining them in one step is fine as a pipeline, but they must be two documented operations in the audit record, because their uncertainties are independent and must be reported separately.

Related References #

• NADCON vs NTv2: choosing the right datum shift — the parent guide that decides when this NADCON engine is the correct tool.
• Core transformation fundamentals and standards — the overarching reference on ISO 19111-compliant, audit-ready transformation pipelines.
• Understanding NTv2 grid shift files in Python — the binary .gsb counterpart to NADCON’s ASCII grids.
• Validating datum alignment with control points — clearing registration error before reading any shift residual.
• Fallback routing strategies for missing grid files — deterministic handling when a coordinate falls outside grid coverage.