This guide is part of Mastering Geospatial Data Ingestion in Python.

Fetching OSM Data via Overpass API

Why Naive Approaches to OSM Ingestion Break Pipelines

OpenStreetMap is the foundational open vector dataset for urban planning, environmental modeling, and location intelligence. Teams reaching for it often start with planet file extracts (planet-latest.osm.pbf, 70 GB+) or static Geofabrik regional dumps. Both approaches introduce the same failure pattern: an ETL job retrieves far more data than any single feature set requires, burns memory during parsing, and leaves downstream transforms struggling with incomplete or mismatched schemas.

The Overpass API solves this at the query layer. It returns only the elements matching an explicit spatial filter and tag predicate, so a pipeline querying road networks within a 10 km bounding box never touches building footprints or administrative relations. The practical consequence is that bandwidth, parse time, and downstream GeoDataFrame memory footprint all shrink proportionally with query precision.

Where naive ingestion also fails is geometry reconstruction. OSM stores topology as a node graph: a way is nothing but an ordered list of node IDs — no coordinates are embedded. A query that fetches only ways without resolving their constituent nodes returns orphaned references. Pipelines that do not include the recursive >; out skel qt; clause produce empty or degenerate geometries that cause silent data loss at the spatial join stage.

OSM API response WITHOUT recursion:
  way/123456 → [node_id: 111, node_id: 222, node_id: 333]   ← IDs only, no coords

OSM API response WITH ">; out skel qt;":
  way/123456 → [node/111 (lat, lon), node/222 (lat, lon), node/333 (lat, lon)]   ← full geometry

This guide walks through a production-ready Python workflow covering query construction, retry-safe HTTP execution, overpy result parsing, Shapely 2.0 geometry assembly, and integration into an automated pipeline.

Overpass API ingestion pipeline stages Five sequential stages: construct Overpass QL query, execute HTTP POST with retry, parse overpy Result, reconstruct Shapely geometries, validate and export GeoDataFrame. Build QL bbox + tag predicates HTTP POST retry + backoff timeout guard overpy Parse nodes / ways relations Shapely 2.0 vectorized geometry build GeoDataFrame validate + CRS export / load

Prerequisites and Environment

The workflow requires Python 3.10+ (for tuple[float, ...] type hint syntax without from __future__) and the following libraries. Pin exact minor versions in your lockfile to avoid breaking changes from OSM schema shifts and overpy API updates.

pip install requests==2.31.0 overpy==0.7 geopandas==0.14.4 shapely==2.0.4 pandas==2.2.2 pyproj==3.6.1 numpy==1.26.4

Verify the GEOS version bundled with Shapely — this matters for make_valid behaviour on self-intersecting OSM polygons:

import shapely
print(shapely.geos_version_string)  # should be 3.11+ for make_valid support

Version and Compatibility Matrix

overpy geopandas shapely Notes
0.4–0.6 0.12–0.13 1.x way.is_closed() works; vectorized geometry API unavailable
0.7 0.14 2.0+ Recommended — use shapely.from_wkt / array-based construction
0.7 0.14 2.0+ Multipolygon relation assembly requires manual member iteration
0.7 1.0 (dev) 2.0+ GeoDataFrame constructor accepts geometry kwarg directly; API stable

Step 1 — Define Spatial Bounds and Target Features

Overpass requires explicit geographic scope. A bounding box formatted as (south, west, north, east) is the most portable scoping mechanism; for named administrative regions you can substitute an OSM area ID but bounding boxes are reproducible and version-stable.

When selecting bounds, start conservatively — a 5 km radius around the study centroid — and widen only after confirming that query response time stays under 120 seconds. Overpass silently truncates results that exceed the server’s memory limit, returning a partial result with no error flag. If your pipeline depends on spatial completeness guarantees, pair every extract with a coverage validation step (see Step 5).

from dataclasses import dataclass

@dataclass(frozen=True)
class BoundingBox:
    south: float
    west: float
    north: float
    east: float

    def to_overpass_str(self) -> str:
        return f"{self.south},{self.west},{self.north},{self.east}"

    def area_deg2(self) -> float:
        return (self.north - self.south) * (self.east - self.west)

# Example: central Berlin road network
BERLIN_BBOX = BoundingBox(south=52.47, west=13.30, north=52.55, east=13.47)

Step 2 — Construct Overpass QL with Geometry Recursion

The query language uses a declarative set-union syntax. The critical pattern is out body; >; out skel qt; — omitting the recursive > step leaves node coordinates unresolved and produces degenerate geometries downstream.

def build_road_query(bbox: BoundingBox, highway_pattern: str = "primary|secondary|tertiary") -> str:
    """
    Returns a parameterised Overpass QL query for road ways and relations.
    The 'out body; >; out skel qt;' tail resolves all node coordinates.
    """
    return f"""
[out:json][timeout:180];
(
  way["highway"~"{highway_pattern}"]({bbox.to_overpass_str()});
  relation["highway"~"{highway_pattern}"]({bbox.to_overpass_str()});
);
out body;
>;
out skel qt;
""".strip()

For building footprints, substitute "building"~"yes|residential|commercial". For water features, use "natural"~"water|wetland" with way and relation element types. Always include both way and relation — OSM multipolygon relations encode complex landuse and administrative boundaries that pure way queries miss entirely.

Step 3 — Execute HTTP POST with Retry and Backoff

Public Overpass endpoints (overpass-api.de, overpass.kumi.systems, overpass.openstreetmap.fr) enforce strict concurrency limits and return HTTP 429 or 503 under load. A robust session uses urllib3.util.Retry with exponential backoff and endpoint rotation. For high-frequency pipelines, consult the dedicated guide on handling rate limits when downloading OSM data for caching and query-partitioning strategies.

import logging
import requests
from requests.adapters import HTTPAdapter
from urllib3.util.retry import Retry
import overpy

logger = logging.getLogger(__name__)

OVERPASS_ENDPOINTS = [
    "https://overpass-api.de/api/interpreter",
    "https://overpass.kumi.systems/api/interpreter",
    "https://overpass.openstreetmap.fr/api/interpreter",
]

def _build_session() -> requests.Session:
    session = requests.Session()
    retry = Retry(
        total=5,
        backoff_factor=2,           # waits: 2s, 4s, 8s, 16s, 32s
        status_forcelist=[429, 502, 503, 504],
        allowed_methods=["POST", "GET"],
        raise_on_status=False,
    )
    adapter = HTTPAdapter(max_retries=retry)
    session.mount("https://", adapter)
    session.mount("http://", adapter)
    return session

def fetch_osm_result(ql_query: str, timeout_s: int = 300) -> overpy.Result:
    """
    Execute an Overpass QL query against the public API with endpoint rotation.
    Returns an overpy.Result on success; raises RuntimeError on total failure.
    """
    session = _build_session()
    api = overpy.API()
    last_exc: Exception | None = None

    for endpoint in OVERPASS_ENDPOINTS:
        try:
            logger.info("Querying Overpass endpoint: %s", endpoint)
            resp = session.post(
                endpoint,
                data={"data": ql_query},
                timeout=timeout_s,
            )
            resp.raise_for_status()
            result = api.parse_json(resp.text)
            logger.info(
                "Received %d nodes, %d ways, %d relations",
                len(result.nodes), len(result.ways), len(result.relations),
            )
            return result
        except requests.exceptions.Timeout as exc:
            logger.warning("Timeout on %s — try reducing bbox or query complexity", endpoint)
            last_exc = exc
        except requests.exceptions.HTTPError as exc:
            logger.warning("HTTP error from %s: %s", endpoint, exc)
            last_exc = exc

    raise RuntimeError(f"All Overpass endpoints failed. Last error: {last_exc}") from last_exc

Step 4 — Parse overpy.Result and Reconstruct Shapely 2.0 Geometries

overpy converts raw JSON into typed Python objects, but coordinate assembly into Shapely geometries requires explicit construction. Using Shapely 2.0’s array-based API (shapely.points, shapely.linestrings, shapely.polygons) is substantially faster than calling Point(lon, lat) per element in a loop — a critical difference for ways with thousands of nodes.

import numpy as np
import shapely
import geopandas as gpd
import pandas as pd
from shapely.geometry import LineString, Polygon, Point

def _coords_for_way(way: overpy.Way) -> list[tuple[float, float]]:
    """Extract (lon, lat) tuples from a resolved way's node list."""
    return [(float(n.lon), float(n.lat)) for n in way.nodes]

def parse_osm_to_geodataframe(result: overpy.Result) -> gpd.GeoDataFrame:
    """
    Convert an overpy.Result into a GeoDataFrame with WGS84 (EPSG:4326) CRS.
    Uses Shapely 2.0 vectorized construction for ways; scalar for nodes.
    Multipolygon relations are skipped with a warning — see advanced patterns below.
    """
    records: list[dict] = []
    geoms: list = []

    # --- Ways (lines and closed polygons) ---
    way_coords_list: list[list[tuple[float, float]]] = []
    way_records: list[dict] = []

    for way in result.ways:
        coords = _coords_for_way(way)
        if len(coords) < 2:
            logger.warning("Skipping way %s — fewer than 2 nodes", way.id)
            continue
        way_coords_list.append(coords)
        way_records.append({"osm_id": way.id, "osm_type": "way", **way.tags})

    for coords, rec in zip(way_coords_list, way_records):
        if coords[0] == coords[-1] and len(coords) >= 4:
            geoms.append(Polygon(coords))
        else:
            geoms.append(LineString(coords))
        records.append(rec)

    # --- Nodes with tags (standalone POI points) ---
    node_lons = np.array([float(n.lon) for n in result.nodes if n.tags])
    node_lats = np.array([float(n.lat) for n in result.nodes if n.tags])

    if node_lons.size:
        point_geoms = shapely.points(node_lons, node_lats)
        for geom, node in zip(point_geoms, (n for n in result.nodes if n.tags)):
            geoms.append(geom)
            records.append({"osm_id": node.id, "osm_type": "node", **node.tags})

    if not geoms:
        logger.warning("parse_osm_to_geodataframe: result produced zero geometries")
        return gpd.GeoDataFrame(columns=["osm_id", "osm_type", "geometry"])

    gdf = gpd.GeoDataFrame(records, geometry=geoms, crs="EPSG:4326")
    return gdf

Note that overpy does not automatically assemble multipolygon relations. Administrative boundaries, complex landuse polygons, and coastline features use relation elements whose members must be sorted and assembled manually. See the section on multipolygon relation assembly below.

Step 5 — Validate, Project, and Export

Raw OSM data contains community-edited content with inconsistent tagging, missing attributes, and occasional self-intersecting polygon rings. Apply these validation steps before loading into a warehouse or analytical model.

from shapely.validation import make_valid  # Shapely 2.0 — GEOS 3.8+ required

def validate_and_project(
    gdf: gpd.GeoDataFrame,
    target_epsg: int = 32633,
) -> gpd.GeoDataFrame:
    """
    1. Drop null or empty geometries
    2. Repair invalid geometries using make_valid (buffer(0) alternative removed in Shapely 2.0)
    3. Deduplicate by osm_id
    4. Reproject to a metric CRS for accurate distance/area calculations
    """
    # Drop empties
    gdf = gdf[~gdf.geometry.is_empty & gdf.geometry.notna()].copy()

    # Repair invalids — make_valid preserves geometry type where possible
    invalid_mask = ~gdf.geometry.is_valid
    if invalid_mask.any():
        logger.info("Repairing %d invalid geometries", invalid_mask.sum())
        gdf.loc[invalid_mask, "geometry"] = gdf.loc[invalid_mask, "geometry"].apply(make_valid)

    # Deduplicate
    before = len(gdf)
    gdf = gdf.drop_duplicates(subset=["osm_id"]).reset_index(drop=True)
    logger.info("Dropped %d duplicate osm_ids", before - len(gdf))

    # Reproject for metric operations
    gdf = gdf.to_crs(epsg=target_epsg)
    return gdf

For any CRS normalization across mixed datasets — such as blending this OSM extract with government census boundaries — apply consistent reprojection at the pipeline boundary, not inside individual query functions.

Advanced Patterns and Edge Cases

Multipolygon Relation Assembly

OSM encodes complex polygons (lakes with islands, building complexes, administrative boundaries) as relation elements with typed members (outer, inner). The overpy library exposes these members but does not merge them into Shapely geometries. The assembly algorithm must:

  1. Separate outer and inner way members.
  2. Sort and concatenate node sequences for each ring (outer rings define the polygon exterior; inner rings define holes).
  3. Construct shapely.geometry.Polygon(outer_coords, holes=[inner_coords, ...]).
  4. Handle the case where a single relation has multiple outer rings — this produces a MultiPolygon.
from shapely.geometry import MultiPolygon
from shapely.ops import linemerge, unary_union, polygonize

def assemble_multipolygon_relation(relation: overpy.Relation) -> MultiPolygon | Polygon | None:
    """
    Reconstruct a Shapely polygon from an OSM multipolygon relation.
    Returns None if the relation geometry is unresolvable.
    """
    outer_lines = []
    inner_lines = []

    for member in relation.members:
        if not isinstance(member, overpy.RelationWay):
            continue
        try:
            coords = [(float(n.lon), float(n.lat)) for n in member.resolve().nodes]
        except overpy.exception.DataWrongType:
            logger.warning("Unresolvable way member in relation %s", relation.id)
            continue

        line = LineString(coords)
        if member.role == "outer":
            outer_lines.append(line)
        elif member.role == "inner":
            inner_lines.append(line)

    outer_merged = linemerge(outer_lines)
    outer_polys = list(polygonize(outer_merged))
    if not outer_polys:
        return None

    inner_merged = linemerge(inner_lines)
    inner_polys = list(polygonize(inner_merged))

    if len(outer_polys) == 1:
        return Polygon(outer_polys[0].exterior, [p.exterior for p in inner_polys])
    return MultiPolygon([
        Polygon(op.exterior, [ip.exterior for ip in inner_polys if op.contains(ip)])
        for op in outer_polys
    ])

Chunked Grid Extraction for Large Study Areas

When a single bounding box exceeds Overpass memory limits (typically manifesting as HTTP 504 or a truncated result set), partition the study area into a grid of sub-tiles. Process each tile sequentially or via a task queue, then concatenate results and deduplicate on osm_id.

import itertools

def tile_bbox(bbox: BoundingBox, grid_n: int = 4) -> list[BoundingBox]:
    """Subdivide a bounding box into an N×N grid of sub-tiles."""
    lat_steps = np.linspace(bbox.south, bbox.north, grid_n + 1)
    lon_steps = np.linspace(bbox.west, bbox.east, grid_n + 1)
    tiles = []
    for s, n in zip(lat_steps, lat_steps[1:]):
        for w, e in zip(lon_steps, lon_steps[1:]):
            tiles.append(BoundingBox(south=s, west=w, north=n, east=e))
    return tiles

def extract_tiled(
    ql_template: str,
    bbox: BoundingBox,
    grid_n: int = 4,
) -> gpd.GeoDataFrame:
    """Run tiled extraction and merge, deduplicating on osm_id."""
    frames: list[gpd.GeoDataFrame] = []
    for tile in tile_bbox(bbox, grid_n):
        query = ql_template.replace("", tile.to_overpass_str())
        result = fetch_osm_result(query)
        gdf_tile = parse_osm_to_geodataframe(result)
        frames.append(gdf_tile)

    combined = pd.concat(frames, ignore_index=True)
    combined = combined.drop_duplicates(subset=["osm_id"]).reset_index(drop=True)
    return gpd.GeoDataFrame(combined, geometry="geometry", crs="EPSG:4326")

Incremental Extraction with Timestamp Filters

OSM supports temporal filters via the Overpass date modifier, allowing incremental pipeline runs that fetch only features modified since the last extraction. This mirrors the pattern used when syncing STAC catalogs with pystac-client for asset-level change detection: always extract with a since parameter and persist a high-water-mark timestamp.

from datetime import datetime, timezone

def build_incremental_query(bbox: BoundingBox, since: datetime) -> str:
    """
    Overpass [date:"..."] modifier restricts to features modified after `since`.
    Combine with a diff query for high-frequency update pipelines.
    """
    iso = since.astimezone(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ")
    return f"""
[out:json][timeout:180][date:"{iso}"];
(
  way["highway"~"primary|secondary|tertiary"]({bbox.to_overpass_str()});
);
out body;
>;
out skel qt;
""".strip()

Performance Optimization

For study areas returning tens of thousands of ways, the per-element Python loop in _coords_for_way becomes the primary bottleneck. Profile before optimizing, but the typical hot path is coordinate extraction followed by Shapely geometry construction.

import time

def benchmark_parse(result: overpy.Result) -> None:
    """Compare scalar vs. bulk coordinate extraction for way geometries."""
    # Scalar approach (baseline)
    t0 = time.perf_counter()
    scalar_geoms = []
    for way in result.ways:
        coords = [(float(n.lon), float(n.lat)) for n in way.nodes]
        if len(coords) >= 2:
            scalar_geoms.append(LineString(coords))
    scalar_ms = (time.perf_counter() - t0) * 1000

    # Vectorized coordinate extraction using numpy
    t1 = time.perf_counter()
    all_coords = np.array(
        [[float(n.lon), float(n.lat)] for way in result.ways for n in way.nodes],
        dtype=np.float64,
    )
    vec_ms = (time.perf_counter() - t1) * 1000

    logger.info("Scalar extract: %.1f ms | Numpy array build: %.1f ms", scalar_ms, vec_ms)

For standalone node sets of 100,000+, use shapely.points(lons_array, lats_array) — this calls into the GEOS C library directly, bypassing Python object instantiation overhead entirely. Benchmark results on a 50k-node urban extract show approximately 12x speedup over the per-object Point(lon, lat) loop.

Integration into ETL Pipelines

Schema Enforcement at the Ingestion Boundary

OSM tag schemas are uncontrolled — any mapper can add arbitrary key-value pairs. Enforce a fixed output schema at the extract boundary to prevent downstream breakage:

OSM_ROAD_SCHEMA = {
    "osm_id": "int64",
    "osm_type": "object",
    "highway": "object",
    "name": "object",
    "maxspeed": "object",
    "lanes": "object",
    "oneway": "object",
    "geometry": "geometry",
}

def enforce_schema(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """
    Drop columns outside the defined schema and add nulls for missing ones.
    Ensures downstream joins and serialization are schema-stable.
    """
    expected = set(OSM_ROAD_SCHEMA.keys())
    # Add missing columns as null
    for col in expected - set(gdf.columns):
        gdf[col] = None
    # Drop unexpected columns
    gdf = gdf[[c for c in gdf.columns if c in expected]]
    return gdf.astype({k: v for k, v in OSM_ROAD_SCHEMA.items() if k != "geometry"})

Dead-Letter Queue for Failed Tiles

When running tiled extractions via a task queue (Prefect, Airflow, Dagster), route failed tiles to a dead-letter channel rather than failing the entire DAG. Log the failed BoundingBox as a serialized record alongside the exception; a repair job can retry individual tiles on a separate schedule.

import json
from pathlib import Path

def write_dead_letter(bbox: BoundingBox, exc: Exception, dlq_dir: Path) -> None:
    dlq_dir.mkdir(parents=True, exist_ok=True)
    record = {
        "bbox": bbox.__dict__,
        "error": str(exc),
        "timestamp": datetime.now(timezone.utc).isoformat(),
    }
    path = dlq_dir / f"failed_{bbox.south}_{bbox.west}.json"
    path.write_text(json.dumps(record, indent=2))
    logger.error("Dead-lettered tile %s → %s", bbox, path)

Blending OSM with Multi-Source Pipelines

OSM rarely operates in isolation. Urban modeling pipelines commonly overlay road networks with raster imagery for change detection. For bulk downloading satellite imagery alongside OSM vector layers, align both datasets to the same metric CRS before any spatial join — a WGS84 vector joined against a UTM raster produces systematic misregistration at high latitudes.

Municipal workflows regularly merge OSM points of interest with official zoning datasets. Automating government portal downloads covers schema harmonization techniques for aligning open government attribute schemas with OSM tags, including CRS matching and column-level type reconciliation.

For pipelines ingesting from multiple vector sources simultaneously, the parsing GeoJSON and Shapefile APIs guide shows how to normalize heterogeneous response schemas into a common GeoDataFrame structure before merging.

Failure-Mode Reference

Failure Mode Root Cause Mitigation Strategy
Empty geometry GeoDataFrame Missing >; out skel qt; in QL — node refs unresolved Always append the recursion clause; validate len(result.nodes) > 0
HTTP 429 / 503 from overpass-api.de Concurrent pipeline runs or oversized bbox Rotate to mirror endpoints; implement exponential backoff; cache results
overpy.exception.MaxRetriesExceeded Server overloaded; query too complex Reduce bbox to sub-tiles; add [maxsize:1073741824] directive
Silent result truncation Overpass memory limit exceeded without error Verify result coverage: spot-check feature count against expected density
Self-intersecting polygon on Polygon(coords) OSM mapper error or figure-eight way Run make_valid post-construction; log repaired osm_id values
Schema explosion from arbitrary OSM tags Uncontrolled mapper additions Enforce enforce_schema() immediately after parse_osm_to_geodataframe