Datum transformations

The 1958 as-built drawings that the district scanned into its GIS in GEOM-330 were tied to survey monuments whose coordinates had been published on NAD27. The district's control is on NAD83(2011), its imagery came with coordinates on WGS 84, and the consultant who surveyed the reservoir in 1995 worked on a 1990s high-accuracy readjustment of NAD83. Putting these together correctly needs datum transformations, which move coordinates from one geodetic datum or frame to another. Geodesy I introduced datums, reference frames and the seven-parameter Helmert transformation. This lesson builds on that foundation with the methods used in practice: grid-based transformations for the distorted older datums, time-dependent Helmert transformations for modern frames, and the accuracy of each.

Conversions, transformations and why the difference matters

Two kinds of operation change coordinates. A conversion changes the way coordinates are expressed within one datum: from latitude and longitude to UTM, from geodetic to ECEF, from meters to feet. Conversions are exact, defined by formulas, and reversible. A transformation changes the datum itself: from NAD27 to NAD83, or from ITRF2020 to a NAD83 realization. Transformations are empirical, derived from observations at points known in both datums, and carry an accuracy that must be stated. The EPSG registry, discussed in lesson 11, keeps the two apart: a conversion has no accuracy, a transformation has a stated accuracy in meters. Changing projection is a conversion; changing datum is a transformation, and a chain from one projected CRS to another on a different datum includes both.

Why old datums need grids

NAD27 was computed in the 1920s and 1930s from triangulation adjusted piece by piece, with its origin at Meades Ranch in Kansas and the Clarke 1866 ellipsoid. Its coordinates contain distortions of up to several meters that vary irregularly across the continent, because each part of the network was fitted to the parts already adjusted. NAD83, adjusted simultaneously with Doppler satellite observations on the GRS80 ellipsoid, moved points by tens of meters to more than a hundred meters. The shift varies smoothly at the large scale but irregularly at the scale of meters, so no set of seven Helmert parameters can reproduce it to better than several meters over a large area.

The solution is a grid of shifts. At each node of a regular grid of latitude and longitude, the shift in latitude and longitude from the old datum to the new is stored, computed from many points known in both. A point's shift is found by interpolating between the four surrounding nodes, usually bilinearly:

where \(x\) and \(y\) are the point's fractional position in the cell. The same method serves for vertical datum changes and geoid models, which Geodesy I described.

In the United States, NGS's NADCON tool and its grids perform the NAD27-to-NAD83 transformation, and its current version, NADCON5, extends the same approach to the later NAD83 realizations: the original NAD83(1986), the HARN or HPGN statewide readjustments of the 1990s, NAD83(CORS96), NAD83(2007) and NAD83(2011). Each step has its own grid and its own stated accuracy. NADCON5 grids are distributed by NGS and are used by GIS software through the PROJ library. Accuracy depends on the region and the step: of the order of tens of centimeters to a meter or two for NAD27 to NAD83(1986), and centimeters for steps between the high-accuracy realizations.

import math
# 1. A grid shift: interpolate latitude and longitude shifts (arc seconds) stored at the
#    nodes of a regular grid, as NADCON grid files do. The node values below are
#    exercise values of realistic size for central Montana, not the published grid.
lat_nodes = [47.25, 47.50, 47.75]                 # degrees, grid spacing 0.25 deg
lon_nodes = [-111.50, -111.25, -111.00]
dlat = [[0.318, 0.315, 0.312],                    # rows: lat_nodes, columns: lon_nodes
        [0.322, 0.319, 0.316],
        [0.326, 0.323, 0.320]]
dlon = [[-3.412, -3.399, -3.386],
        [-3.419, -3.406, -3.393],
        [-3.427, -3.414, -3.401]]

def bilinear(grid, lat, lon):
    i = max(0, min(len(lat_nodes) - 2, int((lat - lat_nodes[0]) // 0.25)))
    j = max(0, min(len(lon_nodes) - 2, int((lon - lon_nodes[0]) // 0.25)))
    y = (lat - lat_nodes[i]) / 0.25; x = (lon - lon_nodes[j]) / 0.25
    return ((1-x)*(1-y)*grid[i][j] + x*(1-y)*grid[i][j+1] + (1-x)*y*grid[i+1][j] + x*y*grid[i+1][j+1])

lat27, lon27 = 47.50591560, -111.29684500          # an old NAD27 monument near CR1 (exercise)
sla, slo = bilinear(dlat, lat27, lon27), bilinear(dlon, lat27, lon27)
lat83, lon83 = lat27 + sla / 3600, lon27 + slo / 3600
m_lat = 111_200 * sla / 3600; m_lon = 111_320 * math.cos(math.radians(lat27)) * slo / 3600
print(f"grid shift: dlat = {sla:+.4f}\" ({m_lat:+.2f} m), dlon = {slo:+.4f}\" ({m_lon:+.2f} m)")
print(f"NAD27 {lat27:.8f}, {lon27:.8f}  ->  NAD83 {lat83:.8f}, {lon83:.8f}")
print(f"total horizontal shift {math.hypot(m_lat, m_lon):.2f} m")

# 2. A time-dependent (14-parameter) Helmert transformation between two frames.
#    Parameters at the reference epoch plus their rates; illustrative values of
#    the size used between ITRF and NAD83, not an official parameter set.
P0 = dict(tx=1.0053, ty=-1.9092, tz=-0.5416, s=-0.05e-9,          # m, m, m, unitless
          rx=26.785e-3, ry=-0.420e-3, rz=10.932e-3)                # arc seconds
RATE = dict(tx=0.0, ty=0.0, tz=0.0, s=0.0,
            rx=0.0532e-3, ry=-0.7423e-3, rz=-0.0316e-3)           # per year
t0 = 2010.0
def helmert14(X, Y, Z, t):
    p = {k: P0[k] + RATE[k] * (t - t0) for k in P0}
    r = {k: math.radians(p[k] / 3600) for k in ("rx", "ry", "rz")}
    s1 = 1 + p["s"]
    Xn = p["tx"] + s1 * (X - r["rz"] * Y + r["ry"] * Z)
    Yn = p["ty"] + s1 * (r["rz"] * X + Y - r["rx"] * Z)
    Zn = p["tz"] + s1 * (-r["ry"] * X + r["rx"] * Y + Z)
    return Xn, Yn, Zn

a = 6378137.0; f = 1 / 298.257222101; e2 = f * (2 - f)
def to_xyz(lat, lon, h):
    p, l = math.radians(lat), math.radians(lon); N = a / math.sqrt(1 - e2 * math.sin(p)**2)
    return (N + h) * math.cos(p) * math.cos(l), (N + h) * math.cos(p) * math.sin(l), (N * (1 - e2) + h) * math.sin(p)
X, Y, Z = to_xyz(47.50600422, -111.29779092, 1004.826)
print("\nframe A -> frame B at three epochs (shift of the same point, m):")
for t in (2010.0, 2020.0, 2026.8):
    Xn, Yn, Zn = helmert14(X, Y, Z, t)
    print(f"  epoch {t:6.1f}: dX = {Xn-X:+.4f}  dY = {Yn-Y:+.4f}  dZ = {Zn-Z:+.4f}  "
          f"|d| = {math.sqrt((Xn-X)**2 + (Yn-Y)**2 + (Zn-Z)**2):.4f}")

grid shift: dlat = +0.3197" (+9.87 m), dlon = -3.4086" (-71.20 m)
NAD27 47.50591560, -111.29684500  ->  NAD83 47.50600439, -111.29779184
total horizontal shift 71.88 m

frame A -> frame B at three epochs (shift of the same point, m):
  epoch 2010.0: dX = +1.2090  dY = -2.5999  dZ = -1.0673  |d| = 3.0595
  epoch 2020.0: dX = +1.0344  dY = -2.6096  dZ = -1.1341  |d| = 3.0276
  epoch 2026.8: dX = +0.9157  dY = -2.6162  dZ = -1.1796  |d| = 3.0123

The interpolated shift at the old monument is +0.320\" in latitude and −3.409\" in longitude, or 9.9 m north and 71.2 m west: the old monument's NAD83 position is 71.9 m from its NAD27 coordinates. The grid values in the code are exercise values chosen to be of the right size for central Montana; for real work, the published grid is used through software. The monument's NAD83 coordinates come out within a few centimeters of CR1, which is why the district could tie the 1958 drawings to its modern control once the grid transformation had been applied, but not before.

Time-dependent Helmert transformations

Modern frames are related by Helmert transformations whose parameters change with time, because frames such as ITRF are realized from stations that move with the tectonic plates, while frames such as NAD83 are attached to the North American plate. The fourteen-parameter form has the seven parameters of Geodesy I, three translations \(T_X, T_Y, T_Z\), a scale change \(s\) and three rotations \(R_X, R_Y, R_Z\), each at a reference epoch \(t_0\), plus a rate of change for each:

The rotation rates between ITRF and a plate-fixed NAD83 are essentially the rotation of the North American plate, about 0.7 milliarcseconds per year about an axis close to the Earth's equatorial plane. The epoch \(t\) at which the transformation is evaluated must be the epoch of the coordinates being transformed. Evaluating the same transformation at the wrong epoch moves points by about 2 cm for every year of error at the district, the speed of the plate there.

The code applies a fourteen-parameter transformation, with parameters of the right size but not an official set, to CR1 at three epochs. The total shift is about 3 m, and its components change by several centimeters per year, about 0.3 m between 2010 and 2026.8 in the X component alone. Those changes are exactly the plate motion that NAD83 removes, which lesson 10 examines.

Sign conventions and parameter sets

Helmert parameters are published in two conventions that differ in the sign of the rotations: the 'position vector' convention used by the IERS and ISO, in which the rotation matrix is as written above, and the 'coordinate frame' convention, which uses its transpose. Applying parameters in the wrong convention reverses the rotations and produces errors of meters. Parameters also come with a direction: from frame A to frame B, or the reverse, and inverting a transformation means more than changing the signs of all the parameters when the rotations and scale are not tiny. EPSG records the convention with each transformation, as 'Position Vector transformation' or 'Coordinate Frame rotation'. Hand-entered parameters should always be checked by transforming a published test point.

Accuracy and the choice of transformation

Between any two datums there may be several published transformations, differing in method, area of use and accuracy. NAD27 to NAD83, for example, can be done with a national grid, a provincial or state grid or a seven-parameter Helmert of much lower accuracy; WGS 84 to NAD83 can be done with a published Helmert, with a chain through ITRF, or with a 'null' transformation that treats them as identical. GIS software chooses among them by rules that are not always visible, and lesson 11 shows how to see and control the choice. The practitioner's rule is simple: choose the transformation whose accuracy and area of use suit the work, record its name and code, and verify it on a point known in both datums.