Resolving Way-Node References Without a Full Node Cache

Turn ways into geometries on a machine that cannot hold every node in memory, by holding only the nodes the ways actually reference.

Prerequisites

Conceptual minimum

A way stores node identifiers, not coordinates. Building its geometry means looking up every identifier, and the naive way to make those lookups fast is a dictionary of every node in the file — which for a country extract is tens of gigabytes.

The saving comes from a simple observation: you almost never need every node. A road-network build references the nodes on highway ways and nothing else, which is around twelve percent of a typical extract. Holding twelve percent of the nodes in a compact array instead of one hundred percent in a Python dictionary is two independent order-of-magnitude wins multiplied together.

A three-pass reference resolution that holds only the nodes it needs A four-stage chain. Pass one reads ways and collects the node identifiers they reference into a sorted int64 array. Pass two reads nodes and keeps only those whose identifier is in that set, tested by binary search. Pass three reads ways again and resolves each reference against the array in logarithmic time, emitting geometry. Memory is proportional to the nodes actually needed rather than to the whole file. Two passes beat one big dictionary pass 1: ways collect the node ids you need a sorted int64 array pass 2: nodes keep only those ids searchsorted membership pass 3: ways again resolve and emit geometry array lookup, O(log n) memory proportional to needed nodes not to the file A road-network extract references perhaps 12 percent of the nodes in the file. Holding only those is an order of magnitude cheaper than holding all of them.
The extra pass costs disk reads, which are cheap and parallelisable. The dictionary it replaces costs memory, which is neither.

The cost is an extra read of the file, and that is a good trade. Reads are sequential, cheap, and get faster with better storage; memory is a hard ceiling that fails abruptly.

Peak memory for five node-reference resolution strategies A bar chart for a 412 million node country extract where highway ways reference 48 million nodes. A Python dictionary of identifier to coordinate tuple needs 41 gigabytes at 104 bytes per node and does not fit. A dictionary holding all nodes more compactly still needs 5 gigabytes. Two sorted numpy arrays need 768 megabytes at 16 bytes per needed node. An osmium dense identifier array needs 3.3 gigabytes at 8 bytes per identifier slot across all identifiers. LMDB on disk needs 90 megabytes resident at about four microseconds per lookup. Memory for the same resolution, five ways 412 M-node country extract, highway ways only (48 M referenced nodes) dict[int, tuple] 41 GB · 104 B/node · does not fit dict, all nodes 5.0 GB · still every node in the file two sorted numpy arrays 768 MB · 16 B/node, needed only osmium dense id-array 3.3 GB · 8 B per id slot, all ids LMDB on disk 90 MB resident · ~4 µs/lookup The numpy option wins on memory because it exploits both facts at once: only the needed nodes, and no per-object Python overhead.
Dense id arrays are sized by the highest identifier in the file, not by how many you need — which is why they lose to a filtered array here.

Runnable solution

python 
#!/usr/bin/env python3
"""Build way geometries holding only the nodes those ways reference."""
from __future__ import annotations

import logging
from pathlib import Path
from typing import Callable, Iterator

import numpy as np
import osmium

logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
logger = logging.getLogger(__name__)

WayFilter = Callable[[osmium.osm.Way], bool]


class NeededNodes(osmium.SimpleHandler):
    """Pass 1 — collect the node ids referenced by the ways we care about."""

    def __init__(self, keep: WayFilter) -> None:
        super().__init__()
        self.keep = keep
        self._chunks: list[np.ndarray] = []
        self._buffer: list[int] = []

    def way(self, w) -> None:
        if not self.keep(w):
            return
        self._buffer.extend(n.ref for n in w.nodes)
        if len(self._buffer) >= 4_000_000:
            self._flush()

    def _flush(self) -> None:
        if self._buffer:
            self._chunks.append(np.fromiter(self._buffer, dtype=np.int64))
            self._buffer.clear()

    def ids(self) -> np.ndarray:
        """Sorted, deduplicated — sorting is what makes pass 3 a binary search."""
        self._flush()
        if not self._chunks:
            return np.empty(0, dtype=np.int64)
        ids = np.unique(np.concatenate(self._chunks))
        logger.info("pass 1: %d distinct node id(s) needed", ids.size)
        return ids


class NodeLocations(osmium.SimpleHandler):
    """Pass 2 — record coordinates for exactly those ids, in id order.

    Coordinates are kept as the raw int32 nanodegree-scaled values osmium exposes,
    which halves the memory against float64 and loses nothing: the scaling is exact.
    """

    def __init__(self, wanted: np.ndarray) -> None:
        super().__init__()
        self.wanted = wanted
        self.lon = np.zeros(wanted.size, dtype=np.int32)
        self.lat = np.zeros(wanted.size, dtype=np.int32)
        self.found = np.zeros(wanted.size, dtype=bool)

    def node(self, n) -> None:
        pos = np.searchsorted(self.wanted, n.id)
        if pos >= self.wanted.size or self.wanted[pos] != n.id:
            return                                    # not one we need
        loc = n.location
        if not loc.valid():
            return
        self.lon[pos] = loc.x                         # already in 1e-7 degree units
        self.lat[pos] = loc.y
        self.found[pos] = True

    def report(self) -> None:
        missing = int((~self.found).sum())
        logger.info("pass 2: %d/%d node(s) located%s", int(self.found.sum()),
                    self.wanted.size,
                    f" — {missing} missing (clipped extract?)" if missing else "")


class GeometryBuilder(osmium.SimpleHandler):
    """Pass 3 — resolve each way's references against the arrays and emit geometry."""

    SCALE = 1e-7

    def __init__(self, keep: WayFilter, ids: np.ndarray,
                 lon: np.ndarray, lat: np.ndarray, found: np.ndarray) -> None:
        super().__init__()
        self.keep, self.ids = keep, ids
        self.lon, self.lat, self.found = lon, lat, found
        self.complete = 0
        self.incomplete = 0
        self.results: list[tuple[int, np.ndarray]] = []

    def way(self, w) -> None:
        if not self.keep(w):
            return
        refs = np.fromiter((n.ref for n in w.nodes), dtype=np.int64)
        pos = np.searchsorted(self.ids, refs)
        # Guard the lookup: a clipped extract can reference nodes it does not contain.
        ok = (pos < self.ids.size) & (self.ids[np.minimum(pos, self.ids.size - 1)] == refs)
        ok &= self.found[np.minimum(pos, self.ids.size - 1)]
        if not ok.all():
            self.incomplete += 1
            return
        coords = np.stack((self.lon[pos] * self.SCALE, self.lat[pos] * self.SCALE), axis=1)
        self.results.append((w.id, coords))
        self.complete += 1


def build_geometries(path: Path, keep: WayFilter) -> list[tuple[int, np.ndarray]]:
    needed = NeededNodes(keep)
    needed.apply_file(str(path))
    ids = needed.ids()

    locations = NodeLocations(ids)
    locations.apply_file(str(path), locations=False)
    locations.report()

    builder = GeometryBuilder(keep, ids, locations.lon, locations.lat, locations.found)
    builder.apply_file(str(path))
    logger.info("pass 3: %d complete way(s), %d incomplete",
                builder.complete, builder.incomplete)
    return builder.results


if __name__ == "__main__":
    highways = lambda w: "highway" in w.tags
    geometries = build_geometries(Path("ireland.osm.pbf"), highways)
Which reference-resolution strategy suits which job A grid of four strategies. Two-pass numpy arrays are best when only a known subset of ways matters and are wasteful when every node is needed anyway. An osmium node cache is best when osmium is already in use and poor when identifiers are sparse and very large. LMDB or another on-disk key-value store is best when memory is the hard limit and poor when lookup latency dominates. PostGIS middle tables are best when the data lands in PostGIS anyway and wrong for a one-off file job. Which strategy fits which shape of job best when avoid when two-pass numpy arrays a known subset of ways you need every node anyway osmium node cache you already use osmium identifiers are sparse and huge LMDB / on-disk KV memory is the hard limit lookup latency dominates PostGIS middle tables the data lands in PostGIS anyway this is a one-off file job The deciding question is whether you need all nodes or a subset. Everything else follows from that.
Subset or everything: that single question eliminates two of the four options before any benchmarking.

Step-by-step walkthrough

NeededNodes buffers into Python lists and flushes into numpy arrays periodically. Appending to a list is fast and appending to a numpy array is not, so the chunked pattern gets the speed of one and the memory of the other. The final np.unique sorts and deduplicates in one operation, and the sort is not incidental — it is what makes every later lookup a binary search.

NodeLocations stores osmium’s raw loc.x and loc.y, which are already integers in units of 100 nanodegrees. Keeping them as int32 rather than converting to float64 halves the memory and loses no precision, because the conversion is an exact multiplication applied only at the end — the encoding described in PBF File Structure Deep Dive.

Note locations=False on the second apply_file. Asking pyosmium to manage locations itself would allocate exactly the node cache this whole approach exists to avoid.

GeometryBuilder guards the lookup twice: once that the identifier is present in the sorted array, and once that a coordinate was actually found for it. Both matter on a clipped extract, where a way can legitimately reference a node the file does not contain — the referential question at the heart of Extract Clipping & Boundary Polygons. Without the guards, searchsorted returns a plausible neighbouring index and the way gets a coordinate belonging to a different node, silently.

Verification

The incomplete-way count is the number to watch:

python 
ratio = builder.incomplete / (builder.complete + builder.incomplete)
logger.info("%.2f%% of ways could not be resolved", 100 * ratio)

On an extract cut with complete_ways or smart this should be zero. Anything above zero means the extract has dangling references, which is a property of how it was cut, not a bug in this code. On a simple-strategy extract a few percent is normal and is the reason that strategy is discouraged.

Spot-check a geometry against an independent source:

python 
way_id, coords = geometries[0]
print(way_id, coords[:3])
# Compare against: https://www.openstreetmap.org/way/<way_id>

And confirm the memory story actually held:

bash 
/usr/bin/time -v python3 build_geometries.py 2>&1 | grep 'Maximum resident'

For a country extract filtered to highways, expect under a gigabyte. Several gigabytes means the filter is matching far more ways than intended — check it before blaming the approach.

Common errors and fixes

Symptom Root cause Fix
Memory still in the tens of GB locations=True on a pass Pass locations=False; manage them yourself
Coordinates land in the wrong place searchsorted hit unguarded Verify ids[pos] == ref before using pos
Many incomplete ways Extract cut with simple Re-cut with complete_ways or smart
Pass 1 is slow Appending to a numpy array per node Buffer in a list, flush in chunks
Coordinates are all zero Node had no valid location Check loc.valid(); track a found mask
Third pass slower than expected searchsorted called per reference Call it once per way, on the whole array

Frequently Asked Questions

Is three passes really faster than one?

Faster in wall-clock on any machine where the one-pass version would swap, and slower on one where it would not. The point is not raw speed but the ceiling: the three-pass version’s memory is set by how many nodes your ways reference, which you control by filtering, while the one-pass version’s memory is set by the file. On a country extract with a highway filter the three-pass version runs comfortably in under a gigabyte; the dictionary version does not run at all.

Why int32 for coordinates rather than float64?

Because the source data is integers. OSM coordinates are stored as scaled integers with seven decimal places of precision, which fits int32 with room to spare, and converting to float only at the point of use loses nothing while halving the array. For 48 million nodes that is 384 megabytes saved for no cost at all.

When should I reach for osmium's own node cache instead?

When you need most of the nodes anyway, or when the pipeline is already built around osmium’s location handlers. The dense id array is sized by the highest node identifier in the file rather than by how many you need, which is why it loses here — but if you genuinely need every node, that overhead disappears and its constant-time lookup beats a binary search.

Can the passes be parallelised?

Pass 2 and pass 3 can, because PBF blocks are independently decodable, as covered in Speed Up OSM Parsing with Multiprocessing in Python. Share the sorted id array read-only across workers — it is a numpy array, so it can be memory-mapped rather than pickled, and each worker writes into a disjoint slice of the coordinate arrays.

Specification reference

pyosmium exposes node coordinates through Node.location, whose x and y attributes are the raw values in units of 100 nanodegrees (1e-7 degrees), and lon/lat as converted floats. apply_file(path, locations=False) disables the built-in location cache, leaving reference resolution to the caller.

Up one level: Node, Way & Relation Data Model.