Writing OSM Features to GeoParquet with PyArrow

Turn a stream of normalised OSM features into a GeoParquet file that GeoPandas, DuckDB and GDAL all open as spatial data, and that a bounding-box filter can read a fraction of.

Prerequisites

Conceptual minimum

GeoParquet is Parquet with a convention. There is no new file format, no new encoding, and no new library — just a geo key in the file’s metadata declaring where the geometry is and what coordinate system it is in.

What turns a Parquet file into a GeoParquet file A four-stage chain. An Arrow table carries typed columns plus a well-known-binary blob for geometry. A geo metadata key declares the version, the primary geometry column, the CRS and the geometry types present. write_table writes it with zstd compression, a row-group size and statistics enabled. Any spatial reader such as GeoPandas, DuckDB or GDAL then finds the geometry. Parquet plus one metadata key is the whole specification Arrow table typed columns + WKB blob geometry as bytes geo metadata version · primary_column crs · geometry_types write_table zstd · row groups statistics on any spatial reader GeoPandas · DuckDB · GDAL finds the geometry Omit the metadata key and the same file is a Parquet table with an unreadable blob column — valid, loadable, and spatially invisible.
The metadata key is the entire difference. Without it the file still loads, still contains every byte of geometry, and no spatial tool will recognise it.

Geometry itself is stored as well-known binary in an ordinary binary column. Everything a spatial reader needs to interpret that column lives in the metadata block, and every interoperability problem people have with GeoParquet is a problem with that block.

Four GeoParquet metadata mistakes and their effect on readers A grid of four mistakes. No geo key at all makes readers treat the file as a plain table with an unreadable blob column. A primary_column naming a missing column causes an error or a silent fallback and a no-geometry-column message. A null CRS on projected data makes readers assume CRS84 degrees so everything plots near zero, zero. Omitting geometry_types is read fine but loses a planning optimisation. Four metadata mistakes and what each reader does with them what readers do how it looks no geo key at all treat it as a plain table blob column, no geometry primary_column names a missing column error, or fall back to none "no geometry column found" crs null but data is projected assume CRS84 degrees everything plots near 0,0 geometry_types omitted read it, no optimisation works, slower planning Only the last is harmless. The third is the dangerous one, because the file loads, plots, and is wrong by the width of a continent.
A null CRS is not "unspecified" — the specification defines it as OGC:CRS84. Writing projected coordinates under it produces a file every reader misreads confidently.

The CRS field deserves particular care. A null CRS is not an omission — the specification defines it to mean OGC:CRS84, longitude then latitude in WGS 84, which is exactly what unprojected OSM data is. That makes null the correct value for a file straight out of an OSM pipeline, and a serious error for one that has been reprojected, which is the trap discussed in Coordinate Reference Systems in OSM.

Runnable solution

python 
#!/usr/bin/env python3
"""Write normalised OSM features as GeoParquet that spatial readers recognise."""
from __future__ import annotations

import json
import logging
from typing import Iterable, Sequence

import h3
import pyarrow as pa
import pyarrow.compute as pc
import pyarrow.parquet as pq
from shapely import to_wkb
from shapely.geometry.base import BaseGeometry

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

ROW_GROUP_ROWS = 250_000
SORT_RESOLUTION = 6          # coarse H3 cell, used only as a locality-preserving sort key


def geo_metadata(geometry_types: Sequence[str], bbox: Sequence[float]) -> bytes:
    """The `geo` key. A null crs means OGC:CRS84 — correct for unprojected OSM."""
    return json.dumps({
        "version": "1.1.0",
        "primary_column": "geometry",
        "columns": {
            "geometry": {
                "encoding": "WKB",
                "geometry_types": list(geometry_types),
                "crs": None,
                "bbox": list(bbox),
            }
        },
    }).encode()


def build_table(features: Iterable[dict]) -> pa.Table:
    """Typed columns for what we query, a map column for everything else."""
    osm_ids, names, kinds, geoms, extra, cells = [], [], [], [], [], []
    for f in features:
        geom: BaseGeometry = f["geometry"]
        tags: dict[str, str] = f["tags"]
        osm_ids.append(f["osm_id"])
        names.append(tags.get("name"))
        kinds.append(tags.get("building"))
        geoms.append(to_wkb(geom, output_dimension=2))
        # Everything not promoted to a column survives in the map, so nothing is lost.
        extra.append([(k, v) for k, v in tags.items() if k not in ("name", "building")])
        centroid = geom.representative_point()
        cells.append(h3.latlng_to_cell(centroid.y, centroid.x, SORT_RESOLUTION))
    return pa.table({
        "osm_id": pa.array(osm_ids, pa.int64()),
        "name": pa.array(names, pa.string()),
        "building": pa.array(kinds, pa.string()),
        "tags": pa.array(extra, pa.map_(pa.string(), pa.string())),
        "h3_r6": pa.array(cells, pa.string()),
        "geometry": pa.array(geoms, pa.binary()),
    })


def write_geoparquet(table: pa.Table, path: str, geometry_types: Sequence[str],
                     bbox: Sequence[float]) -> None:
    """Sort for locality, attach the geo metadata, write with statistics on."""
    ordered = table.sort_by([("h3_r6", "ascending"), ("osm_id", "ascending")])
    meta = dict(ordered.schema.metadata or {})
    meta[b"geo"] = geo_metadata(geometry_types, bbox)
    ordered = ordered.replace_schema_metadata(meta)
    pq.write_table(
        ordered, path,
        compression="zstd", compression_level=3,
        row_group_size=ROW_GROUP_ROWS,
        write_statistics=True,
    )
    groups = -(-ordered.num_rows // ROW_GROUP_ROWS)
    logger.info("wrote %s: %d rows in %d row group(s)", path, ordered.num_rows, groups)


def verify(path: str) -> None:
    """Assert the file is readable as spatial data before anything downstream sees it."""
    pf = pq.ParquetFile(path)
    meta = pf.schema_arrow.metadata or {}
    if b"geo" not in meta:
        raise ValueError(f"{path}: no `geo` metadata — readers will not see the geometry")
    geo = json.loads(meta[b"geo"])
    primary = geo["primary_column"]
    if primary not in pf.schema_arrow.names:
        raise ValueError(f"{path}: primary_column {primary!r} is not in the schema")
    logger.info("%s: %d row group(s), primary column %r, crs %r",
                path, pf.num_row_groups, primary, geo["columns"][primary]["crs"])

Step-by-step walkthrough

build_table implements the hybrid schema the parent topic recommends: name and building become real typed columns because they are what queries filter on, and everything else goes into a map<string,string> so the long tail of OSM tagging survives without inflating the schema to thousands of mostly-null columns.

The h3_r6 column is not there to be queried — it is a sort key. A coarse H3 cell groups geographically nearby features into adjacent rows, which is what makes the per-row-group min/max statistics selective. The same job can be done with a quadkey or a Hilbert index; the property that matters is locality preservation, discussed in Spatial Index Selection.

write_geoparquet sorts, then attaches metadata, then writes. The order matters: replace_schema_metadata returns a new table, and sorting after attaching would carry the metadata forward anyway, but writing before sorting would produce a file whose statistics are useless.

verify reads the file back and asserts the two things that actually break readers. It costs milliseconds and it catches the failure that otherwise appears as “this file has no geometry” in someone else’s notebook a week later.

Bytes read by a spatial filter at four row-group sizes A bar chart of bytes read for a bounding-box query selecting about 0.4 percent of 14.2 million sorted building features. A single row group reads all 980 megabytes in 2.10 seconds. One million rows per group reads 172 megabytes in 0.41 seconds. A quarter of a million rows per group reads 61 megabytes in 0.17 seconds. Fifty thousand rows per group reads 38 megabytes in 0.19 seconds but adds seven percent to file size. The same quarter-million grouping on unsorted rows reads 812 megabytes in 1.78 seconds. Row-group size against a bounding-box query 14.2 M building features sorted by H3 cell, filter selecting ~0.4% of rows one row group (no pruning) 980 MB read · 2.10 s 1 M rows per group 172 MB read · 0.41 s 250 k rows per group 61 MB read · 0.17 s 50 k rows per group 38 MB read · 0.19 s · +7% file size unsorted, 250 k rows 812 MB read · 1.78 s The last row is the same row-group size with the rows in arrival order — sorting is what makes the statistics selective, and without it the groups barely prune at all.
Sorting and row-group size only work together. Small groups over unsorted rows prune almost nothing, because every group spans the whole extent.

Verification

Open the file with a reader that was not involved in writing it:

python 
import geopandas as gpd
gdf = gpd.read_parquet("buildings.parquet")
print(gdf.crs, len(gdf), gdf.geometry.geom_type.value_counts().to_dict())

Three things should be true. The CRS prints as EPSG:4326 (readers resolve a null CRS to CRS84, which is equivalent). The row count matches what was written. And the geometry types match what the metadata declared — a mismatch means the declaration was copied rather than derived.

Then check the pruning actually works:

python 
import duckdb
duckdb.sql("INSTALL spatial; LOAD spatial;")
duckdb.sql("""
  SELECT count(*) FROM 'buildings.parquet'
  WHERE h3_r6 BETWEEN '861f8d4ffffffff' AND '861f8d5ffffffff'
""").show()

Compare the bytes read for that query against a full scan. A well-sorted file with quarter-million-row groups should read single-digit percentages of the file for a small area.

Common errors and fixes

Symptom Root cause Fix
GeoPandas: “no geometry column” geo metadata missing Attach it before writing
Everything plots off West Africa Projected coords with a null CRS Set the PROJJSON CRS, or write CRS84
Filter reads the whole file Rows unsorted, or one row group Sort by a locality key; set row_group_size
ArrowInvalid on the tags column Mixed types in the tag map Coerce every tag value to str first
File much larger than expected compression=None, the default Pass compression="zstd"
Geometry column is all nulls to_wkb given a Shapely 1.x geometry Upgrade to Shapely 2.0+, or use geom.wkb

Frequently Asked Questions

Should I write one file or many?

Many, partitioned by something with a few hundred distinct values — a country, a region, a coarse cell. One enormous file cannot be read in parallel by row-range and forces every consumer through one metadata footer; millions of tiny files make listing the dataset slower than reading it. Aim for partitions in the low hundreds of megabytes.

Is WKB the only geometry encoding allowed?

WKB is the encoding every reader supports and the safe default. GeoParquet 1.1 also permits a native Arrow geometry encoding which is faster to read because it avoids parsing WKB per row, but support is not yet universal. Write WKB unless you control every consumer.

Do I need to write the bbox in the metadata?

It is optional and worth including. Readers use it to skip a file entirely when it cannot intersect a query, which for a partitioned dataset means most files are never opened. Compute it from the data rather than from the boundary you cut with, since the two differ wherever a clipping strategy kept features that spill past the line.

How do I add a column later without rewriting everything?

You cannot, within a file — Parquet is immutable. What you can do is write the new column as a separate dataset keyed on osm_id and join at read time, or accept the rewrite for the partitions that need it. This is the practical argument for keeping the tag map: promoting a tag to a column later becomes a schema change rather than a re-extract from the PBF.

Writing incrementally

Building the whole table in memory before writing works up to a few million features and stops working somewhere below a country-sized layer. pq.ParquetWriter takes the schema up front and accepts batches, which keeps peak memory at one batch rather than one layer:

python 
writer = pq.ParquetWriter(path, schema, compression="zstd",
                          compression_level=3, write_statistics=True)
for batch in batches:
    writer.write_table(batch, row_group_size=ROW_GROUP_ROWS)
writer.close()

There is a real trade here, and it is worth stating plainly. Streaming batches means you cannot sort globally, because sorting needs every row at once — so the locality that makes row-group statistics selective is lost unless the batches already arrive in spatial order. Two ways out: sort each batch and accept coarse, per-batch locality, which recovers most of the pruning when batches are large; or write unsorted and run a compaction pass afterwards that reads, sorts and rewrites, which costs one extra full pass and gives the ideal layout.

Which is right depends on how often the file is read. For a dataset written nightly and queried thousands of times a day, the compaction pass pays for itself before breakfast. For an intermediate artefact read once by the next pipeline stage, per-batch sorting is enough and the extra pass is waste.

The schema must be identical across every batch — a column that is all-null in one batch and typed in another will raise on the second write, which is the usual reason an incremental writer fails halfway through a long run. Build the schema once from the promoted-column list rather than inferring it per batch.

Up one level: Exporting OSM to GeoParquet & PostGIS.