Converting Lat/Lon to Slippy-Map Tile Numbers

To find the tile covering a point, compute x = floor((lon + 180) / 360 × 2^z) and derive y from the Web Mercator latitude formula y = floor((1 − asinh(tan(lat)) / π) / 2 × 2^z) — both using radians for the trig and both flooring to an integer. That pair (z, x, y) is exactly the XYZ address in a {z}/{x}/{y}.pbf request. This page gives the forward transform, its inverse (tile → bounding box), the TMS flip, and a copy-pasteable Python module that ties them together against a real coordinate.

When to Use This

Reach for this conversion whenever you need to go from a geographic point to a specific tile rather than letting a map client do it implicitly:

Which tile covers this point — locating the single tile a POI, incident, or address falls in.
Pre-warming a CDN or cache — enumerating the {z}/{x}/{y} set over a bounding box to prime edge caches before launch.
Testing a specific tile — pulling one .pbf to inspect its contents while debugging a build, instead of clicking around a live map.
Validating a pipeline — asserting that a known landmark lands in the tile you expect after regenerating a tileset.

The concepts behind these formulas — the EPSG:3857 grid, the extent space, and the y-axis conventions — are covered in the parent tile coordinate systems guide.

Specification Detail

Forward: lon/lat → tile x/y (XYZ, north origin)

At zoom z the grid is n = 2^z tiles per axis:

text

x = floor( (lon + 180) / 360 · n )
y = floor( (1 − asinh(tan(lat_rad)) / π) / 2 · n )

x = floor( (lon + 180) / 360 · n )
y = floor( (1 − asinh(tan(lat_rad)) / π) / 2 · n )

asinh(tan(lat)) is mathematically identical to the more familiar ln(tan(lat) + sec(lat)) Mercator form and to atanh(sin(lat)); all three require lat in radians. Latitude must be clamped to the Mercator limit ±85.05112878° before the trig, or tan blows up toward the poles. The x term needs no trig at all — longitude maps linearly across the grid, which is why only the y axis carries the Mercator nonlinearity.

Two edge cases fall out of the floor: a coordinate at exactly lon = +180° or the southern latitude limit scales to n rather than n − 1, one past the last valid index, so both results are clamped with min(index, n − 1).

Inverse: tile x/y → lon/lat bounding box

The north-west corner of tile (x, y, z):

text

lon = x / n · 360 − 180
lat = degrees( atan( sinh( π · (1 − 2·y / n) ) ) )

lon = x / n · 360 − 180
lat = degrees( atan( sinh( π · (1 − 2·y / n) ) ) )

Evaluating this at (x, y) gives the tile’s NW corner and at (x+1, y+1) its SE corner — together the full geographic bbox.

TMS flip

text

y_tms = n − 1 − y_xyz        (self-inverse)

y_tms = n − 1 − y_xyz        (self-inverse)

XYZ is what URLs and PMTiles use; TMS is what MBTiles rows and some GDAL tooling use.

Production Command Block

A self-contained module with the forward transform, the inverse, a bbox helper, and the TMS flip, plus a worked example and a live-tile check:

python

#!/usr/bin/env python3
"""slippytiles.py — lon/lat <-> slippy-map tile numbers (Web Mercator, EPSG:3857)."""
import math

MAX_LAT = 85.05112878  # Web Mercator latitude limit

def deg2num(lon: float, lat: float, z: int) -> tuple[int, int]:
    """Longitude/latitude (EPSG:4326) -> XYZ tile (x, y) at zoom z."""
    lat = max(min(lat, MAX_LAT), -MAX_LAT)          # clamp before trig
    lat_rad = math.radians(lat)
    n = 2 ** z
    x = int((lon + 180.0) / 360.0 * n)
    y = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
    return min(x, n - 1), min(y, n - 1)             # clamp east/south edges

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

def tile_bbox(x: int, y: int, z: int) -> tuple[float, float, float, float]:
    """XYZ tile -> (west, south, east, north) geographic bbox."""
    west, north = num2deg(x, y, z)
    east, south = num2deg(x + 1, y + 1, z)
    return west, south, east, north

def xyz_to_tms(x: int, y: int, z: int) -> tuple[int, int]:
    """Flip the y-axis between XYZ (north origin) and TMS (south origin)."""
    return x, (2 ** z - 1) - y

if __name__ == "__main__":
    # Worked example: Berlin, Germany at zoom 12
    lon, lat, z = 13.404954, 52.520008, 12
    x, y = deg2num(lon, lat, z)
    print(f"XYZ tile: z={z} x={x} y={y}")            # z=12 x=2200 y=1343
    print(f"TMS  y  : {xyz_to_tms(x, y, z)[1]}")     # 2752
    print(f"bbox    : {tile_bbox(x, y, z)}")
    # -> (13.359375, 52.4827..., 13.4472..., 52.5327...)  — Berlin falls inside

#!/usr/bin/env python3
"""slippytiles.py — lon/lat <-> slippy-map tile numbers (Web Mercator, EPSG:3857)."""
import math

MAX_LAT = 85.05112878  # Web Mercator latitude limit

def deg2num(lon: float, lat: float, z: int) -> tuple[int, int]:
    """Longitude/latitude (EPSG:4326) -> XYZ tile (x, y) at zoom z."""
    lat = max(min(lat, MAX_LAT), -MAX_LAT)          # clamp before trig
    lat_rad = math.radians(lat)
    n = 2 ** z
    x = int((lon + 180.0) / 360.0 * n)
    y = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
    return min(x, n - 1), min(y, n - 1)             # clamp east/south edges

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

def tile_bbox(x: int, y: int, z: int) -> tuple[float, float, float, float]:
    """XYZ tile -> (west, south, east, north) geographic bbox."""
    west, north = num2deg(x, y, z)
    east, south = num2deg(x + 1, y + 1, z)
    return west, south, east, north

def xyz_to_tms(x: int, y: int, z: int) -> tuple[int, int]:
    """Flip the y-axis between XYZ (north origin) and TMS (south origin)."""
    return x, (2 ** z - 1) - y

if __name__ == "__main__":
    # Worked example: Berlin, Germany at zoom 12
    lon, lat, z = 13.404954, 52.520008, 12
    x, y = deg2num(lon, lat, z)
    print(f"XYZ tile: z={z} x={x} y={y}")            # z=12 x=2200 y=1343
    print(f"TMS  y  : {xyz_to_tms(x, y, z)[1]}")     # 2752
    print(f"bbox    : {tile_bbox(x, y, z)}")
    # -> (13.359375, 52.4827..., 13.4472..., 52.5327...)  — Berlin falls inside

Then confirm the tile actually exists on a running server (Web Mercator tiles are served in XYZ):

bash

# Fetch the exact tile the module computed and confirm a non-empty MVT body
curl -s -o berlin.pbf "https://tiles.example.com/basemap/12/2200/1343.pbf"
ls -l berlin.pbf     # a populated urban tile is typically tens to hundreds of KB

# Fetch the exact tile the module computed and confirm a non-empty MVT body
curl -s -o berlin.pbf "https://tiles.example.com/basemap/12/2200/1343.pbf"
ls -l berlin.pbf     # a populated urban tile is typically tens to hundreds of KB

Interaction Effects

The (z, x, y) this produces is the direct input to a tile fetch — drop x and y straight into a curl of /{z}/{x}/{y}.pbf as shown above to pull one tile for inspection. From there, feed the downloaded .pbf into a decoder to read its features: the decoding MVT tiles to GeoJSON in Python walkthrough turns that byte payload back into inspectable geometry, closing the loop from coordinate to tile to feature.

The choice of z interacts directly with how deep your tileset was built. Requesting a z above the tileset’s --maximum-zoom returns nothing unless overzooming is configured, so pair this conversion with calculating the optimal max zoom for urban datasets when deciding which zoom to test at. And if the backend is MBTiles rather than PMTiles, apply xyz_to_tms() before querying the tiles table — MBTiles stores tile_row in south-origin TMS order.

To enumerate every tile over a region — the core of cache pre-warming — convert the two bbox corners and iterate the inclusive index range between them:

python

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

# Warm a small metro extent at z14 before launch
for z, x, y in tiles_in_bbox(13.09, 52.34, 13.76, 52.68, 14):
    ...  # curl or HEAD each /{z}/{x}/{y}.pbf to prime the CDN edge cache

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

# Warm a small metro extent at z14 before launch
for z, x, y in tiles_in_bbox(13.09, 52.34, 13.76, 52.68, 14):
    ...  # curl or HEAD each /{z}/{x}/{y}.pbf to prime the CDN edge cache

Note that the NW corner yields the smallest y because XYZ counts rows southward from the north — reversing the corners silently produces an empty range.

Performance Impact

The math is trivial — a handful of floating-point operations per point, effectively free next to any network or disk cost. Two directional notes:

Bulk conversion: for millions of points (pre-warming a continental bbox, indexing a point dataset), vectorize with NumPy — replace the scalar math calls with numpy.radians, numpy.arctanh(numpy.sin(...)), and numpy.floor, and process whole arrays at once for a 50–100× speedup over a Python loop.
Enumeration cost: the number of tiles over a bbox grows as 4^z, so pre-warming a city at z=16 can be tens of thousands of tiles. Enumerate lazily and cap the zoom to what you actually serve rather than materializing the full list.
Caching the constant: 2 ** z and math.pi are the only per-call constants; when converting many points at one zoom, hoist n = 2 ** z out of the loop so it is computed once rather than per point.

For reference, the widely used mercantile library implements exactly these formulas; dropping in mercantile.tile(lon, lat, z) gives the same (x, y) as deg2num above and is worth adopting once you need adjacent helpers like tile parents, children, and Web Mercator bounds.

Common Mistakes

1. Feeding degrees where radians are required

math.tan() and math.asinh() expect radians. Passing lat in degrees produces a y that is wrong by a large, silent margin — the point lands in a plausible-looking but incorrect tile.

python

lat_rad = math.radians(lat)          # correct
# y = ... math.tan(lat) ...          # WRONG: lat still in degrees

lat_rad = math.radians(lat)          # correct
# y = ... math.tan(lat) ...          # WRONG: lat still in degrees

2. Forgetting to floor (or flooring after the flip)

The tile index is int(...) of the scaled value; skipping the floor yields a fractional “tile” that breaks any URL or lookup. Always floor to get the integer index, and apply the TMS flip only to the already-floored integer, never to the float.

3. XYZ vs TMS confusion, or latitude past the Mercator limit

Querying an MBTiles tiles table with an XYZ y returns the wrong row (or none) — apply xyz_to_tms() first. Separately, a latitude beyond ±85.0511° sends tan(lat_rad) toward infinity and yields a nonsensical y; the MAX_LAT clamp in deg2num prevents both the overflow and an out-of-range tile.

bash

# Symptom of a missed flip: the tile exists in XYZ but the MBTiles query is empty
sqlite3 basemap.mbtiles \
  "SELECT length(tile_data) FROM tiles
   WHERE zoom_level=12 AND tile_column=2200 AND tile_row=2752;"   # TMS row, not 1343

# Symptom of a missed flip: the tile exists in XYZ but the MBTiles query is empty
sqlite3 basemap.mbtiles \
  "SELECT length(tile_data) FROM tiles
   WHERE zoom_level=12 AND tile_column=2200 AND tile_row=2752;"   # TMS row, not 1343

Parent: Tile Coordinate Systems & the Slippy Map Grid

Decoding MVT Tiles to GeoJSON in Python — take the .pbf this conversion locates and turn its protobuf body back into inspectable GeoJSON features.
Calculating Optimal Max Zoom for Urban Datasets — decide which z is worth requesting so the tile numbers you compute actually resolve to generated tiles.