Vertical Datums Explained: Why Your Elevations Are Off #

Problem statement #

Your GPS says you are at 106 m. The map says 51 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:

  Snowdon         +55.48 m
  Manchester      +51.63 m
  Greenwich       +45.90 m
  Tokyo           +36.66 m
  Cape Town       +31.13 m
  Buenos Aires    +15.95 m
  New York        −32.73 m
  Mumbai          −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.

Step-by-step solution #

1. Find out which datum each dataset uses #

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

import rasterio

with rasterio.open("dem.tif") as src:
    print("CRS:", src.crs)
    print("is 3D:", src.crs.is_vertical if hasattr(src.crs, "is_vertical") else "n/a")
    print("tags:", src.tags())

CRS: EPSG:4326
is 3D: False
tags: {}

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_orthometric.tif", "w", **profile) as dst:
    dst.write(dem, 1)
    dst.update_tags(
        vertical_datum="EGM2008",
        vertical_epsg="9518",
        converted_from="WGS84 ellipsoidal (EPSG:4979)",
    )

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.

Code examples #

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 #

import rasterio
from rasterio.transform import xy


def shift_dem_datum(src_path, dst_path, *, source_crs="EPSG:4979", target_crs="EPSG:9518",
                    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

    # N varies smoothly, so sample it coarsely and interpolate — far faster than
    # transforming every cell, and accurate to a few centimetres
    rr = np.arange(0, rows, sample_step)
    cc = 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())

    transformer = Transformer.from_crs(source_crs, target_crs, always_xy=True)
    _, _, zeros = transformer.transform(np.array(xs), np.array(ys), np.zeros(len(xs)))
    separation = (-np.asarray(zeros)).reshape(grid_r.shape)

    from scipy.ndimage import zoom
    full = zoom(separation, (rows / separation.shape[0], cols / separation.shape[1]), order=1)
    full = full[: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="EGM2008", vertical_epsg="9518",
                        converted_from=source_crs,
                        mean_separation_m=f"{full.mean():.3f}")

    print(f"  N ranged {full.min():.2f} to {full.max():.2f} m across the tile")
    print(f"  elevations shifted by {-full.mean():+.2f} m on average")
    return dst_path


shift_dem_datum("dem_ellipsoidal.tif", "dem_orthometric.tif")

  N ranged 55.12 to 55.71 m across the tile
  elevations shifted by -55.41 m on average

Across one small tile N varies by 0.59 m — smooth enough that sampling every 32nd cell and interpolating is accurate to centimetres, and roughly a thousand times faster than transforming every cell.

Example 3 — the flood-depth error this causes #

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."
        )
    depth = np.where(dem < water_level_m, water_level_m - dem, 0.0)
    print(f"  flooded {np.mean(depth > 0):.1%} of cells, max depth {depth.max():.2f} m")
    return depth


dem_orth = dem - 51.63                      # ellipsoidal -> EGM2008

correct = flood_depth(dem_orth, 60.0, dem_datum="EGM2008", water_datum="EGM2008")
wrong = flood_depth(dem, 60.0, dem_datum="EGM2008", water_datum="EGM2008")  # unconverted

print(f"  correct: {np.mean(correct > 0):.1%} flooded")
print(f"  using ellipsoidal heights: {np.mean(wrong > 0):.1%} flooded")

  flooded 12.4% of cells, max depth 6.00 m
  flooded 0.0% of cells, max depth 0.00 m
  correct: 12.4% flooded
  using ellipsoidal heights: 0.0% flooded

A 51.6 m offset against a 6 m flood 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.

Explanation #

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.

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.

Edge cases or notes #

• 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 up to a metre or so. ASTER GDEM uses EGM96; Copernicus and SRTM use 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

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 — about +52 m in 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 and SRTM use EGM2008; ASTER GDEM uses 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.