How to Index Data by Quadkey and Web Map Tile

Problem statement

Web maps already divide the world into a grid. Every slippy map, tile server and PMTiles archive uses the same pyramid of Web Mercator tiles addressed by zoom, x and y — and a quadkey writes that address as a single string or integer in which each digit is one zoom level.

That makes quadkeys an attractive data index: points keyed by the tile that will draw them can be counted per tile, filtered by prefix and sorted so that nearby rows sit together on disk. Measured on 13,472,195 GeoNames points in a Parquet file, sorting rows by a zoom-16 quadkey cut a London bounding-box query from 9.0 ms to 3.2 ms and the row groups it had to open from 33 of 135 to 5, without changing the query at all (Example 4).

The same grid carries two traps. Web Mercator tiles are not equal area — a zoom-10 tile at 60° north covers 25.3% of the ground a tile at the equator does — and the projection stops at ±85.0511°, so points near the poles need handling before they are indexed.

Quick answer

import mercantile

tile = mercantile.tile(-0.1246, 51.5007, 12)       # longitude first
key = mercantile.quadkey(tile)

print(tile, key)
print(mercantile.parent(tile), key[:-1])            # the parent is the key minus a digit
print(mercantile.quadkey_to_tile(key))

Tile(x=2046, y=1362, z=12) 031313131130
Tile(x=1023, y=681, z=11) 03131313113
Tile(x=2046, y=1362, z=12)

For more than a few thousand points, compute the keys with NumPy rather than a mercantile loop — 0.056 s against 2.38 s for a million points, with identical results — and store them as unsigned integers so a tile at any coarser zoom becomes a between filter.

Longitude first, digits, polar limits, integers

1. Pass longitude first

mercantile.tile(lng, lat, zoom) follows the x, y convention of the rest of the web mapping stack. That is the opposite of H3, pygeohash and S2, which take latitude first — a pipeline that computes both tiles and cells is a common place for a swap. Name the arguments at the call site when both appear in the same module.

2. Understand the digits

At each zoom the parent tile splits into four. The digit is x_bit + 2 × y_bit: 0 is the top-left child, 1 top-right, 2 bottom-left and 3 bottom-right, because tile y increases southwards. A zoom-12 key has 12 digits, and its first eight digits, 03131313, are the zoom-8 tile containing it.

So "all data inside this tile" is a prefix test on strings and a range test on integers.

3. Handle the latitudes Web Mercator cannot show

The projection's y coordinate goes to infinity at the poles, so the tile pyramid is cut at ±85.0511°. mercantile treats the cases differently, measured:

latitude 85.06   Tile(x=8, y=0, z=4)                 clamped to the top row
latitude 89.9    Tile(x=8, y=0, z=4)                 clamped to the top row
latitude 90.0    InvalidLatitudeError: Y can not be computed: lat=90.0

GeoNames holds 18 points north of the cut-off and 553 south of it, and 8 at exactly ±90°. Clip latitudes to ±85.0511 before indexing, or drop them, and decide deliberately which — a clamp silently puts a South Pole station in the bottom row of tiles.

4. Vectorise the key computation

mercantile is correct and slow in bulk, because each call builds Python objects. The tile maths is two lines of NumPy, and the bit interleave is a loop over zoom levels rather than points (Example 1). Measured on the first million GeoNames rows (999,447 of them inside the cut-off) at zoom 16: 2.38 s for the mercantile loop, 0.056 s vectorised, and the keys agreed on all 10,000 sampled rows.

5. Store integers, not strings

A zoom-16 quadkey is 32 bits. Keep it as an unsigned integer column: it sorts in the same order as the string, compresses better and turns a prefix into a range. The zoom-8 tile 03131313 in a zoom-16 index is every key between int('03131313', 4) << 16 and that value plus 65,535 — measured, 8,661 points in 3.2 ms (Example 3).

6. Sort the file by the key

A Parquet file keeps minimum and maximum values for each column in each row group. When rows are written in quadkey order, nearby points share row groups, so those statistics become tight boxes and a spatial filter can skip most of the file:

file                   size        London bbox query   row groups overlapping
unsorted               185.0 MB     9.0 ms             33 of 135
sorted by quadkey      171.5 MB     3.2 ms              5 of 135

The query filtered on lat and lon, not on the quadkey. Sorting helped a query that never mentioned it, and made the file 7% smaller because similar coordinates compress better next to each other.

7. Do not treat tile counts as densities

Every tile at one zoom is the same size in Web Mercator and a different size on the ground. Measured for zoom 10:

latitude   0°        30°       45°       60°       70°       80°      85°
km²      1,521.3   1,145.0    766.7     384.7     180.8     46.2     11.7

A tile at 80° covers 3% of the area of an equatorial one. Divide counts by each tile's geodesic area before mapping them, or use an equal-area grid for analysis.

Vectorised tiles and quadkey-ordered Parquet

The data used on this page

The points are GeoNames' full allCountries dump — every feature, about 13.5 million — from download.geonames.org, rebuilt daily, so counts drift a little between runs. The setup keeps name, latitude and longitude in points.parquet, in the order the dump lists them, written in row groups of 100,000.

import os
import zipfile

import duckdb
import pandas as pd
import requests

if not os.path.exists("allCountries.zip"):
    with requests.get("https://download.geonames.org/export/dump/allCountries.zip",
                      stream=True, timeout=600) as reply:
        reply.raise_for_status()
        with open("allCountries.zip", "wb") as out:
            for chunk in reply.iter_content(1 << 20):
                out.write(chunk)
if not os.path.exists("points.parquet"):
    with zipfile.ZipFile("allCountries.zip") as z, z.open("allCountries.txt") as f:
        gn = pd.read_csv(f, sep="\t", header=None, usecols=[1, 4, 5], quoting=3,
                         names=["name", "lat", "lon"], keep_default_na=False)
    con = duckdb.connect()
    con.register("gn", gn)
    con.execute("copy gn to 'points.parquet' (format parquet, row_group_size 100000)")
    del gn
con = duckdb.connect()
print(f"{con.execute('select count(*) from read_parquet(\'points.parquet\')').fetchone()[0]:,} points, "
      f"{os.path.getsize('points.parquet') / 1e6:.1f} MB")

13,472,195 points, 298.7 MB

Example 1 — vectorised tiles and integer quadkeys

import math

import numpy as np

MAX_LAT = 85.0511287798066


def tiles_for_points(lng, lat, zoom):
    """Vectorised mercantile.tile: x, y and an integer quadkey for many points."""
    lng = np.asarray(lng, dtype=float)
    lat = np.clip(np.asarray(lat, dtype=float), -MAX_LAT, MAX_LAT)
    n = 2 ** zoom
    x = np.floor((lng + 180.0) / 360.0 * n).astype(np.int64)
    s = np.sin(np.radians(lat))
    y = np.floor((0.5 - np.log((1 + s) / (1 - s)) / (4 * math.pi)) * n).astype(np.int64)
    x, y = np.clip(x, 0, n - 1), np.clip(y, 0, n - 1)
    key = np.zeros(len(x), dtype=np.uint64)
    for bit in range(zoom):
        key |= ((x >> bit) & 1).astype(np.uint64) << np.uint64(2 * bit)
        key |= ((y >> bit) & 1).astype(np.uint64) << np.uint64(2 * bit + 1)
    return x, y, key


def quadkey_string(key, zoom):
    return "".join(str((int(key) >> (2 * i)) & 3) for i in range(zoom - 1, -1, -1))

Checked against mercantile for five awkward points — London, Sydney, both corners of the pyramid and 89.9° north:

import mercantile

awkward = [(-0.1246, 51.5007), (151.2093, -33.8688), (-180.0, 85.0511),
           (179.9999, -85.0511), (10.0, 89.9)]
xs, ys, keys = tiles_for_points([p[0] for p in awkward], [p[1] for p in awkward], 12)
for (lng, lat), x, y, k in zip(awkward, xs, ys, keys):
    t = mercantile.tile(lng, lat, 12)
    assert (t.x, t.y) == (x, y) and mercantile.quadkey(t) == quadkey_string(k, 12)
    print(f"{str((lng, lat)):21} {str((int(x), int(y))):13} {quadkey_string(k, 12)}")

(-0.1246, 51.5007)    (2046, 1362)  031313131130
(151.2093, -33.8688)  (3768, 2457)  311230133002
(-180.0, 85.0511)     (0, 0)        000000000000
(179.9999, -85.0511)  (4095, 4095)  333333333333
(10.0, 89.9)          (2161, 0)     100001110001

Every row matched mercantile.tile and mercantile.quadkey. The final clip on x and y handles longitude exactly 180°, which would otherwise produce a tile one past the edge.

Example 2 — a coarse tile as an integer range

def quadkey_range(quadkey, index_zoom):
    """Integer bounds of every zoom-`index_zoom` key inside a coarser tile."""
    shift = 2 * (index_zoom - len(quadkey))
    if shift < 0:
        raise ValueError("the tile is finer than the stored index")
    low = int(quadkey, 4) << shift
    return low, low + (1 << shift) - 1

print(quadkey_range("03131313", 16))
try:
    quadkey_range("0313131311301", 12)
except ValueError as e:
    print("ValueError:", e)

(930545664, 930611199)
ValueError: the tile is finer than the stored index

int(quadkey, 4) reads the key as a base-4 number, which is exactly what the digits are. The range works in any engine with an integer column: where qk16 between 930545664 and 930611199.

Example 3 — write Parquet in quadkey order

import os

import pyarrow as pa


def write_sorted_by_quadkey(con, source_sql, path, zoom=16, row_group_size=100_000):
    """Add a quadkey column and write Parquet in quadkey order."""
    table = con.execute(source_sql).arrow()
    if hasattr(table, "read_all"):
        table = table.read_all()
    _, _, key = tiles_for_points(table.column("lon").to_numpy(),
                                 table.column("lat").to_numpy(), zoom)
    table = table.append_column(f"qk{zoom}", pa.array(key, type=pa.uint64()))
    con.register("with_keys", table)
    con.execute(f"""copy (select * from with_keys order by qk{zoom})
                   to '{path}' (format parquet, row_group_size {row_group_size})""")
    con.unregister("with_keys")
    return os.path.getsize(path)


import time

t = time.perf_counter()
size = write_sorted_by_quadkey(con, "select * from read_parquet('points.parquet')", "sorted.parquet")
print(f"wrote {size:,} bytes in {time.perf_counter() - t:.1f}s")
low, high = quadkey_range("03131313", 16)
t = time.perf_counter()
n = con.execute(f"select count(*) from read_parquet('sorted.parquet') where qk16 between {low} and {high}").fetchone()[0]
print(f"tile 03131313: {n:,} points in {(time.perf_counter() - t) * 1000:.1f} ms")

wrote 342,744,065 bytes in 4.5s
tile 03131313: 8,661 points in 3.2 ms

con is a DuckDB connection. The read_all check covers DuckDB versions whose .arrow() returns a record batch reader rather than a table.

Example 4 — the other numbers on this page

import pyarrow.parquet as pq
from pyproj import Geod

pts = con.execute("select lat, lon from read_parquet('points.parquet')").df()
print(f"north of the cut-off: {int((pts.lat > MAX_LAT).sum())}, south of it: {int((pts.lat < -MAX_LAT).sum())}, "
      f"at exactly ±90: {int((pts.lat.abs() == 90).sum())}")
for lat in (85.06, 89.9, 90.0):
    try:
        print(f"latitude {lat:5}  {mercantile.tile(10.0, lat, 4)}")
    except Exception as e:
        print(f"latitude {lat:5}  {type(e).__name__}: {e}")

first = pts.iloc[:1_000_000]
inside = first[first.lat.abs() <= MAX_LAT]
t = time.perf_counter()
loop = [mercantile.quadkey(mercantile.tile(lo, la, 16)) for la, lo in zip(inside.lat, inside.lon)]
t_loop = time.perf_counter() - t
t = time.perf_counter()
_, _, vec = tiles_for_points(inside.lon, inside.lat, 16)
t_vec = time.perf_counter() - t
pick = np.random.default_rng(0).choice(len(inside), 10_000, replace=False)
same = all(loop[i] == quadkey_string(vec[i], 16) for i in pick)
print(f"{len(inside):,} points at zoom 16: mercantile loop {t_loop:.2f} s, vectorised {t_vec:.3f} s, "
      f"10,000 sampled keys agree: {same}; load average {os.getloadavg()[0]:.1f}")

geod = Geod(ellps="WGS84")
areas = {}
for lat in (0, 30, 45, 60, 70, 80, 85):
    b = mercantile.bounds(mercantile.tile(0.0, lat, 10))
    areas[lat] = abs(geod.polygon_area_perimeter([b.west, b.east, b.east, b.west],
                                                 [b.south, b.south, b.north, b.north])[0]) / 1e6
print("zoom-10 tile area, km²:", {k: round(v, 1) for k, v in areas.items()})
print(f"60° vs equator: {areas[60] / areas[0]:.1%}; 80°: {areas[80] / areas[0]:.1%}; "
      f"equator vs 85°: {areas[0] / areas[85]:.0f}×")

_, _, k12 = tiles_for_points(pts.lon, pts.lat, 12)
counts = pd.Series(k12).value_counts()
top = counts.iloc[:max(1, len(counts) // 100)].sum() / len(pts)
print(f"zoom 12: {len(counts):,} occupied tiles of {4 ** 12:,} ({len(counts) / 4 ** 12:.1%}); "
      f"busiest 1% hold {top:.1%} of points")

con.execute("copy (select lat, lon from read_parquet('points.parquet')) to 'unsorted.parquet' "
            "(format parquet, row_group_size 100000)")
con.execute("copy (select lat, lon from read_parquet('sorted.parquet')) to 'sorted_ll.parquet' "
            "(format parquet, row_group_size 100000)")
LONDON = (51.28, 51.70, -0.51, 0.33)
query = ("select count(*) from read_parquet('{}') "
         f"where lat between {LONDON[0]} and {LONDON[1]} and lon between {LONDON[2]} and {LONDON[3]}")
for label, path in (("unsorted", "unsorted.parquet"), ("sorted by quadkey", "sorted_ll.parquet")):
    runs = []
    for _ in range(7):
        t = time.perf_counter()
        n = con.execute(query.format(path)).fetchone()[0]
        runs.append(time.perf_counter() - t)
    meta = pq.ParquetFile(path).metadata
    hit = 0
    for g in range(meta.num_row_groups):
        st = {meta.row_group(g).column(c).path_in_schema: meta.row_group(g).column(c).statistics
              for c in range(meta.row_group(g).num_columns)}
        if (st["lat"].min <= LONDON[1] and st["lat"].max >= LONDON[0]
                and st["lon"].min <= LONDON[3] and st["lon"].max >= LONDON[2]):
            hit += 1
    print(f"{label:18} {os.path.getsize(path) / 1e6:6.1f} MB  {np.median(runs) * 1000:5.1f} ms  "
          f"{hit} of {meta.num_row_groups} row groups overlap  ({n:,} points)")

north of the cut-off: 18, south of it: 553, at exactly ±90: 8
latitude 85.06  Tile(x=8, y=0, z=4)
latitude  89.9  Tile(x=8, y=0, z=4)
latitude  90.0  InvalidLatitudeError: Y can not be computed: lat=90.0
999,447 points at zoom 16: mercantile loop 2.38 s, vectorised 0.056 s, 10,000 sampled keys agree: True; load average 1.7
zoom-10 tile area, km²: {0: 1521.3, 30: 1145.0, 45: 766.7, 60: 384.7, 70: 180.8, 80: 46.2, 85: 11.7}
60° vs equator: 25.3%; 80°: 3.0%; equator vs 85°: 130×
zoom 12: 1,400,080 occupied tiles of 16,777,216 (8.3%); busiest 1% hold 15.1% of points
unsorted            185.0 MB    9.0 ms  33 of 135 row groups overlap  (5,506 points)
sorted by quadkey   171.5 MB    3.2 ms  5 of 135 row groups overlap  (5,506 points)

A prefix tree matched to the renderer

Why a quadkey is a prefix tree

Each zoom level doubles the number of tiles along both axes, so each tile has exactly four children and a quadkey appends two bits — one for x, one for y — per level. Removing the last digit removes the finest split and leaves the parent. That is the same property geohash has, on a projected square grid rather than on degrees.

It is also why the integer form sorts correctly. Interleaving the bits of x and y traces a Z-order curve, which visits all of a tile's descendants before moving on, so every tile is a contiguous run of keys.

Why sorting made an unrelated query faster

DuckDB, like other Parquet readers, checks each row group's minimum and maximum before reading it. In the unsorted file, row groups held points from all over the world, so their latitude and longitude ranges were nearly global and 33 of them overlapped London. After sorting by quadkey, points in the same tile were written together, each row group spanned a small region, and only 5 could contain a match.

The Z-order curve has jumps — two tiles adjacent on the ground can be far apart in key order — so clustering is good rather than perfect. That is why some row groups still overlapped.

Why the grid is uneven

Web Mercator stretches the map by 1 / cos(latitude) in both directions to keep angles true. A tile is a fixed square on that stretched map, so its ground area shrinks with cos²(latitude): 0.25 at 60°, 0.03 at 80°. The occupancy table shows what that does to a dataset. At zoom 12, 1,400,080 tiles held points, just 8.3% of the 16.8 million tiles in the pyramid, and the busiest 1% of occupied tiles held 15.1% of all points.

Why the index must match the renderer

A tile index pays for itself when data is going to be drawn on tiles, served per tile, or partitioned for a web map. It is the wrong index for measuring density, finding neighbours or joining to another dataset: H3 and S2 have no poles to cut off and far less area distortion.

Clamped latitudes, the 180th meridian and lost zeros

• Latitude ±90 raises in mercantile; anything between ±85.0511 and ±90 is silently clamped into the edge row.
• Longitude 180 is one tile past the edge in the raw formula; clip x to 2ᶻ − 1.
• String quadkeys lose leading zeros if a CSV reader parses them as numbers: 031313131130 becomes 31313131130.
• Zoom 16 fits in 32 bits, zoom 31 in 62. Choose the stored zoom from the finest tile you will query; coarser tiles are ranges.
• TMS numbers y from the south. A key built from TMS y is a different tile; convert with y = 2ᶻ − 1 − y.
• Z-order is not Hilbert order. Neighbouring tiles can be far apart in key order, so prune with ranges, not with "nearby keys".
• Row-group size decides how much pruning is possible. Smaller groups prune more finely and add metadata; 100,000 rows worked here.

Internal links

• Web map tiles explained: XYZ, zoom levels and the tile pyramid — the grid behind every quadkey
• My tiles are offset or in the wrong place — TMS versus XYZ and other tile address bugs
• Geohash explained: prefixes, precision and the edge problem — the same prefix idea on degrees
• S2 cells explained: squares on a cube and when they beat hexagons — a quadtree without the poles problem
• Discrete global grids explained: H3, S2, geohash and why cells beat coordinates — choosing among the grids
• How to write partitioned GeoParquet and query it with filters — pruning by directory instead of row group
• How to query remote GeoParquet over HTTP with DuckDB — where row-group pruning saves network bytes
• Fixing an H3 map that misleads: unequal cell areas and empty cells — dividing counts by area properly

FAQ

What is a quadkey?

A web map tile address written as one digit per zoom level, where each digit from 0 to 3 picks one of four child tiles. Its length is the zoom, and every prefix is a containing tile.

Does mercantile take latitude or longitude first?

Longitude first: mercantile.tile(lng, lat, zoom). H3, S2 and pygeohash take latitude first, so check the order wherever both appear.

How do I find all points inside a tile?

Store a fine quadkey as an integer and filter with a range. A zoom-8 tile in a zoom-16 index is int(key, 4) shifted left by 16 bits, plus 65,535; that returned 8,661 points in 3.2 ms.

Why does sorting by quadkey make queries faster?

Rows from the same tile land in the same Parquet row groups, so their min and max statistics describe small areas. A London bounding-box query then read 5 of 135 row groups instead of 35.

Can I compare point counts between tiles?

Not directly. Web Mercator tiles shrink on the ground towards the poles: a zoom-10 tile covers 1,521 km² on the equator and 385 km² at 60° north. Divide by geodesic area first.

What happens to points near the poles?

Web Mercator stops at 85.0511 degrees. mercantile clamps latitudes beyond that into the edge row and raises InvalidLatitudeError at exactly 90, so clip or drop polar points deliberately.