State plane coordinate systems

In the United States, each state has its own official plane coordinate system, defined in state law and designed by the National Geodetic Survey (NGS). The State Plane Coordinate System (SPCS) has been the framework for engineering, surveying and land records in most states for nearly a century, and it is about to change for the first time since 1983. This lesson describes how SPCS was designed, the zones in Montana and its neighbors, the units problem that has caused more errors than any projection formula, and the State Plane Coordinate System of 2022 (SPCS2022) that will accompany the modernized National Spatial Reference System.

Origins: SPCS 27

The system began in 1933, when the US Coast and Geodetic Survey designed a plane coordinate system for North Carolina at the request of a highway engineer who wanted surveys tied to the national control network without the burden of geodetic computation. The idea spread quickly, and by the late 1930s every state had a system. The goal was a maximum scale error of 1 part in 10,000 (100 ppm), small enough that surveyors could use plane trigonometry with a single correction. To achieve it, states long in the north-south direction were given transverse Mercator zones and states long in the east-west direction Lambert conformal conic zones, and larger states were divided into several zones. These SPCS 27 zones were defined on NAD27 in US survey feet, and many old plats, deeds and highway plans still carry their coordinates.

SPCS 83

With the adoption of NAD83, NGS redefined the zones as SPCS 83, with defining parameters in meters. Most zones kept their SPCS 27 geometry with new false origins, so that SPCS 83 coordinates are clearly distinguishable from SPCS 27 ones, and some states took the opportunity to change their zones. Montana, which had three Lambert zones under SPCS 27 (North, Central and South), adopted a single statewide Lambert zone, identified as zone 2500, with standard parallels at 45° N and 49° N, an origin at 44°15' N and 109°30' W, a false easting of 600,000 m and a false northing of zero. A single zone simplified statewide data and land records at the cost of distortion up to about 600 ppm in the middle of the state, larger than the 1933 design goal. Idaho kept three transverse Mercator zones and Wyoming four; North Dakota and South Dakota kept two Lambert zones each.

The foot problem

SPCS 83 was defined in meters, but most users work in feet, and there were two feet. The US survey foot, defined in 1893 as 1200/3937 m, was used for geodetic surveying; the international foot of exactly 0.3048 m, adopted in 1959 for all other purposes, is shorter by 2 parts per million. Each state chose by law which foot to use with SPCS 83, and a few, Montana among them, chose the international foot. Two ppm sounds negligible, but state plane coordinates are millions of feet: at 1.5 million feet the difference is 3 ft, enough to put a property corner in the neighbor's yard or a culvert in the wrong place. Software that converts a Montana file using the wrong foot moves every point by several feet in a direction away from the false origin. To end the confusion, NIST and NGS deprecated the US survey foot at the end of 2022, and SPCS2022 uses the international foot wherever feet are used.

import math
a = 6378137.0; f = 1 / 298.257222101; e = math.sqrt(f * (2 - f))   # GRS80

def _m(p): return math.cos(p) / math.sqrt(1 - (e * math.sin(p)) ** 2)
def _t(p): return math.tan(math.pi / 4 - p / 2) / ((1 - e * math.sin(p)) / (1 + e * math.sin(p))) ** (e / 2)

class LCC:
    """Lambert conformal conic with two standard parallels (degrees, m)."""
    def __init__(self, lat1, lat2, lat0, lon0, fe=0.0, fn=0.0):
        p1, p2, p0 = map(math.radians, (lat1, lat2, lat0))
        self.n = (math.log(_m(p1)) - math.log(_m(p2))) / (math.log(_t(p1)) - math.log(_t(p2)))
        self.F = _m(p1) / (self.n * _t(p1) ** self.n)
        self.r0 = a * self.F * _t(p0) ** self.n
        self.lon0, self.fe, self.fn = lon0, fe, fn
    def forward(self, lat, lon):
        p = math.radians(lat); r = a * self.F * _t(p) ** self.n
        th = self.n * math.radians(lon - self.lon0)
        E = self.fe + r * math.sin(th); N = self.fn + self.r0 - r * math.cos(th)
        k = r * self.n / (a * _m(p))                  # point scale factor
        return E, N, k, math.degrees(th)               # th is the grid convergence
    def inverse(self, E, N):
        dx, dy = E - self.fe, self.r0 - (N - self.fn)
        r = math.copysign(math.hypot(dx, dy), self.n); th = math.atan2(dx, dy)
        t = (r / (a * self.F)) ** (1 / self.n); p = math.pi / 2 - 2 * math.atan(t)
        for _ in range(30):
            pn = math.pi / 2 - 2 * math.atan(t * ((1 - e * math.sin(p)) / (1 + e * math.sin(p))) ** (e / 2))
            if abs(pn - p) < 1e-15: break
            p = pn
        return math.degrees(p), self.lon0 + math.degrees(th / self.n)
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}\""

US_FT = 1200 / 3937          # the US survey foot, in m (retired after 2022)
INT_FT = 0.3048              # the international foot, in m

lat, lon = 47.50600422, -111.29779092          # pillar CR1, Great Falls
mt = LCC(49, 45, 44.25, -109.5, 600000, 0)     # SPCS 83 Montana (2500), a Lambert zone
E, N, k, g = mt.forward(lat, lon)
print(f"SPCS 83 Montana:  E = {E:.3f} m  N = {N:.3f} m  k = {k:.7f}  convergence {g:+.4f} deg")
print(f"  international feet: E = {E/INT_FT:.3f}  N = {N/INT_FT:.3f}")
print(f"  US survey feet:     E = {E/US_FT:.3f}  N = {N/US_FT:.3f}")
print(f"  using the wrong foot moves the point {(E/INT_FT - E/US_FT):.3f} ft east and "
      f"{(N/INT_FT - N/US_FT):.3f} ft north (2 parts per million of the coordinate)")

# A transverse Mercator state plane zone for comparison: Idaho Central (SPCS 83)
iE, iN, ik, ig = tm_forward(lat, lon, -114.0, 1 - 1/19000, 500000.0, 0.0)
print(f"\nThe same point on a TM zone with k0 = 1 - 1/19000 and central meridian 114 W:")
print(f"  E = {iE:.3f} m  k = {ik:.7f}  (an example of a zone used outside its area)")

print("\nMontana zone scale at the latitude of CR1 and in the district:")
for L in (47.3, 47.4, 47.5, 47.6, 47.7):
    print(f"  latitude {L}: k = {mt.forward(L, -111.3)[2]:.7f}  ({(mt.forward(L, -111.3)[2]-1)*1e6:+7.1f} ppm)")

SPCS 83 Montana:  E = 464652.589 m  N = 363352.692 m  k = 0.9994297  convergence -1.3151 deg
  international feet: E = 1524450.750  N = 1192102.008
  US survey feet:     E = 1524447.701  N = 1192099.624
  using the wrong foot moves the point 3.049 ft east and 2.384 ft north (2 parts per million of the coordinate)

The same point on a TM zone with k0 = 1 - 1/19000 and central meridian 114 W:
  E = 703553.631 m  k = 1.0004564  (an example of a zone used outside its area)

Montana zone scale at the latitude of CR1 and in the district:
  latitude 47.3: k = 0.9994052  ( -594.8 ppm)
  latitude 47.4: k = 0.9994155  ( -584.5 ppm)
  latitude 47.5: k = 0.9994288  ( -571.2 ppm)
  latitude 47.6: k = 0.9994452  ( -554.8 ppm)
  latitude 47.7: k = 0.9994646  ( -535.4 ppm)

CR1 has Montana state plane coordinates of E 464,652.589 m and N 363,352.692 m, or 1,524,450.750 ft and 1,192,102.008 ft in international feet. The same meters expressed in US survey feet are 3.05 ft smaller in easting and 2.38 ft smaller in northing. The point scale factor at CR1 is 0.9994297, 570 ppm short, and across the district the scale runs from −595 to −535 ppm, so a 1,000 ft grid distance corresponds to about 1,000.57 ft on the ellipsoid before height is considered. Projected on a transverse Mercator zone meant for central Idaho, the same point has a scale factor of 1.00046, an example of how quickly distortion grows when a zone is used outside its area.

SPCS2022

NGS is replacing NAD83 and NAVD 88 with the modernized National Spatial Reference System, including the North American Terrestrial Reference Frame of 2022 (NATRF2022) described in lesson 10, and SPCS2022 is the plane coordinate system that goes with it. Its design differs from SPCS 83 in several ways. Each state may have a statewide zone, for those who want one system for the whole state, and in many states a layer of smaller low-distortion zones designed with local surveyors, engineers and GIS users. Zones are designed to minimize distortion at the topographic surface rather than on the ellipsoid, so that grid distances in populated areas are close to ground distances, as lesson 8 explains. All zones are defined in meters, and feet, where used, are international feet. Zone parameters were developed with the states and published by NGS on its development websites; SPCS2022 becomes official with the rest of the modernized NSRS, after a beta period in which NGS released the new tools for public testing.

For a practitioner in Montana, this means that coordinates in the coming years will exist in at least three forms: SPCS 27 on NAD27 in the oldest records, SPCS 83 on one of several NAD83 realizations, and SPCS2022 on NATRF2022. The numbers in each will be different, and a coordinate without its system, datum, realization, epoch and unit is not usable for precise work. The district's data standard therefore records all of these for every coordinate set.

Choosing a system for the district

Leah's recommendation followed the uses. The GIS master data remained on SPCS 83 Montana in international feet on NAD83(2011), shared with the county's parcel data and the state's cadastral layer. Survey control and legal surveys would also be reported on SPCS 83, with grid distances converted to ground distances by the combined factor of lesson 8, until the state adopts SPCS2022. Engineering design for the canal rehabilitation, where contractors stake from coordinates with total stations and expect grid distances to agree with measured distances, would use a low-distortion projection designed for the district, which lesson 8 also describes, with a documented transformation from state plane. Federal data arriving in UTM zone 12 would be transformed on load. Each choice was written into the district's data standard with its EPSG code where one exists.