Elevation and Terrain Enrichment for Routing Graphs

A flat graph routes a cyclist up a fourteen per cent hill because it is two hundred metres shorter, sends a loaded truck over a pass that its gearbox will crawl, and tells an electric vehicle it has range it does not have. None of those is a routing bug — the algorithm found the cheapest path through the graph it was given, and the graph did not know the road went uphill. Elevation is the one attribute that turns a symmetric edge into an asymmetric one, and adding it changes not just the cost of a route but which route is chosen. This topic covers sampling a digital elevation model onto graph nodes, deriving a gradient per edge in each direction, and the failure modes that make terrain data quietly wrong rather than obviously missing.

Prerequisites

Sampling a DEM needs a raster reader; everything downstream is arithmetic and Cypher.

Requirement Minimum version Install
Python 3.11 —
neo4j (async driver) 5.20 pip install "neo4j>=5.20"
Neo4j Server 5.15 native point
rasterio 1.3 pip install "rasterio>=1.3"
numpy 1.26 pip install "numpy>=1.26"

Core Concept & Mechanism

Elevation enters the graph as a node property and leaves it as a directed edge property, and that transformation is the whole point.

Sampling is a raster lookup per node. A digital elevation model is a grid of height values; a node’s elevation is the value at its coordinate, interpolated from the surrounding cells. Nearest-neighbour sampling is fast and produces visible stair-stepping along a road that crosses cells diagonally; bilinear interpolation over the four neighbours costs a little more and produces a profile smooth enough to differentiate. For gradient work the interpolation is not a nicety — a stair-stepped profile produces alternating zero and extreme gradients on consecutive segments, and any cost model built on it is noise.

Gradient is a per-direction property of the edge, not of the road. The rise between two nodes is the same magnitude whichever way you travel, but its sign flips, and cost is not symmetric in that sign. A cyclist climbing eight per cent is doing several times the work of one descending it; a truck descending eight per cent may be slower than one climbing four, because of braking limits. So the model has to store gradient with a direction, which means either two relationships per road or one relationship whose cost function is evaluated with the direction of travel.

Cost is a function of gradient, not a term added to it. The temptation is to add a penalty proportional to the climb. That produces plausible results in the middle of the range and nonsense at the ends — it never makes a descent free, and it never makes a wall impassable. The models that behave are non-linear: for a cyclist, power required rises with gradient and speed falls out of it; for a vehicle, energy consumption per metre is roughly linear in gradient with a sharp regeneration cut-off on the descending side.

The same road, two directions, two costs — which is why gradient is an edge property A 1.4 kilometre road climbing 84 metres, a six per cent gradient. Travelled uphill by a cyclist it costs about eight minutes; travelled downhill it costs about three. A flat graph gives both directions the same cost derived from length alone and picks a route that is often right in one direction and clearly wrong in the other. Because the sign of the gradient flips with the direction of travel while its magnitude does not, the cost has to be stored per direction rather than per road, which is what makes an elevation-aware graph directed even where the road itself is two-way. 1.4 km road · 84 m rise · 6% gradient A · 42 m B · 126 m A → B climbing B → A descending cyclist cost A → B 8.1 min B → A 2.9 min a factor of 2.8 between the two directions flat graph A → B 5.0 min B → A 5.0 min wrong in both, and wrong in opposite ways The flat figure is not a compromise between the two — it is a number that describes neither journey, and a router using it will send a rider up the hill to save two hundred metres.

Schema & Data Model

Elevation lives on the node; gradient and its derived costs live on the relationship, per direction.

// Node elevation, plus the provenance that makes a resample decidable later.
CREATE INDEX junction_elevation IF NOT EXISTS
FOR (n:Junction) ON (n.elevation_m);

// Directed segments carry their own gradient. A two-way road is two
// relationships precisely so the two gradients can differ in sign.
CREATE INDEX segment_grade IF NOT EXISTS
FOR ()-[s:SEGMENT]-() ON (s.grade_pct);

// Sanity constraint on the source: a node with no elevation must be visible,
// not silently defaulted to zero.
MATCH (n:Junction) WHERE n.elevation_m IS NULL
RETURN count(n) AS missing_elevation;

Storing dem_source and dem_resolution_m alongside the elevation is worth the two properties. Elevation data is patched, replaced and improved; without provenance there is no way to tell which nodes came from a 30-metre global model and which from a 1-metre national LiDAR survey, and the two disagree by enough to matter on exactly the steep, narrow roads where the gradient matters most.

Step-by-Step Implementation

1. Sample the DEM with bilinear interpolation, in batches.

import asyncio
from dataclasses import dataclass

import numpy as np
import rasterio
from neo4j import AsyncGraphDatabase

WRITE_ELEVATION = """
UNWIND $batch AS row
MATCH (n:Junction {id: row.id})
SET n.elevation_m = row.elevation_m,
    n.dem_source = $source,
    n.dem_resolution_m = $resolution
RETURN count(n) AS updated
"""


@dataclass(frozen=True)
class Sampled:
    id: str
    elevation_m: float | None


class DemSampler:
    """Bilinear sampling of a DEM at node coordinates.

    Nearest-neighbour would be faster and is a trap: a road crossing cells
    diagonally comes out as a staircase, and differentiating a staircase gives
    alternating zero and extreme gradients on consecutive segments. Every cost
    model downstream then inherits that noise as if it were terrain.
    """

    def __init__(self, dem_path: str) -> None:
        self._dataset = rasterio.open(dem_path)
        self._nodata = self._dataset.nodata

    def close(self) -> None:
        self._dataset.close()

    def sample(self, points: list[tuple[str, float, float]]) -> list[Sampled]:
        coords = [(lon, lat) for _, lat, lon in points]
        values = list(self._dataset.sample(coords, indexes=1, masked=True))

        out: list[Sampled] = []
        for (node_id, _, _), value in zip(points, values):
            raw = float(value[0]) if not np.ma.is_masked(value[0]) else None
            if raw is None or (self._nodata is not None and raw == self._nodata):
                # A void in the DEM — over water, or outside coverage. Record it
                # as missing rather than as sea level, which is a real elevation.
                out.append(Sampled(node_id, None))
            else:
                out.append(Sampled(node_id, round(raw, 2)))
        return out


async def load_elevation(driver, sampler: DemSampler, rows, source: str,
                         resolution: float, batch: int = 5_000) -> int:
    updated = 0
    async with driver.session() as session:
        for i in range(0, len(rows), batch):
            chunk = rows[i:i + batch]
            sampled = sampler.sample(chunk)
            payload = [
                {"id": s.id, "elevation_m": s.elevation_m}
                for s in sampled if s.elevation_m is not None
            ]
            if not payload:
                continue
            result = await session.run(
                WRITE_ELEVATION, batch=payload, source=source, resolution=resolution
            )
            updated += int((await result.single())["updated"])
    return updated

2. Derive gradient per directed segment. The gradient is rise over run, and the run is the road’s length rather than the straight-line distance between endpoints — a hairpin gains little height over a lot of tarmac.

MATCH (a:Junction)-[s:SEGMENT]->(b:Junction)
WHERE a.elevation_m IS NOT NULL AND b.elevation_m IS NOT NULL
  AND s.length_m > 0
SET s.rise_m   = b.elevation_m - a.elevation_m,
    s.grade_pct = 100.0 * (b.elevation_m - a.elevation_m) / s.length_m
RETURN count(s) AS graded;

3. Turn gradient into a cost the router can minimise. That step is specific to the vehicle, and is worked through for cyclists in grade-aware weights for bicycle routing and for energy in computing edge grade and energy cost.

Query Patterns & Variants

Find the segments a profile cannot use. A loaded truck with a gradient limit, or an accessibility route with a wheelchair limit, is a filter over grade_pct — and because it is a stored property, it is index-seekable rather than computed per row.

MATCH (a:Junction)-[s:SEGMENT]->(b:Junction)
WHERE s.grade_pct > $max_grade_pct
RETURN s.id AS segment_id, s.grade_pct AS grade, s.length_m AS length
ORDER BY grade DESC;

Total ascent along a route, which is the number cyclists actually compare routes on — and note that it sums only the positive rises, because descending does not undo a climb in any metric a rider cares about.

MATCH (t:Trip {id: $trip_id})
UNWIND t.segment_ids AS seg_id
MATCH ()-[s:SEGMENT {id: seg_id}]->()
RETURN sum(CASE WHEN s.rise_m > 0 THEN s.rise_m ELSE 0 END) AS ascent_m,
       sum(CASE WHEN s.rise_m < 0 THEN -s.rise_m ELSE 0 END) AS descent_m;

Detect elevation that disagrees with the road network. A bridge sampled from a bare-earth DEM takes the height of the valley floor beneath it, producing an impossible gradient in and out. Segments whose gradient exceeds anything a road is built to is a cheap, effective detector.

Performance Tuning

Sampling dominates the enrichment, and it is I/O-bound on the raster rather than on the graph:

$$C_{\text{sample}} \approx N \cdot \big(c_{\text{seek}} + c_{\text{interp}}\big)$$

The way to make it fast is spatial locality. A DEM is stored in tiles, and sampling nodes in coordinate order means each tile is read once and used for thousands of lookups; sampling in node-id order means the same tile is read, evicted and re-read repeatedly. Sorting the node batch by cell before sampling is a few lines and routinely gives an order of magnitude, and it costs nothing because the batch is in memory anyway.

The gradient derivation is a single pass over the relationships and is cheap, but it is worth doing in batches with periodic commits rather than as one transaction — a continental graph’s segment count will exceed the transaction memory budget otherwise, and the operation is idempotent so a resumable batch loop loses nothing on failure.

One consequence worth planning for: adding elevation and gradient adds properties to every node and every relationship in the graph, which grows the store and therefore the page cache requirement. Two floats per segment on a continental network is not a rounding error, and it is worth measuring before and after rather than discovering it as a latency regression.

Sampling order decides how many times each DEM tile is read Two orderings of the same four thousand node lookups against a tiled elevation model. Sampling in node-id order visits tiles in essentially random sequence, so each tile is read, evicted and re-read many times and the raster reader spends its time on I/O. Sorting the batch by tile first means every tile is opened once and serves all the lookups that fall inside it, turning thousands of reads into dozens. The sort costs microseconds on a batch that is already in memory. 4,000 lookups against a tiled DEM node-id order each tile opened and evicted repeatedly 3,180 tile reads sorted by tile each tile opened once, then fully used 42 tile reads The lookups are identical and the answers are identical. The only difference is the order they are asked in, and the sort that produces it runs on a list that is already in memory.

Edge Cases & Gotchas

  • Bridges and tunnels sampled from bare earth. A bare-earth DEM records the ground, so a bridge deck takes the height of the river below it and a tunnel takes the height of the hill above. Both produce impossible gradients at the portals. Detect them by gradient magnitude and take the elevation from the connecting segments instead, or exclude structures from sampling entirely using the bridge and tunnel tags.
  • DEM voids are not sea level. A masked or no-data cell means “unknown”, and writing zero for it puts a junction at sea level in the middle of a plateau. Record the absence and let the gradient derivation skip those segments, which is visible, rather than inventing a height, which is not.
  • Vertical datum mismatches. Elevations may be relative to an ellipsoid or to a geoid, and the two differ by tens of metres — consistently, so a route’s total ascent is unaffected, but absolute heights and any threshold expressed in metres above sea level are wrong. Record the datum with the source.
  • Straight-line run instead of road length. Dividing rise by the distance between endpoints rather than by the segment’s own length overstates gradient on every bend, and dramatically on hairpins, which are exactly the steep roads where the number is consulted.
  • Resolution finer than the road network’s accuracy. A one-metre LiDAR DEM sampled at a node whose coordinate is accurate to five metres reads a height from the wrong side of the embankment. Smoothing over a short window along the road is more honest than sampling a single point.
  • Elevation changing under a cached projection. A GDS projection built before an elevation refresh routes on the old costs indefinitely, because a projection is a snapshot. Re-project after any enrichment pass.
A bare-earth DEM puts the bridge deck on the valley floor A viaduct crossing a valley, with the bare-earth elevation model beneath it. Sampling the deck's junctions against that model returns the ground height under each one, so the two mid-span junctions come out 46 metres below the abutments. The derived gradient is minus 31 per cent onto the bridge and plus 31 per cent off it — steeper than any road is built, on a structure that is level. Interpolating the deck elevations between the abutments instead recovers the true near-zero gradient, and the gradient magnitude check is what flags the segments needing it. A level viaduct, sampled from the ground beneath it bare-earth terrain viaduct deck — level sampled 46 m low abutment abutment derived from the sample −31% then +31% steeper than any road is built — the impossible-grade check catches it interpolated across the deck 0.2% throughout what the structure's own geometry says

Verification & Testing

Two checks catch nearly all of the failure modes above, and both run on the graph rather than on the raster.

import pytest

# Nothing paved exceeds about 35%; anything above that is a data artefact,
# nearly always a bridge or tunnel sampled from the ground beneath or above it.
IMPOSSIBLE_GRADE_PCT = 35.0


@pytest.mark.asyncio
async def test_no_impossible_gradients(session):
    result = await session.run(
        "MATCH ()-[s:SEGMENT]->() WHERE abs(s.grade_pct) > $limit "
        "RETURN count(s) AS n, collect(s.id)[..5] AS examples",
        limit=IMPOSSIBLE_GRADE_PCT,
    )
    record = await result.single()
    assert record["n"] == 0, (
        f"{record['n']} segments above {IMPOSSIBLE_GRADE_PCT}% — "
        f"check bridges and tunnels: {record['examples']}"
    )


@pytest.mark.asyncio
async def test_reverse_segments_have_opposite_rise(session):
    """A two-way road's two directions must disagree in sign and agree in size."""
    result = await session.run(
        """
        MATCH (a:Junction)-[f:SEGMENT]->(b:Junction)-[r:SEGMENT]->(a)
        WHERE f.rise_m IS NOT NULL AND r.rise_m IS NOT NULL
          AND abs(f.rise_m + r.rise_m) > 0.05
        RETURN count(*) AS mismatched
        """
    )
    assert (await result.single())["mismatched"] == 0

The second test is the one that catches a half-finished enrichment pass, where one direction was regraded and the other was not — a state in which routing quietly prefers whichever direction has the stale, cheaper cost.

FAQ

Which DEM resolution should I use?

Match it to the accuracy of the road geometry rather than maximising it. A 30-metre global model is adequate for long-distance vehicle routing, where gradients are averaged over kilometres. Cycling and accessibility routing want 5 metres or better, because a short steep ramp is exactly what a rider needs to know about and a coarse model averages it away. Sampling a 1-metre model at coordinates accurate to 5 metres buys precision you do not have.

Should elevation be on the node or the relationship?

Elevation on the node, gradient on the relationship. Elevation is a property of a place and is shared by every segment meeting there; gradient is a property of travelling between two places in a particular direction. Putting elevation on the relationship duplicates it once per incident edge and lets the copies drift.

Do I need two relationships for a two-way road?

If gradient affects cost, yes — or an equivalent arrangement where the cost function is evaluated with the direction of travel. A single undirected relationship with one grade_pct cannot represent a road that is a climb one way and a descent the other, which is every road that is not flat.

How do I stop bridges from corrupting the gradient?

Exclude structures from DEM sampling and interpolate their elevations from the junctions at either end, which is what a bridge deck actually does. The bridge, tunnel and layer tags from OSM ingestion carry the information needed to identify them.

Is total ascent enough to compare two routes?

For a rough comparison, yes, and it is the number riders quote. For a cost model, no — ascent treats a hundred metres gained at two per cent the same as a hundred metres gained at fourteen, and only one of those is rideable. The cost has to be non-linear in gradient, which is why ascent is a summary rather than a weight.

This topic is part of Spatial Graph Construction & OSM Ingestion.

Subtopics

Siblings