Grade is the attribute that separates a plausible ETA from an accurate one on anything but flat terrain, and it is the one attribute OSM almost never carries. A handful of ways have incline tags; the rest have to be derived by sampling a digital elevation model along their geometry. This page covers that derivation end to end, and it feeds directly into the cost work described in speed profile calibration for heavy vehicles within the wider OSM Graph Architecture & Network Modeling pipeline.

The naive version of this — subtract the elevation at one end of an edge from the other, divide by length — is quick to write and wrong in three separate ways. It averages away the pitches that actually slow a vehicle, it produces nonsense on bridges and in tunnels, and it silently treats DEM voids as sea level. Each of those failures is easy to fix once you know it exists, and each is invisible until someone notices ETAs are wrong in the hills.

Prerequisites

Libraries. rasterio for DEM access, shapely for geometry densification, numpy for the vectorised profile maths, and pyproj if the DEM and graph disagree on CRS.

pip install rasterio>=1.3 shapely>=2.0 numpy>=1.24 pyproj>=3.6 geopandas>=0.14

Elevation data. Copernicus GLO-30 is the current default for most of the world: 30 metre posting, void-filled, and freely redistributable. SRTM 1-arcsecond remains widely used and covers 60°N to 56°S. National products — LIDAR-derived DTMs published by mapping agencies — are far better where they exist and are worth the extra handling for cycling or pedestrian work.

Storage. A national extract’s DEM coverage runs to tens of gigabytes as individual tiles. Building a single mosaicked Cloud-Optimised GeoTIFF up front costs disk but turns every subsequent sample into one windowed read rather than a tile lookup.

A graph with geometry. Edge geometries, not just endpoints. A graph built with only node coordinates cannot be densified, and retrofitting geometry after the fact means re-parsing the PBF.

Conceptual architecture

The pipeline is a sampling problem wrapped in two correction layers. Sampling produces a raw elevation profile per edge; the corrections remove the artefacts that make raw profiles unusable.

From DEM tiles to a directional grade attribute DEM tiles and edge geometries feed a densification and sampling stage that produces a raw elevation profile per edge. Two correction layers follow — structure suppression for bridges and tunnels, and void handling for missing DEM cells — before the profile is aggregated into signed per-direction grade values. Sampling is the easy part; the two correction layers are where accuracy comes from DEM mosaic edge geometries densify + sample every 25–50 m structure suppression bridge · tunnel void handling never fill with zero grade_fwd grade_rev = −grade_fwd Skipping structure suppression puts a false 9 % descent on every valley-crossing viaduct in the extract. Skipping void handling puts a false 400 m cliff wherever the DEM has a gap, because the fill value reads as sea level. Both artefacts survive every downstream check, because a steep grade is a perfectly valid number. Carry a confidence flag alongside the value so a corrected edge can be weighted differently from a directly sampled one.

The critical design decision is the sampling interval, and it is a trade between fidelity and cost. Sampling every 25 metres on a national extract produces hundreds of millions of raster reads; sampling every 200 metres misses the pitches that matter for anything heavier than a car. Vehicle profiles tolerate 50 metres comfortably; cycling and pedestrian profiles want 20 to 25.

Step-by-step implementation

1. Mosaic the DEM once

Individual tiles work but make every sample a lookup-plus-open. A single mosaicked, internally tiled COG turns the whole thing into windowed reads against one file handle.

# Build one internally tiled, compressed mosaic from a directory of DEM tiles
gdalbuildvrt dem.vrt dem_tiles/*.tif

gdal_translate dem.vrt dem_mosaic.tif \
  -of COG \
  -co COMPRESS=DEFLATE \
  -co BLOCKSIZE=512 \
  -co OVERVIEWS=NONE

Skip overviews. They exist to accelerate zoomed-out rendering; every read this pipeline performs is at full resolution, so overviews cost storage and buy nothing.

2. Densify edge geometries

Sampling only at existing vertices inherits OSM’s own vertex spacing, which is driven by how the mapper drew the road rather than by terrain. Interpolate a fixed interval instead.

# requires: shapely, numpy (pip install shapely numpy)
import numpy as np
from shapely.geometry import LineString
from shapely.ops import transform


def densify(line: LineString, interval_m: float, project) -> np.ndarray:
    """Return sample points at a fixed spacing along a line, in the DEM's CRS."""
    metric = transform(project, line)
    n = max(int(np.ceil(metric.length / interval_m)), 1)
    # Include both endpoints so the profile covers the full edge
    fractions = np.linspace(0.0, 1.0, n + 1)
    pts = [metric.interpolate(f, normalized=True) for f in fractions]
    return np.array([(p.x, p.y) for p in pts]), metric.length

Densifying in a projected CRS matters for the same reason it always does: an interval expressed in degrees is a different distance at every latitude, so the effective sampling density would vary across the extract.

3. Sample elevation along the profile

rasterio.sample reads a list of coordinates in one call, which is dramatically faster than per-point reads.

# requires: rasterio, numpy (pip install rasterio numpy)
import rasterio
import numpy as np


def sample_profile(dem: rasterio.DatasetReader, xy: np.ndarray) -> np.ndarray:
    """Elevation at each sample point, with the DEM's nodata mapped to NaN."""
    vals = np.fromiter(
        (v[0] for v in dem.sample(xy)), dtype="float64", count=len(xy)
    )
    nodata = dem.nodata
    if nodata is not None:
        # NaN propagates visibly; a nodata sentinel silently becomes a cliff
        vals = np.where(vals == nodata, np.nan, vals)
    return vals

Mapping nodata to NaN rather than leaving the sentinel is the single most important line here. DEM nodata values are typically large negatives such as −32768; left in place they produce grades of several thousand percent, and any clamp applied downstream turns them into a plausible-looking maximum rather than an obvious error.

4. Aggregate the profile into a directional grade

A profile compresses to a grade in more than one way, and the right choice depends on what the cost function does with it.

# requires: numpy (pip install numpy)
import numpy as np


def grade_from_profile(
    elev: np.ndarray, length_m: float, *, mode: str = "mean_abs_weighted"
) -> float:
    """Signed grade in percent, aggregated from an elevation profile."""
    if np.isnan(elev).all() or length_m <= 0:
        return 0.0

    # Linearly bridge isolated voids so one missing cell does not kill the edge
    idx = np.arange(len(elev))
    good = ~np.isnan(elev)
    if good.sum() < 2:
        return 0.0
    elev = np.interp(idx, idx[good], elev[good])

    rises = np.diff(elev)
    if mode == "net":
        # Net rise over the edge — understates rolling terrain
        return 100.0 * (elev[-1] - elev[0]) / length_m
    # Signed by net direction, magnitude from total climbing effort
    sign = np.sign(elev[-1] - elev[0]) or 1.0
    return float(100.0 * sign * np.abs(rises).sum() / length_m)

The net mode is the naive calculation and is appropriate only where terrain is genuinely monotonic across an edge. The default weights by total vertical movement, which is much closer to what a loaded vehicle actually experiences: a rolling edge costs energy on every rise regardless of where it ends up.

5. Suppress structures and flag corrections

# requires: numpy (pip install numpy)
def apply_structure_correction(elev: np.ndarray, tags: dict) -> tuple[np.ndarray, bool]:
    """Replace sampled terrain with a linear ramp across bridges and tunnels."""
    if tags.get("bridge") in (None, "no") and tags.get("tunnel") in (None, "no"):
        return elev, False
    # A structure's deck runs between its endpoints, not over the terrain beneath
    ramp = np.linspace(elev[0], elev[-1], len(elev))
    return ramp, True

Using the sampled endpoints as the ramp’s ends is deliberate: the structure meets the terrain at each end, so those two samples are valid even when everything between them is not. The boolean travels with the edge so a downstream consumer can distinguish a measured grade from an inferred one.

Configuration reference

Parameter Recommended value Notes
Sampling interval 50 m vehicle, 25 m cycling Halving it roughly doubles raster reads for a modest fidelity gain on vehicle profiles
DEM product Copernicus GLO-30 Void-filled and redistributable; national LIDAR products are better where available
DEM resolution 30 m Finer data mostly resolves detail a vehicle cost function averages away again
Aggregation mode weighted by total climb net only where terrain is monotonic across a typical edge
Void fill linear interpolation, flagged Never zero; a zero fill reads as sea level and produces cliffs
Structure handling linear ramp between endpoints Applies to bridge and tunnel; both suffer the same terrain-versus-deck mismatch
Grade storage signed, per direction grade_rev must equal −grade_fwd; assert it
Clamp range ±25 % Real roads exceed 20 % only rarely; beyond that suspect an artefact, not a hill

Production optimization and scaling

Batch reads by tile, not by edge. Sorting the sample points by DEM tile before reading turns a random-access pattern into a sequential one. On a national extract this is frequently a five- to tenfold difference in wall-clock time, entirely from page-cache behaviour.

Choosing the sampling interval Raster read count rises steeply as the sampling interval shortens, while grade fidelity flattens out. The vehicle recommendation of 50 metres sits where fidelity has largely converged; the cycling recommendation of 25 metres buys the remaining detail at four times the read cost. Read cost and fidelity against sampling interval 200 m 100 m 50 m 25 m raster reads grade fidelity vehicle cycling Vehicle cost functions average grade over hundreds of metres anyway, so the detail below 50 m is discarded downstream. Cycling and cargo-bike profiles respond to short pitches, which is the one case where the extra reads change a route.

Compute once, store on the edge. Grade does not change between graph rebuilds unless the geometry does. Persist it as an edge attribute in the same Parquet sidecar that carries the other computed weights, keyed by edge id, and recompute only for edges whose geometry changed — the same incremental logic described in incremental graph rebuilds without full reprocessing.

Parallelise by spatial block. Sampling is embarrassingly parallel across disjoint areas, and blocking by DEM tile keeps each worker’s reads confined to one file region. Avoid parallelising by edge id, which scatters reads across the whole mosaic in every worker.

Keep the profile, not just the aggregate. Storing the sampled profile costs more than storing a single grade value, but it lets you re-derive the aggregate under a different mode without re-sampling. For a graph that feeds both vehicle and cycling profiles, that saves a full re-run every time the cycling cost model is retuned.

Validation and testing

The useful checks are the ones that catch artefacts rather than errors, because artefacts produce valid-looking numbers.

# requires: numpy, pandas (pip install numpy pandas)
import numpy as np
import pandas as pd


def audit_grades(edges: pd.DataFrame) -> dict:
    """Surface the signatures of the three common grade artefacts."""
    g = edges["grade_fwd"].to_numpy()
    return {
        "edges": len(edges),
        "abs_grade_over_20pct": int((np.abs(g) > 20).sum()),
        "abs_grade_over_40pct": int((np.abs(g) > 40).sum()),
        "exactly_zero": int((g == 0).sum()),
        "reverse_mismatch": int(
            (~np.isclose(edges["grade_rev"], -edges["grade_fwd"])).sum()
        ),
        "corrected_share": round(edges["grade_inferred"].mean(), 4),
    }

A non-zero reverse_mismatch means the sign convention broke somewhere, which is the defect that most reliably survives review. A large abs_grade_over_40pct count points at unfilled voids. A suspiciously high exactly_zero count means edges are falling through to the fallback rather than being sampled — usually a CRS mismatch putting the sample points outside the DEM’s footprint.

Terrain profile against the actual road deck on a viaduct Sampling a DEM along a viaduct follows the valley floor, producing a steep descent and climb across a road that is physically level. Replacing the sampled section with a linear ramp between the structure's endpoints recovers the true deck profile. A 600 m viaduct: what the DEM sees, and what the vehicle drives 0 m 300 m 600 m elevation DEM samples — the valley floor deck — linear ramp between endpoints Sampled naively this edge reports a 9 % descent followed by a 9 % climb — on a road that is level to within a metre. The endpoints remain valid because the structure meets the terrain there; only the span between them needs replacing.

Troubleshooting

Grades of several hundred percent appear on scattered edges

Root cause: DEM nodata sentinels are being treated as elevations. A value of −32768 next to a real elevation of 120 produces an apparent drop of hundreds of metres over a few tens of metres of road.

Fix: Map the raster’s declared nodata to NaN immediately after sampling, then interpolate across isolated gaps. Never substitute zero — it reads as sea level and produces the same artefact with a smaller magnitude, which is harder to spot.

Every edge reports a grade of exactly zero

Root cause: The sample points and the DEM are in different coordinate reference systems, so every read falls outside the raster footprint and returns nodata.

Fix: Reproject the sample coordinates into the DEM’s CRS before sampling, and assert that at least one sample per edge is valid. A silent all-nodata result is indistinguishable from flat terrain in the output.

Reverse edges appear to climb in both directions

Root cause: A single grade value was copied to both directed edges rather than negated for the reverse.

Fix: Store grade_fwd and derive grade_rev = -grade_fwd at edge-insertion time, then assert the two sum to zero across the whole graph. This is the defect most likely to survive code review, because both values look individually reasonable.

ETAs are accurate on motorways and poor on rural roads

Root cause: The sampling interval is coarser than the terrain variation on smaller roads. Motorways are engineered to gentle gradients that survive coarse sampling; a rural lane’s pitches do not.

Fix: Reduce the interval for minor road classes rather than globally. Sampling highway=unclassified and below at 25 metres while leaving trunk roads at 50 recovers most of the accuracy for a fraction of the additional reads.

Cycling routes avoid a flat riverside path

Root cause: The path runs under bridges, and the DEM records the bridge decks above it rather than the path surface.

Fix: This is the inverse of the usual structure problem and needs the same treatment: where a way is tagged layer below zero or runs beneath a mapped structure, prefer interpolation over sampling. Flag these edges so their grade is treated as inferred.