How to Build Vector Tiles from a GeoDataFrame in Python
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 59,391 real building polygons:
one zoom-15 tile, 182 buildings
MVT 5.80 kB gzipped 4.54 kB
MVT + 2 attributes 8.02 kB gzipped 5.73 kB
the same as GeoJSON 95.47 kB gzipped 14.72 kB
generalisation at zoom 10
simplify to 1 pixel 57.9% of vertices kept
drop below 2 pixels 3 of 59,391 features kept
Generalisation is not an optimisation you add later. At low zoom it is the difference between 3 features and 59,391.
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=()):
bounds = mercantile.bounds(tile)
inside = gdf.cx[bounds.west:bounds.east, bounds.south:bounds.north]
if inside.empty:
return None
def to_tile(x, y, z=None):
return ((np.asarray(x) - bounds.west) / (bounds.east - bounds.west) * extent,
(np.asarray(y) - bounds.south) / (bounds.north - bounds.south) * 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}], extents=extent)
Step-by-step solution
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) # manyAt zoom 16 there are 625 tiles over a city extent. Simplifying inside the tile loop does the same work 625 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
sindex = gdf.sindex
candidates = list(sindex.intersection(tile_bounds))gdf.cx[...] is convenient and rebuilds its filter each call. With 625 tiles and 59,391 features, using the index directly is the difference between minutes and seconds.
3. Buffer the selection past the tile edge
pad_x = (bounds.east - bounds.west) * 64 / 4096
inside = gdf.cx[bounds.west - pad_x:bounds.east + pad_x, ...]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, extents=4096), 6)Store .pbf gzipped. Every client expects Content-Encoding: gzip, and it roughly halves the transfer — 5.80 kB to 4.54 kB on the measured tile.
Code examples
Example 1 — a complete pyramid builder
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."""
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(4326)
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.bounds(tile)
pad_x = (tb.east - tb.west) * buffer_units / extent
pad_y = (tb.north - tb.south) * buffer_units / extent
idx = list(sindex.intersection(
(tb.west - pad_x, tb.south - pad_y,
tb.east + pad_x, tb.north + pad_y)))
if not idx:
continue
rows = subset.iloc[idx]
def to_tile(x, y, z=None, tb=tb):
return ((np.asarray(x) - tb.west) / (tb.east - tb.west) * extent,
(np.asarray(y) - tb.south) / (tb.north - tb.south) * 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}],
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 writtenReprojecting 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.
layer_specs: {"buildings": (gdf, ("building",)), "roads": (gdf2, ("highway",))}
"""
bounds = mercantile.bounds(tile)
layers = []
for name, (gdf, properties) in layer_specs.items():
inside = gdf.cx[bounds.west:bounds.east, bounds.south:bounds.north]
if inside.empty:
continue
def to_tile(x, y, z=None):
return ((np.asarray(x) - bounds.west) /
(bounds.east - bounds.west) * extent,
(np.asarray(y) - bounds.south) /
(bounds.north - bounds.south) * 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, 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 encodedLayer 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
import os
from concurrent.futures import ProcessPoolExecutor
import mercantile
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)
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(build_zoom, tasks):
total += count
print(f" {total:,} tiles across {len(list(zooms))} zoom levels")
return totalParallelising 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.
Explanation
Why generalisation belongs outside the tile loop
Two reasons, one of speed and one of correctness.
Speed: with 625 tiles at zoom 16, simplifying inside the loop repeats the work 625 times 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 that is 625 queries against 59,391 features.
gdf.cx[...] builds its mask each time, so the cost is linear in the feature count per tile — 37 million comparisons across the level. Using sindex.intersection is an R-tree lookup: logarithmic per query.
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 5.80 kB raw and 4.54 kB gzipped — 22%. 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.
Edge cases or notes
- 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 measured tile of 182 buildings with two attributes was 5.73 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 5.80 kB to 4.54 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.