Converting OSM Geometries to Quadkeys

Attach a tile-aligned cell key to OSM features so containment, roll-up and range queries become plain string operations — and know what the latitude distortion costs.

Prerequisites

Conceptual minimum

A quadkey is the slippy-map tile address written as a base-4 string. Zoom level 3 tile (5, 3) becomes "213" — three digits, one per zoom level, each naming which quadrant of its parent the tile occupies.

Deriving a quadkey from a longitude and latitude A four-stage chain. Longitude and latitude in EPSG:4326 are the input. They become tile x and y at a zoom level through the Web Mercator tile grid, taking the floor of a projection. The x and y bits are interleaved into a Morton code. The code is read as base-4 digits, one per zoom level, giving the quadkey string. Interleaving is the whole trick lon, lat degrees, EPSG:4326 the input tile x, y at zoom z Web Mercator tile grid floor of a projection interleave the bits x and y, alternating a Morton code base-4 digits one per zoom level the quadkey string Because each digit is one zoom level, a quadkey is its own prefix tree: truncating the string zooms out, and a string comparison is a containment test.
One digit per zoom level is what makes the string a prefix tree — truncate it to zoom out, compare prefixes to test containment.

That one-digit-per-level structure is the whole reason to use them. Truncating a quadkey zooms out; testing containment is child.startswith(parent); a range query over a prefix is a BETWEEN on a string column. None of that needs a spatial library, an index type, or a database extension.

Quadkeys compared with H3 cells A grid of five properties. A quadkey maps exactly to a map tile while an H3 cell has no tile relationship. Quadkey containment is a string prefix test while H3 needs its library. Quadkey ground area shrinks with the cosine of latitude while H3 area is near-uniform worldwide. Quadkey neighbour distances differ between diagonal and orthogonal directions while hexagons are uniform. Both support range queries on one column. What a quadkey gives you that a cell id does not quadkey H3 cell maps to a map tile exactly — it is the tile no relationship containment by string prefix yes — startswith() no — needs the H3 library equal ground area no — shrinks as cos(lat) near-uniform worldwide neighbour distance uniform no — diagonal ≠ orthogonal yes — hexagons range query on one column yes — BETWEEN on a string yes, at one resolution Quadkeys win when the answer has to line up with tiles or be queried with plain string operations. H3 wins when the cells have to be comparable to each other.
The prefix property is the reason to reach for quadkeys: containment, roll-up and range queries all become plain string operations in any database.

The cost is that a quadkey cell is a Web Mercator tile, and Web Mercator tiles shrink toward the poles.

Ground width of a zoom-14 quadkey cell by latitude A bar chart of the ground width of one zoom-14 quadkey cell. At the equator and in Nairobi at one degree south it is 2445 metres. In Cairo at 30 degrees north it is 2118 metres, 87 percent of the equatorial width. In Berlin at 52 degrees north it is 1505 metres, 62 percent. In Reykjavík at 64 degrees north it is 1072 metres, 44 percent. Cell size at the latitudes that matter ground width of one zoom-14 quadkey cell equator (0°) 2 445 m Nairobi (1°S) 2 445 m Cairo (30°N) 2 118 m · 87% of the equator Berlin (52°N) 1 505 m · 62% Reykjavík (64°N) 1 072 m · 44% A density computed per quadkey cell is not comparable between Berlin and Nairobi, because the Berlin cell covers well under half the ground.
This is the cost of tile alignment. If cells are being compared to each other rather than to tiles, it is the wrong scheme.

Runnable solution

python 
#!/usr/bin/env python3
"""Attach quadkeys to OSM geometries: point cells, and covers for extents."""
from __future__ import annotations

import logging
import math
from typing import Iterator

import numpy as np
from shapely.geometry.base import BaseGeometry

logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
logger = logging.getLogger(__name__)

MAX_LATITUDE = 85.05112878          # where Web Mercator is clipped


def lonlat_to_tile(lon: float, lat: float, zoom: int) -> tuple[int, int]:
    """Web Mercator tile containing a point. Latitude is clamped, not wrapped."""
    lat = max(-MAX_LATITUDE, min(MAX_LATITUDE, lat))
    n = 1 << zoom
    x = int((lon + 180.0) / 360.0 * n)
    sin_lat = math.sin(math.radians(lat))
    y = int((0.5 - math.log((1 + sin_lat) / (1 - sin_lat)) / (4 * math.pi)) * n)
    # A point exactly on the eastern or northern edge lands one tile out.
    return min(x, n - 1), min(y, n - 1)


def tile_to_quadkey(x: int, y: int, zoom: int) -> str:
    """Interleave the x and y bits, most significant first, as base-4 digits."""
    digits: list[str] = []
    for level in range(zoom, 0, -1):
        mask = 1 << (level - 1)
        digit = 0
        if x & mask:
            digit += 1
        if y & mask:
            digit += 2
        digits.append(str(digit))
    return "".join(digits)


def quadkey_to_tile(quadkey: str) -> tuple[int, int, int]:
    """The inverse — useful for turning a stored key back into a tile URL."""
    x = y = 0
    zoom = len(quadkey)
    for index, char in enumerate(quadkey):
        mask = 1 << (zoom - index - 1)
        digit = int(char)
        if digit & 1:
            x |= mask
        if digit & 2:
            y |= mask
    return x, y, zoom


def point_quadkey(lon: float, lat: float, zoom: int) -> str:
    return tile_to_quadkey(*lonlat_to_tile(lon, lat, zoom), zoom)


def cover(geom: BaseGeometry, zoom: int, max_cells: int = 4096) -> list[str]:
    """Every quadkey whose tile intersects the geometry's bounds.

    Bounds, not the geometry: this is a coarse cover, meant as a candidate filter.
    Refining it against the real geometry is the caller's job.
    """
    west, south, east, north = geom.bounds
    x0, y1 = lonlat_to_tile(west, south, zoom)      # south → larger y
    x1, y0 = lonlat_to_tile(east, north, zoom)
    count = (x1 - x0 + 1) * (y1 - y0 + 1)
    if count > max_cells:
        raise ValueError(
            f"cover would be {count} cells at zoom {zoom}; "
            f"use a coarser zoom or cover the parts separately")
    return [tile_to_quadkey(x, y, zoom)
            for x in range(x0, x1 + 1)
            for y in range(y0, y1 + 1)]


def parent(quadkey: str, levels: int = 1) -> str:
    """Zoom out by truncation — the property the whole scheme exists for."""
    if levels >= len(quadkey):
        return ""
    return quadkey[:-levels]


def contains(parent_key: str, child_key: str) -> bool:
    """Containment as a string operation, no geometry involved."""
    return child_key.startswith(parent_key)


def quadkeys_for_points(lons: np.ndarray, lats: np.ndarray, zoom: int) -> list[str]:
    """Vectorised tile arithmetic; only the digit assembly stays in Python."""
    lats = np.clip(lats, -MAX_LATITUDE, MAX_LATITUDE)
    n = 1 << zoom
    xs = np.clip(((lons + 180.0) / 360.0 * n).astype(np.int64), 0, n - 1)
    sin_lat = np.sin(np.radians(lats))
    ys = np.clip(
        ((0.5 - np.log((1 + sin_lat) / (1 - sin_lat)) / (4 * np.pi)) * n).astype(np.int64),
        0, n - 1)
    return [tile_to_quadkey(int(x), int(y), zoom) for x, y in zip(xs, ys)]

Using it as a partition and filter key:

sql 
-- Everything inside a zoom-8 area, with no spatial index at all.
SELECT count(*) FROM features WHERE quadkey LIKE '02313021%';

-- The same as a range scan, which an ordinary B-tree serves.
SELECT count(*) FROM features
WHERE quadkey >= '02313021' AND quadkey < '02313022';

-- Roll up to zoom 10 for a density surface.
SELECT left(quadkey, 10) AS cell, count(*) FROM features GROUP BY 1;

Step-by-step walkthrough

lonlat_to_tile clamps latitude to the Web Mercator limit rather than letting the logarithm run away. Beyond about 85.05 degrees the projection goes to infinity, and OSM does contain nodes above that — research stations, sea-ice features — so an unclamped implementation raises or produces a nonsensical tile on real data.

Both coordinates are also clamped to n - 1. A point exactly on the eastern or northern edge of the world computes a tile index one past the last valid one, which is a one-in-a-million input that appears reliably in a planet-scale run.

tile_to_quadkey walks bits from most significant to least, which is what makes the resulting string a prefix tree: the first digit describes the coarsest subdivision, so shorter strings are ancestors of longer ones.

cover refuses to produce an unbounded number of cells. A bounding box covering a country at zoom 16 is tens of millions of tiles, and returning that list is never what the caller wanted — the guard turns a memory exhaustion into a clear error naming the fix.

contains and parent are string operations with no geometry, which is the payoff. Compare with the same operations on H3, which need the library on both sides of the query.

quadkeys_for_points vectorises the arithmetic but assembles digits in a Python loop, because base-4 string construction does not vectorise usefully in numpy. For very large batches, storing the interleaved integer instead of the string keeps everything in numpy — the string is only needed where prefix operations are.

Verification

Check the round trip and the prefix property, which together validate the bit interleaving:

python 
def test_round_trip():
    for lon, lat, z in ((13.405, 52.520, 14), (-74.006, 40.713, 18), (0.0, 0.0, 1)):
        qk = point_quadkey(lon, lat, z)
        x, y, zoom = quadkey_to_tile(qk)
        assert (x, y) == lonlat_to_tile(lon, lat, z) and zoom == z

def test_prefix_is_containment():
    child = point_quadkey(13.405, 52.520, 16)
    for levels in range(1, 8):
        assert contains(parent(child, levels), child)

def test_poles_are_clamped():
    assert point_quadkey(0.0, 89.9, 10) == point_quadkey(0.0, MAX_LATITUDE, 10)

Then verify against an external reference — a quadkey is a Bing Maps tile address, so it can be checked visually:

python 
print(point_quadkey(13.405, 52.520, 14))   # → 12022332303121
# https://www.bing.com/maps?... or any quadkey→tile viewer

Finally, sanity-check the cell size at your working latitude before committing to a zoom level:

python 
def cell_width_m(lat: float, zoom: int) -> float:
    return 40_075_016.686 * math.cos(math.radians(lat)) / (1 << zoom)

for lat in (0, 30, 52, 64):
    print(f"{lat:>3}°  {cell_width_m(lat, 14):7.0f} m")

Common errors and fixes

Symptom Root cause Fix
math domain error near the poles Latitude not clamped Clamp to ±85.05112878
Tile index one past the maximum Point exactly on the world edge Clamp x and y to n - 1
Prefix test fails for a known parent Bits interleaved least-significant first Walk levels from zoom down to 1
Cover returns millions of cells Zoom too fine for the extent Guard the count; use a coarser zoom
Densities differ between cities Cell area shrinks with latitude Use an equal-area scheme for comparison
Keys sort oddly in the database Stored as an integer, losing leading zeros Store as a fixed-length string

Frequently Asked Questions

Quadkey or H3?

Quadkey when the cells must line up with map tiles, or when you want containment and roll-up as string operations in a database with no extensions. H3 when the cells are compared with each other — densities, coverage percentages, anything aggregated across latitudes — because equal area is exactly what quadkeys do not offer.

Should I store the string or the integer?

The string, if you intend to use prefix operations, and pad it to a fixed length so lexical order matches spatial order. The interleaved integer is more compact and supports range queries just as well, but LIKE 'prefix%' is far more readable than the equivalent bit arithmetic, and readability wins in a column other people will query.

What zoom level should I use?

Pick from cell size at your working latitude, not from the zoom number. Zoom 14 is roughly 2.4 km at the equator and 1.5 km in Berlin; zoom 16 is roughly 600 m and 380 m. Choose the coarsest level at which the smallest thing you need to distinguish still lands in its own cell.

Can a quadkey index a polygon rather than a point?

Only as a cover — a set of keys whose tiles intersect it. That is a legitimate coarse filter, and it is exactly the filter half of the filter-then-refine pattern in Building an R-tree Index over OSM Geometries. Storing a cover means one row per cell per polygon, which grows quickly; storing a single key for the polygon’s centroid does not answer containment questions at all.

Specification reference

A quadkey of length z addresses a Web Mercator tile at zoom z. Digit i (from the left) encodes the tile’s quadrant at level i + 1: bit 0 of the digit is the x bit, bit 1 the y bit, giving values 0 through 3. Tile (x, y, z) is derived from longitude and latitude with the standard slippy-map formulae, with latitude clamped to ±85.05112878 degrees.

Up one level: Spatial Index Selection: R-tree, H3 or Quadkey.