Tile Coordinate Systems & the Slippy Map Grid

Almost every misaligned, offset, or upside-down tileset traces back to a coordinate-system mistake made before the first tile was ever rendered. The pipeline juggles at least four distinct coordinate spaces — geographic degrees, projected meters, integer tile indices, and tile-local grid units — and confusing any two of them produces features that land in the ocean, a basemap flipped north-for-south, or an off-by-one seam at every tile edge. This guide fixes the vocabulary and the exact transforms so that data entering the tile generation pipeline stays spatially correct from source projection all the way to the {z}/{x}/{y} request the browser makes.

Prerequisites

Input geometry must already be in EPSG:4326 before it reaches the grid math. If your source is in a national grid or UTM zone, reproject first; the same EPSG:4326 enforcement step used for GeoParquet inputs applies to any format feeding a tile build.

Core Concept: Four Coordinate Systems, One Pipeline

The tiling scheme is not a single coordinate system — it is a chain of four, and every rendering bug is really a confusion between two adjacent links in that chain.

Web Mercator and the ±85.0511° limit

EPSG:3857 projects the sphere onto a square. Because the Mercator projection stretches vertically without bound as latitude approaches the poles, the tiling scheme truncates the world at latitude ±85.0511° — the exact value where the projected map becomes square (y reaches the same ±20 037 508.34 m extent as x). The precise constant is lat = ±(2·atan(e^π) − π/2)·180/π ≈ ±85.05112878°. Anything beyond that latitude has no tile and simply does not exist in a Web Mercator tileset; this is why polar data never appears on a standard slippy map.

A second consequence of the projection is distortion: Web Mercator preserves angles (it is conformal) but grossly inflates area toward the poles. A tile at 60° latitude represents half the ground width and a quarter the ground area of a tile at the equator, even though both render as the same on-screen square. This is why the ground-resolution figures later in this guide are quoted at the equator and must be scaled by cos(latitude) for any other row — a detail that matters when budgeting how many zoom levels a high-latitude city actually needs.

XYZ versus TMS: the y-axis flip

The tile grid has a settled x convention — column 0 is the westernmost, increasing east — but two competing conventions for y:

The conversion between them at zoom z is a single reflection:

y_tms = 2^z − 1 − y_xyz

The flip is its own inverse, so applying it twice returns the original. This one formula is behind the overwhelming majority of “my whole basemap is upside-down” reports.

Quadkeys

A quadkey collapses (z, x, y) into a single base-4 string of length z, where each character (0–3) selects a quadrant at successive zoom levels — the encoding Bing Maps popularized and the conceptual basis for the Hilbert-curve ordering used inside PMTiles archives. Interleaving the bits of x and y produces the digit at each level, which is why quadkeys sort spatially and make good cache keys. The digit at zoom level i (counting from the most significant bit) is bit_i(x) + 2·bit_i(y):

def tile_to_quadkey(x: int, y: int, z: int) -> str:
    digits = []
    for i in range(z, 0, -1):
        digit = 0
        mask = 1 << (i - 1)
        if x & mask: digit += 1
        if y & mask: digit += 2
        digits.append(str(digit))
    return "".join(digits)   # tile z2/1/1 -> "03"

Because the string prefix of a deep tile is exactly the quadkey of its parent, quadkeys make range scans over a subtree trivial — a property object stores exploit for locality.

Tile-local extent coordinates

Once a tile is selected, MVT does not store geographic coordinates inside it. Each tile carries its own integer grid, 0..extent (default 4096), where (0,0) is the tile’s north-west corner and (4096,4096) is just past its south-east corner. Converting a geographic vertex into these units is the last transform in the chain and is covered in detail in the MVT encoding internals reference.

The extent grid is deliberately decoupled from screen pixels: a tile is 4096 grid units wide regardless of whether it renders at 256 or 512 CSS pixels, so the same tile serves retina and non-retina clients without re-encoding. Values slightly outside 0..4096 are legal and expected — geometry is buffered a few units past each edge so that lines and polygon fills join seamlessly across the tile boundary rather than showing a hairline seam.

Step-by-Step: Walking the Coordinate Chain

Step 1 — Project lon/lat to Web Mercator

Transform geographic degrees into normalized Mercator space in [0, 1], clamping latitude to the projection limit first.

import math

MAX_LAT = 85.05112878  # Web Mercator latitude limit

def lonlat_to_mercator_norm(lon: float, lat: float) -> tuple[float, float]:
    lat = max(min(lat, MAX_LAT), -MAX_LAT)          # clamp to the Mercator limit
    x = (lon + 180.0) / 360.0                        # 0 at −180°, 1 at +180°
    siny = math.sin(math.radians(lat))
    y = 0.5 - math.log((1 + siny) / (1 - siny)) / (4 * math.pi)
    return x, y                                      # both in [0, 1], y=0 at north

Verify: lonlat_to_mercator_norm(0, 0) returns (0.5, 0.5) — the equator/prime-meridian origin sits at the center of the normalized square.

Step 2 — Compute the tile x/y at a zoom

Scale the normalized position by the grid size 2^z and floor.

def tile_xy(lon: float, lat: float, z: int) -> tuple[int, int]:
    nx, ny = lonlat_to_mercator_norm(lon, lat)
    n = 2 ** z
    x = min(int(nx * n), n - 1)   # clamp the +180° / south edge into the last tile
    y = min(int(ny * n), n - 1)
    return x, y                    # XYZ (north-origin) tile indices

Verify: for Berlin (13.405, 52.520) at z=12 this returns (2200, 1343) in XYZ.

Step 3 — Convert XYZ ↔ TMS

def xyz_to_tms(x: int, y: int, z: int) -> tuple[int, int]:
    return x, (2 ** z - 1) - y     # same function converts back the other way

Apply this exactly once, at the boundary where an XYZ-addressed request meets TMS-indexed storage — never inside the grid math itself.

Step 4 — Map a pixel within a tile to a coordinate

To place a tile-local extent point back into geography, first recover the tile’s bounding box, then interpolate the 0..extent grid across it.

def num2deg(x: int, y: int, z: int) -> tuple[float, float]:
    """North-west corner (lon, lat) of an XYZ tile."""
    n = 2 ** z
    lon = x / n * 360.0 - 180.0
    lat = math.degrees(math.atan(math.sinh(math.pi * (1 - 2 * y / n))))
    return lon, lat

def extent_point_to_lonlat(px: int, py: int, x: int, y: int, z: int,
                           extent: int = 4096) -> tuple[float, float]:
    lon_nw, lat_nw = num2deg(x, y, z)
    lon_se, lat_se = num2deg(x + 1, y + 1, z)
    lon = lon_nw + (lon_se - lon_nw) * (px / extent)
    lat = lat_nw + (lat_se - lat_nw) * (py / extent)   # linear approx within a tile
    return lon, lat

Note the latitude interpolation is linear here for illustration; because Mercator y is nonlinear in latitude, a fully correct inverse interpolates in normalized Mercator space and then applies num2deg. Within a single tile the error is sub-pixel and rarely matters, but at low zooms it does.

Optimization Knobs

Resolution per zoom (Web Mercator, 256 px tiles, at the equator)

Resolution shrinks with cos(latitude), so a tile at 60°N covers half the ground width of one at the equator — a key input when deciding how deep to build.

Integration with Adjacent Pipeline Stages

The coordinate chain threads through the entire build. Tippecanoe and ogr2ogr require EPSG:4326 input, then internally reproject to the EPSG:3857 grid during encoding — which is why the GeoParquet ingestion path enforces 4326 before a single feature is written. Inside each tile, geometry lives in the 0..4096 extent space described by the MVT encoding internals.

Storage is where the y-axis convention becomes a landmine. MBTiles stores its tile_row in TMS order (south origin), while PMTiles and virtually every {z}/{x}/{y} URL a browser requests use XYZ (north origin). A tile server reading MBTiles must apply y_tms = 2^z − 1 − y on the way out; a server that forgets serves a vertically mirrored map. When debugging, always confirm which convention each end of the connection speaks before touching the grid math.

The grid also drives batch operations that never touch a browser. Cache pre-warming enumerates the tile range covering a bounding box — convert the bbox corners with the forward transform, then iterate x from the west tile to the east tile and y from the north tile to the south tile:

def tiles_in_bbox(west, south, east, north, z):
    x0, y0 = tile_xy(west, north, z)   # NW corner -> smallest x, smallest y
    x1, y1 = tile_xy(east, south, z)   # SE corner -> largest x, largest y
    for x in range(x0, x1 + 1):
        for y in range(y0, y1 + 1):
            yield z, x, y

Because the tile count over a region scales as 4^z, enumerating deep zooms over a continent produces millions of addresses; cap the zoom to what you actually serve and generate the list lazily. The exact forward and inverse formulas this section leans on are given as a standalone module in the lat/lon to tile-number conversion guide.

Troubleshooting

1. Entire basemap is upside-down (north and south swapped)

The classic XYZ/TMS mismatch — a client requesting XYZ is being served TMS-ordered rows, or vice versa.

# Compare a known tile against a reference. Berlin at z12 is XYZ y=1343, TMS y=2752.
python3 -c "z=12; y=1343; print('tms_y =', (2**z - 1) - y)"   # -> 2752

Fix by applying y_tms = 2^z − 1 − y at the MBTiles boundary, or set "scheme": "tms" on the MapLibre source if the origin is genuinely south-west.

2. Features shifted thousands of kilometres or clustered near 0°,0°

EPSG:4326 coordinates are being treated as if they were already EPSG:3857 meters (or the reverse). Degrees fed into a meters-based grid collapse the whole dataset toward the origin.

ogrinfo -al -so input.geojson | grep -Ei "srs|geogcrs|projcrs"
# Then force the source CRS before tiling:
ogr2ogr -t_srs EPSG:4326 -makevalid fixed_4326.geojson input.geojson

3. Nothing renders above ~85° latitude

This is correct behavior, not a bug. Web Mercator truncates at ±85.0511°; polar features have no tile. If you must show them, the dataset needs a polar projection and a non-Mercator tiling scheme — a different pipeline entirely.

4. Off-by-one tile at the far east or south edge

At exactly lon = +180° or the southern latitude limit, int(norm * 2^z) yields 2^z, which is one past the last valid index 2^z − 1.

n = 2 ** z
x = min(int(nx * n), n - 1)   # the min() clamp prevents the off-by-one

Always clamp the floored index to n − 1 so edge coordinates fall into the last real tile rather than a phantom column.

Further Reading

Converting Lat/Lon to Slippy-Map Tile Numbers — the exact forward and inverse formulas as a self-contained Python module (deg2num, num2deg, tile bbox, xyz_to_tms), with a worked city example and a check against a live tile URL.

Parent: Vector Tile Architecture & Format Fundamentals

Related

• MVT Encoding Internals: The Protobuf Tile Structure — how the 0..4096 tile-local extent grid quantizes geometry, and the delta-encoded command sequences that store it.
• MBTiles Architecture & Limits — why MBTiles stores its tile_row in TMS order and where the y-flip must be applied on read.
• Zoom Level Optimization Strategies — turning the 2^z grid-size math into a concrete minimum and maximum zoom for a dataset.