Handling Vertical Datums and Geoid Separation in 3D Geospatial Data

This page shows how to reconcile the two height systems a digital twin mixes — ellipsoidal heights on WGS84 (EPSG:4979) and orthometric heights on a geoid model (NAVD88, EPSG:5703; or EGM2008, EPSG:3855) — by computing the geoid separation N explicitly with pyproj compound-CRS transforms, and how a datum mismatch shows up as a vertical step at tile seams. The core operation is a compound-CRS transform such as EPSG:32618+5703 → EPSG:4979 that carries a metric, orthometric point into the ellipsoidal frame Cesium renders in, applying the geoid grid per point rather than a constant.

You hit this whenever two datasets that “look” co-registered disagree in Z: a GNSS-derived cloud on ellipsoidal heights sits tens of metres above a DEM on NAVD88, or two adjacent tiles authored against different geoid models leave a cliff at their shared edge. The offset is the geoid separation, and it is not a bug in the geometry — it is an undeclared vertical datum. This guide isolates and quantifies that separation. The broader datum-management strategy lives in coordinate reference systems for 3D assets, and the horizontal EPSG:4326 → UTM path is covered in converting WGS84 to local projected coordinates; here the focus is strictly the vertical.

Geoid separation and a Z step at a tile seam Two adjacent terrain tiles share an edge; the left tile is referenced to the ellipsoid and the right to the geoid, so their surfaces meet at a vertical step equal to the geoid separation N, the gap between the ellipsoid and geoid surfaces. ellipsoid (EPSG:4979) geoid (NAVD88 / EGM2008) N Z step = N at the seam tile A: ellipsoidal Z tile B: orthometric Z
The geoid separation N is the gap between ellipsoid and geoid; author two tiles against different vertical datums and N reappears as a vertical step at their shared edge.

Prerequisites

  • Python 3.9+ with pyproj>=3.6 (PROJ 9.x) and numpy>=1.24: pip install "pyproj>=3.6" numpy.
  • PROJ vertical grids reachable at runtime — us_noaa_g2018u0.tif (GEOID18, behind NAVD88) and us_nga_egm2008_1.tif (EGM2008). Enable the PROJ CDN or install proj-data; without the grid, PROJ silently returns ellipsoidal heights and every separation reads as zero.
  • Data whose vertical datum you can state. The examples use a horizontal EPSG:32618 (UTM 18N) paired with orthometric NAVD88 (EPSG:5703) or EGM2008 (EPSG:3855), and the geographic-3D ellipsoidal frame EPSG:4979.
  • A published benchmark or a control point with a known orthometric height for validation.

Enable network grids first so a missing geoid fails loudly:

python
import pyproj
pyproj.network.set_network_enabled(active=True)
print("PROJ", pyproj.proj_version_str, "| network:", pyproj.network.is_network_enabled())

Step-by-Step

1. Compute the geoid separation N per point

N is the signed gap between ellipsoidal height h and orthometric height H, with H = h − N. Recover it directly: transform your orthometric compound CRS (EPSG:32618+5703) into the ellipsoidal frame (EPSG:4979) and subtract. N varies continuously across the extent, which is exactly why a constant will not do.

python
import numpy as np
from pyproj import Transformer

# Orthometric NAVD88 (EPSG:32618+5703) -> ellipsoidal WGS84 3D (EPSG:4979)
to_ellip = Transformer.from_crs("EPSG:32618+5703", "EPSG:4979", always_xy=True)

easting  = np.array([583120.4, 583402.9, 584010.7])   # metres, EPSG:32618
northing = np.array([4507880.1, 4508110.7, 4508640.2])
ortho_H  = np.array([14.62, 11.08, 9.55])              # NAVD88 orthometric, metres

lon, lat, ellip_h = to_ellip.transform(easting, northing, ortho_H)
separation = ellip_h - ortho_H                          # N, metres
print("geoid separation N (m):", np.round(separation, 3))

2. Convert orthometric heights to ellipsoidal for Cesium delivery

Web runtimes place geometry on the WGS84 ellipsoid (EPSG:4978 ECEF, reached via EPSG:4979). So the delivery boundary needs ellipsoidal Z, not the orthometric Z the twin analyses in. Do the conversion with the same compound transform, keeping the internal store orthometric.

python
from pyproj import Transformer

to_ecef = Transformer.from_crs("EPSG:4979", "EPSG:4978", always_xy=True)

# lon, lat, ellip_h came from step 1 — now push to geocentric ECEF for the tileset.
x, y, z = to_ecef.transform(lon, lat, ellip_h)
radii = np.linalg.norm(np.column_stack([x, y, z]), axis=1)
assert (6.3e6 < radii).all() and (radii < 6.5e6).all(), "not on the ellipsoid — datum chain wrong"
print("ECEF radius (m):", np.round(radii, 1))

3. Quantify the Z step two datums produce at a seam

If tile A is authored on EGM2008 and tile B on NAVD88, points that are physically identical carry different Z. Transform one benchmark through both compound CRSs into the shared ellipsoidal frame and the difference is the seam step you would see in the viewer.

python
from pyproj import Transformer

navd88  = Transformer.from_crs("EPSG:32618+5703", "EPSG:4979", always_xy=True)
egm2008 = Transformer.from_crs("EPSG:32618+3855", "EPSG:4979", always_xy=True)

e, n, h_ortho = 583120.4, 4507880.1, 14.62     # same ground point, one number per datum
_, _, ellip_from_navd88  = navd88.transform(e, n, h_ortho)
_, _, ellip_from_egm2008 = egm2008.transform(e, n, h_ortho)

seam_step = ellip_from_navd88 - ellip_from_egm2008
print(f"seam Z step from mixed vertical datums: {seam_step*100:.1f} cm")

4. Assert one vertical datum across a dataset before tiling

The durable fix is a gate: refuse any tile whose declared vertical sub-CRS differs from the project datum. Resolve the compound CRS and compare the vertical component by EPSG code.

python
from pyproj import CRS

PROJECT_VERTICAL = 5703        # NAVD88 for the whole twin

def assert_vertical_datum(crs_string: str) -> None:
    crs = CRS.from_user_input(crs_string)
    if not crs.is_compound:
        raise ValueError(f"{crs_string} has no vertical datum; Z is ambiguous")
    vert = crs.sub_crs_list[1]
    if vert.to_epsg() != PROJECT_VERTICAL:
        raise ValueError(f"vertical datum EPSG:{vert.to_epsg()} != project EPSG:{PROJECT_VERTICAL}")
    print(f"{crs_string}: vertical datum OK ({vert.name})")

assert_vertical_datum("EPSG:32618+5703")

Expected Output & Verification

Running steps 1–3 over the sample control points around New York, where the NAVD88 geoid separation is roughly −33 m, prints:

text
geoid separation N (m): [-33.412 -33.398 -33.371]
ECEF radius (m): [6368521.7 6368512.4 6368498.9]
seam Z step from mixed vertical datums: 27.4 cm

Verify against a benchmark whose published NAVD88 height you trust — the residual is your real vertical accuracy, and a residual near a round 30–100 m means the geoid grid never loaded:

python
known_ortho = 38.512                       # published NAVD88 benchmark height, m
_, _, ellip_bm = to_ellip.transform(583500.0, 4508900.0, known_ortho)

# The benchmark's ellipsoidal height must equal H + N (N ≈ -33.4 m here).
n_at_benchmark = ellip_bm - known_ortho
assert -34.5 < n_at_benchmark < -32.5, f"N={n_at_benchmark:.2f} m off — check the geoid grid"
print("benchmark separation N:", round(n_at_benchmark, 2), "m")

Two signatures matter. A separation that comes back exactly 0.0 means PROJ fell back to ellipsoidal heights because the grid was missing — the transform ran but did nothing. A seam step in step 3 on the order of 0.2–0.5 m is the real difference between GEOID18/NAVD88 and EGM2008 in the mid-latitudes, and it is precisely the cliff you would otherwise ship into a tileset.

Common Errors

Geoid separation reads exactly 0.0 for every point. PROJ could not find the vertical grid and silently returned ellipsoidal heights. Call pyproj.network.set_network_enabled(active=True) before building the transformer, pre-download the grid with projsync --file us_noaa_g2018u0.tif, and assert transformer.get_grids_used() reports available=True.

pyproj.exceptions.CRSError: Invalid projection: EPSG:32618+3855. The compound + syntax needs pyproj>=3.0 and PROJ 7+, and EGM2008 as a height CRS is EPSG:3855 (not a geographic code). Upgrade pyproj, or build the compound explicitly with CRS.from_epsg(32618) + CRS.from_epsg(3855).

Heights are off by a clean 33 m only in one region of a tiled twin. One batch of tiles was authored on ellipsoidal Z and the rest on orthometric, so the geoid separation appears as a step exactly where the batches meet. Gate every tile through the vertical-datum check in step 4 and re-transform the offending batch through the correct compound CRS rather than shifting it by a constant, which is wrong wherever N changes.

Back to Coordinate Reference Systems for 3D Assets.