Reading elevation along a few hundred thousand edge geometries is a problem where the naive implementation is not slightly slow but catastrophically slow, and the fix has nothing to do with the raster library. This page covers the mechanics of doing it efficiently and correctly, as the sampling stage of elevation and terrain data integration within the OSM Graph Architecture & Network Modeling pipeline. The output feeds the aggregation described in computing edge grade from elevation profiles.

Two things dominate the outcome: the order in which points are read, and whether nodata is recognised as nodata. Everything else is detail.

When to use this approach

Use a batched, block-ordered sampling pass when:

  • The extract is regional or larger. Below a few thousand edges, ordering makes no measurable difference and a straightforward loop is easier to read.
  • Sampling runs on every graph rebuild. Anything in the rebuild critical path is worth optimising once; a one-off analysis is not.
  • The DEM spans multiple source tiles. Tile-boundary crossings are where unordered access hurts most, because each crossing may reopen a file.
  • You need the raw profile, not just a grade. Storing the profile lets the aggregation mode change later without re-reading the raster — which is worth a great deal when tuning cost models.

A single small city with a DEM that fits comfortably in page cache needs none of this.

It is worth being precise about why the naive version is slow, because the instinct is to blame the raster library. It is not the library. A modern DEM is stored as a grid of compressed blocks, and reading any single cell requires decompressing the whole block containing it. Reading points in edge order means walking the road network, which wanders across the raster in a pattern that has nothing to do with block layout — so the same block is decompressed, evicted, and decompressed again, sometimes dozens of times. Reading in block order decompresses each block exactly once. The values are identical; only the number of decompressions changes, and on a national extract that number differs by more than an order of magnitude.

Read order decides whether sampling is fast The same set of sample points is read twice. In edge order the reads jump across the raster, touching many blocks repeatedly. Sorted by raster block, each block is loaded once and all its points are read together. Same points, same values, very different read cost edge order — random access 6 blocks touched 14 times block order — sequential access 6 blocks touched 6 times The values are identical either way — this is purely about how many times each raster block is faulted in and evicted again. Scatter the results back to their edges afterwards; keeping an index array makes that a single vectorised assignment.

Implementation

The whole technique is: flatten every edge’s points into one array, remember which edge each came from, sort by block, read, then scatter back.

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


def sample_edges(dem_path: str, edge_points: dict[int, np.ndarray],
                 src_crs: str = "EPSG:4326") -> dict[int, np.ndarray]:
    """Sample elevation for many edges with block-ordered reads.

    edge_points maps edge_id -> (n, 2) array of lon/lat sample coordinates.
    """
    edge_ids = list(edge_points)
    counts = np.array([len(edge_points[e]) for e in edge_ids])
    flat = np.vstack([edge_points[e] for e in edge_ids])
    # Index recording which edge each flattened point belongs to
    owner = np.repeat(np.arange(len(edge_ids)), counts)

    with rasterio.open(dem_path) as dem:
        transformer = Transformer.from_crs(src_crs, dem.crs, always_xy=True)
        xs, ys = transformer.transform(flat[:, 0], flat[:, 1])

        # Convert to row/col so points can be ordered by raster block
        rows, cols = dem.index(xs, ys)
        rows = np.asarray(rows)
        cols = np.asarray(cols)
        bh, bw = dem.block_shapes[0]
        block_key = (rows // bh).astype(np.int64) * (1 << 20) + (cols // bw)
        order = np.argsort(block_key, kind="stable")

        vals = np.full(len(flat), np.nan)
        ordered = np.column_stack([xs[order], ys[order]])
        read = np.fromiter(
            (v[0] for v in dem.sample(ordered)), dtype="float64", count=len(ordered)
        )
        if dem.nodata is not None:
            read = np.where(read == dem.nodata, np.nan, read)
        vals[order] = read

    # Scatter back into per-edge profiles
    out, start = {}, 0
    for i, eid in enumerate(edge_ids):
        n = counts[i]
        out[eid] = vals[start:start + n]
        start += n
    return out

The stable sort matters: it preserves the original order within a block, so points from the same edge stay adjacent, which helps the page cache a second time. The block_key packing assumes fewer than a million block columns, which holds comfortably for any DEM that fits on a disk.

Reprojecting once for the whole array rather than per edge is the other significant saving. pyproj transformations are vectorised and the per-call overhead dominates for small inputs, so a hundred thousand two-point calls cost far more than one two-hundred-thousand-point call.

Key parameters and tuning

Parameter Recommended value Notes
DEM product Copernicus GLO-30 Void-filled, better high-latitude coverage, clear redistribution terms
Mosaic format COG, BLOCKSIZE=512 Larger blocks reduce block count but read more bytes per touch
Compression DEFLATE Elevation compresses well; the CPU cost is far below the I/O saved
Overviews none Every read here is full resolution; overviews are pure storage cost
Resampling nearest Bilinear only helps when the interval approaches the cell size
Batch size whole extract if memory allows A 200 M-point float64 array is 1.6 GB; chunk by partition above that
Sort stability kind="stable" Keeps same-edge points adjacent within a block
Nodata handling map to NaN immediately The single highest-value line in the whole routine

Integration points

Feeding grade aggregation. The per-edge profile arrays are exactly what computing edge grade from elevation profiles consumes. Persisting the profiles rather than only the aggregate lets the cost model be retuned without re-reading the raster, which on a national extract is the difference between a minute and an hour.

Handing off to correction. Edges tagged bridge or tunnel still need sampling, because the correction described in handling DEM voids and bridge-tunnel artefacts uses the endpoint samples to build its ramp. Skipping them at sample time removes the information the correction needs.

Scoping to changed geometry. Sampling is expensive enough to be worth restricting to edges whose geometry actually changed, using the same partition-scoped logic as incremental graph rebuilds without full reprocessing. Elevation does not change between rebuilds; only the geometry sampling it does.

Where the sampling time actually goes Sampling a national extract takes over four hours with a per-point loop and per-edge reprojection. Batching the reprojection, then batching the reads, then ordering by raster block reduces it to about eleven minutes, with block ordering contributing the largest single saving. National extract, 47 M sample points — cumulative effect of each change per-point loop 4 h 20 min + batched reprojection 2 h 26 min + batched reads 1 h 13 min + block ordering 11 min — the largest single saving Block ordering changes no values and adds one sort. It is almost always the highest-return change available in this stage. Beyond this point the run is I/O bound on the mosaic, so further gains come from faster storage rather than from code. A CRS mismatch produces flat terrain, not an error With matched coordinate systems the sample points fall inside the DEM footprint and return elevations. With a mismatch they land far outside it, every read returns nodata, and the resulting profile is indistinguishable from genuinely flat ground. Where the sample points land when the CRS is wrong DEM footprint matched CRS — real elevations degrees read as metres — points land here every read returns nodata grade computes to 0.0 % on every edge — looks like flat terrain Assert that at least one sample per edge is valid. A silent all-nodata result is the one failure mode that produces no symptom at all. Sampling a peak of known elevation catches it immediately and costs one test.

Validation checklist

  • Sample count matches expectation. For each edge, the profile length should equal ceil(length / interval) + 1. A mismatch means densification and sampling disagree about the interval.
  • Nodata never survives as a number. Assert that no returned value equals the raster’s declared nodata. One survivor produces a grade in the thousands of percent.
  • A known summit reads correctly. Sample a peak with a published elevation and confirm the value is within the DEM’s stated vertical error. This catches CRS mismatches that otherwise produce plausible but wrong numbers.
  • Block ordering changes nothing but time. Run once ordered and once unordered on a subset and assert the profiles are identical. If they differ, the scatter-back index is wrong.
  • Coverage is complete. Count edges whose profile is entirely NaN. A non-trivial count means the DEM footprint does not cover the extract, usually at a border.