NAD83(CSRS) and epochs

Grant's crew observed pillar CR1 with GNSS for 24 hours in October 2026 and processed the data with Natural Resources Canada's CSRS-PPP service. The result agreed with the published NAD83(CSRS) coordinates to 3 mm in easting and 2 mm in northing, but only after the software was told which epoch to report. Asked for coordinates at the epoch of observation, it gave values 13 mm east, 17 mm south and 11 mm lower than the published ones. Nothing was wrong with the pillar or the survey: the published coordinates referred to epoch 2010.0, and the ground at Lethbridge had moved in sixteen years. This lesson explains NAD83(CSRS), its versions and epochs, and how coordinates are moved between epochs.

From NAD83 to NAD83(CSRS)

The original NAD83 of 1986 was a continental adjustment of terrestrial and Doppler observations, accurate at the metre level and distorted at the decimetre level. Canada built on it with the Canadian Spatial Reference System, realized through the Canadian Active Control System of continuously operating GNSS stations, the Canadian Base Network of high-accuracy monuments and provincial networks. NAD83(CSRS) is NAD83 redefined by a fourteen-parameter transformation from the International Terrestrial Reference Frame: it keeps NAD83's origin, orientation and scale as they were in 1986 and stays attached to the stable interior of the North American plate, while being as accurate and consistent as ITRF itself.

Because the transformation from ITRF is redefined each time a new ITRF is released, NAD83(CSRS) has versions. Each version is tied to a particular ITRF and has a reference epoch, the date to which its published coordinates refer. Natural Resources Canada adopted version 8 in 2024, based on ITRF2020 and with a reference epoch of 2010.0; earlier versions were tied to earlier ITRFs, and some provinces still publish their control in an older version or at another reference epoch, such as 1997.0 or 2002.0. The differences between versions at the same epoch are typically millimetres to a centimetre or two; the differences between epochs can be several centimetres. A coordinate stated as 'NAD83(CSRS)' without a version and epoch is therefore ambiguous at the centimetre level.

Why coordinates have epochs

The North American plate moves relative to the ITRF at about 2 cm per year, rotating about a pole near the equator off South America. NAD83(CSRS) removes that rotation, so stations in the plate's interior have nearly constant NAD83(CSRS) coordinates. 'Nearly' is important. Glacial isostatic adjustment, the slow rebound of the crust after the ice sheets melted, lifts the land around Hudson Bay by about a centimetre per year and causes slight subsidence farther out; the hinge between uplift and subsidence crosses the Prairies. Western British Columbia moves northeast by several millimetres per year as the Juan de Fuca plate pushes beneath it, and coastal regions shift suddenly in large earthquakes. Over decades, these motions add up to centimetres, enough to matter for precise surveys, deformation monitoring and the integration of data collected years apart.

The velocity grid and epoch propagation

To move coordinates between epochs, Natural Resources Canada publishes a velocity grid for NAD83(CSRS), derived from GNSS station velocities and a model of glacial isostatic adjustment. The current grid, used with version 8, gives east, north and up velocities at regular intervals across Canada. A coordinate at epoch \(t_1\) is moved to epoch \(t_2\) by adding the velocity times the time difference:

where \(M\) and \(N\) are the radii of curvature. NRCan's TRX tool performs these epoch transformations together with transformations between NAD83(CSRS) and ITRF realizations, using the velocity grid or user-supplied velocities, and converts the results to UTM, MTM and the provincial projections. CSRS-PPP reports NAD83(CSRS) coordinates at the epoch the user chooses, applying the same velocities. The rule for practice is to bring all coordinates in a project to one stated version and epoch, normally the one in which the provincial control is published, before combining them.

import math
a = 6378137.0; f = 1 / 298.257222101; e2 = f * (2 - f)
MAS = math.radians(1 / 3.6e6)                      # one milliarcsecond in radians

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)]
def enu(lat, lon, v):
    p, l = math.radians(lat), math.radians(lon)
    e = -math.sin(l) * v[0] + math.cos(l) * v[1]
    n = -math.sin(p) * math.cos(l) * v[0] - math.sin(p) * math.sin(l) * v[1] + math.cos(p) * v[2]
    u = math.cos(p) * math.cos(l) * v[0] + math.cos(p) * math.sin(l) * v[1] + math.sin(p) * v[2]
    return e, n, u

lat, lon, h = 49.69673058, -112.82708789, 893.412      # CR1 in NAD83(CSRS) at epoch 2010.0
X = to_xyz(lat, lon, h)

# 1. Why ITRF coordinates of CR1 keep changing: rotation of the North American plate.
#    Euler vector of approximately the published size (mas per year).
w = [0.03 * MAS, -0.69 * MAS, -0.02 * MAS]
v = [w[1] * X[2] - w[2] * X[1], w[2] * X[0] - w[0] * X[2], w[0] * X[1] - w[1] * X[0]]
ve, vn, vu = enu(lat, lon, v)
print(f"plate motion at CR1 in ITRF: east {ve*1000:+.1f} mm/yr, north {vn*1000:+.1f} mm/yr, "
      f"speed {math.hypot(ve, vn)*1000:.1f} mm/yr toward azimuth {math.degrees(math.atan2(ve, vn)) % 360:.0f} deg")
print(f"  over 2010.0 -> 2026.8 that is {math.hypot(ve, vn)*16.8:.3f} m")

# 2. Inside NAD83(CSRS) the plate motion is removed; what is left is small.
#    Velocities interpolated from a velocity grid (exercise values, mm/yr):
vg_e, vg_n, vg_u = 0.6, -0.9, -0.8
t_obs, t_ref = 2026.8, 2010.0
# A PPP result for CR1 observed at epoch 2026.8 (exercise: equal to the 2010.0 value
# moved forward by the grid velocities, plus 4 mm of noise)
M = a * (1 - e2) / (1 - e2 * math.sin(math.radians(lat))**2) ** 1.5
Nr = a / math.sqrt(1 - e2 * math.sin(math.radians(lat))**2)
def move(lat, lon, h, de, dn, du):
    return (lat + math.degrees(dn / (M + h)), lon + math.degrees(de / ((Nr + h) * math.cos(math.radians(lat)))), h + du)
obs = move(lat, lon, h, vg_e/1000*(t_obs-t_ref) + 0.003, vg_n/1000*(t_obs-t_ref) - 0.002, vg_u/1000*(t_obs-t_ref) + 0.002)
back = move(*obs, -vg_e/1000*(t_obs-t_ref), -vg_n/1000*(t_obs-t_ref), -vg_u/1000*(t_obs-t_ref))
for name, p in (("published, epoch 2010.0", (lat, lon, h)), ("PPP, epoch 2026.8", obs), ("PPP moved to 2010.0", back)):
    de = math.radians(p[1] - lon) * (Nr + h) * math.cos(math.radians(lat)) * 1000
    dn = math.radians(p[0] - lat) * (M + h) * 1000
    print(f"  {name:24s} {p[0]:.9f} {p[1]:.9f} {p[2]:9.4f}   dE {de:+6.1f} mm  dN {dn:+6.1f} mm  dU {(p[2]-h)*1000:+6.1f} mm")

plate motion at CR1 in ITRF: east -15.1 mm/yr, north -9.1 mm/yr, speed 17.6 mm/yr toward azimuth 239 deg
  over 2010.0 -> 2026.8 that is 0.296 m
  published, epoch 2010.0  49.696730580 -112.827087890  893.4120   dE   +0.0 mm  dN   +0.0 mm  dU   +0.0 mm
  PPP, epoch 2026.8        49.696730426 -112.827087709  893.4006   dE  +13.1 mm  dN  -17.1 mm  dU  -11.4 mm
  PPP moved to 2010.0      49.696730562 -112.827087848  893.4140   dE   +3.0 mm  dN   -2.0 mm  dU   +2.0 mm

The first part of the code computes the motion of CR1 in ITRF from an Euler vector of approximately the published size for the North American plate: about 15 mm per year west and 9 mm per year south, 17.6 mm per year toward an azimuth of 239°, or 0.30 m between 2010.0 and 2026.8. That is the motion NAD83(CSRS) removes. The second part uses exercise velocities of the kind the grid gives in southern Alberta, a few tenths of a millimetre to about a millimetre per year, to move a simulated PPP result at epoch 2026.8 back to 2010.0. At the observation epoch the PPP coordinates differ from the published ones by 13 mm, −17 mm and −11 mm; moved to 2010.0 they agree to 3 mm, −2 mm and +2 mm, the noise added in the simulation. Without the epoch correction a careful 24-hour survey would appear to disagree with control by 2 cm.

Projected coordinates and epochs

Projections are applied after the datum and epoch are settled: UTM, 3TM or a low-distortion projection on NAD83(CSRS) v8 at epoch 2010.0 is a fully defined system, while 'UTM on NAD83(CSRS)' is not, at the centimetre level. The EPSG registry has entries for NAD83(CSRS) versions, and newer software records the coordinate epoch alongside the coordinates of a data set. For most GIS and mapping the epoch can be ignored, because centimetres are below the data's accuracy. For control surveys, deformation monitoring, machine control and any work that combines GNSS observations made years apart, it cannot.

Looking ahead: NATRF2022 and Canada

The US National Geodetic Survey is replacing NAD83 with the North American Terrestrial Reference Frame of 2022 (NATRF2022), and Natural Resources Canada works closely with it. NATRF2022 is also plate-fixed, but it is aligned with ITRF2020's origin and scale rather than with the 1986 NAD83, so coordinates in it differ from NAD83(CSRS) by about a metre or more. NRCan has published a beta version of TRX for testing transformations to NATRF2022, and the timing and form of any Canadian adoption will be announced by NRCan and the provinces. In the meantime, surveyors near the border, like Coulee Geomatics, meet both systems and must label each carefully.

Writing a complete coordinate statement

The district's data standard requires every coordinate set delivered by a surveyor or contractor to carry a statement of its reference system. Leah's template had six lines, and she explained each to the district's staff with CR1 as the example. The statement looks pedantic until the first time two surveys of the same marks disagree by 2 cm and nobody can tell whether the marks moved, the epochs differ or one survey is wrong. With the statement, the question is settled in minutes.

Heights deserve special mention. Vertical velocities from glacial isostatic adjustment are larger than horizontal ones over much of Canada, reaching about a centimetre per year near Hudson Bay, and ellipsoidal heights therefore change with epoch more than positions do. Orthometric heights in CGVD2013 change too, because the land moves while the geoid stays nearly fixed, so a benchmark's published height is also tied to the epoch of its determination. In southern Alberta the vertical rates are small, a millimetre per year or less, and over a decade the effect is about the size of a careful levelling's error; for a monitoring survey of the canal's embankments, which looks for settlement of a few millimetres, it is not negligible and must be modelled.