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(CSRS), its imagery came with coordinates on a recent ITRF realization, and the consultant who surveyed the reservoir in 1995 worked on the original NAD83 adjustment. 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 metres 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 metres. 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 metres 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 metres to more than a hundred metres. The shift varies smoothly at the large scale but irregularly at the scale of metres, so no set of seven Helmert parameters can reproduce it to better than several metres 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 Canada the national NAD27-to-NAD83 grid is distributed in the NTv2 format (National Transformation version 2), designed by the Geodetic Survey of Canada in the 1990s and since adopted by many other countries. NTv2 files hold nested sub-grids, so that a coarse national grid can be refined where control is dense. Several provinces produced their own NTv2 grids to account for local distortions, and some, including Alberta, have used grids for the change from the original NAD83 adjustment to NAD83(CSRS). Accuracy depends on the grid and the region: of the order of tens of centimetres to a metre or two for NAD27 to NAD83 over most of the country, and centimetres for provincial grids between NAD83 realizations.

import math
# 1. A grid shift: interpolate latitude and longitude shifts (arc seconds) stored at the
#    nodes of a regular grid, as NTv2 grid files do. The node values below are
#    exercise values of realistic size for southern Alberta, not the published grid.
lat_nodes = [49.50, 49.75, 50.00]                 # degrees, grid spacing 0.25 deg
lon_nodes = [-113.00, -112.75, -112.50]
dlat = [[0.512, 0.509, 0.506],                    # rows: lat_nodes, columns: lon_nodes
        [0.516, 0.513, 0.510],
        [0.520, 0.517, 0.514]]
dlon = [[-3.951, -3.938, -3.925],
        [-3.958, -3.945, -3.932],
        [-3.966, -3.953, -3.940]]

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 = 49.69658830, -112.82598900          # 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(49.69673058, -112.82708789, 893.412)
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.5131" (+15.85 m), dlon = -3.9475" (-78.96 m)
NAD27 49.69658830, -112.82598900  ->  NAD83 49.69673082, -112.82708552
total horizontal shift 80.53 m

frame A -> frame B at three epochs (shift of the same point, m):
  epoch 2010.0: dX = +1.1975  dY = -2.6227  dZ = -1.0399  |d| = 3.0650
  epoch 2020.0: dX = +1.0174  dY = -2.6328  dZ = -1.1075  |d| = 3.0320
  epoch 2026.8: dX = +0.8949  dY = -2.6396  dZ = -1.1534  |d| = 3.0164

The interpolated shift at the old monument is +0.513\" in latitude and −3.948\" in longitude, or 15.9 m north and 79.0 m west: the old monument's NAD83 position is 80.5 m from its NAD27 coordinates. The grid values in the code are exercise values chosen to be of the right size for southern Alberta; for real work, the published grid is used through software. The monument's NAD83 coordinates come out within a few centimetres 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 centimetres 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 metres. 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.