Vertical Datums Explained: Why Your Elevations Are Off

Problem statement

Your GPS says you are at 106 m. The map says 55 m. Neither is wrong.

import pyproj
from pyproj import Transformer

pyproj.network.set_network_enabled(True)
to_orthometric = Transformer.from_crs("EPSG:4979", "EPSG:9518", always_xy=True)

for lon, lat, name in [(-2.2426, 53.4808, "Manchester"), (72.80, 19.00, "Mumbai")]:
    _, _, h = to_orthometric.transform(lon, lat, 0.0)
    print(f"{name:12} ellipsoidal 0 m -> orthometric {h:+.2f} m")

Manchester   ellipsoidal 0 m -> orthometric -51.63 m
Mumbai       ellipsoidal 0 m -> orthometric +67.75 m

A height of zero above the WGS84 ellipsoid is 51.6 m below sea level in Manchester and 67.8 m above it in Mumbai. The same number means two things 119 m apart.

This is the vertical datum problem. Horizontal CRS errors put things in the wrong place and are obvious on a map. Vertical datum errors shift everything by a smooth, plausible amount and are invisible.

Quick answer

There are two zeros in common use, and they are not the same surface:

The difference is the geoid separation N:

orthometric height (above sea level)  =  ellipsoidal height  −  N

N is not small and not constant:

PLACES = {"Snowdon": (-4.0763, 53.0685), "Manchester": (-2.2426, 53.4808),
          "Greenwich": (0.0, 51.4779), "Tokyo": (139.6917, 35.6895),
          "Cape Town": (18.4241, -33.9249), "Buenos Aires": (-58.3816, -34.6037),
          "New York": (-74.006, 40.7128), "Mumbai": (72.80, 19.00)}
for name, (lon, lat) in PLACES.items():
    _, _, h = to_orthometric.transform(lon, lat, 0.0)
    print(f"  {name:13} N = {-h:+7.2f} m")

  Snowdon       N =  +55.48 m
  Manchester    N =  +51.63 m
  Greenwich     N =  +45.90 m
  Tokyo         N =  +36.72 m
  Cape Town     N =  +31.13 m
  Buenos Aires  N =  +15.96 m
  New York      N =  -32.73 m
  Mumbai        N =  -67.75 m

A 123 m range across those eight places. If two datasets disagree by a suspiciously smooth few tens of metres, this is why.

Which datum, when it matters, converting with the real grid

1. Find out which datum each dataset uses

This is documentation work, not code — the information is usually not in the file:

Here is the Copernicus DEM of Snowdonia, straight from its public bucket:

import rasterio
from pyproj import CRS

with rasterio.open("dem.tif") as src:
    crs = CRS.from_user_input(src.crs)
    print("CRS:", src.crs, "| axes:", len(crs.axis_info), "| has a vertical part:",
          any(sub.is_vertical for sub in crs.sub_crs_list) or crs.is_vertical)
    print("tags:", src.tags())

CRS: EPSG:4326 | axes: 2 | has a vertical part: False
tags: {'AREA_OR_POINT': 'Point'}

EPSG:4326 is a 2D CRS. It says nothing at all about the vertical reference, and most DEMs ship exactly like this — the vertical datum lives in a PDF.

What the common products use:

2. Recognise when it matters

It matters whenever the absolute height is used:

• flood depth — water level minus ground level
• bridge or cable clearance
• "how far above sea level is this"
• combining GPS-surveyed points with a DEM

It does not matter when only differences are used:

• slope and aspect
• hillshade
• contour shapes (the labels shift, the lines do not)
• height above ground from DSM − DTM

The rule: a constant offset cancels in a subtraction. If your calculation subtracts two heights from the same dataset, you are safe. If it compares heights from two sources, or against an absolute figure, you are not.

3. Convert with pyproj, using the real grid

import pyproj
from pyproj import Transformer

pyproj.network.set_network_enabled(True)      # fetch the geoid grid from the PROJ CDN

to_orthometric = Transformer.from_crs("EPSG:4979", "EPSG:9518", always_xy=True)
lon, lat, ellipsoidal = -2.2426, 53.4808, 106.4

_, _, orthometric = to_orthometric.transform(lon, lat, ellipsoidal)
print(f"{ellipsoidal:.1f} m ellipsoidal -> {orthometric:.1f} m above sea level")

106.4 m ellipsoidal -> 54.8 m above sea level

EPSG:4979 is WGS84 3D (ellipsoidal height). EPSG:9518 is WGS84 with EGM2008 orthometric height. The transformation needs the EGM2008 grid file; set_network_enabled(True) downloads and caches it.

Without the grid, PROJ silently returns the input unchanged. That is the trap — a transform that appears to work and does nothing:

pyproj.network.set_network_enabled(False)
_, _, h = Transformer.from_crs("EPSG:4979", "EPSG:9518", always_xy=True).transform(lon, lat, 0.0)
print(f"without the grid: {h:+.2f} m")     # should be -51.63

without the grid: +0.00 m

4. Always verify against a known point

def check_datum(transformer, lon, lat, expected_n, tolerance=1.0):
    _, _, h = transformer.transform(lon, lat, 0.0)
    n = -h
    ok = abs(n - expected_n) < tolerance
    print(f"  N = {n:+.2f} m (expected {expected_n:+.2f}) — {'ok' if ok else 'GRID MISSING'}")
    return ok


check_datum(to_orthometric, -2.2426, 53.4808, expected_n=51.63)

  N = +51.63 m (expected +51.63) — ok

One assertion against a value you looked up once. It costs nothing and it catches the silent-no-op failure, which is otherwise indistinguishable from success.

5. Record the datum with the data

with rasterio.open("dem.tif") as src:
    profile, dem = src.profile, src.read(1)
with rasterio.open("dem_tagged.tif", "w", **profile) as dst:
    dst.write(dem, 1)
    dst.update_tags(vertical_datum="EGM2008 height", vertical_epsg="3855",
                    source="Copernicus DEM GLO-30 product handbook")
with rasterio.open("dem_tagged.tif") as src:
    print(src.tags())

{'source': 'Copernicus DEM GLO-30 product handbook', 'vertical_datum': 'EGM2008 height', 'vertical_epsg': '3855', 'AREA_OR_POINT': 'Area'}

A GeoTIFF carries the horizontal CRS automatically and the vertical datum not at all. Writing it into the tags is the only thing that stops the next person guessing — and tags only travel if you copy them: the source's AREA_OR_POINT=Point did not survive the profile copy, and the new file claims Area.

A converter that refuses to do nothing

The data used on this page

The heights on this page come from PROJ's own copy of the EGM2008 geoid grid, fetched from the PROJ CDN on first use. The DEM is the Copernicus GLO-30 model of Snowdonia and Caernarfon Bay, read from its public bucket; the window crosses the 4°W tile boundary, so two tiles are merged.

import numpy as np
import pyproj
import rasterio
from pyproj import Transformer
from rasterio.merge import merge

pyproj.network.set_network_enabled(True)
BBOX = (-4.40, 53.02, -4.00, 53.12)
COP = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{0}_00_DEM/"
       "Copernicus_DSM_COG_10_N53_00_{0}_00_DEM.tif")
tiles = [COP.format("W005"), COP.format("W004")]
mosaic, transform = merge(tiles, bounds=BBOX)
with rasterio.open(tiles[0]) as src:
    profile = src.profile | {"driver": "GTiff", "width": mosaic.shape[2], "height": mosaic.shape[1],
                             "transform": transform, "compress": "deflate"}
    tags = src.tags()
with rasterio.open("dem.tif", "w", **profile) as dst:
    dst.write(mosaic)
    dst.update_tags(**tags)
print(f"dem.tif {mosaic.shape[2]} × {mosaic.shape[1]}, {mosaic.min():.1f} to {mosaic.max():.1f} m")

dem.tif 960 × 360, -0.6 to 1074.5 m

Example 1 — a converter that cannot silently do nothing

import numpy as np
import pyproj
from pyproj import Transformer

pyproj.network.set_network_enabled(True)

KNOWN_SEPARATIONS = {          # lon, lat -> EGM2008 N, for verification
    "manchester": (-2.2426, 53.4808, 51.63),
    "new_york": (-74.006, 40.7128, -32.73),
    "mumbai": (72.80, 19.00, -67.75),
}


class DatumGridMissing(RuntimeError):
    pass


def ellipsoidal_to_orthometric(lon, lat, height, *, target="EPSG:9518", check="manchester"):
    """Convert ellipsoidal heights to orthometric, verifying the grid is actually present."""
    transformer = Transformer.from_crs("EPSG:4979", target, always_xy=True)

    check_lon, check_lat, expected_n = KNOWN_SEPARATIONS[check]
    _, _, probe = transformer.transform(check_lon, check_lat, 0.0)
    if abs(-probe - expected_n) > 1.0:
        raise DatumGridMissing(
            f"transform returned N={-probe:.2f} at {check}, expected {expected_n:+.2f}. "
            f"The geoid grid is not available — enable pyproj network access or "
            f"install proj-data."
        )

    lon, lat, height = np.atleast_1d(lon), np.atleast_1d(lat), np.atleast_1d(height)
    _, _, out = transformer.transform(lon, lat, height)
    print(f"  converted {len(out)} height(s); mean shift {np.mean(out - height):+.2f} m")
    return out


gps_heights = np.array([106.4, 98.2, 121.7])
lons = np.array([-2.2426, -2.2500, -2.2350])
lats = np.array([53.4808, 53.4770, 53.4850])

above_sea_level = ellipsoidal_to_orthometric(lons, lats, gps_heights)
print(np.round(above_sea_level, 1))

  converted 3 height(s); mean shift -51.63 m
[54.8 46.6 70.1]

The probe is the whole point. Without it, a missing grid produces heights that are 51.6 m too high and entirely plausible.

Example 2 — converting a whole DEM between vertical datums

The common real case runs the other way: the Copernicus DEM is in EGM2008 heights, and a GNSS survey, an ICESat-2 track or a photogrammetry model you want to compare it with is in ellipsoidal heights. So convert the DEM:

from rasterio.transform import xy
from scipy.ndimage import zoom


def shift_dem_datum(src_path, dst_path, *, source_crs="EPSG:9518", target_crs="EPSG:4979",
                    sample_step=32):
    """Apply a geoid separation surface to every cell of a DEM."""
    with rasterio.open(src_path) as src:
        dem = src.read(1).astype("float64")
        profile = src.profile.copy()
        transform, nodata = src.transform, src.nodata
    rows, cols = dem.shape
    transformer = Transformer.from_crs(source_crs, target_crs, always_xy=True)

    # N varies smoothly, so sample it coarsely and interpolate
    rr, cc = np.arange(0, rows, sample_step), np.arange(0, cols, sample_step)
    grid_r, grid_c = np.meshgrid(rr, cc, indexing="ij")
    xs, ys = xy(transform, grid_r.ravel(), grid_c.ravel())
    _, _, zeros = transformer.transform(np.array(xs), np.array(ys), np.zeros(len(xs)))
    offset = np.asarray(zeros).reshape(grid_r.shape)
    full = zoom(offset, (rows / offset.shape[0], cols / offset.shape[1]), order=1)[:rows, :cols]

    shifted = dem + full
    if nodata is not None:
        shifted = np.where(dem == nodata, nodata, shifted)
    profile.update(dtype="float32")
    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(shifted.astype("float32"), 1)
        dst.update_tags(vertical_datum="WGS 84 ellipsoidal height", vertical_epsg="4979",
                        converted_from=source_crs, mean_offset_m=f"{full.mean():.3f}")
    print(f"  offset ranged {full.min():+.2f} to {full.max():+.2f} m across the tile")
    return full, transformer


sampled, transformer = shift_dem_datum("dem.tif", "dem_ellipsoidal.tif")

with rasterio.open("dem.tif") as src:            # check the shortcut against every cell
    rows, cols = np.indices((src.height, src.width))
    xs, ys = xy(src.transform, rows.ravel(), cols.ravel())
_, _, exact = transformer.transform(np.array(xs), np.array(ys), np.zeros(rows.size))
error = np.abs(np.asarray(exact).reshape(rows.shape) - sampled)

print(f"  sampled every 32nd cell: max error {error.max() * 100:.1f} cm "
      f"against transforming all {rows.size:,} cells")

  offset ranged +55.16 to +55.52 m across the tile
  sampled every 32nd cell: max error 1.3 cm against transforming all 345,600 cells

Across this tile N varies by well under a metre — smooth enough that sampling every 32nd cell and interpolating is accurate to about a centimetre. On a tile this size transforming every cell is also quick; the sampling earns its keep on continental DEMs with billions of cells, where it cuts the PROJ calls by a factor of about a thousand.

Example 3 — the flood-depth error this causes

A still-water level of 4 m above sea level along the Caernarfon Bay coast — a round number for illustration, not a forecast. Compare it with the DEM in the right datum, and with the ellipsoidal copy from Example 2:

def flood_depth(dem, water_level_m, *, dem_datum, water_datum):
    if dem_datum != water_datum:
        raise ValueError(
            f"DEM is {dem_datum}, water level is {water_datum} — convert one first. "
            f"Ignoring this produces a depth error of tens of metres."
        )
    land = np.isfinite(dem)
    depth = np.where(land & (dem < water_level_m), water_level_m - dem, 0.0)
    print(f"  flooded {np.mean(depth[land] > 0):.2%} of land cells, max depth {depth.max():.2f} m")
    return depth


def read_land(path):
    with rasterio.open(path) as src, rasterio.open("dem.tif") as ref:
        z = src.read(1).astype("float64")
        return np.where(ref.read(1) == 0, np.nan, z)     # Copernicus: open sea is exactly 0

orthometric, ellipsoidal = read_land("dem.tif"), read_land("dem_ellipsoidal.tif")
correct = flood_depth(orthometric, 4.0, dem_datum="EGM2008", water_datum="EGM2008")
wrong = flood_depth(ellipsoidal, 4.0, dem_datum="EGM2008", water_datum="EGM2008")  # mislabelled
try:
    flood_depth(ellipsoidal, 4.0, dem_datum="WGS84 ellipsoid", water_datum="EGM2008")
except ValueError as err:
    print(f"  ValueError: {err}")

  flooded 3.30% of land cells, max depth 4.57 m
  flooded 0.00% of land cells, max depth 0.00 m
  ValueError: DEM is WGS84 ellipsoid, water level is EGM2008 — convert one first. Ignoring this produces a depth error of tens of metres.

A 55 m offset against a 4 m water level is not a small error — it is the difference between a flood and no flood. The raise in the function is doing more work than the arithmetic, and it only works if the datum labels are true.

Geoid, ellipsoid and national datums

Why the geoid is not a sphere or an ellipsoid

The geoid is an equipotential surface of Earth's gravity field: the shape a global ocean would take with no tides, currents or wind. Water finds the level where gravitational potential is constant, so the geoid bulges where mass is concentrated — over mountain ranges and dense mantle — and dips over ocean trenches and lighter crust.

The result is a lumpy surface deviating from the best-fit ellipsoid by roughly −105 m to +85 m worldwide. It has no formula. It is represented as a grid of separation values, which is why the conversion needs a data file rather than arithmetic, and why a missing file breaks it.

"Height above sea level" means height above the geoid, because that is what a spirit level and a tide gauge measure.

Why GNSS gives ellipsoidal heights

A GNSS receiver solves for its position from satellite ranges, and the natural output is Cartesian coordinates in an Earth-centred frame — converted to latitude, longitude and ellipsoidal height, because the ellipsoid is a formula the receiver can evaluate.

Getting orthometric height requires the geoid model, which is a large data file. Handheld receivers and phones usually apply one silently, often EGM96 rather than EGM2008. Survey-grade receivers give you the choice and record which they used.

So "my GPS says 106 m" is ambiguous until you know whether a geoid model was applied and which one.

Why national datums are a third layer

Most countries have their own vertical datum, defined by a historic tide gauge: Ordnance Datum Newlyn in Britain, NAVD88 in the United States, NAP in the Netherlands. These predate satellites and differ from EGM2008 by decimetres to a couple of metres — small compared with the geoid separation, large compared with survey tolerances.

PROJ can put numbers on both gaps. Here are EGM96 and Ordnance Datum Newlyn, each against EGM2008, at a few points:

egm2008 = Transformer.from_crs("EPSG:4979", "EPSG:9518", always_xy=True)
egm96 = Transformer.from_crs("EPSG:4979", "EPSG:9707", always_xy=True)     # WGS 84 + EGM96 height
odn = Transformer.from_crs("EPSG:4979", "EPSG:7405", always_xy=True)       # OSGB36 + ODN height (OSGM15)
for name, (lon, lat) in {"Snowdon": (-4.0763, 53.0685), "Manchester": (-2.2426, 53.4808),
                         "Newlyn": (-5.54, 50.10), "Edinburgh": (-3.19, 55.95),
                         "Mumbai": (72.80, 19.00), "New York": (-74.006, 40.7128)}.items():
    h08 = egm2008.transform(lon, lat, 0.0)[2]
    line = f"  {name:11} EGM96 − EGM2008 {egm96.transform(lon, lat, 0.0)[2] - h08:+.2f} m"
    if lat > 49.5 and lon < 2:
        line += f"   ODN − EGM2008 {odn.transform(lon, lat, 0.0)[2] - h08:+.2f} m"
    print(line)

  Snowdon     EGM96 − EGM2008 +0.90 m   ODN − EGM2008 +0.59 m
  Manchester  EGM96 − EGM2008 +0.23 m   ODN − EGM2008 +0.62 m
  Newlyn      EGM96 − EGM2008 +0.07 m   ODN − EGM2008 +0.57 m
  Edinburgh   EGM96 − EGM2008 +0.22 m   ODN − EGM2008 +0.62 m
  Mumbai      EGM96 − EGM2008 +0.78 m
  New York    EGM96 − EGM2008 +0.04 m

So there are three surfaces in play: the ellipsoid, the global geoid, and the national datum. Converting between the first two needs a global grid; between the second and third needs a national one, which PROJ can also fetch when network access is enabled.

For most GIS work, EGM2008 is close enough to any national datum. For engineering, it is not.

Why the silent no-op is the real danger

If PROJ cannot find the geoid grid, it does not raise. It performs the transformation it can — a null vertical shift — and returns the input height unchanged.

That behaviour is defensible (a horizontal transform still succeeded) and catastrophic in practice, because the failure looks exactly like success. Your heights are wrong by the local geoid separation, smoothly, with no error anywhere.

The probe in Example 1 is the defence: transform one point whose separation you know, and refuse to proceed if the answer is zero.

2D CRSs, missing grids and EGM96 against EGM2008

• EPSG:4326 is 2D. It carries no vertical datum. Use EPSG:4979 when you mean WGS84 with ellipsoidal height.
• A missing grid returns the input unchanged. Always verify against a known separation.
• pyproj.network.set_network_enabled(True) fetches grids from the PROJ CDN and caches them. In an offline environment, install proj-data instead.
• EGM96 and EGM2008 differ, by a few centimetres to about a metre at the places in the explanation's comparison, and by more in parts of the world with sparse gravity data. SRTM and ASTER GDEM use EGM96; Copernicus uses EGM2008.
• Phone GPS altitude is unreliable regardless of datum — vertical GNSS accuracy is typically two to three times worse than horizontal.
• Slope, aspect, hillshade and contour geometry are unaffected by a constant vertical offset. Only contour labels change.
• Bathymetry uses yet another zero, often lowest astronomical tide, so depths and elevations do not simply concatenate across a coastline.
• Record the datum in file tags. GeoTIFF carries the horizontal CRS automatically and the vertical datum not at all.

Internal links

• Digital elevation models explained: DEM, DSM and DTM — which surface, and which datum, your DEM uses
• Projected vs geographic CRS explained — the horizontal counterpart
• Datum shifts and transformations in Python — the horizontal version of this problem
• How to calculate slope and aspect from a DEM in Python — unaffected by vertical datum
• How to extract an elevation profile along a line in Python — where absolute heights start to matter
• The raster data model explained — what a GeoTIFF does and does not record
• How to set up a Python GIS environment that actually works — installing PROJ grids
• Coordinate precision explained — how many decimals a height deserves
• Extruded buildings float above or sink into the ground — the geoid separation in practice

FAQ

Why does my GPS altitude differ from the map?

Your receiver is probably reporting height above the WGS84 ellipsoid while the map reports height above the geoid. The difference is the geoid separation — between about +46 m and +56 m across Britain, −68 m in Mumbai.

What is the geoid separation?

The vertical distance between the ellipsoid and the geoid at a location, usually written N. It ranges from about −105 m to +85 m worldwide and varies smoothly.

Does the vertical datum affect slope or hillshade?

No. Those depend on differences between neighbouring cells, and a smooth offset cancels. It affects absolute heights: flood depth, clearance, height above sea level.

How do I convert in Python?

pyproj.Transformer.from_crs("EPSG:4979", "EPSG:9518") with network access enabled so the EGM2008 grid can be downloaded. Verify against a known separation afterwards.

Why did my conversion do nothing?

The geoid grid was not available, and PROJ returns the input unchanged rather than raising. Probe one point with a known separation before trusting any result.

Which vertical datum does my DEM use?

Check the product documentation — it is not in the file. Copernicus DEM uses EGM2008; SRTM and ASTER GDEM use EGM96; national products use national datums.

Is EGM2008 the same as my national datum?

Close but not identical — typically within a metre. Fine for most GIS work, not for engineering or survey.