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. Here is the free global default, the Copernicus DEM, cut to Snowdonia and the coast of Caernarfon Bay:

import rasterio

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

960 360 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 last 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:

Which surface, which CRS, which NoData, which datum

1. Establish what surface you have

The metadata rarely says "DSM" outright:

import numpy as np

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

tags: {'AREA_OR_POINT': 'Point'}
range -0.6 to 1074.5

Nothing there says surface or terrain. The one tag present, AREA_OR_POINT=Point, says each value is a sample at the cell's reference point rather than an average over the cell — a detail that matters when you line this grid up with another one. Three checks that do settle the surface question:

• The product name. Copernicus DEM, SRTM and ASTER are all DSMs — the Copernicus files even say so in their names, Copernicus_DSM_COG_10_…. Anything called "terrain" or "bare earth" is a DTM.
• Look at a forest edge. A DSM shows a step at the treeline; a DTM does not. Example 4 measures the step.
• 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

with rasterio.open("elevation.tif") as src:
    print(src.crs, src.crs.is_geographic)
    print(src.res)
    bounds, res, declared = src.bounds, src.res, src.nodata

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 widens its longitude spacing by latitude band so that cells stay roughly square on the ground — but only roughly. At this latitude:

import math

lat = (bounds.bottom + bounds.top) / 2
x_m = res[0] * 111_320 * math.cos(math.radians(lat))
y_m = res[1] * 110_574
print(f"at {lat:.2f}°N a cell is {x_m:.1f} m x {y_m:.1f} m")

land = np.where(dem > 0, dem, np.nan)
def slope_deg(z, dx, dy):
    gy, gx = np.gradient(z, dy, dx)
    return np.degrees(np.arctan(np.hypot(gx, gy)))
right, square = slope_deg(land, x_m, y_m), slope_deg(land, y_m, y_m)
print(f"treating the cells as {y_m:.1f} m squares changes slope by a median "

      f"{np.nanmedian(np.abs(square - right) / right):.1%}, and by up to {np.nanmax(np.abs(square - right)):.1f}° in one cell")

at 53.07°N a cell is 27.9 m x 30.7 m
treating the cells as 30.7 m squares changes slope by a median 4.4%, and by up to 2.8° in one cell

Not 30 × 30, and not square. Assuming square pixels introduces an error into slope that is small enough to survive review and large enough to matter on the steep ground where slope is usually the point.

3. Find out how NoData is encoded

print("declared nodata:", declared)
print(f"cells exactly 0 m: {(dem == 0).sum():,}   below 0 m: {(dem < 0).sum()}   "
      f"below -100 m: {(dem < -100).sum()}")
print(f"mean with the sea: {dem.mean():.1f} m   without it: {dem[dem != 0].mean():.1f} m")

declared nodata: None
cells exactly 0 m: 48,983   below 0 m: 80   below -100 m: 0
mean with the sea: 237.8 m   without it: 277.0 m

A declared nodata of None does not mean there is nothing to mask. Copernicus has no voids, but it sets the open sea to exactly 0, and in this window that is thousands of cells of Caernarfon Bay pulling every statistic towards sea level. The handful of slightly negative cells are real — coastal land below mean sea level is common. Other products encode missing data as -32768 or -9999, which is worse: a hole 32 km deep dominates every statistic and every colour ramp. Decide per product, and write it down:

dem = np.where(dem == 0, np.nan, dem)     # Copernicus: open sea is exactly 0

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 from place to place — 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")

-0.6 m to 1074.5 m

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

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

Copernicus and Environment Agency surfaces, loaded and compared

The data used on this page

Three open sources. The Copernicus DEM (GLO-30) for Snowdonia and Caernarfon Bay, read from its public bucket; the window crosses the 4°W tile boundary, so two tiles are merged. The Environment Agency's lidar surface and terrain models (DSM and DTM) for 1.5 km² of central Leeds at 1 m, from its Web Coverage Services. And, for Example 4, the lidar terrain model of Epping Forest at 5 m.

import numpy as np
import rasterio
import requests
from rasterio.merge import merge

BBOX = (-4.40, 53.02, -4.00, 53.12)          # Snowdon to Caernarfon Bay, WGS84
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("elevation.tif", "w", **profile) as dst:
    dst.write(mosaic)
    dst.update_tags(**tags)

EA = "https://environment.data.gov.uk/spatialdata/{}/wcs"
LIDAR = {
    "dtm": ("lidar-composite-digital-terrain-model-dtm-1m",
            "13787b9a-26a4-4775-8523-806d13af58fc__Lidar_Composite_Elevation_DTM_1m"),
    "dsm": ("lidar-composite-digital-surface-model-first-return-dsm-1m",
            "df4e3ec3-315e-48aa-aaaf-b5ae74d7b2bb__Lidar_Composite_Elevation_FZ_DSM_1m"),
}
def lidar(kind, box, path, scale=None):
    slug, coverage = LIDAR[kind]
    params = [("service", "WCS"), ("version", "2.0.1"), ("request", "GetCoverage"),
              ("CoverageId", coverage), ("format", "image/tiff"),
              ("subset", f"E({box[0]},{box[2]})"), ("subset", f"N({box[1]},{box[3]})")]
    if scale:
        params.append(("SCALEFACTOR", str(scale)))
    r = requests.get(EA.format(slug), params=params, timeout=300)
    r.raise_for_status()
    with open(path, "wb") as f:
        f.write(r.content)

LEEDS = (429000, 432500, 430500, 434000)        # British National Grid
lidar("dsm", LEEDS, "leeds_dsm.tif")
lidar("dtm", LEEDS, "leeds_dtm.tif")
lidar("dtm", (540000, 196000, 545000, 201000), "epping_dtm.tif", scale=0.2)

for name in ("elevation.tif", "leeds_dsm.tif", "leeds_dtm.tif", "epping_dtm.tif"):
    with rasterio.open(name) as src:
        print(f"{name:15} {src.width:>5} × {src.height:<5} {src.crs}  {src.dtypes[0]}")

elevation.tif     960 × 360   EPSG:4326  float32
leeds_dsm.tif    1500 × 1500  EPSG:27700  float32
leeds_dtm.tif    1500 × 1500  EPSG:27700  float32
epping_dtm.tif   1000 × 1000  EPSG:27700  float32

Example 1 — loading a DEM with every assumption checked

import math


def load_dem(path, *, nodata_values=(), 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

    print(path)
    if declared is not None:
        dem = np.where(dem == declared, np.nan, dem)
    for value in nodata_values:                     # undeclared codes you know about
        hidden = dem == value
        print(f"  {hidden.sum():,} cells equal to {value} — treated as NoData")
        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]:.6f} x {res[1]:.6f})")
    else:
        cell_x, cell_y = res

        print(f"  projected CRS: cell is {cell_x:.1f} x {cell_y:.1f} {crs.linear_units}")

    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("elevation.tif", nodata_values=[0], expect_max=1085)

elevation.tif
  48,983 cells equal to 0 — treated as NoData
  geographic CRS: cell is 27.9 m x 30.7 m at 53.07°N (NOT 0.000417 x 0.000278)
  (360, 960) · 85.8% valid · -0.6 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

The Environment Agency flies lidar and publishes both surfaces from the same survey, on the same 1 m grid:

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, min {np.nanmin(height):.1f} m")
    print(f"  {standing.mean():.1%} of cells have something over {min_height} m tall")
    return height


heights = object_heights("leeds_dsm.tif", "leeds_dtm.tif")
with rasterio.open("leeds_dsm.tif") as src:
    row, col = np.unravel_index(np.nanargmax(heights), heights.shape)
    print(f"  tallest cell at {src.xy(row, col)[0]:.0f} E, {src.xy(row, col)[1]:.0f} N")

leeds_dsm.tif
  projected CRS: cell is 1.0 x 1.0 metre
  (1500, 1500) · 100.0% valid · 23.7 m to 135.1 m
leeds_dtm.tif
  projected CRS: cell is 1.0 x 1.0 metre
  (1500, 1500) · 100.0% valid · 23.5 m to 58.3 m
  height above ground: median 0.32 m, max 107.0 m, min -6.2 m
  41.4% of cells have something over 2.0 m tall
  tallest cell at 429882 E, 432936 N

The tallest cell sits on Bridgewater Place, the office tower by Leeds station that is 112 m to its highest point; 107 m of it shows in a 1 m lidar surface. Two in five cells of a city centre have something over 2 m standing on them. And the minimum is negative: in a few cells the surface model lies below the terrain model, which is impossible physically and normal in practice — two products interpolated separately disagree at the edges of buildings and under bridges. Clip at zero before reporting heights.

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 raster version of misaligned layers.

Example 3 — finding the tiles in a STAC catalogue

The setup above built the tile URLs by hand. A STAC catalogue finds them for you:

import os
from pystac_client import Client

os.environ["AWS_NO_SIGN_REQUEST"] = "YES"      # the catalogue gives s3:// links to a public bucket
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: {sorted(i.id for i in items)}")

found, _ = merge([i.assets["data"].href for i in items], bounds=BBOX)
print(f"{found.shape[1:]}, {found.min():.1f} m to {found.max():.1f} m, "
      f"same as elevation.tif: {np.array_equal(found, mosaic)}")

2 tiles: ['Copernicus_DSM_COG_10_N53_00_W004_00_DEM', 'Copernicus_DSM_COG_10_N53_00_W005_00_DEM']
(360, 960), -0.6 m to 1074.5 m, same as elevation.tif: True

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 one tile with a window over the whole box would silently return only the part that tile covers. merge(..., bounds=...) reads just the window from each Cloud-Optimised GeoTIFF rather than the whole 2400 × 3600 tiles. See searching and downloading with STAC for the mechanics.

Example 4 — how far above the ground a global DSM sits in a forest

Epping Forest has lidar ground under its trees. Put the Copernicus DSM and ESA's WorldCover land-cover map onto that grid and compare them by class:

from rasterio.warp import reproject, Resampling

with rasterio.open("epping_dtm.tif") as src:
    ground = src.read(1, masked=True).filled(np.nan)
    grid = {"dst_transform": src.transform, "dst_crs": src.crs}

def onto_grid(url, resampling):
    out = np.full(ground.shape, np.nan, dtype="float32")
    with rasterio.open(url) as src:
        reproject(rasterio.band(src, 1), out, resampling=resampling, **grid)
    return out

copernicus = onto_grid("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N51_00_E000_00_DEM/"
                       "Copernicus_DSM_COG_10_N51_00_E000_00_DEM.tif", Resampling.bilinear)
cover = onto_grid("https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/"
                  "ESA_WorldCover_10m_2021_v200_N51E000_Map.tif", Resampling.nearest)

above = copernicus - ground
for code, name in [(10, "trees"), (30, "grassland"), (50, "built-up")]:
    cells = (cover == code) & np.isfinite(above)
    print(f"{name:10} {cells.sum():>8,} cells   Copernicus above lidar ground: "
          f"median {np.median(above[cells]):5.1f} m, 90th percentile {np.percentile(above[cells], 90):5.1f} m")

trees       693,657 cells   Copernicus above lidar ground: median  10.3 m, 90th percentile  16.2 m
grassland   147,575 cells   Copernicus above lidar ground: median  -0.1 m, 90th percentile   2.2 m
built-up    120,672 cells   Copernicus above lidar ground: median   1.1 m, 90th percentile   2.5 m

Over grassland the two agree to within about ten centimetres at the median. Over the forest, the global model sits about ten metres above the ground, because the radar that built it saw the canopy. Nothing in the file tells you that; only a comparison with a bare-earth model does.

Why global DEMs are surface models

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 too high across the whole wooded area — a median of about ten metres over Epping Forest in Example 4, more in taller or denser woodland.

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.

Sea as zero, tile edges, dtypes and negative heights

• A nodata of None does not mean there is nothing to mask. Copernicus uses exactly 0 for open sea; other products use -32768 or -9999. Check each product's documentation and the values.

• 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 — down to −6 m in a few cells of central Leeds — from the two products being interpolated separately. 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
• DEM hydrology explained — what a DEM says about where water goes
• Bathymetry explained: depths, datums and grids — the same ideas below sea level
• SAR geometry explained — why radar needs a DEM to be mapped at all

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.