How to Export OpenStreetMap Data to GeoPackage or Shapefile in Python

On this page

Problem statement

Getting OpenStreetMap data into a GeoDataFrame is half the job; the other half is saving it in a form that QGIS, ArcGIS or a colleague can open. OSM data resists the usual formats. One query returns points, lines and polygons together. Tag keys contain colons, such as addr:street and socket:type1_combo, and number in the hundreds. Values are all strings. A shapefile cannot hold mixed geometry types, truncates field names to ten characters and caps a file at 2 GB. A GeoPackage handles all of it, if the layers are organised sensibly.

This page reads every amenity in the Washington, D.C. OpenStreetMap extract of 1 January 2026 with all its tags. It writes the data to a GeoPackage, one layer per geometry type, and shows what goes wrong when the same data is written to a shapefile.

Quick answer

Split by geometry type and write one GeoPackage layer each:

fp = (osmium.FileProcessor(str(PBF)).with_areas()
      .with_filter(osmium.filter.KeyFilter("amenity"))
      .with_filter(osmium.filter.GeoInterfaceFilter()))
A = gpd.GeoDataFrame.from_features(fp, crs=4326)
A = A[A.geom_type != "LineString"]                              # line copies of closed ways
out = Path("data/dc_amenities.gpkg"); out.unlink(missing_ok=True)
for gtype, part in A.groupby(A.geom_type):
    keep = part.dropna(axis=1, how="all")                       # columns this layer actually uses
    keep.to_file(out, layer=f"amenity_{gtype.lower()}", driver="GPKG")
print(f"{len(A):,} features, {A.shape[1] - 1} tag columns in all")
print({l: (n, c) for l, n, c in [(lay, len(gpd.read_file(out, layer=lay)), gpd.read_file(out, layer=lay).shape[1] - 1) for lay in gpd.list_layers(out).name]})

16,802 features, 615 tag columns in all
{'amenity_multipolygon': (8054, 374), 'amenity_point': (8748, 473)}

The 16,802 amenity points and polygons carry 615 different tag keys. Split by type, the polygon layer needs 374 columns and the point layer 473.

GeoPackage or shapefile?

Why one layer is awkward, the shapefile, the columns to keep

1. Why a single layer is awkward

Write everything to one layer and check what the file records as its geometry type:

mixed = Path("data/dc_amenities_mixed.gpkg"); mixed.unlink(missing_ok=True)
A.to_file(mixed, layer="amenity", driver="GPKG")
print(gpd.list_layers(mixed).to_string(index=False))
print(f"geometry types in the layer: {gpd.read_file(mixed).geom_type.value_counts().to_dict()}")

   name geometry_type
amenity       Unknown
geometry types in the layer: {'Point': 8748, 'MultiPolygon': 8054}

Written as one layer, the GeoPackage stores 8,748 points and 8,054 multipolygons under a geometry type GDAL reports as Unknown. Some desktop tools then guess, or show only one type.

What a shapefile does to tag names.

2. The shapefile

Try the same data as a shapefile:

import warnings
try:
    A.to_file("data/dc_amenities.shp")
    print("mixed shapefile written")
except Exception as e:
    print(f"mixed geometry types: {type(e).__name__}: {str(e)[:120]}")
pts = A[A.geom_type == "Point"].dropna(axis=1, how="all")
with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    pts.to_file("data/dc_amenity_points.shp")
names = gpd.read_file("data/dc_amenity_points.shp").columns
import re
renamed = [n for n in names if re.search(r"_\d+$", n)]
print(f"points only: {pts.shape[1] - 1} columns written, {len(names) - 1} read back; examples {[n for n in names if ':' in n or len(n) == 10][:5]}")
print(f"names that collided after truncation and were renamed: {len(renamed)}, e.g. {renamed[:4]}")
print(f"warnings raised: {len(w)}")

mixed geometry types: FeatureError: Could not add feature to layer at index 8748: Attempt to write non-point (MULTIPOLYGON) geometry to point shapefile.
points only: 473 columns written, 473 read back; examples ['addr_house', 'addr_postc', 'addr_state', 'addr_stree', 'contact_fa']
names that collided after truncation and were renamed: 69, e.g. ['healthca_1', 'communit_1', 'operator_1', 'operator_2']
warnings raised: 299

The shapefile driver refuses the mixed data at the first polygon, feature 8,748. Even the points alone come back altered: 473 columns survive, but with laundered, truncated names such as addr_house and contact_fa. 69 of them collided after truncation and were renamed, among them operator_1 and operator_2, and the write raised 299 warnings.

3. Keep the columns that matter

Hundreds of mostly empty tag columns make any file awkward. Keep the tags present on at least 1% of features and put the rest in one JSON column:

share = A.drop(columns="geometry").notna().mean()
common = share[share >= 0.01].index.tolist()
rest = [c for c in A.columns if c not in common + ["geometry"]]
slim = A[common + ["geometry"]].copy()
slim["other_tags"] = A[rest].apply(lambda r: json.dumps({k: v for k, v in r.items() if pd.notna(v)}) if r.notna().any() else None, axis=1)
slim.to_file("data/dc_amenities_slim.gpkg", layer="amenity", driver="GPKG")
print(f"columns kept: {len(common)} of {A.shape[1] - 1}; features with other tags: {slim.other_tags.notna().sum():,}")
print(f"GeoPackage sizes: per-type {Path('data/dc_amenities.gpkg').stat().st_size / 1e6:.1f} MB, slim {Path('data/dc_amenities_slim.gpkg').stat().st_size / 1e6:.1f} MB")

columns kept: 60 of 615; features with other tags: 3,750
GeoPackage sizes: per-type 13.7 MB, slim 6.9 MB

Keeping the 60 tags present on at least 1% of features, and folding the rest into other_tags for 3,750 features, halves the GeoPackage from 13.7 MB to 6.9 MB without losing any tag.

Reading the export back

The data used on this page

The Geofabrik OpenStreetMap extract for the District of Columbia dated 1 January 2026 (© OpenStreetMap contributors, ODbL). Exported files must keep the attribution "© OpenStreetMap contributors" and the ODbL terms.

import json
from pathlib import Path

import geopandas as gpd
import osmium
import pandas as pd
import requests

Path("data").mkdir(exist_ok=True)
PBF = Path("data/district-of-columbia-260101.osm.pbf")
if not PBF.exists():
    r = requests.get("https://download.geofabrik.de/north-america/us/district-of-columbia-260101.osm.pbf", timeout=600)
    r.raise_for_status()
    PBF.write_bytes(r.content)
print(f"{PBF.name} ready")

district-of-columbia-260101.osm.pbf ready

The dated D.C. extract is ready.

Example 1: read it back

Open the GeoPackage the way another tool would, layer by layer, and check a tag with a colon survived:

for lay in gpd.list_layers(out).name:
    g = gpd.read_file(out, layer=lay)
    print(f"{lay:>22}: {len(g):,} features; has 'addr:street': {'addr:street' in g.columns}; CRS {g.crs.to_epsg()}")

  amenity_multipolygon: 8,054 features; has 'addr:street': True; CRS 4326
         amenity_point: 8,748 features; has 'addr:street': True; CRS 4326

Both layers read back intact, with their colon-named columns such as addr:street and their CRS, EPSG:4326.

Halving the file without losing tags.

A GeoPackage holds many layers

OSM to GeoPackage.

GeoPackage is a single SQLite file that holds any number of layers, each with its own geometry type, CRS and columns, with field names of any length and characters. That fits OSM data: one layer per geometry type keeps desktop GIS happy, colons in tag keys survive, and there is no 2 GB limit. A layer can be declared with the generic type GEOMETRY and hold points and polygons together, as GeoPandas does when it is given mixed data, but many tools then style and query it poorly, so splitting by type is the practical choice.

Shapefile is the format OSM data most often breaks. Each file holds one geometry type. Field names are limited to ten characters, so tag keys lose their colons and their ends, and keys that share their first ten characters, such as operator:wikidata and operator:wikipedia, end up with the same name, so the driver renames them with numeric suffixes. The DBF format limits text fields to 254 characters and the file to 2 GB. Whatever the format, the width of OSM's tag table is the real problem: one rare tag on one feature becomes a column for all of them. Keeping the common tags as columns and the rest in a JSON column, as GDAL's own OSM driver does with other_tags, gives a file that is both complete and usable.

Lists, field types and attribution

Before sharing an export.

Lists and dicts: if a reader produces list- or dict-valued columns (OSMnx can), convert them to strings before writing, or the write fails.
Field types: OSM values are strings; convert numbers such as capacity explicitly if you need them numeric.
GeoParquet: for large extracts and analysis in Python or DuckDB, to_parquet is faster and smaller than GeoPackage.
Attribution: OSM-derived files are covered by the ODbL; keep the attribution and licence with the data.

Internal links

How to read an OpenStreetMap PBF file in Python — the input side
Fixing OSM features counted twice — why the line copies are dropped
OpenStreetMap ODbL licence explained — what exporting the data obliges

FAQ

How do I save OpenStreetMap data to a GeoPackage in Python?

Read it into a GeoDataFrame, split it by geometry type, and write each part with to_file(path, layer=..., driver="GPKG").

Why can't I save OSM data as a shapefile?

Shapefiles hold one geometry type, truncate field names to ten characters, and cannot handle OSM's mixed geometries and long tag names.

How do I handle hundreds of OSM tag columns?

Keep the common tags as columns and store the rest in a single JSON text column.

Yes; OSM data is under the ODbL and requires attribution to OpenStreetMap contributors.