How to Build Vector Tiles from a GeoDataFrame in Python
On this page
Problem statement
Building a vector tile pyramid in Python means doing three things per tile: select the features, transform them into tile-local coordinates, and encode them as MVT — with generalisation applied per zoom level before any of that.
The measurements that shape the design, on Leiden's 58,879 OpenStreetMap buildings (measured on the format page and on the generalisation page):
one zoom-15 tile of central Leiden
MVT 69.61 kB gzipped 45.83 kB
MVT + 2 attributes 86.07 kB gzipped 48.12 kB
the same as GeoJSON 928.48 kB gzipped 126.63 kB
generalisation at zoom 10
simplify to 1 pixel 54.5% of vertices kept
drop below 2 pixels 1 of 58,879 features kept
Generalisation is not an optimisation you add later. At low zoom it is the difference between one feature and 58,879.
Quick answer
import mercantile
import mapbox_vector_tile as mvt
import numpy as np
from shapely.ops import transform
def encode_tile(gdf, tile, layer="layer", extent=4096, properties=()):
"""gdf in EPSG:3857 — the tile grid is Web Mercator."""
b = mercantile.xy_bounds(tile)
inside = gdf.cx[b.left:b.right, b.bottom:b.top]
if inside.empty:
return None
def to_tile(x, y, z=None):
return ((np.asarray(x) - b.left) / (b.right - b.left) * extent,
(np.asarray(y) - b.bottom) / (b.top - b.bottom) * extent)
features = [{"geometry": transform(to_tile, row.geometry),
"properties": {k: str(row[k]) for k in properties}}
for _, row in inside.iterrows() if row.geometry is not None]
return mvt.encode([{"name": layer, "features": features}],
default_options={"extents": extent})import geopandas as gpd
merc = gpd.read_parquet("buildings.parquet").to_crs(3857)
tile = mercantile.tile(4.4938, 52.1585, 15) # central Leiden
print(f"{len(encode_tile(merc, tile, 'buildings', properties=('building',))) / 1e3:.2f} kB")83.19 kB
Generalise per zoom, index once, buffer the selection
1. Generalise once per zoom, not once per tile
for zoom in zooms:
subset = generalise(gdf, zoom) # once
for tile in mercantile.tiles(*bounds, zooms=[zoom]):
encode_tile(subset, tile) # manyBOUNDS = tuple(gpd.read_parquet("buildings.parquet").total_bounds)
for zoom in (10, 14, 16):
print(f"zoom {zoom}: {len(list(mercantile.tiles(*BOUNDS, zooms=[zoom]))):,} tiles over Leiden")zoom 10: 1 tiles over Leiden
zoom 14: 25 tiles over Leiden
zoom 16: 306 tiles over Leiden
At zoom 16 there are 306 tiles over Leiden alone. Simplifying inside the tile loop does the same work 306 times, and it produces inconsistent geometry at tile boundaries because each tile simplifies its own clipped copy.
2. Build a spatial index before the tile loop
import time
tiles16 = list(mercantile.tiles(*BOUNDS, zooms=[16]))
boxes = [tuple(mercantile.xy_bounds(t)) for t in tiles16]
start = time.perf_counter()
via_cx = sum(len(merc.cx[x0:x1, y0:y1]) for x0, y0, x1, y1 in boxes)
cx_s = time.perf_counter() - start
sindex = merc.sindex
start = time.perf_counter()
via_index = sum(len(sindex.intersection(bounds)) for bounds in boxes)
index_s = time.perf_counter() - start
print(f"{len(boxes):,} tiles: .cx {cx_s:.2f} s, sindex {index_s:.3f} s ({via_cx:,} and {via_index:,} hits)")306 tiles: .cx 0.46 s, sindex 0.002 s (61,637 and 61,668 hits)
gdf.cx[...] is convenient and rebuilds its filter each call; the index is an R-tree lookup. The hit counts differ slightly because .cx tests the geometries while the index returns every bounding-box candidate.
3. Buffer the selection past the tile edge
pad = (b.right - b.left) * 64 / 4096
inside = gdf.cx[b.left - pad:b.right + pad, b.bottom - pad:b.top + pad]Without a buffer, a polygon outline is cut exactly at the seam and the stroke is visibly clipped. Sixty-four units at extent=4096 is 1.6% of the tile, enough for typical line widths.
4. Transform to tile coordinates
Tile coordinates are integers from 0 to extent relative to the tile's own corner, with y increasing north in the MVT convention while the tile grid numbers rows southward. The encoder handles the flip; you supply coordinates in the tile's local space.
5. Encode, compress and write
payload = gzip.compress(mvt.encode(layers, default_options={"extents": 4096}), 6)Store .pbf gzipped. Every client expects Content-Encoding: gzip, and it cut the measured tile by a third — 69.61 kB to 45.83 kB.
A pyramid builder and multi-layer tiles
The data used on this page
The buildings are every building in Leiden in OpenStreetMap — about 59,000 footprints with their building type and, where mapped, a name, height and number of levels — fetched with OSMnx through an Overpass mirror and saved as GeoParquet in EPSG:4326. The multi-layer example adds Leiden's roads from the same source. The OpenStreetMap answers are read as the Overpass mirror gave them on 9 October 2026.
import io
import os
import requests
import zipfile
import geopandas as gpd
import osmnx as ox
# OpenStreetMap as Overpass returned it on 9 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
reply = requests.get("https://data.spatialworkflow.io/osm-cache-build-vector-tiles-python-2026-10-09.zip", timeout=300)
reply.raise_for_status()
with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
archive.extractall("cache")
ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api" # overpass-api.de mirror
ox.settings.overpass_rate_limit = False
if not os.path.exists("buildings.parquet"):
found = ox.features_from_place("Leiden, Netherlands", {"building": True}).reset_index()
found = found[found.geom_type.isin(["Polygon", "MultiPolygon"])]
keep = [c for c in ("osmid", "building", "name", "height", "building:levels") if c in found.columns]
found[keep + ["geometry"]].rename(columns={"building:levels": "levels"}).to_parquet("buildings.parquet")
buildings = gpd.read_parquet("buildings.parquet")
print(f"{len(buildings):,} OSM buildings in Leiden, EPSG:{buildings.crs.to_epsg()}; "
f"columns {[c for c in buildings.columns if c != 'geometry']}")
if not os.path.exists("roads.parquet"):
roads = ox.features_from_place("Leiden, Netherlands", {"highway": True}).reset_index()
roads = roads[roads.geom_type == "LineString"][["highway", "geometry"]]
roads.to_parquet("roads.parquet")
print(f"{len(gpd.read_parquet('roads.parquet')):,} road segments")58,879 OSM buildings in Leiden, EPSG:4326; columns ['building', 'name', 'height', 'levels']
15,201 road segments
Example 1 — a complete pyramid builder
The builder is a module, so Example 3 can hand it to worker processes.
%%writefile tilebuilder.py
import gzip
import math
import os
import mercantile
import numpy as np
import mapbox_vector_tile as mvt
from shapely.ops import transform
def build_tiles(gdf, bounds, out_dir, zooms=range(10, 17), layer="layer",
properties=(), extent=4096, buffer_units=64,
min_visible_px=2.0, compress=True):
"""A tile pyramid with per-zoom generalisation. gdf in EPSG:4326."""
os.makedirs(out_dir, exist_ok=True)
lat = (bounds[1] + bounds[3]) / 2
utm_crs = gdf.estimate_utm_crs()
written, total_bytes = 0, 0
for zoom in zooms:
m_per_px = 156543.03392 * math.cos(math.radians(lat)) / (2 ** zoom)
utm = gdf.to_crs(utm_crs)
keep = utm.geometry.area >= (min_visible_px * m_per_px) ** 2
utm = utm[keep].copy()
utm["geometry"] = utm.geometry.simplify(m_per_px,
preserve_topology=True)
subset = utm.to_crs(3857) # the tile grid's own CRS
sindex = subset.sindex
print(f" z{zoom}: {len(subset):,}/{len(gdf):,} features "
f"({keep.mean():.1%}) at {m_per_px:.1f} m/px")
for tile in mercantile.tiles(*bounds, zooms=[zoom]):
tb = mercantile.xy_bounds(tile)
pad = (tb.right - tb.left) * buffer_units / extent
idx = list(sindex.intersection(
(tb.left - pad, tb.bottom - pad, tb.right + pad, tb.top + pad)))
if not idx:
continue
rows = subset.iloc[idx]
def to_tile(x, y, z=None, tb=tb):
return ((np.asarray(x) - tb.left) / (tb.right - tb.left) * extent,
(np.asarray(y) - tb.bottom) / (tb.top - tb.bottom) * extent)
features = []
for _, row in rows.iterrows():
if row.geometry is None or row.geometry.is_empty:
continue
features.append({
"geometry": transform(to_tile, row.geometry),
"properties": {k: ("" if row.get(k) is None else str(row[k]))
for k in properties},
})
if not features:
continue
payload = mvt.encode([{"name": layer, "features": features}],
default_options={"extents": extent})
if compress:
payload = gzip.compress(payload, 6)
directory = os.path.join(out_dir, str(tile.z), str(tile.x))
os.makedirs(directory, exist_ok=True)
with open(os.path.join(directory, f"{tile.y}.pbf"), "wb") as f:
f.write(payload)
written += 1
total_bytes += len(payload)
print(f" {written:,} tiles, {total_bytes / 1e6:.2f} MB, "
f"mean {total_bytes / max(written, 1) / 1024:.1f} kB per tile")
return writtenfrom tilebuilder import build_tiles
buildings = gpd.read_parquet("buildings.parquet")
start = time.perf_counter()
build_tiles(buildings, BOUNDS, "tiles", zooms=range(10, 17), layer="buildings", properties=("building",))
print(f" {time.perf_counter() - start:.1f} s") z10: 1/58,879 features (0.0%) at 93.8 m/px
z11: 11/58,879 features (0.0%) at 46.9 m/px
z12: 211/58,879 features (0.4%) at 23.4 m/px
z13: 947/58,879 features (1.6%) at 11.7 m/px
z14: 3,603/58,879 features (6.1%) at 5.9 m/px
z15: 33,438/58,879 features (56.8%) at 2.9 m/px
z16: 43,751/58,879 features (74.3%) at 1.5 m/px
280 tiles, 1.30 MB, mean 4.5 kB per tile
11.9 s
Reprojecting to UTM once per zoom rather than once per tile is the other large saving. Reprojection is expensive and its result is the same for every tile at that zoom.
Example 2 — several layers in one tile
import mapbox_vector_tile as mvt
def encode_multilayer(layer_specs, tile, extent=4096):
"""One tile carrying several named layers, each GeoDataFrame in EPSG:3857.
layer_specs: {"buildings": (gdf, ("building",)), "roads": (gdf2, ("highway",))}
"""
b = mercantile.xy_bounds(tile)
layers = []
for name, (gdf, properties) in layer_specs.items():
inside = gdf.cx[b.left:b.right, b.bottom:b.top]
if inside.empty:
continue
def to_tile(x, y, z=None):
return ((np.asarray(x) - b.left) / (b.right - b.left) * extent,
(np.asarray(y) - b.bottom) / (b.top - b.bottom) * extent)
layers.append({"name": name, "features": [
{"geometry": transform(to_tile, row.geometry),
"properties": {k: str(row[k]) for k in properties}}
for _, row in inside.iterrows() if row.geometry is not None]})
if not layers:
return None
encoded = mvt.encode(layers, default_options={"extents": extent})
sizes = {l["name"]: len(l["features"]) for l in layers}
print(f" z{tile.z}/{tile.x}/{tile.y}: {sizes}, {len(encoded) / 1024:.2f} kB")
return encodedroads = gpd.read_parquet("roads.parquet").to_crs(3857)
tile_bytes = encode_multilayer({"buildings": (merc, ("building",)), "roads": (roads, ("highway",))}, tile)
print(sorted(mvt.decode(tile_bytes))) z15/16793/10800: {'buildings': 2699, 'roads': 704}, 97.92 kB
['buildings', 'roads']
Layer names are the contract with every style that consumes the tiles. Changing one silently breaks every style, so treat them as a public interface.
Example 3 — parallel generation
%%writefile -a tilebuilder.py
def build_zoom(args):
gdf, bounds, zoom, out_dir, layer, properties = args
return build_tiles(gdf, bounds, out_dir, zooms=[zoom],
layer=layer, properties=properties)from concurrent.futures import ProcessPoolExecutor
import importlib
import tilebuilder
importlib.reload(tilebuilder)
def build_parallel(gdf, bounds, out_dir, zooms=range(10, 17), workers=4,
layer="layer", properties=()):
"""One process per zoom level — they share nothing."""
tasks = [(gdf, bounds, z, out_dir, layer, properties) for z in zooms]
total = 0
with ProcessPoolExecutor(workers) as pool:
for count in pool.map(tilebuilder.build_zoom, tasks):
total += count
print(f" {total:,} tiles across {len(list(zooms))} zoom levels")
return total
start = time.perf_counter()
build_parallel(buildings, BOUNDS, "tiles_parallel", workers=4, layer="buildings", properties=("building",))
print(f" {time.perf_counter() - start:.1f} s") 280 tiles across 7 zoom levels
7.0 s
Parallelising by zoom is the simplest split: each level does its own generalisation and writes to its own directory, so there is no coordination at all.
For a very large dataset, splitting by tile range within a zoom is better balanced — the high zooms have thousands of times more tiles than the low ones.
Generalisation belongs outside the tile loop
Why generalisation belongs outside the tile loop
Two reasons, one of speed and one of correctness.
Speed: with 306 tiles at zoom 16, simplifying inside the loop repeats the work once per tile on overlapping subsets.
Correctness: simplification is not local. Simplifying a polygon that has been clipped to a tile gives a different result from simplifying it whole and then clipping. Neighbouring tiles then disagree along their shared edge, producing visible cracks and slivers.
Generalise the whole layer once per zoom, then clip.
Why the spatial index matters so much
Selecting features for a tile is a bounding-box query, run once per tile. At zoom 16 over Leiden that is 306 queries against 58,879 features.
gdf.cx[...] builds its mask each time, so the cost is linear in the feature count per tile — 18 million comparisons across the level. Using sindex.intersection is an R-tree lookup: logarithmic per query. Measured, the whole level took 0.49 s with .cx and 0.002 s with the index.
The index is built once and reused, so the saving grows with the tile count.
Why the buffer is not optional
MVT clips geometry to the tile. A polygon outline drawn with a 2-pixel stroke, clipped exactly at the tile edge, shows the stroke ending abruptly — and the neighbouring tile shows the same, so the seam is a visible line across the map.
Buffering the selection means the geometry continues past the edge, the encoder clips it slightly outside the visible area, and the stroke crosses the seam cleanly.
Sixty-four units at extent=4096 is 1.6% of the tile, which covers typical stroke widths and label placement.
Why to write gzipped
The measured tile is 69.61 kB raw and 45.83 kB gzipped — a third smaller. On a pyramid of thousands of tiles that is a real saving in storage and transfer.
Every vector-tile client expects Content-Encoding: gzip, and storing the tiles already compressed means the server does not recompress on every request.
The one thing to get right is telling the client: serve the header, or record tile_compression in the archive metadata.
Indexes, buffers and extents
- Generalise per zoom, outside the tile loop — for speed and for seam consistency.
- Build the spatial index once and reuse it across tiles.
- Buffer the selection by 8–64 tile units.
extent=4096is the convention.- Keep layer names stable; they are the contract with every style.
- Write
.pbfgzipped and record the compression. - Skip empty tiles rather than writing empty archives.
- Parallelise by zoom for simplicity, by tile range for balance.
Internal links
- Vector tiles explained: MVT, layers and why they are not GeoJSON — the format being produced
- Generalisation for zoom levels explained — the per-zoom step
- Web map tiles explained: XYZ, zoom levels and the tile pyramid — the tile grid
- PMTiles and MBTiles explained — packaging the output
- How to serve a web map from PMTiles with MapLibre — displaying them
- How to prepare a GeoDataFrame for the web — the preparation before tiling
- How to use the spatial index directly in GeoPandas (sindex) — the index used here
- My tiles are offset or in the wrong place — when the grid is wrong
FAQ
How do I build vector tiles in Python?
Generalise the layer once per zoom, enumerate tiles with mercantile, select features per tile with a spatial index, transform to tile-local coordinates and encode with mapbox_vector_tile.
Should I simplify inside or outside the tile loop?
Outside, once per zoom. Inside is both slower and wrong: simplifying clipped copies makes neighbouring tiles disagree at the seam.
Why do I see lines along my tile boundaries?
No buffer on the selection. Geometry clipped exactly at the edge has its stroke cut; buffer by 8–64 tile units.
How big should a tile be?
Aim for a few tens of kilobytes. The densest tile of central Leiden, 2,699 buildings with two attributes, was 48.12 kB gzipped; tiles above about 500 kB are a sign of missing generalisation.
Do I need to gzip the tiles?
Yes. Every client expects it, and it took the measured tile from 69.61 kB to 45.83 kB.
How do I put several layers in one tile?
mvt.encode takes a list of layer dictionaries, each with a name and features. Styles refer to layers by name, so keep them stable.
How do I speed up generation?
Build the spatial index once, generalise once per zoom, and parallelise across zoom levels or tile ranges.