โš™๏ธ How-To H3 & global grids Web maps & tiles Performance & scale H3 & Global Grids ยท step 14 of 20 Updated

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,464,117 GeoNames points in a Parquet file, sorting rows by a zoom-16 quadkey cut a London bounding-box query from 10.5 ms to 2.8 ms and the row groups it had to open from 35 of 135 to 5, without changing the query at all.

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.083 s against 2.77 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.

Diagram of a quadkey: each zoom level splits a tile into four quadrants numbered 0 to 3, and digits are appended from coarse to fine.
A quadkey's length is its zoom, and every prefix is the tile that contains it.

Step-by-step solution

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ยฐ. Seven of those raised inside the first million rows. 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 1,000,000 points at zoom 16: 2.77 s for the mercantile loop, 0.083 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 1.4 ms, matching a direct x/y check exactly.

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               298.3 MB    10.5 ms             35 of 135
sorted by quadkey      277.8 MB     2.8 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.

Bar chart of the ground area of one zoom-10 web map tile at latitudes from 0 to 85 degrees.
The same tile address covers 130 times more ground on the equator than at 85ยฐ north.

Code examples

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:

(-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
>>> quadkey_range("03131313", 16)
(930545664, 930611199)
>>> quadkey_range("0313131311301", 12)
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)
wrote 399,644,572 bytes in 5.3s
tile 03131313: 8,661 points in 1.6 ms

con is a DuckDB connection. This run also carried a name column, hence the larger file than step 6's. The read_all check covers DuckDB versions whose .arrow() returns a record batch reader rather than a table.

Explanation

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 35 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,399,925 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.

Grid comparing an unsorted and a quadkey-sorted Parquet file on size, London query time and row groups read.
The sort cost one pass at write time; every later spatial filter gets it for free.

Edge cases or notes

  • 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.

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 1.4 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.