How do I subtract one DEM from another to find elevation change?

First co-register the two DEMs to an identical grid (same CRS, cell size, extent, alignment), then subtract them cell by cell with gdal_calc.py or the QGIS raster calculator to produce a DEM of Difference (DoD) where each cell is new minus old elevation.

How do I align two DEMs with different resolution and extent before comparing?

Use gdalwarp to reproject both to a common CRS, resample to a common cell size with -tr, and snap to a shared grid with -te and -tap so cell centres line up exactly. Resample the finer DEM to the coarser resolution rather than the reverse.

How do I tell real elevation change from noise in a DEM difference?

Estimate the combined vertical uncertainty from each DEM's stated accuracy, set a level-of-detection threshold (often the root-sum-square of the two errors, sometimes at 1.96 sigma for 95% confidence), and mask out differences below it before interpreting change.

Do the two DEMs need the same vertical datum to be compared?

Yes. A mismatch between ellipsoidal and orthometric heights, or between two geoid models, introduces a systematic offset of metres to tens of metres that masquerades as uniform elevation change. Reconcile vertical datums before differencing.

How to Compare Two DEMs

The short answer

Comparing two DEMs means producing a DEM of Difference (DoD) — a raster where each cell is one surface's elevation minus the other's — and then deciding which of those differences are real. The work is 80% preparation: the two DEMs must share an identical CRS, cell size, extent, and grid alignment, and the same vertical datum, before a single subtraction is valid. The interpretation step is to threshold the DoD against the combined vertical uncertainty so that you map genuine change (erosion, deposition, subsidence, construction) and not sensor noise or processing artefacts.

Step 1 — Inspect both DEMs before touching them

Read the metadata, do not assume:

gdalinfo dem_2019.tif
gdalinfo dem_2024.tif

Record for each: horizontal CRS, vertical units and vertical datum, pixel size, origin/extent, NoData value, and data type. Mismatches in any of these will silently corrupt a difference. The two most dangerous and least visible are vertical datum (ellipsoidal vs orthometric, or different geoid models) and NoData handling (a NoData treated as 0 produces a cliff of fake change at the data edge).

Confirm both are genuine elevation rasters (float or integer height values), not 8-bit rendered hillshades — a surprisingly common mistake that makes everything downstream meaningless.

Step 2 — Reconcile the vertical datum

A DEM derived from GNSS/lidar may store ellipsoidal heights, while a national DEM stores orthometric heights (above a geoid such as EGM2008 or a national vertical datum). The difference is the geoid undulation, which ranges from roughly −100 m to +85 m globally and varies smoothly across a region. Difference two DEMs on different vertical datums and you get a broad, gently varying offset that looks like regional uplift but is an artefact.

If both DEMs are on the same vertical datum, proceed. If not, transform one — PROJ supports vertical transformations with geoid grids (for example cs2cs or gdalwarp with a compound CRS such as EPSG:4979 to EPSG:4326+5773). When you cannot fully reconcile the datum, at least estimate and document the residual offset by comparing stable, flat reference areas.

Step 3 — Co-register to a common grid

Both rasters must occupy the same grid: identical CRS, cell size, extent, and cell-centre alignment. Pick the coarser resolution as the target — upsampling the finer DEM invents detail and biases the comparison.

# reproject + resample both to a shared 5 m grid, snapped and clipped to a common extent
gdalwarp -t_srs EPSG:32631 -tr 5 5 -tap \
  -te 448000 5410000 452000 5414000 \
  -r bilinear dem_2019.tif dem_2019_aligned.tif

gdalwarp -t_srs EPSG:32631 -tr 5 5 -tap \
  -te 448000 5410000 452000 5414000 \
  -r bilinear dem_2024.tif dem_2024_aligned.tif

Notes on the flags: -tr sets the target resolution; -tap (target-aligned pixels) forces both outputs onto the same grid lines; -te clips to a shared bounding box so the rasters are pixel-for-pixel coincident. Use bilinear for continuous elevation (cubic can overshoot at breaklines; nearest-neighbour is wrong for continuous data). For change detection at the limit of resolution, consider a horizontal co-registration step (the Nuth & Kääb approach aligns DEMs by minimising the relationship between elevation residuals and terrain aspect/slope) to remove sub-pixel horizontal shifts that otherwise dominate the DoD on steep slopes.

Step 4 — Compute the DEM of Difference

With aligned grids, subtract:

gdal_calc.py -A dem_2024_aligned.tif -B dem_2019_aligned.tif \
  --outfile=dod.tif --calc="A-B" --NoData=-9999

Positive cells gained elevation (deposition, vegetation growth, construction); negative cells lost it (erosion, subsidence, excavation). In QGIS the Raster Calculator does the same with the expression "dem_2024@1" - "dem_2019@1", but verify both layers report identical extent and resolution first, or QGIS will resample on the fly with defaults you did not choose.

Step 5 — Separate signal from noise

A raw DoD is never zero on stable ground; it scatters around zero because of each DEM's vertical error. To map real change, estimate the level of detection (LoD).

If the two DEMs have vertical RMSE values σ₁ and σ₂, the uncertainty of their difference is the root-sum-square:

σ_dod = sqrt(σ₁² + σ₂²)

A common minimum LoD is σ_dod itself; for a 95% confidence threshold use 1.96 × σ_dod. For example, two airborne lidar DEMs each with σ ≈ 0.15 m give σ_dod ≈ 0.21 m, and a 95% LoD ≈ 0.42 m. Mask the DoD so cells with |Δz| < LoD are set to NoData:

gdal_calc.py -A dod.tif --outfile=dod_thresh.tif \
  --calc="A*(abs(A)>0.42)" --NoData=0

Calibrate σ empirically by sampling the DoD over stable, flat areas (paved surfaces, bedrock) that should not have changed — the standard deviation there is your real combined uncertainty, often larger than the nominal spec. Then quantify volumes: multiply each surviving cell's Δz by the cell area and sum the positive and negative parts separately to get gross deposition and gross erosion, and their algebraic sum for net volume change.

Common pitfalls and why they happen

Differencing un-aligned grids. Without -tap/-te, cell centres are offset by a fraction of a pixel, producing a checkerboard or aspect-correlated false signal. It happens because the source DEMs rarely share an origin.
Vertical datum mismatch. Produces a smooth regional offset mistaken for tectonics or settlement; invisible because basemaps don't show height.
Upsampling the coarse DEM. Fabricates resolution and biases volumes; people do it to "keep detail" from the finer DEM.
Ignoring horizontal misregistration on steep terrain. A 1 m horizontal shift on a 30° slope creates ~0.6 m of false vertical difference, concentrated on slopes facing one direction.
Not masking by LoD. Reporting every cell as change inflates volumes with noise that nearly cancels but biases totals.

QA and validation

Histogram the DoD over stable areas — it should be centred near zero; a non-zero mean signals a datum or registration offset to correct.
Confirm both inputs share CRS, cell size, extent, and alignment with gdalinfo before differencing.
Cross-check a few change cells against independent evidence (imagery dates, survey, known works).
Report volumes with their LoD and the stable-area-derived uncertainty, not as exact figures.

Bathyl perspective

We treat a DoD as a measurement with error bars, never a finished map. The deliverable states each DEM's source, date, vertical datum, the alignment method, the level of detection, and the volumes split into gross gain, gross loss, and net — so a reviewer can judge how much of the change is real.