Transverse Mercator and UTM

Most of the coordinates that pass through Coulee Geomatics' office are transverse Mercator coordinates: UTM from GNSS receivers, provincial orthophotos and the district's GIS, Alberta's 3TM from the provincial survey control database, and 10TM from the forestry and resource layers. This lesson develops the projection that underlies all of them. The transverse Mercator (TM) is the conformal cylinder of lesson 3, turned so that it touches the ellipsoid along a meridian. Distortion depends almost entirely on distance east or west of that central meridian, which makes it ideal for regions long in the north-south direction, and its formulas, once a matter of laborious series, are now accurate to fractions of a millimetre anywhere they are used.

From Lambert and Gauss to Krüger

Lambert described the transverse Mercator for the sphere in 1772. Gauss developed the ellipsoidal version in the 1820s for the survey of Hanover, and Louis Krüger published practical series for it in 1912, which is why it is called the Gauss-Krüger projection in much of Europe. For most of the twentieth century surveyors computed TM coordinates with power series in the longitude difference, of which the best known in the English-speaking world are Thomas's and Redfearn's of the 1940s and 1950s. Those series are accurate to a millimetre within about 3° to 4° of the central meridian and degrade quickly beyond it, which was one reason for limiting zones to a few degrees.

Krüger's other series, written in powers of the third flattening \(n = f/(2 - f)\) instead of the longitude difference, was rediscovered and extended to higher order in the early 2000s, notably by Charles Karney, who showed in 2011 that a sixth-order version is accurate to about 5 nanometres within 3,900 km of the central meridian. Modern software, including the PROJ library behind most GIS programs, uses this form. It works through the conformal latitude and two intermediate coordinates \(\xi'\) and \(\eta'\) of a transverse spherical Mercator, followed by trigonometric series that convert them to the ellipsoid:

\(A\) is the rectifying radius, which turns the meridian arc into a length, and the coefficients \(\alpha_j\) are polynomials in \(n\). The inverse uses a second set of coefficients \(\beta_j\) and a short iteration for latitude. The code in this lesson implements both directions, together with the point scale factor and the grid convergence, in about sixty lines of plain Python. It was checked against an independent Redfearn-series computation and agrees to better than a millimetre at the district.

Scale and convergence on a TM zone

Near the central meridian the scale factor of a TM projection with central scale \(k_0\) is approximately

where \(E'\) is the distance east or west of the central meridian and \(R_m\) the mean radius of curvature. Scale grows with the square of the distance from the central meridian and is almost independent of latitude, so lines of equal scale run roughly parallel to it. With \(k_0 = 1\) the projection is tangent and all scales are greater than one. Reducing \(k_0\) below one makes it secant: two lines of exact scale appear on either side of the central meridian, scales between them are less than one and scales outside are greater than one.

The grid convergence \(\gamma\), the angle between grid north and the meridian at a point, is approximately

where \(\Delta\lambda\) is the longitude difference from the central meridian. West of the central meridian \(\gamma\) is negative in the convention used here: true north points east of grid north, by about 1.4° at CR1. A geodetic azimuth \(\alpha\) measured from true north becomes a grid azimuth \(\alpha - \gamma\), plus the small arc-to-chord correction of lesson 2 for long lines. Confusing the sign of the convergence is one of the commonest errors in survey computations, and the safest check is a sketch showing the meridian and the grid line at the point.

The Universal Transverse Mercator system

UTM, adopted by the United States Army in the late 1940s and later by most of the world's mapping agencies, divides the Earth between 80° S and 84° N into 60 zones, each 6° of longitude wide, numbered eastward from zone 1 at 180° to 174° W. Each zone is a transverse Mercator with these parameters:

The scale factor of 0.9996 on the central meridian was chosen so that the scale is exactly one along two lines about 180 km east and west of it, and reaches about 1.0010 at the zone boundary on the equator, 333 km from the central meridian. Farther north the boundaries are closer to the central meridian, so the edge scale is smaller: at 50° N the boundary is 215 km out and its scale about 1.0002. Over the whole of a zone the scale stays within about −400 to +1,000 ppm, roughly 1 part in 2,500 at worst, adequate for topographic mapping at 1:50,000 and for GIS, not for construction. There are exceptions to the standard zone pattern in south-western Norway and around Svalbard, a reminder that zone numbers should never be computed blindly from longitude where such exceptions apply.

import math
A_GRS80 = 6378137.0
F_GRS80 = 1 / 298.257222101

def tm_setup(a=A_GRS80, f=F_GRS80):
    n = f / (2 - f)
    A = a / (1 + n) * (1 + n**2 / 4 + n**4 / 64 + n**6 / 256)
    al = [n/2 - 2*n**2/3 + 5*n**3/16 + 41*n**4/180 - 127*n**5/288 + 7891*n**6/37800,
          13*n**2/48 - 3*n**3/5 + 557*n**4/1440 + 281*n**5/630 - 1983433*n**6/1935360,
          61*n**3/240 - 103*n**4/140 + 15061*n**5/26880 + 167603*n**6/181440,
          49561*n**4/161280 - 179*n**5/168 + 6601661*n**6/7257600,
          34729*n**5/80640 - 3418889*n**6/1995840,
          212378941*n**6/319334400]
    be = [n/2 - 2*n**2/3 + 37*n**3/96 - n**4/360 - 81*n**5/512 + 96199*n**6/604800,
          n**2/48 + n**3/15 - 437*n**4/1440 + 46*n**5/105 - 1118711*n**6/3870720,
          17*n**3/480 - 37*n**4/840 - 209*n**5/4480 + 5569*n**6/90720,
          4397*n**4/161280 - 11*n**5/504 - 830251*n**6/7257600,
          4583*n**5/161280 - 108847*n**6/3991680,
          20648693*n**6/638668800]
    e = math.sqrt(f * (2 - f))
    return n, A, al, be, e

def tm_forward(lat, lon, lon0, k0, fe=0.0, fn=0.0, a=A_GRS80, f=F_GRS80):
    """Geodetic latitude and longitude (degrees) to TM easting, northing,
    point scale factor k and grid convergence gamma (degrees)."""
    n, A, al, be, e = tm_setup(a, f)
    phi = math.radians(lat); dl = math.radians(lon - lon0)
    t = math.sinh(math.atanh(math.sin(phi)) - e * math.atanh(e * math.sin(phi)))
    xi1 = math.atan2(t, math.cos(dl)); eta1 = math.atanh(math.sin(dl) / math.sqrt(1 + t * t))
    xi, eta = xi1, eta1; p, q = 1.0, 0.0
    for j in range(1, 7):
        xi += al[j-1] * math.sin(2*j*xi1) * math.cosh(2*j*eta1)
        eta += al[j-1] * math.cos(2*j*xi1) * math.sinh(2*j*eta1)
        p += 2*j*al[j-1] * math.cos(2*j*xi1) * math.cosh(2*j*eta1)
        q += 2*j*al[j-1] * math.sin(2*j*xi1) * math.sinh(2*j*eta1)
    E = fe + k0 * A * eta; N = fn + k0 * A * xi
    k = k0 * A / a * math.sqrt(1 + ((1 - n) / (1 + n) * math.tan(phi))**2) * math.sqrt(p*p + q*q) / math.sqrt(t*t + math.cos(dl)**2)
    g1 = math.atan2(t * math.tan(dl), math.sqrt(1 + t*t))
    g2 = math.atan2(q, p)
    return E, N, k, math.degrees(g1 + g2)

def tm_inverse(E, N, lon0, k0, fe=0.0, fn=0.0, a=A_GRS80, f=F_GRS80):
    """TM easting and northing back to latitude and longitude (degrees)."""
    n, A, al, be, e = tm_setup(a, f)
    xi = (N - fn) / (k0 * A); eta = (E - fe) / (k0 * A)
    xi1, eta1 = xi, eta
    for j in range(1, 7):
        xi1 -= be[j-1] * math.sin(2*j*xi) * math.cosh(2*j*eta)
        eta1 -= be[j-1] * math.cos(2*j*xi) * math.sinh(2*j*eta)
    chi = math.asin(math.sin(xi1) / math.cosh(eta1))
    lam = math.atan2(math.sinh(eta1), math.cos(xi1))
    t0 = math.tan(chi); t = t0
    for _ in range(10):
        sig = math.sinh(e * math.atanh(e * t / math.sqrt(1 + t*t)))
        tp = t * math.sqrt(1 + sig*sig) - sig * math.sqrt(1 + t*t)
        dt = (t0 - tp) / math.sqrt(1 + tp*tp) * (1 + (1 - e*e) * t*t) / ((1 - e*e) * math.sqrt(1 + t*t))
        t += dt
        if abs(dt) < 1e-14: break
    return math.degrees(math.atan(t)), lon0 + math.degrees(lam)

def dms(x):
    s = '-' if x < 0 else ''; x = abs(x); d = int(x); m = int((x - d) * 60); sec = ((x - d) * 60 - m) * 60
    if sec >= 59.99995: sec -= 60; m += 1
    if m == 60: m = 0; d += 1
    return f"{s}{d}°{m:02d}'{sec:07.4f}\""

def utm_zone(lon): return int((lon + 180) // 6) + 1
def utm(lat, lon, zone=None):
    z = zone or utm_zone(lon); lon0 = -183 + 6 * z
    E, N, k, g = tm_forward(lat, lon, lon0, 0.9996, 500000.0, 0.0 if lat >= 0 else 10000000.0)
    return z, E, N, k, g

lat, lon = 49.69673058, -112.82708789          # pillar CR1
z, E, N, k, g = utm(lat, lon)
print(f"CR1  {dms(lat)} N  {dms(-lon)} W")
print(f"UTM zone {z}N (central meridian {-183 + 6*z} deg)")
print(f"  E = {E:12.3f} m   N = {N:12.3f} m")
print(f"  point scale factor k = {k:.8f}   grid convergence = {g:+.6f} deg = {dms(g)}")
la, lo = tm_inverse(E, N, -111, 0.9996, 500000.0)
print(f"  inverse: {dms(la)} N  {dms(-lo)} W   (error {abs(la-lat)*111.2e6:.3f} mm, {abs(lo-lon)*72e6:.3f} mm)")

print("\nScale factor across UTM zone 12 at the latitude of CR1")
for dlon in (0, 1, 2, 3, 4):
    _, E1, _, k1, _ = utm(lat, -111 - dlon, 12)
    print(f"  {dlon} deg from the central meridian: E = {E1:10.1f}  k = {k1:.6f}  ({(k1-1)*1e6:+7.1f} ppm)")

# The same point in the neighbouring zone 11 (overlap use near a zone boundary)
z2, E2, N2, k2, g2 = utm(lat, lon, 11)
print(f"\nCR1 in zone 11 (central meridian -117): E = {E2:.3f}  N = {N2:.3f}  k = {k2:.6f}")

CR1  49°41'48.2301" N  112°49'37.5164" W
UTM zone 12N (central meridian -111 deg)
  E =   368239.741 m   N =  5506515.119 m
  point scale factor k = 0.99981324   grid convergence = -1.393594 deg = -1°23'36.9381"
  inverse: 49°41'48.2301" N  112°49'37.5164" W   (error 0.000 mm, 0.000 mm)

Scale factor across UTM zone 12 at the latitude of CR1
  0 deg from the central meridian: E =   500000.0  k = 0.999600  ( -400.0 ppm)
  1 deg from the central meridian: E =   427883.7  k = 0.999664  ( -336.1 ppm)
  2 deg from the central meridian: E =   355771.0  k = 0.999856  ( -144.5 ppm)
  3 deg from the central meridian: E =   283665.4  k = 1.000175  ( +174.9 ppm)
  4 deg from the central meridian: E =   211570.6  k = 1.000622  ( +622.0 ppm)

CR1 in zone 11 (central meridian -117): E = 800894.074  N = 5513275.846  k = 1.000712

CR1 lies in zone 12 at E 368,239.741 m, N 5,506,515.119 m. It is 1.83° of longitude west of the central meridian, so its easting is 131.8 km less than the false easting. Its point scale factor is 0.99981324, which means a short line at CR1 is about 187 ppm shorter on the grid than on the ellipsoid, and its grid convergence is −1.3936°, or −1°23'36.9\". The inverse returns the original latitude and longitude to well under a millimetre. Across the zone, at the latitude of CR1, scale runs from −400 ppm on the central meridian to −145 ppm 2° away and +175 ppm 3° away, at the zone edge. Projected into zone 11 instead, CR1 has an easting of 800,894 m and a scale factor of 1.000712, which shows why a point should be expressed in its own zone except where a project that straddles a boundary needs one consistent grid.

Zone boundaries, overlaps and grid references

A feature that crosses a zone boundary has coordinates in two zones that cannot be combined directly: the same point has completely different eastings in the two. Mapping agencies handle this by printing an overlap strip, by extending zones a little past their boundaries for projects that straddle them, or by choosing a single zone for a whole project even if part of it lies outside. Extending a TM zone by a degree or so costs a few hundred ppm of extra scale error, which may be acceptable for mapping but should be checked for engineering. Software that converts between zones must use the full projection formulas, inverse then forward; adding a constant shift to the coordinates is wrong.

For grid references, the Military Grid Reference System (MGRS) and the US National Grid divide each UTM zone into latitude bands and 100 km squares labelled by two letters, so that a point can be given as, for example, '12U' followed by the square letters and a short easting and northing within the square. CR1 lies in grid zone 12U. Canada's topographic maps of the National Topographic System show UTM grid lines, as do the provincial orthophoto mosaics, and searchers, emergency services and the military read positions from them. Precision is set by the number of digits: four digits in each of the easting and northing within the square give a position to 10 m, five digits to 1 m.

Choosing k₀ and the width of a TM zone

UTM is one point in a family of choices. Narrowing the zone and raising \(k_0\) toward one reduces distortion. A 3° zone with \(k_0 = 0.9999\) has scale within about −100 to +100 ppm at middle latitudes, the choice made by Alberta's 3TM and by Canada's modified transverse Mercator (MTM) zones of lesson 7. Widening the zone and lowering \(k_0\) covers more ground with one system at the cost of larger errors: Alberta's 10TM, with \(k_0 = 0.9992\), covers the whole province in one zone. State plane TM zones in the United States are mostly about 2° to 3° wide with \(k_0\) between 0.9999 and 0.99996. And a low-distortion projection for a single project may use a TM with \(k_0\) slightly greater than one, chosen to match the height of the ground, which lesson 8 develops.