Understanding OSM multipolygon relations for GIS Jump to heading#

Take an OpenStreetMap type=multipolygon relation and reconstruct it into an OGC-valid polygon — building exterior rings from outer members, subtracting inner members as holes, and repairing winding order — before it reaches PostGIS, so a lake island or a country enclave is never silently swallowed or inverted.

Prerequisites Jump to heading#

Confirm each item before running the code below; an unmet prerequisite is the usual cause of a relation that “parses” but produces geometry that fails ST_IsValid only after ingestion.

• pyosmium ≥ 3.6.0 installed (pip install "osmium>=3.6") — it wraps libosmium and resolves node locations when apply_file(..., locations=True) is set.
• Shapely ≥ 2.0 installed (pip install "shapely>=2.0") for LinearRing, polygonize, make_valid, and the vectorized GEOS backend.
• pyproj ≥ 3.4 installed for the metric area sanity check (pip install "pyproj>=3.4").
• An .osm.pbf extract whose complete_ways are present at clip edges, so every outer/inner member fully resolves (re-clip with osmium extract --strategy=complete_ways if not).
• A target projected CRS chosen for the area check (EPSG:3857 for a quick metric proxy; EPSG:3035 for European equal-area work).
• A node-location index strategy selected (flex_mem for country extracts; sparse_file_array on an NVMe scratch volume for planet-scale runs).

Conceptual minimum Jump to heading#

OpenStreetMap encodes complex areal features through multipolygon relations: a structural primitive that aggregates several linear way members into one topological area. Each member carries a role of outer or inner, and that role — not the geometric winding direction — is the authoritative signal for how the ring participates. As the parent Node-Way-Relation Data Model explains, OSM does not guarantee ring orientation, so inferring exterior versus hole from segment direction is unreliable; you must trust the role, build outer rings first, then subtract every inner ring that falls inside them. The assembled signed area is the sum of outer areas minus the sum of inner areas:

$$\text{Area} = \sum_{o \in \text{outer}} |A_o| \; - \sum_{i \in \text{inner}} |A_i|$$

Topological validity adds three hard constraints: rings must be closed and non-self-intersecting, they may share nodes only at explicit boundary intersections, and an inner ring must lie wholly inside exactly one outer ring. Overlapping interiors or an unclosed ring violate the OGC Simple Features specification and abort geometry construction in standard GIS engines. Tag authority follows the same role discipline — the relation’s key-value pairs are canonical, and member-way tags apply only when a way is used standalone — so a strict key allowlist drawn from Tag Taxonomy & Key-Value Standards prevents inner-ring attributes from bleeding onto the assembled feature.

Runnable solution Jump to heading#

This two-pass pyosmium handler collects way coordinates, then assembles each type=multipolygon relation: it validates member resolution, builds rings by role, repairs validity with make_valid, runs a projected area sanity check, enforces canonical winding, and routes defects to a log instead of crashing the stream. It targets pyosmium>=3.6.0 and Shapely>=2.0.

import logging
import osmium
import shapely.geometry as geom
from shapely.geometry.polygon import orient
from shapely.validation import make_valid
from shapely.ops import transform as shp_transform, polygonize, unary_union
from pyproj import Transformer

logger = logging.getLogger("osm.multipolygon")

# Build ONE transformer per process; the constructor queries the PROJ
# operation database, so never rebuild it inside the relation loop.
_transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)


def _project(g):
    return shp_transform(_transformer.transform, g)


class MultipolygonETL(osmium.SimpleHandler):
    """Two-pass multipolygon handler: collect way coordinates, then assemble relations.

    Call ``apply_file(path, locations=True)`` so pyosmium resolves node
    coordinates and exposes them on each NodeRef in ``w.nodes``.
    """

    def __init__(self) -> None:
        super().__init__()
        self.way_coords: dict[int, list[tuple[float, float]]] = {}
        self.defect_log: list[str] = []
        self.relation_count = 0

    def way(self, w) -> None:
        # With locations=True each NodeRef carries a valid .location.
        coords = [
            (nr.location.lon, nr.location.lat)
            for nr in w.nodes
            if nr.location.valid()
        ]
        if coords:
            self.way_coords[w.id] = coords

    def relation(self, r) -> None:
        if r.tags.get("type") != "multipolygon":
            return

        self.relation_count += 1
        outer_rings: list[geom.LinearRing] = []
        inner_rings: list[geom.LinearRing] = []
        missing_ways: list[int] = []

        for member in r.members:
            if member.type != "w":
                continue  # multipolygon members are ways; skip stray node/relation refs
            coords = self.way_coords.get(member.ref)
            if not coords:
                missing_ways.append(member.ref)
                continue
            if len(coords) < 3:
                self.defect_log.append(
                    f"Relation {r.id}: way {member.ref} has fewer than 3 nodes"
                )
                continue

            ring = geom.LinearRing(coords)
            if member.role == "outer":
                outer_rings.append(ring)
            elif member.role == "inner":
                inner_rings.append(ring)

        if missing_ways:
            self.defect_log.append(f"Relation {r.id}: unresolved ways {missing_ways}")
            return
        if not outer_rings:
            self.defect_log.append(f"Relation {r.id}: no outer rings defined")
            return

        try:
            # One outer ring: pair it with all inners directly.
            # Multiple outer rings: let polygonize associate holes to the right shell.
            if len(outer_rings) == 1:
                poly = geom.Polygon(outer_rings[0], inner_rings)
            else:
                polys = list(polygonize(outer_rings + inner_rings))
                poly = unary_union(polys) if polys else geom.Polygon()

            valid_poly = make_valid(poly)
            if valid_poly.is_empty:
                self.defect_log.append(f"Relation {r.id}: empty geometry after validation")
                return

            # Project to a metric CRS for an area sanity check (< 1 m² is degenerate).
            if _project(valid_poly).area < 1.0:
                self.defect_log.append(f"Relation {r.id}: degenerate projected area")
                return

            # Enforce canonical winding (CCW outer / CW inner) before ingestion.
            if isinstance(valid_poly, geom.Polygon):
                valid_poly = orient(valid_poly, sign=1.0)

            self._ingest_to_postgis(r.id, valid_poly, dict(r.tags))

        except Exception as exc:  # GEOS topology failures surface here
            self.defect_log.append(f"Relation {r.id}: geometry construction failed: {exc}")

    def _ingest_to_postgis(self, rel_id: int, poly, tags: dict) -> None:
        # Placeholder for production DB insertion with ST_GeomFromWKB.
        pass


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
    handler = MultipolygonETL()
    handler.apply_file("extract.osm.pbf", locations=True, idx="flex_mem")
    logger.info("processed %d multipolygon relations", handler.relation_count)
    for defect in handler.defect_log[:20]:
        logger.warning("DEFECT %s", defect)

Step-by-step walkthrough Jump to heading#

1. Transformer caching — _transformer is built once at module scope. Its constructor queries the PROJ database, so rebuilding it per relation would dominate runtime; always_xy=True pins argument order to (lon, lat) to match OSM’s storage.
2. First pass (way) — with locations=True, every NodeRef in w.nodes carries a resolved .location, so the handler caches (lon, lat) arrays keyed by way id for later assembly. For planet runs, swap this dict for an on-disk LMDB/SQLite store.
3. Relation filter — only relations tagged type=multipolygon proceed; everything else returns immediately so the stream stays cheap.
4. Role-driven sorting — each way member is looked up; resolved rings are partitioned into outer_rings and inner_rings strictly by their role attribute, never by winding direction.
5. Reference closure — any unresolved member (missing_ways) aborts that relation with a logged defect rather than emitting a partial, deceptively valid shape — the same closure discipline covered in Error Handling in Large OSM Extracts.
6. Assembly — a single outer ring is combined with all inners via Polygon(shell, holes); multiple outer rings are handed to shapely.ops.polygonize, which associates each hole with its containing shell, then unary_union merges the parts into a MultiPolygon.
7. Validity repair — make_valid resolves self-touches and bowties into an OGC-valid geometry; an empty result is logged and dropped.
8. Area sanity check — the geometry is projected to EPSG:3857 and rejected if its area is under 1 m², catching collapsed or degenerate rings before they reach the database.
9. Winding repair — orient(poly, sign=1.0) rewrites the shell counter-clockwise and holes clockwise, the canonical order PostGIS and most renderers expect.

Verification Jump to heading#

Confirm the reconstruction is correct before wiring it into the next stage:

• Relation count. The processed N multipolygon relations log line should match osmium fileinfo --extended extract.osm.pbf relation counts filtered to type=multipolygon.
• Defect ratio. A healthy country extract logs defects for well under 1% of relations; a spike points to a clip that dropped complete_ways at the boundary.
• Hole presence. For a known feature with a hole (a lake with an island, an enclave), assert poly.interiors is non-empty — a missing interior means an inner member was misclassified or unresolved.
• Validity gate. After ingestion, SELECT count(*) FROM features WHERE NOT ST_IsValid(geom) must return 0; a non-zero count means a ring slipped past make_valid.
• Winding. ST_IsPolygonCCW(ST_ExteriorRing(geom)) should be true for every shell after orient(..., sign=1.0).

Common errors and fixes Jump to heading#

Specification reference Jump to heading#

A type=multipolygon relation builds its area from outer and inner member roles, not from way winding order; an inner ring must be fully contained within a single outer ring, and rings must be closed and non-self-intersecting. See the OSM Wiki on Relation:multipolygon for role rules and the OGC Simple Features access spec for the validity constraints GEOS enforces. OSM stores all coordinates in WGS 84 (EPSG:4326); reconstruct rings first, then project — the PBF File Structure Deep Dive covers the granularity/offset decode that yields those coordinates.

For loading into the database, osm2pgsql --slim --flat-nodes --hstore preserves multipolygon tags; follow it with ST_MakeValid and ST_CollectionExtract(geom, 3) to guarantee polygon output, then re-run ST_IsValid after each diff merge so community edits cannot quietly reintroduce topology regressions. Projected geometry from this stage feeds metric work most often through Spatial Indexing for OSM Extracts.

Related Jump to heading#

• Node-Way-Relation Data Model — how members and roles resolve into geometry across all relation types.
• Coordinate Reference Systems in OSM — why you reconstruct in WGS 84 and project only afterward.
• Spatial Indexing for OSM Extracts — R-tree and H3 structures for the assembled polygons.
• Tag Taxonomy & Key-Value Standards — the key allowlist that prevents inner-ring tag bleed.
• Error Handling in Large OSM Extracts — quarantine patterns for the relations this handler rejects.
• PBF File Structure Deep Dive — the decode that delivers node coordinates to your handler.

Up one level: Node-Way-Relation Data Model.