Digital Elevation Models Explained: DEM, DSM and DTM #

Problem statement #

You download "a DEM", compute building heights from it, and the buildings come out 0 m tall. Or you compute a flood extent and the water flows through a forest as if the trees were solid ground.

Both failures come from the same misunderstanding. "DEM" is a category, and the two products inside it measure different surfaces:

import rasterio

with rasterio.open("elevation.tif") as src:
    print(src.width, src.height, src.crs, src.dtypes[0], src.nodata)
    print(src.res)

2400 3600 EPSG:4326 float32 None
(0.00041666666666666664, 0.0002777777777777778)

Four things in that output decide whether your analysis works, and none of them is the elevation. The CRS is geographic — so the "cell size" is in degrees. The two resolutions differ — the pixels are not square. nodata is None — so sea and voids may be encoded as something you have to discover. And nothing at all says whether this surface includes trees and buildings.

Quick answer #

Three terms, and the middle one is the ambiguous one:

# height above ground = surface minus terrain
building_height = dsm - dtm

That subtraction is the whole reason the distinction matters. Get a DSM when you wanted a DTM and hydrology runs over treetops; get a DTM when you wanted a DSM and every building is 0 m tall.

The global products you will actually meet:

Step-by-step solution #

1. Establish what surface you have #

The metadata rarely says "DSM" outright. Three checks that do settle it:

import numpy as np

with rasterio.open("elevation.tif") as src:
    dem = src.read(1)
    print("tags:", src.tags())
    print(f"range {np.nanmin(dem):.1f} to {np.nanmax(dem):.1f}")

• The product name. Copernicus DEM, SRTM and ASTER are all DSMs. Anything called "terrain" or "bare earth" is a DTM.
• Look at a forest edge. A DSM shows a 20 m step at the treeline; a DTM does not.
• Look at a city. A DSM has blocky rectangles where buildings are.

Getting this wrong is not subtle in its consequences and is completely silent at load time.

2. Check the CRS before you compute anything from it #

print(src.crs, src.crs.is_geographic)
print(src.res)

EPSG:4326 True
(0.00041666666666666664, 0.0002777777777777778)

Most global DEMs ship in EPSG:4326. That means the cell size is in degrees, and any calculation involving horizontal distance — slope, aspect, hillshade, volume — is wrong unless you convert.

Worse, the two resolutions differ. Copernicus DEM varies its longitude spacing by latitude band so that cells stay roughly square on the ground. At 53°N:

import math
lat = 53.065
x_m = src.res[0] * 111_320 * math.cos(math.radians(lat))
y_m = src.res[1] * 110_574
print(f"cell {x_m:.1f} m x {y_m:.1f} m")

cell 27.9 m x 30.7 m

Not 30 × 30, and not square. Assuming square pixels here introduces a 4.5% error into slope — small enough to survive review and large enough to matter.

3. Find out how NoData is encoded #

print("declared nodata:", src.nodata)
print("suspicious values:", np.unique(dem[dem < -100])[:5])

declared nodata: None
suspicious values: [-32767. -9999.]

A declared nodata of None does not mean there is none. Common encodings are -32768, -9999, and 0 for sea. Treating -32767 as an elevation gives you a hole 32 km deep, which then dominates every statistic and every colour ramp.

dem = np.where(dem < -100, np.nan, dem)

4. Understand the vertical reference #

Elevation is measured from somewhere, and global DEMs and GPS use different zeros. Copernicus DEM heights are above the EGM2008 geoid; a raw GPS fix is above the WGS84 ellipsoid. The two differ by tens of metres and the difference varies geographically — see vertical datums explained.

If you are only computing slope or drawing contours, this does not matter: those depend on differences in height, and a constant offset cancels. If you are computing flood depth, clearance or absolute height above sea level, it matters enormously.

5. Sanity-check against something you know #

print(f"{np.nanmin(dem):.1f} m to {np.nanmax(dem):.1f} m")

54.0 m to 1074.5 m

Snowdon's summit is 1,085 m. A maximum of 1,074.5 m from a 30 m grid — where no cell centre falls exactly on the peak — is right. A maximum of 3,000 m or −32,767 would not be.

One known elevation is enough to catch unit errors, datum errors and NoData contamination at once.

Code examples #

Example 1 — loading a DEM with every assumption checked #

import math

import numpy as np
import rasterio


def load_dem(path, *, nodata_below=-100, expect_max=None):
    """Read a DEM, normalise NoData, and report the cell size in metres."""
    with rasterio.open(path) as src:
        dem = src.read(1).astype("float64")
        crs, transform, res = src.crs, src.transform, src.res
        declared = src.nodata
        bounds = src.bounds

    if declared is not None:
        dem = np.where(dem == declared, np.nan, dem)
    hidden = np.isfinite(dem) & (dem < nodata_below)
    if hidden.any():
        print(f"  {hidden.sum():,} cells below {nodata_below} m — undeclared NoData "
              f"(values {np.unique(dem[hidden])[:3]})")
        dem = np.where(hidden, np.nan, dem)

    if crs.is_geographic:
        lat = (bounds.bottom + bounds.top) / 2
        cell_x = res[0] * 111_320 * math.cos(math.radians(lat))
        cell_y = res[1] * 110_574
        print(f"  geographic CRS: cell is {cell_x:.1f} m x {cell_y:.1f} m at {lat:.2f}°N "
              f"(NOT {res[0]} x {res[1]})")
    else:
        cell_x, cell_y = res
        print(f"  projected CRS: cell is {cell_x:.1f} x {cell_y:.1f} "
              f"{crs.axis_info[0].unit_name}")

    valid = np.isfinite(dem)
    print(f"  {dem.shape} · {valid.mean():.1%} valid · "
          f"{np.nanmin(dem):.1f} m to {np.nanmax(dem):.1f} m")

    if expect_max is not None:
        gap = abs(np.nanmax(dem) - expect_max)
        flag = "ok" if gap < expect_max * 0.05 else "SUSPECT"
        print(f"  max vs expected {expect_max} m: off by {gap:.1f} m — {flag}")

    return dem, {"transform": transform, "crs": crs,
                 "cell_x": cell_x, "cell_y": cell_y}


dem, meta = load_dem("snowdonia_glo30.tif", expect_max=1085)

  geographic CRS: cell is 27.9 m x 30.7 m at 53.07°N (NOT 0.0004166666666666666 x 0.0002777777777777778)
  (252, 239) · 100.0% valid · 54.0 m to 1074.5 m
  max vs expected 1085 m: off by 10.5 m — ok

The cell_x/cell_y values are what every downstream terrain calculation needs. Returning them alongside the array is what stops the next function guessing.

Example 2 — measuring object heights from a DSM and a DTM #

def object_heights(dsm_path, dtm_path, *, min_height=2.0):
    """Height above ground — the reason the DSM/DTM distinction exists."""
    dsm, dsm_meta = load_dem(dsm_path)
    dtm, dtm_meta = load_dem(dtm_path)

    if dsm.shape != dtm.shape:
        raise ValueError(f"grids differ: {dsm.shape} vs {dtm.shape} — resample one first")
    if dsm_meta["transform"] != dtm_meta["transform"]:
        raise ValueError("transforms differ — the grids are not aligned")

    height = dsm - dtm
    standing = height > min_height

    print(f"  height above ground: median {np.nanmedian(height):.2f} m, "
          f"max {np.nanmax(height):.1f} m")
    print(f"  {standing.mean():.1%} of cells have something over {min_height} m tall")
    return height


heights = object_heights("area_dsm.tif", "area_dtm.tif")

  height above ground: median 0.04 m, max 31.2 m
  8.3% of cells have something over 2.0 m tall

The transform check is not optional. Two grids of the same shape can still be offset by half a cell, and subtracting them then produces a height map full of edge artefacts that look like real structures — the vector-raster version of misaligned layers.

A median of 0.04 m is the reassuring number: over most of the area the two surfaces agree, as they should where there is nothing standing.

Example 3 — fetching a DEM from a STAC catalogue #

import os

os.environ["GDAL_DISABLE_READDIR_ON_OPEN"] = "EMPTY_DIR"
os.environ["AWS_NO_SIGN_REQUEST"] = "YES"

import rasterio
from pystac_client import Client
from rasterio.windows import from_bounds

BBOX = [-4.10, 53.03, -4.00, 53.10]          # Snowdonia, WGS84

catalog = Client.open("https://earth-search.aws.element84.com/v1")
items = list(catalog.search(collections=["cop-dem-glo-30"], bbox=BBOX).items())
print(f"{len(items)} tiles: {[i.id for i in items]}")

with rasterio.open(items[0].assets["data"].href) as src:
    window = from_bounds(*BBOX, src.transform)
    dem = src.read(1, window=window)
    transform = src.window_transform(window)

print(f"{dem.shape}, {dem.min():.1f} m to {dem.max():.1f} m")

2 tiles: ['Copernicus_DSM_COG_10_N53_00_W005_00_DEM', 'Copernicus_DSM_COG_10_N53_00_W004_00_DEM']
(252, 239), 54.0 m to 1074.5 m

Two tiles for a small bounding box, because the extent crosses a one-degree tile boundary — the normal case, and the reason mosaicking is part of most terrain workflows.

Reading a window rather than the whole 2400 × 3600 tile costs a fraction of the bandwidth. See searching and downloading with STAC for the mechanics.

Explanation #

Why global DEMs are almost all DSMs #

The global products are built from radar (SRTM, Copernicus/TanDEM-X) or optical stereo (ASTER). Both measure whatever reflects the signal first — canopy, roof, ground where nothing stands on it.

Producing a DTM requires classifying returns as ground or not-ground and interpolating the ground surface underneath everything else. That is straightforward with LiDAR, which records multiple returns per pulse and gets some energy through gaps in canopy, and very hard with radar or stereo imagery.

So: global coverage means DSM. Bare-earth DTMs exist where someone has flown LiDAR, which is mostly national programmes over populated areas.

Why forest is the worst case #

Over a city a DSM is obviously a DSM — the buildings are visible as blocks. Over forest it is much subtler. The surface is smooth, plausible and 15–25 m too high, uniformly, across the whole wooded area.

Every downstream product inherits the error:

• Flood modelling — water flows over the canopy surface, so valleys under trees route incorrectly.
• Viewsheds — you can see over the treetops.
• Slope — near-correct in the forest interior, badly wrong at every forest edge, where a 20 m step becomes a cliff.

The forest edge artefact is the most reliable way to tell a DSM from a DTM by eye.

Why resolution is not accuracy #

A 30 m DEM has 30 m cells. It does not have 30 m vertical accuracy, and it does not resolve 30 m features.

Copernicus GLO-30 specifies about 4 m absolute vertical accuracy, better in flat terrain and worse on steep slopes where a small horizontal error produces a large vertical one. And a feature needs several cells across it to be represented at all — a 30 m grid does not show a 30 m gully, it shows a hint of one.

The practical consequences: do not interpolate a 30 m DEM to 5 m and expect 5 m detail, and do not report elevations to more decimal places than the accuracy supports.

Why voids exist and what to do about them #

SRTM has holes — radar shadow in steep terrain, and no coverage above 60°N. ASTER has cloud artefacts. Copernicus DEM is void-filled and is the reason it replaced SRTM as the default.

Voids are not always flagged as NoData. Sometimes they are filled with an interpolated surface that is smoother than the real terrain, which shows up as an unnaturally flat patch in a mountain range. If a slope map has a suspiciously smooth region, check whether it is a filled void rather than a plateau.

Edge cases or notes #

• A nodata of None does not mean there is none. Check for -32768, -9999 and 0-as-sea explicitly.
• Copernicus DEM tiles are one degree square, so most study areas need two or more mosaicked.
• Longitude spacing changes with latitude band in Copernicus DEM, to keep ground cells roughly square. Never assume the two resolutions are equal.
• Elevation dtype is usually float32, so int arithmetic on it silently truncates.
• Sea level is not always zero. Some products use a small positive value over water, others NoData, others real bathymetry.
• DSM minus DTM can be negative in small amounts from independent errors in the two products. Clip at zero before reporting heights.
• The vertical datum matters for absolute heights, not for slope, aspect or contours, which depend only on differences.
• LiDAR DTMs may be interpolated under buildings, so "ground" beneath a large building is a guess.

Internal links #

• Vertical datums explained: why your elevations are off — the zero the heights are measured from
• Slope and aspect explained — the first thing most people compute from a DEM
• How to calculate slope and aspect from a DEM in Python — with the cell size handled properly
• The raster data model explained — dtype, NoData and the affine transform
• How to search and download satellite imagery with STAC in Python — where the DEM above came from
• How to merge and mosaic rasters with rasterio — joining one-degree tiles
• Raster and vector do not line up in Python — misaligned grids in a DSM/DTM subtraction
• Raster resampling explained — before matching two elevation grids

FAQ #

What is the difference between a DSM and a DTM? #

A DSM is the top surface — canopy, roofs, ground where nothing stands. A DTM is the bare ground with everything removed. Their difference is the height of objects.

Which one is a "DEM"? #

Either. The term is generic, so check the product documentation. Almost every free global product — Copernicus, SRTM, ASTER — is a DSM.

Which global DEM should I use? #

Copernicus DEM GLO-30 unless you have a reason not to. It is current, void-filled and free, and it superseded SRTM as the practical default.

Why is my cell size a fraction of a degree? #

Because the DEM is in EPSG:4326. Convert to metres using the latitude before any slope, aspect or volume calculation.

Why are the x and y resolutions different? #

Copernicus DEM varies longitude spacing by latitude band so cells stay roughly square on the ground. At 53°N a GLO-30 cell is about 27.9 m by 30.7 m.

What are the huge negative values in my DEM? #

Undeclared NoData, usually −32768 or −9999. Mask anything below a sensible floor before computing statistics.

Does the vertical datum matter? #

For absolute heights — flood depth, clearance, height above sea level — yes, and the difference is tens of metres. For slope, aspect and contours, no: a constant offset cancels.