Decoding MVT Tiles to GeoJSON in Python

To turn a .pbf vector tile back into GeoJSON in Python, decode the protobuf with the mapbox_vector_tile package and then apply an affine transform that maps the tile’s local 0..extent grid onto the EPSG:4326 bounding box implied by its (z, x, y) address. The decode is one function call; the coordinate rescale is the part that trips people up, because the tile bytes carry no geographic coordinates of their own.

When to Use This

Reach for a Python decoder when you need to know what actually shipped in a tile rather than what your source data or Tippecanoe logs claim. Two situations dominate:

  • Debugging production tiles. A layer is missing, an attribute is absent, or a feature is clipped. Decoding the exact .pbf the CDN served answers “is the data in the tile?” definitively, isolating generation bugs from styling bugs.
  • CI assertions on tile contents. After a nightly Tippecanoe build, a test decodes a handful of representative tiles and asserts that expected layers exist, feature counts are within tolerance, and required attribute keys survived filtering — catching regressions before they reach clients.

For a broad structural walkthrough of the fields you are decoding, read the MVT encoding internals guide first; this page is the runnable Python counterpart.

Specification Detail

mapbox_vector_tile.decode(raw_bytes) returns a dict keyed by layer name. Each layer is a dict with extent, version, and a features list; each feature has geometry (already GeoJSON-shaped, but in tile-local coordinates) and properties (the keys[]/values[] tables already joined for you).

Element Shape after decode() Notes
top level { "roads": {...}, "water": {...} } one entry per Layer.name
layer["extent"] int, usually 4096 the divisor for rescaling — never hard-code it
layer["features"] list[dict] each has geometry, properties, id, type
feature["geometry"] GeoJSON geometry coordinates are 0..extent, y-axis up by default
feature["properties"] dict decoded typed attributes

The affine transform from tile-local coordinates to lon/lat has two stages. First, the tile’s Web Mercator bounds come from its address: at zoom z there are 2**z tiles per axis, so tile (x, y) spans a known slice of the mercator square. Second, within the tile, a local coordinate cx in 0..extent is a simple linear interpolation across that slice. Because MVT’s local grid has its origin at the tile’s top-left with y increasing downward, but mapbox_vector_tile yields geometry with y increasing upward, you flip y during rescale. The longitude is linear in mercator x; the latitude requires the inverse Web Mercator (atan(sinh(...))) because latitude is non-linear in the projection.

One consequence of the linear-per-tile interpolation is worth internalizing: the latitude edges (north, south) are computed once per tile with the non-linear inverse, and everything inside the tile is then interpolated linearly between them. That approximation is exact enough at the scale of a single tile — a tile spans a tiny fraction of a degree at anything past zoom 8 — that the error is well below MVT’s own quantization grid. Do not be tempted to run the inverse-mercator per vertex; it is both slower and no more accurate than interpolating within the pre-computed tile bounds.

The properties dict also preserves attribute types faithfully. MVT’s values[] table stores typed scalars — string, int, uint, sint, float, double, bool — and mapbox_vector_tile restores them to the matching Python types. So an assertion like assert feat["properties"]["population"] > 0 works directly against an int, and a value that unexpectedly arrives as a str is a strong signal that mixed source types forced a string coercion during generation.

Production Command

The block below fetches one tile over HTTP, handles gzip, decodes it, rescales every coordinate to EPSG:4326, and writes one GeoJSON FeatureCollection per layer. It depends only on httpx and mapbox-vector-tile.

python
import gzip
import json
import math
from pathlib import Path

import httpx
import mapbox_vector_tile


def tile_lonlat_bounds(z: int, x: int, y: int) -> tuple[float, float, float, float]:
    """Return (west, south, east, north) in EPSG:4326 for an XYZ tile."""
    n = 2.0 ** z
    west = x / n * 360.0 - 180.0
    east = (x + 1) / n * 360.0 - 180.0
    north = math.degrees(math.atan(math.sinh(math.pi * (1 - 2 * y / n))))
    south = math.degrees(math.atan(math.sinh(math.pi * (1 - 2 * (y + 1) / n))))
    return west, south, east, north


def rescale(coords, extent, bounds):
    """Recursively rescale tile-local (0..extent) coords to lon/lat."""
    west, south, east, north = bounds
    if isinstance(coords[0], (list, tuple)):
        return [rescale(c, extent, bounds) for c in coords]
    cx, cy = coords
    lon = west + (cx / extent) * (east - west)
    lat = south + (cy / extent) * (north - south)  # y is already up after decode
    return [lon, lat]


def fetch_tile(url: str) -> bytes:
    resp = httpx.get(url, headers={"Accept-Encoding": "gzip"})
    resp.raise_for_status()
    body = resp.content
    # httpx transparently gunzips Content-Encoding; guard against double-gzip
    # (tiles stored pre-gzipped and served without the header) via magic bytes.
    if body[:2] == b"\x1f\x8b":
        body = gzip.decompress(body)
    return body


def tile_to_geojson(url: str, z: int, x: int, y: int) -> dict[str, dict]:
    raw = fetch_tile(url)
    decoded = mapbox_vector_tile.decode(raw)
    bounds = tile_lonlat_bounds(z, x, y)
    collections: dict[str, dict] = {}
    for layer_name, layer in decoded.items():
        extent = layer["extent"]
        features = []
        for feat in layer["features"]:
            geom = feat["geometry"]
            geom = {
                **geom,
                "coordinates": rescale(geom["coordinates"], extent, bounds),
            }
            features.append({
                "type": "Feature",
                "geometry": geom,
                "properties": feat.get("properties", {}),
                "id": feat.get("id"),
            })
        collections[layer_name] = {
            "type": "FeatureCollection",
            "features": features,
        }
    return collections


if __name__ == "__main__":
    Z, X, Y = 10, 512, 384
    url = f"https://tiles.example.com/v2/{Z}/{X}/{Y}.pbf"
    for name, fc in tile_to_geojson(url, Z, X, Y).items():
        out = Path(f"{name}_{Z}_{X}_{Y}.geojson")
        out.write_text(json.dumps(fc))
        print(f"{name}: {len(fc['features'])} features -> {out}")

Run it against any XYZ endpoint. For a (z, x, y) you extracted from a local archive, serve the archive first (pmtiles serve roads.pmtiles --port 8080) and point the URL at http://localhost:8080.

Interaction Effects

Gzip and Content-Encoding. mapbox_vector_tile.decode() expects raw protobuf. If bytes still carry the gzip magic header 1f 8b, decode raises immediately. httpx and requests decompress transparently when the server sets Content-Encoding: gzip, but tiles are frequently stored pre-gzipped on disk and served with Content-Type: application/x-protobuf and no encoding header — so the client never decompresses. The magic-byte guard in fetch_tile() handles both cases. This is the same double-compression hazard that surfaces on the delivery side; see the parent MVT encoding internals troubleshooting notes.

Layer name equals style source-layer. The keys of the decode() result are the exact Layer.name values a MapLibre style must reference as source-layer. A quick print(list(decoded.keys())) is the fastest way to confirm the contract that keeps styles and tiles in sync.

Decoding straight from an archive. When the tiles live in a PMTiles or MBTiles container rather than behind an HTTP endpoint, skip the network entirely. For MBTiles, remember the tiles table is TMS-indexed, so flip y before decoding:

python
import sqlite3, gzip, mapbox_vector_tile

con = sqlite3.connect("roads.mbtiles")
z, x, y = 10, 512, 384
tms_y = (2 ** z - 1) - y  # MBTiles stores rows in TMS order
row = con.execute(
    "SELECT tile_data FROM tiles WHERE zoom_level=? AND tile_column=? AND tile_row=?",
    (z, x, tms_y),
).fetchone()
blob = row[0]
if blob[:2] == b"\x1f\x8b":
    blob = gzip.decompress(blob)
decoded = mapbox_vector_tile.decode(blob)

Feeding these bytes into the same rescale() and tile_lonlat_bounds() helpers — using the XYZ y, not tms_y — produces identical GeoJSON without a running tile server.

Cross-checking with tippecanoe-decode. For a second opinion, decode the same tile with the Tippecanoe tool and diff feature counts:

bash
tippecanoe-decode tile.pbf 10 512 384 | python3 -c "import sys,json; d=json.load(sys.stdin); print(sum(len(l['features']) for l in d['layers']))"

If the Python and tippecanoe-decode counts disagree, you are almost certainly decoding stale or double-compressed bytes.

Performance Impact

Decoding is cheap per tile — a typical 50–150 KB tile decodes in single-digit milliseconds — but it is pure-Python object construction, so cost scales with feature count, not byte count. A dense z14 building tile with 20,000 features can take 50–100 ms and allocate several megabytes of dicts. The rescale recursion adds a modest constant factor per vertex.

The practical guidance is to sample, not sweep. A CI job that decodes every tile in a continental archive will run for hours and dwarf the Tippecanoe build itself. Instead, decode a fixed set of representative tiles — one dense urban tile, one sparse rural tile, one coastline tile per zoom band — which catches the regressions that matter (missing layers, dropped attributes, clipped coastlines) in seconds. When you do need bulk decoding, stream tiles out of the archive with pmtiles or a direct SQLite read rather than fetching each over HTTP, which removes network latency from the loop entirely.

Common Mistakes

Forgetting to gunzip. The error is blunt:

text
google.protobuf.message.DecodeError: Error parsing message

or a mapbox_vector_tile failure on the first field. It means gzip bytes reached the decoder. Check the first two bytes; if they are 1f 8b, gzip.decompress() before decoding. This is the single most common failure and the reason the guard above exists.

Not rescaling — coordinates stuck in 0..4096. If your output GeoJSON has coordinates like [2048, 3711] instead of [-73.98, 40.71], you skipped the affine transform. The geometry is correct in tile space; it just has not been mapped to lon/lat. Apply rescale() with the tile’s bounds, and read extent from the layer rather than assuming 4096 — a tile generated with a custom --extent will otherwise be scaled wrong.

Y-flip and TMS confusion. Two independent y-axis conventions collide here. First, if the archive is TMS-indexed (y origin at the south, common in older MBTiles), convert to XYZ with xyz_y = 2**z - 1 - tms_y before computing bounds, or every tile lands in the wrong hemisphere. Second, MVT’s internal grid has y increasing downward while mapbox_vector_tile returns y increasing upward — the code above already accounts for the latter, so do not flip a second time. The full coordinate derivation lives in converting lat/lon to slippy map tile numbers, the inverse of the bounds math used here.


Parent: MVT Encoding Internals: The Protobuf Tile Structure