Clipping an OSM Extract with a .poly Boundary

Cut a region out of a larger .osm.pbf using a real boundary rather than a bounding box, and end up with a file whose features are complete and whose extent you have actually checked.

Prerequisites

Conceptual minimum

Clipping asks one question of every node — is it inside the polygon — and then asks a much harder question about every way and relation that references nodes on both sides. The --strategy flag answers the second question, and the Extract Clipping & Boundary Polygons topic compares the four available answers in detail. For this procedure, use smart: it keeps relations intact, which is what makes multipolygon buildings and route relations survive the cut.

The boundary format is Osmosis .poly, a plain-text description of one or more rings. It is worth knowing precisely because it has no error detection at all.

The three parts of the Osmosis .poly boundary format Three panels. Structure: a name line, then one or more ring sections each opened by an identifier line, followed by coordinate pairs and closed by END, with a final END closing the file. Coordinates: longitude first then latitude, whitespace-separated, decimal or scientific notation, always WGS 84 degrees, with the first and last pair matching. Holes: a ring identifier prefixed with an exclamation mark marks an enclave subtracted from the outer ring; order does not matter and deeper nesting is undefined. The .poly format in full — there is not much of it Structure Line 1: a name, any text Then one or more ring sections Each: an id line, then coordinates Each ring closed by END File closed by a final END Coordinates Longitude first, then latitude Whitespace-separated, any amount Decimal or scientific notation WGS 84 degrees, always First and last pair should match Holes Ring id prefixed with ! e.g. `!2` for an enclave Subtracted from the outer ring Order does not matter Nesting deeper than one is undefined There is no CRS declaration and no version marker, so a file in the wrong axis order is indistinguishable from a valid one until you look at the output.
The format carries no CRS and no version, which is why a file written latitude-first is accepted, cut, and produces an empty extract without complaint.

Two properties of that format cause almost every problem people have with it. Coordinates are longitude first — the opposite of the order most people say aloud, and the same axis-order trap discussed in Coordinate Reference Systems in OSM. And there is no validation: a syntactically fine file describing a polygon in the wrong hemisphere cuts cleanly and produces nothing.

The five steps from a boundary to a verified regional extract A five-stage chain: obtain a boundary from an administrative relation or your own GeoJSON; simplify it at about a hundred metres tolerance, taking a typical 180 thousand vertices down to two thousand; write it as a .poly file with longitude then latitude, one pair per line, terminated by END twice; run osmium extract with the polygon and strategy flags in one pass; and verify with osmium fileinfo that the counts are non-zero and a bounding box is present. Five steps, and only one of them is the extract command get a boundary admin relation or GeoJSON from OSM or your own simplify ~100 m tolerance 180 k verts → 2 k write .poly lon lat, one pair per line END, END osmium extract --polygon --strategy one pass verify fileinfo counts non-zero, has a bbox The two steps most often skipped are simplification and verification, and they are the two that cause the slow runs and the silently empty outputs.
Simplification and verification are the two steps teams skip, and they are exactly the ones that prevent an hour-long run and a silently empty file.

Runnable solution

The whole procedure, from a GeoJSON boundary to a verified extract:

python 
#!/usr/bin/env python3
"""Clip a regional extract from a parent .osm.pbf using a simplified .poly boundary."""
from __future__ import annotations

import json
import logging
import subprocess
from pathlib import Path

from shapely.geometry import shape, MultiPolygon, Polygon
from shapely.ops import unary_union

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

SIMPLIFY_DEGREES = 0.001   # ~100 m at mid latitudes — plenty for a clip boundary


def load_boundary(geojson_path: Path) -> MultiPolygon:
    """Read a GeoJSON Feature/FeatureCollection into one simplified MultiPolygon."""
    data = json.loads(geojson_path.read_text())
    features = data["features"] if data.get("type") == "FeatureCollection" else [data]
    geoms = [shape(f["geometry"] if "geometry" in f else f) for f in features]
    merged = unary_union(geoms)
    simplified = merged.simplify(SIMPLIFY_DEGREES, preserve_topology=True)
    if isinstance(simplified, Polygon):
        simplified = MultiPolygon([simplified])
    before = sum(len(p.exterior.coords) for p in merged.geoms) if hasattr(merged, "geoms") else len(merged.exterior.coords)
    after = sum(len(p.exterior.coords) for p in simplified.geoms)
    logger.info("boundary simplified: %d → %d exterior vertices", before, after)
    return simplified


def write_poly(mp: MultiPolygon, name: str, out: Path) -> Path:
    """Write a MultiPolygon as an Osmosis .poly file — longitude first, then latitude."""
    lines: list[str] = [name]
    ring_id = 0
    for poly in mp.geoms:
        ring_id += 1
        lines.append(str(ring_id))
        for lon, lat in poly.exterior.coords:
            lines.append(f"   {lon:.7E}   {lat:.7E}")
        lines.append("END")
        for interior in poly.interiors:          # holes are marked with a leading !
            ring_id += 1
            lines.append(f"!{ring_id}")
            for lon, lat in interior.coords:
                lines.append(f"   {lon:.7E}   {lat:.7E}")
            lines.append("END")
    lines.append("END")
    out.write_text("\n".join(lines) + "\n")
    logger.info("wrote %s (%d rings)", out, ring_id)
    return out


def clip(parent: Path, poly: Path, output: Path, base_url: str | None = None) -> None:
    """Run osmium extract with the smart strategy, carrying replication metadata."""
    cmd = [
        "osmium", "extract",
        "--polygon", str(poly),
        "--strategy", "smart",
        "--overwrite",
        "-o", str(output),
    ]
    if base_url:
        cmd[-3:-3] = ["--output-header", f"osmosis_replication_base_url={base_url}"]
    cmd.append(str(parent))
    logger.info("running: %s", " ".join(cmd))
    subprocess.run(cmd, check=True)


def verify(path: Path, boundary: MultiPolygon, tolerance_deg: float = 0.5) -> None:
    """Assert the extract is non-empty and its bbox overlaps the boundary we asked for."""
    info = json.loads(subprocess.run(
        ["osmium", "fileinfo", "--extended", "--json", str(path)],
        capture_output=True, text=True, check=True,
    ).stdout)
    counts = info["data"]["count"]
    if counts["nodes"] == 0:
        raise ValueError(f"{path}: zero nodes — the boundary is probably inverted")
    bbox = info["data"]["bbox"]                       # left, bottom, right, top
    want = boundary.bounds                             # minx, miny, maxx, maxy
    for got, expected, label in zip(bbox, want, ("left", "bottom", "right", "top")):
        if abs(got - expected) > tolerance_deg + 1.0:  # +1° slack for complete_ways spill
            raise ValueError(f"{path}: {label} edge at {got:.3f}, expected near {expected:.3f}")
    logger.info("%s verified: %d nodes, %d ways, %d relations",
                path, counts["nodes"], counts["ways"], counts["relations"])


if __name__ == "__main__":
    boundary = load_boundary(Path("boundaries/ireland.geojson"))
    poly = write_poly(boundary, "ireland", Path("boundaries/ireland.poly"))
    clip(Path("europe-latest.osm.pbf"), poly, Path("extracts/ireland.osm.pbf"),
         base_url="https://planet.osm.org/replication/minute/")
    verify(Path("extracts/ireland.osm.pbf"), boundary)

Step-by-step walkthrough

load_boundary merges every feature into a single geometry and simplifies it. The tolerance of 0.001 degrees is roughly a hundred metres, which is far finer than a clipping boundary needs and still removes the great majority of vertices from an administrative relation. preserve_topology=True guarantees the simplified ring stays valid and does not self-intersect — the defect class covered in Detecting Self-Intersecting OSM Polygons with Shapely.

write_poly emits the format exactly. Note the coordinate order in the loop: Shapely stores coordinates as (x, y), which for geographic data is (lon, lat), and that is also what .poly wants — so this is one of the few places where no swap is needed. Interior rings are written with a ! prefix so they become holes rather than additional outer rings.

clip builds the command. The --output-header insertion is what keeps the extract catchable by the diff stream later; without it, the workflow in Catching Up a Stale OSM Extract with pyosmium has no anchor to start from.

verify is the step that turns a hopeful run into a checked one. It reads the counts and the output bounding box from osmium fileinfo and asserts both. The one-degree slack allows for the smart strategy legitimately spilling past the boundary where a long way crosses it.

Four verification steps and the failure each one catches A grid of four checks. osmium fileinfo with the extended flag proves the file parses and gives counts, catching an empty or truncated cut. Comparing the output bounding box against the requested boundary proves the cut landed where asked and catches inverted axis order. Comparing counts against the parent file proves a plausible share and catches a boundary that is an order of magnitude wrong. Rendering a sample proves features look complete and catches a simple-strategy shred. What each verification command actually proves proves catches osmium fileinfo -e the file parses; counts an empty or truncated cut compare bbox to the boundary the cut landed where asked inverted axis order count vs the parent a plausible share of the parent a boundary an order of magnitude wrong render a sample features look complete a simple-strategy shred The bbox comparison is the one to automate: it is a four-number assertion and it catches the failure that produces a valid, empty, exit-zero file.
Automate the bounding-box comparison. It is four numbers and it catches the one failure that exits zero with a valid file containing nothing.

Verification

Run the script and expect three log lines: a vertex reduction of at least an order of magnitude, a ring count matching the number of separate landmasses in your boundary, and non-zero counts at the end. Then check the numbers by hand once:

bash 
osmium fileinfo --extended extracts/ireland.osm.pbf | grep -E 'Bounding box|Number of'

The bounding box should sit within about a degree of your boundary’s bounds. Node counts on the order of tens of millions for a country and hundreds of thousands for a city are the right magnitude; a count in the hundreds means the boundary is wrong, not that the region is empty.

Common errors and fixes

Message or symptom Root cause Fix
Zero nodes from the verifier Latitude and longitude swapped in the .poly Write longitude first; the script above already does
Open file ... exists osmium refuses to overwrite Pass --overwrite, or remove the output first
Run takes over an hour on a small region Boundary has tens of thousands of vertices Simplify before writing the .poly
Unknown file format on the boundary .poly written with a missing final END Every ring ends with END, and the file ends with one more
Output far larger than expected Boundary rings not closed, so the polygon is degenerate Ensure the first and last coordinate pair match
Extract has no replication anchor Parent had none and no header was set Pass --output-header with the base URL

Specification reference

The .poly format places the polygon name on the first line, then one section per ring: a section identifier, one coordinate pair per line as longitude then latitude, and END. A section identifier prefixed with ! marks the ring as a hole. A final END closes the file. Coordinates are WGS 84 degrees; the format carries no coordinate-system declaration.

Frequently Asked Questions

Can I use a GeoJSON boundary directly instead of converting to .poly?

Yes, on osmium-tool 1.11 and later — pass the GeoJSON file to --polygon and it is detected by extension. The conversion step above exists for two other reasons. It gives you a place to simplify, which is where the runtime saving comes from, and .poly is what most other OSM tooling accepts, so keeping the canonical boundary in that format avoids maintaining two. If your boundaries are generated fresh each run from a GIS source, skipping the conversion and simplifying in memory is perfectly reasonable.

How much should I simplify the boundary?

Enough that the vertex count drops by an order of magnitude, and no more. A tolerance of about a hundred metres takes a typical administrative relation from six figures of vertices to low thousands and moves the boundary by less than the width of the roads that cross it. Going coarser starts to cut across features you meant to include; going finer buys accuracy that a clipping operation cannot use, because the strategy already spills past the line wherever a way crosses it.

Why does the output bounding box extend past my boundary?

Because smart and complete_ways keep entire ways when any node of the way is inside. A motorway entering the region at its edge brings all of its nodes, including the ones far outside, and the header bounding box is computed from the nodes actually present. This is correct behaviour and the reason the verifier above allows a degree of slack. If you need output strictly bounded, clip the geometries in your own processing after loading them.

Should the boundary come from OSM or from an official source?

For clipping, from whichever is more stable. An OSM administrative relation is free and current but can be edited between runs, which means an extract cut last month and one cut today may not cover the same area — a difference that shows up as unexplained row-count drift. Pinning a boundary file in version control, whatever its origin, makes the cut reproducible, which usually matters more than the boundary being canonical.

Up one level: Extract Clipping & Boundary Polygons.