Encoding OSM Geometry into MVT with Python

Turn a handful of OSM features into a valid vector tile from first principles, so that when a generator produces something unexpected you know exactly which stage to look at.

Prerequisites

Conceptual minimum

Encoding is four transformations applied in a fixed order, and getting the order wrong is the source of most defects.

Project geographic degrees to Web Mercator metres. Locate the tile: at zoom \(z\) the world spans \(2^{z}\) tiles in each direction, so a tile’s origin in projected space follows from its x and y. Quantise the geometry into the tile’s local grid by subtracting that origin and scaling by the extent. Encode the resulting integers as command and parameter arrays.

The reason order matters is that simplification and validity must happen on the quantised coordinates. A polygon simplified in degrees and then rounded onto a coarse grid can acquire zero-length segments; a ring whose orientation is checked before rounding can flip afterwards if two vertices collapse.

The four encoding stages and what each one must not be moved past Four stacked stages in order. Projection converts geographic degrees into the Web Mercator plane and must happen before anything tile-specific. Tile location computes the tile's origin and span in projected units from its zoom and coordinates. Quantisation subtracts that origin, scales by the extent and rounds to integers, which is the only lossy step. Encoding writes command integers and zigzag deltas, and both validity checking and winding enforcement must happen after quantisation rather than before it. Four stages, and two checks that belong at the end Project Degrees to Web Mercator metres before anything tile-local Locate Tile origin and span from z, x, y pure arithmetic Quantise Subtract, scale, round to integers the only lossy step Encode Commands, counts and zigzag deltas validate here, not earlier Checking validity and winding before quantisation tests geometry that will not be the geometry actually stored in the tile.
Every correctness check belongs on the rounded integers, because those are the coordinates a client will read.

Runnable solution

python 
from __future__ import annotations

import logging
import math
from dataclasses import dataclass
from typing import Any, Iterable

from shapely.geometry import LineString, Point, Polygon, mapping
from shapely.geometry.base import BaseGeometry
from shapely.ops import transform

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
logger = logging.getLogger("osm.mvt.encode")

EXTENT = 4096
EARTH_CIRCUMFERENCE = 40_075_016.686
ORIGIN_SHIFT = EARTH_CIRCUMFERENCE / 2.0

MOVE_TO, LINE_TO, CLOSE_PATH = 1, 2, 7


@dataclass(frozen=True)
class TileId:
    z: int
    x: int
    y: int

    def span(self) -> float:
        """Width of this tile in projected metres."""
        return EARTH_CIRCUMFERENCE / (2 ** self.z)

    def origin(self) -> tuple[float, float]:
        """Top-left corner of the tile in projected metres."""
        s = self.span()
        return (-ORIGIN_SHIFT + self.x * s, ORIGIN_SHIFT - self.y * s)


def to_mercator(lon: float, lat: float) -> tuple[float, float]:
    """WGS 84 degrees -> Web Mercator metres. Latitude is clamped at the poles."""
    lat = max(min(lat, 85.05112878), -85.05112878)
    x = lon * ORIGIN_SHIFT / 180.0
    y = math.log(math.tan((90.0 + lat) * math.pi / 360.0)) / (math.pi / 180.0)
    return x, y * ORIGIN_SHIFT / 180.0


def quantiser(tile: TileId, extent: int = EXTENT):
    """Return a function mapping projected metres onto the tile's integer grid."""
    ox, oy = tile.origin()
    span = tile.span()
    scale = extent / span

    def q(x: float, y: float) -> tuple[int, int]:
        # y increases downward in tile space, so it is subtracted, not added.
        return (round((x - ox) * scale), round((oy - y) * scale))

    return q


def command(cid: int, count: int) -> int:
    return (cid & 0x7) | (count << 3)


def zigzag(value: int) -> int:
    return (value << 1) ^ (value >> 31)


def encode_ring(points: list[tuple[int, int]], cursor: list[int],
                close: bool) -> list[int]:
    """Encode one ring or line, advancing the SHARED cursor."""
    out: list[int] = [command(MOVE_TO, 1)]
    dx, dy = points[0][0] - cursor[0], points[0][1] - cursor[1]
    out += [zigzag(dx), zigzag(dy)]
    cursor[0], cursor[1] = points[0]

    rest = points[1:-1] if close else points[1:]
    if rest:
        out.append(command(LINE_TO, len(rest)))
        for px, py in rest:
            out += [zigzag(px - cursor[0]), zigzag(py - cursor[1])]
            cursor[0], cursor[1] = px, py
    if close:
        out.append(command(CLOSE_PATH, 1))
    return out


def signed_area(points: list[tuple[int, int]]) -> float:
    """Shoelace on GRID coordinates; positive means clockwise in screen space."""
    total = 0.0
    for i in range(len(points) - 1):
        x1, y1 = points[i]
        x2, y2 = points[i + 1]
        total += (x2 - x1) * (y2 + y1)
    return total / 2.0


def encode_polygon(geom: Polygon, q) -> list[int]:
    cursor = [0, 0]
    out: list[int] = []
    rings = [list(geom.exterior.coords)] + [list(r.coords) for r in geom.interiors]
    for index, ring in enumerate(rings):
        grid = [q(x, y) for x, y in ring]
        # Drop consecutive duplicates created by rounding, or the ring is invalid.
        deduped = [grid[0]]
        for pt in grid[1:]:
            if pt != deduped[-1]:
                deduped.append(pt)
        if len(deduped) < 4:
            logger.warning("ring collapsed to %d point(s) at this zoom", len(deduped))
            continue
        if deduped[0] != deduped[-1]:
            deduped.append(deduped[0])
        # Exterior clockwise, interior counter-clockwise — checked AFTER rounding.
        want_clockwise = index == 0
        if (signed_area(deduped) > 0) != want_clockwise:
            deduped.reverse()
        out += encode_ring(deduped, cursor, close=True)
    return out


def encode_feature(geom: BaseGeometry, tile: TileId) -> list[int]:
    q = quantiser(tile)
    if isinstance(geom, Point):
        cursor = [0, 0]
        gx, gy = q(geom.x, geom.y)
        return [command(MOVE_TO, 1), zigzag(gx), zigzag(gy)]
    if isinstance(geom, LineString):
        cursor = [0, 0]
        return encode_ring([q(x, y) for x, y in geom.coords], cursor, close=False)
    if isinstance(geom, Polygon):
        return encode_polygon(geom, q)
    raise TypeError(f"unsupported geometry type {geom.geom_type}")


def build_layer(name: str, features: Iterable[tuple[BaseGeometry, dict[str, Any]]],
                tile: TileId) -> dict[str, Any]:
    """Assemble a layer with deduplicated key and value tables."""
    keys: list[str] = []
    values: list[Any] = []
    key_index: dict[str, int] = {}
    value_index: dict[Any, int] = {}
    encoded: list[dict[str, Any]] = []

    for geom, attrs in features:
        merc = transform(to_mercator, geom)
        pairs: list[int] = []
        for key, value in attrs.items():
            if key not in key_index:
                key_index[key] = len(keys)
                keys.append(key)
            if value not in value_index:
                value_index[value] = len(values)
                values.append(value)
            pairs += [key_index[key], value_index[value]]
        encoded.append({"geometry": encode_feature(merc, tile), "tags": pairs})

    logger.info("layer %s: %d feature(s), %d key(s), %d distinct value(s)",
                name, len(encoded), len(keys), len(values))
    return {"name": name, "extent": EXTENT, "version": 2,
            "keys": keys, "values": values, "features": encoded}


if __name__ == "__main__":
    tile = TileId(z=14, x=9111, y=5455)
    layer = build_layer("pois", [
        (Point(19.9373, 50.0617), {"class": "pharmacy", "rank": 2}),
    ], tile)
    logger.info("first feature geometry array: %s", layer["features"][0]["geometry"])

Step-by-step walkthrough

  1. Clamp latitude before projecting. Web Mercator is undefined at the poles; clamping at roughly 85.05 degrees is the standard convention and avoids an infinity in the logarithm.
  2. Compute the tile origin arithmetically. The origin and span follow from the zoom and the tile coordinates alone, with no lookup table, which makes the encoder trivially testable.
  3. Invert the y axis during quantisation. Tile space has its origin at the top-left with y increasing downwards, the opposite of projected space. Missing this mirrors every tile vertically, which looks plausible for a symmetric shape and obviously wrong for a coastline.
  4. Share one cursor across the whole feature. encode_ring takes a mutable cursor and advances it. Resetting it per ring is the classic bug that scatters holes to the tile origin.
  5. Drop duplicates created by rounding. Two source vertices a few centimetres apart round to the same grid point at low zoom; leaving both produces a zero-length segment that makes the ring invalid.
  6. Reject collapsed rings explicitly. A ring reduced below four points no longer encloses anything and is logged and skipped rather than emitted as broken geometry.
  7. Enforce winding after rounding. The shoelace test runs on grid coordinates, because that is where the orientation a client sees is decided.
  8. Deduplicate attributes into tables. Keys and values are interned and features reference them by index, which is what keeps a layer with one repeated class value cheap regardless of feature count.
Three encoder bugs, their symptom on the rendered map, and the check that catches each Three panels. A mirrored tile comes from forgetting that tile space has its y axis pointing downwards, and renders as geography flipped vertically within each tile, caught by comparing a known asymmetric coastline against a reference. Scattered holes come from resetting the geometry cursor between rings, and render as stray shapes near the tile corner, caught by decoding a polygon with a hole and checking the ring start positions. Inverted polygons come from checking ring orientation before rounding, and render as missing landmasses, caught by asserting signed area on the quantised ring. Three bugs that render as something, just not the right thing Mirrored tile Cause: y axis not inverted Renders flipped per tile Symmetric shapes look fine Check: an asymmetric coastline Scattered holes Cause: cursor reset per ring Holes appear at the corner Exterior ring looks correct Check: decode ring starts Inverted polygon Cause: winding checked early Landmass renders as a hole Rounding flipped the area sign Check: shoelace on the grid All three produce output a decoder accepts without complaint, which is why each needs a positive assertion rather than an absence of errors.
A tile that decodes is not a tile that is correct; the decoder validates structure, not geography.
What each geometry type contributes to a tile and what to watch on each A grid of three geometry types against three concerns. A point contributes a single move command and two parameters, has no winding or validity concern, and is the cheapest feature possible. A line contributes a move followed by one line command and two parameters per vertex, has no winding concern, and is watched for vertex ordering because unordered vertices inflate the deltas. A polygon contributes a move, a line and a close command per ring, has strict winding requirements, and is watched for rings collapsing to fewer than four points after rounding. Three geometry types, three different things to watch Commands emitted Winding Watch for Point one move none nothing Line move plus line none vertex ordering Polygon move, line, close strict collapsed rings Multi-part repeat per part per ring shared cursor The bottom row is where the shared cursor matters most: each additional part continues from where the previous one ended.
Points are nearly free, lines are cheap if ordered, and polygons carry every rule the format has.

Verification

  • Round-trip a known point. Encode a coordinate whose tile and grid position you computed by hand, decode it, and confirm the integers match.
  • A polygon with a hole decodes with two rings. The interior ring must start near the exterior, not at the tile origin.
  • Signed areas have opposite signs. Exterior and interior rings must wind oppositely on the quantised coordinates.
  • An asymmetric shape is not mirrored. Encode a recognisable coastline and compare against a reference rendering; a flipped y axis is invisible on symmetric test data.
  • The value table is small. For a layer with a closed vocabulary, distinct value count should be in the tens regardless of feature count; if it tracks feature count, something unique is being encoded per feature.

Common errors and fixes

Symptom Root cause One-line fix
Geometry flipped vertically y axis not inverted for tile space Subtract from the tile origin rather than adding
Holes near the tile corner Cursor reset between rings Share one cursor across the whole feature
Ring reported invalid Duplicate points after rounding Drop consecutive duplicates before encoding
Landmass renders as a hole Winding checked before quantisation Run the shoelace test on grid coordinates
Enormous geometry arrays Vertices not in spatial order Order coordinates along the geometry before encoding
Huge value table A unique attribute per feature Remove identifiers from tile attributes
Decoder rejects parameters Signed integers written without zigzag Apply zigzag encoding to every delta

Specification reference

Vector tile geometry is encoded as a sequence of unsigned integers holding commands and parameters. A command integer encodes a command identifier in its three least significant bits and a repeat count in the remaining bits; MoveTo and LineTo are followed by two parameters per repetition, and ClosePath by none. Parameters are zigzag-encoded deltas relative to the previous position, with the cursor initialised at the tile origin for each feature. See the Mapbox Vector Tile specification for the command values, the zigzag definition and the polygon winding rules.

Frequently Asked Questions

Why does my tile render upside down?

Because tile coordinate space has its origin at the top-left with y increasing downwards, while projected coordinate space has y increasing northwards. Quantising by adding to the origin rather than subtracting from it mirrors every tile vertically. The bug is easy to miss because a symmetric test shape looks identical either way — use an asymmetric feature such as a real coastline to catch it.

Should the geometry cursor reset between rings?

No. The cursor is initialised once per feature at the tile origin and advances continuously through every ring of that feature. Resetting it between rings means each ring’s first delta is measured from the origin instead of from the previous ring’s last point, which places holes and subsequent parts of a multi-part geometry near the tile corner. The exterior ring still looks right, which is why the bug survives a quick visual check.

Why check winding after quantisation rather than before?

Because rounding can change it. A thin sliver whose vertices collapse onto the grid can end up with a signed area of the opposite sign, and a ring that was correctly wound in source coordinates is then wrongly wound in the tile. Since the client only ever sees the quantised coordinates, that is where the check has to happen — along with the duplicate-point removal that makes the ring valid in the first place.

How do I keep attribute tables small?

Encode a closed vocabulary rather than raw values, and never put a per-feature unique identifier in the tile attributes. Keys and values are interned per layer, so a class attribute with six possible values costs six entries no matter how many features carry it, while a name or an OSM identifier adds one entry per feature and can dominate the tile. If a client genuinely needs identifiers, weigh that cost explicitly rather than including them by habit.

Up one level: The Mapbox Vector Tile Spec & Tile Geometry.