S2 Cells Explained: Squares on a Cube and When They Beat Hexagons
Problem statement
H3 gets most of the attention in Python GIS, but a lot of production geospatial infrastructure runs on a different grid. BigQuery's S2_CELLIDFROMPOINT and S2_COVERINGCELLIDS return S2 cells, and so do many location services and databases built on Google's S2 geometry library. If you receive data keyed by a 64-bit integer or a short token such as 487604d, it is probably an S2 cell.
S2 makes the opposite trade to H3. Its cells are quadrilaterals, so a cell has neighbours at two different distances. In exchange, its children tile the parent exactly. Measured over 500 random level-10 cells, the four children's areas matched the parent to within 1.29 × 10⁻¹³ — floating-point noise. H3's seven children leave 7.14% of their area outside the parent.
That single property decides most of the "S2 or H3?" question. This guide explains how S2 is built, what the measurements say about area and neighbours, and when squares on a cube are the better index.
Quick answer
S2 projects the sphere onto the six faces of a cube, divides each face into a quadtree, and numbers the cells along a Hilbert curve:
import s2sphere as s2
point = s2.LatLng.from_degrees(51.5007, -0.1246) # latitude first
leaf = s2.CellId.from_lat_lng(point) # level 30, about 0.74 cm²
cell = leaf.parent(12) # level 12, about 4 km²
print(cell.to_token(), cell.id(), cell.level())
print(cell.range_min().id(), cell.range_max().id())
487604d 5221366109382377472 12
5221366040662900737 5221366178101854207
The last line is the property that matters. Every S2 cell is a contiguous range of leaf ids, so points stored once as sorted leaf ids can be counted for any cell at any level with two binary searches. Measured on 1,000,000 points, a level-8 cell count took 5.3 microseconds.
Choose S2 when you need exact hierarchy, range scans or coverings of arbitrary regions. Choose H3 when you need uniform neighbours for smoothing, flows or hexagon maps.
Step-by-step solution
1. Start from the cube, not the map
S2 does not bisect latitude and longitude, which is what makes geohash cells shrink towards the poles. It places the sphere inside a cube and projects each point onto the nearest face. A non-linear transform then evens out the distortion that a straight gnomonic projection would leave near the face corners.
The six faces are level 0. London falls on face 2, the face centred on the north pole.
2. Quarter every cell at each level
Each level splits every cell into four, from level 0 down to level 30. The number of cells at level L is 6 × 4ᴸ, so the average area is the Earth's surface divided by that:
level cells average area
0 6 85,010,980 km²
7 98,304 5,188.7 km²
12 100,663,296 5.07 km²
30 6.9 × 10¹⁸ 0.74 cm²
The level-7 row was checked by enumerating all 98,304 cells: their mean exact area was 5,188.7 km², matching the formula.
3. Read the id as a path
A cell id is 64 bits: three for the face, then two per level describing which child was taken, then a single trailing 1 bit that marks where the path ends. A coarser cell has its trailing bit further left.
That is why range_min() and range_max() exist. Every leaf inside a cell shares the cell's path prefix, so all their ids fall between two integers. The token is the id in hexadecimal with trailing zeros removed — 487604d for London at level 12, and 487604c4 to 487604dc for its four children.
4. Expect unequal areas, but bounded ones
The cube projection keeps area variation within a fixed band everywhere on Earth, not one that grows with latitude. Measured across every cell at comparable sizes:
grid level cells min km² max km² max/min
S2 7 98,304 3,175.4 6,529.1 2.056
H3 4 288,122 896.6 2,136.0 2.382 (1.970 without pentagons)
The spread is similar to H3's, and it does not follow latitude. Three level-12 cells measured 4.303 km² on the equator, 3.804 km² in London and 5.607 km² at 70° north — cell size depends on position within a cube face rather than on distance from the equator.
5. Accept that neighbours are not equidistant
A square has four edge neighbours and four corner neighbours. Measured around London:
S2 level 13: 8 neighbours, centre distances 813–1,425 m (ratio 1.752)
H3 res 8: 6 neighbours, centre distances 836–911 m (ratio 1.090)
For a moving-window smoother, a flow model or a diffusion process, that 75% difference between near and far neighbours is a real bias. For containment and range queries it is irrelevant.
6. Use coverings for arbitrary regions
S2's RegionCoverer approximates a shape with a mix of cell levels: large cells for the interior, small ones along the edge. You set the budget with max_cells. Measured on Great Britain's bounding rectangle, 771,126 km²:
max_cells cells used levels area covered
4 4 3–4 4.921 × the rectangle
8 8 4–5 2.720 ×
20 14 4–10 1.602 ×
100 43 4–12 1.174 ×
500 365 4–13 1.022 ×
A covering always contains the region; the budget controls how much extra it includes. At 500 cells it was within 2.2%.
7. Store leaf ids once and query any level
Because a cell is a range, one sorted column of leaf ids serves every level. On 109,179 points in Great Britain, counting the points in London's cells took 37 µs at level 6 and 4–7 µs at levels 8 to 12. Nothing is re-indexed when the question changes resolution.
Code examples
Example 1 — the facts about a cell worth logging
import s2sphere as s2
EARTH_RADIUS_KM = 6371.0088
def s2_cell(lat, lng, level):
"""The S2 cell containing a point, with the facts worth logging."""
cell_id = s2.CellId.from_lat_lng(s2.LatLng.from_degrees(lat, lng)).parent(level)
cell = s2.Cell(cell_id)
centre = s2.LatLng.from_point(cell.get_center())
return {
"token": cell_id.to_token(),
"id": cell_id.id(),
"level": cell_id.level(),
"area_km2": round(cell.exact_area() * EARTH_RADIUS_KM ** 2, 3),
"centre": (round(centre.lat().degrees, 6), round(centre.lng().degrees, 6)),
"range": (cell_id.range_min().id(), cell_id.range_max().id()),
}
>>> s2_cell(51.5007, -0.1246, 12)
{'token': '487604d', 'id': 5221366109382377472, 'level': 12, 'area_km2': 3.804,
'centre': (51.504425, -0.129155), 'range': (5221366040662900737, 5221366178101854207)}
exact_area() returns steradians; multiplying by the radius squared gives km². Log the token rather than the integer when a person will read it.
Example 2 — one sorted column, every level
import numpy as np
import s2sphere as s2
class LeafIndex:
"""Points stored once as sorted leaf ids; any cell at any level is a range scan."""
def __init__(self, lats, lngs):
ids = [s2.CellId.from_lat_lng(s2.LatLng.from_degrees(a, b)).id()
for a, b in zip(lats, lngs)]
self.order = np.argsort(np.array(ids, dtype=np.uint64))
self.ids = np.array(ids, dtype=np.uint64)[self.order]
def rows_in(self, cell_id):
lo = np.searchsorted(self.ids, np.uint64(cell_id.range_min().id()), side="left")
hi = np.searchsorted(self.ids, np.uint64(cell_id.range_max().id()), side="right")
return self.order[lo:hi]
def count_in(self, cell_id):
return len(self.rows_in(cell_id))
indexed 109,179 GB points in 0.46s
level 6: 11,534 points, 37 us; area 16,017.1 km2
level 8: 3,954 points, 7 us; area 975.4 km2
level 10: 1,030 points, 7 us; area 60.8 km2
level 12: 202 points, 4 us; area 3.8 km2
rows_in returns positional indices into the original arrays, so the same object answers "which rows" as well as "how many". The same idea works in any database: store the leaf id as an unsigned 64-bit column, sort or index it, and filter with between range_min and range_max.
Example 3 — covering a region with a cell budget
import math
import s2sphere as s2
EARTH_RADIUS_KM = 6371.0088
def cover_rectangle(south, west, north, east, max_cells=20, max_level=30):
"""Approximate a lat/lng rectangle with S2 cells and report the excess area."""
rect = s2.LatLngRect.from_point_pair(s2.LatLng.from_degrees(south, west),
s2.LatLng.from_degrees(north, east))
coverer = s2.RegionCoverer()
coverer.min_level, coverer.max_level, coverer.max_cells = 0, max_level, max_cells
cells = coverer.get_covering(rect)
covered = sum(s2.Cell(c).exact_area() for c in cells) * EARTH_RADIUS_KM ** 2
exact = (EARTH_RADIUS_KM ** 2 * math.radians(east - west)
* (math.sin(math.radians(north)) - math.sin(math.radians(south))))
print(f"max_cells={max_cells}: {len(cells)} cells at levels "
f"{min(c.level() for c in cells)}-{max(c.level() for c in cells)}, "
f"covering {covered / exact:.3f}x the rectangle")
return cells
max_cells=8: 8 cells at levels 4-5, covering 2.720x the rectangle
max_cells=20: 14 cells at levels 4-10, covering 1.602x the rectangle
max_cells=100: 43 cells at levels 4-12, covering 1.174x the rectangle
A covering turns "points inside this region" into a handful of integer range filters — the pattern behind S2-based region queries in databases. Filter by the covering first, then test the exact geometry on the survivors if the excess matters.
Explanation
Why squares nest and hexagons cannot
A square divides into four smaller squares with no remainder. A hexagon cannot be divided into smaller hexagons at all: seven hexagons make a rough flower shape whose outline zigzags across the parent's edges. H3 accepts that and rotates each resolution slightly, which is why its children spill 7.14% of their area outside the parent.
S2's children are exactly their parent, measured to 1.29 × 10⁻¹³. So a count at level 12 always sums to the count at level 11, and rolling a result up the hierarchy never moves a point between cells.
Why the Hilbert curve makes ranges
The cell id orders leaves along a Hilbert curve, a path that visits every cell of a quadrant before leaving it, at every level. A cell is therefore one unbroken stretch of the curve, and one unbroken stretch of integers.
A Z-order curve, which geohash follows, has the same range property but makes long jumps between consecutive cells. The Hilbert curve keeps consecutive ids adjacent on the ground, so points close in id order are also close in space, which helps any storage engine that sorts and compresses by key.
H3 has a weaker version of this. At a fixed stored resolution, the resolution-9 ids under one resolution-5 parent formed a contiguous integer range in all 42 cases tested. That works for the hierarchy's parents, which are not the same as the containing cells.
Why the neighbour distances are uneven
Every quadrilateral grid has diagonal neighbours. On a square grid their centres are √2, or 1.41 times, further than the edge neighbours. S2 cells are distorted squares, so the measured ratio around London was 1.752. H3's hexagons have six edge neighbours and no corner neighbours, and their measured ratio was 1.090.
That difference is invisible in a containment query and dominant in anything that treats "neighbour" as "equally near": spatial weights, k-ring smoothing, movement between cells.
Why S2 is slower in Python
s2sphere 0.2.5 is a pure-Python port of the C++ library. Indexing 1,000,000 points took 4.7 s, against 0.80 s for h3-py, which calls compiled C. The gap is the binding, not the grid; the C++ library itself is fast. For bulk indexing in Python, compute leaf ids where a compiled implementation exists — BigQuery, a database extension — or accept the cost once and store the result.
Edge cases or notes
LatLng.from_degreestakes latitude first, like H3 and unlike shapely. Swapping silently indexes a different place.- Level 30 is a leaf of about 0.74 cm². Store leaves and derive coarser cells with
parent(level); a leaf id loses nothing. - Tokens and ids are interchangeable.
CellId.from_token('487604d')returns the level-12 cell; use tokens in logs and URLs, integers in storage. - Signed 64-bit columns overflow. Ids above 2⁶³ appear negative in databases without an unsigned type, so range filters must use the same representation on both sides.
- Cells are geodesic quadrilaterals. Their edges are great-circle arcs, so drawing a cell by joining its four corners in a flat CRS is an approximation that grows with cell size.
- A covering is a superset. Test the exact geometry after the range filter when the excess area matters.
- Area is not monotonic with latitude. Normalise counts by
exact_area()rather than assuming cells near the poles are smaller.
Internal links
- Discrete global grids explained: H3, S2, geohash and why cells beat coordinates — where S2 sits among the grids
- H3 explained: how a hexagonal grid indexes the whole Earth — the hexagonal alternative
- H3 hierarchy explained: parents, children and why they do not nest exactly — the 7.14% overhang in detail
- Geohash explained: prefixes, precision and the edge problem — the latitude–longitude grid S2 improves on
- How to index data by quadkey and web map tile — another quadtree, on Web Mercator
- How to compact an H3 cell set to store a region cheaply — H3's version of a mixed-level covering
- Spatial indexes explained: R-trees and why spatial joins are fast — the tree-based alternative to cell ids
- How to use H3 neighbours for k-ring smoothing and buffers — where equidistant neighbours matter
FAQ
What is an S2 cell?
A region of the Earth's surface defined by projecting the sphere onto a cube, splitting each face into a quadtree and numbering the cells along a Hilbert curve. Each cell has a 64-bit id and a short hexadecimal token.
Is S2 better than H3?
For containment, it is. S2 children tile their parent exactly, and every cell is a range of leaf ids. H3 is better for neighbourhood work, because its six neighbours are almost equidistant: a centre-distance ratio of 1.09 against 1.75 for S2.
How big is an S2 cell at each level?
The average is the Earth's area divided by 6 × 4 to the power of the level: about 5,189 km² at level 7, 5.07 km² at level 12 and 0.74 square centimetres at level 30. Individual cells vary by about a factor of two either way.
Why do S2 cells vary in area if they are not built on latitude?
The cube projection stretches cells differently depending on where they sit on a face. Measured at level 12, a cell on the equator was 4.30 km², one in London 3.80 km² and one at 70° north 5.61 km².
What is S2 used for?
Region coverings and fast range queries: BigQuery exposes S2 cell functions, and S2-style indexes let a database answer "points in this area" with a few integer range filters.
Is there a fast S2 library for Python?
s2sphere is easy to install but pure Python, and it indexed a million points in 4.7 seconds against 0.80 for h3-py. The C++ library is fast; use a compiled implementation for bulk work, or index once and store the ids.