S2 Cells Explained: Squares on a Cube and When They Beat Hexagons

On this page

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 8.7 × 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 2.6 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.

Each level quarters every cell, so a cell's four children share its boundary exactly — the thing hexagons cannot do.

A cube, quartered cells and ids as paths

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 on London's meridian measured 4.316 km² on the equator, 3.804 km² in London and 5.168 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 the bounding rectangle of GeoNames' United Kingdom points, 905,396 km² (Example 3):

max_cells   cells used   levels   area covered
    4             4       3–4        4.191 × the rectangle
    8             8       3–5        2.268 ×
   20            17       4–6        1.478 ×
  100            58       4–9        1.150 ×
  500           314       4–13       1.031 ×

A covering always contains the region; the budget controls how much extra it includes. At 500 cells it was within 3.1%.

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 GeoNames points for the United Kingdom, counting the points in London's cells took 34 µs at level 6 and 4–6 µs at levels 8 to 12. Nothing is re-indexed when the question changes resolution.

Each grid wins two rows. Which two matter depends on whether the work is containment or neighbourhoods.

Logging cells and covering regions over Britain

The data used on this page

The points are GeoNames: every entry in the Great Britain file (GB.zip), and for the speed tests a 1,000,000-point random sample of the full allCountries dump. Both come from download.geonames.org and are rebuilt daily, so counts drift a little between runs.

import os
import zipfile

import numpy as np
import pandas as pd
import requests


def fetch(name):
    if not os.path.exists(name):
        with requests.get(f"https://download.geonames.org/export/dump/{name}", stream=True, timeout=600) as r:
            r.raise_for_status()
            with open(name, "wb") as out:
                for chunk in r.iter_content(1 << 20):
                    out.write(chunk)


def read_points(name):
    with zipfile.ZipFile(name) as z, z.open(name.replace(".zip", ".txt")) as f:
        return pd.read_csv(f, sep="\t", header=None, usecols=[4, 5], quoting=3, names=["lat", "lon"])


fetch("GB.zip")
fetch("allCountries.zip")
gb = read_points("GB.zip")
if not os.path.exists("million.parquet"):
    read_points("allCountries.zip").sample(1_000_000, random_state=0).to_parquet("million.parquet")
million = pd.read_parquet("million.parquet")
print(f"Great Britain: {len(gb):,} points; speed sample: {len(million):,}")

Great Britain: 109,179 points; speed sample: 1,000,000

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()),
    }

print(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))


import time

t = time.perf_counter()
index = LeafIndex(gb.lat, gb.lon)
print(f"indexed {len(gb):,} GB points in {time.perf_counter() - t:.2f}s")
london = s2.CellId.from_lat_lng(s2.LatLng.from_degrees(51.5007, -0.1246))
for level in (6, 8, 10, 12):
    cell = london.parent(level)
    t = time.perf_counter()
    n = index.count_in(cell)
    us = (time.perf_counter() - t) * 1e6
    print(f"level {level}: {n:,} points, {us:.0f} us; area {s2.Cell(cell).exact_area() * EARTH_RADIUS_KM ** 2:.1f} km2")

indexed 109,179 GB points in 0.39s
level 6: 11,534 points, 33 us; area 16017.1 km2
level 8: 3,954 points, 6 us; area 975.4 km2
level 10: 1,030 points, 4 us; area 60.8 km2
level 12: 202 points, 6 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


mainland = gb[gb.lon.between(-9, 2.5) & gb.lat.between(49.5, 61)]   # the file also lists overseas territories
south, west, north, east = mainland.lat.min(), mainland.lon.min(), mainland.lat.max(), mainland.lon.max()
area = (EARTH_RADIUS_KM ** 2 * math.radians(east - west)
        * (math.sin(math.radians(north)) - math.sin(math.radians(south))))
print(f"Great Britain's bounding rectangle: {south:.2f} to {north:.2f} N, {west:.2f} to {east:.2f} E, {area:,.0f} km²")
for budget in (4, 8, 20, 100, 500):
    cover_rectangle(south, west, north, east, max_cells=budget)

Great Britain's bounding rectangle: 49.70 to 60.96 N, -8.97 to 2.49 E, 905,396 km²
max_cells=4: 4 cells at levels 3-4, covering 4.191x the rectangle
max_cells=8: 8 cells at levels 3-5, covering 2.268x the rectangle
max_cells=20: 17 cells at levels 4-6, covering 1.478x the rectangle
max_cells=100: 58 cells at levels 4-9, covering 1.150x the rectangle
max_cells=500: 314 cells at levels 4-13, covering 1.031x 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.

Example 4 — the measurements quoted above

import h3

rng = np.random.default_rng(0)
worst = 0.0
for _ in range(500):
    leaf = s2.CellId.from_lat_lng(s2.LatLng.from_degrees(rng.uniform(-80, 80), rng.uniform(-180, 180)))
    parent = leaf.parent(10)
    kids = sum(s2.Cell(c).exact_area() for c in parent.children())
    worst = max(worst, abs(kids - s2.Cell(parent).exact_area()) / s2.Cell(parent).exact_area())
print(f"500 level-10 cells: children's area differs from the parent's by at most {worst:.3g}")

earth = 4 * math.pi * EARTH_RADIUS_KM ** 2
for level in (0, 7, 12, 30):
    print(f"level {level:2}: {6 * 4 ** level:,} cells, average {earth / (6 * 4 ** level):,.6g} km²")
level7 = []
for face in range(6):
    start = s2.CellId.from_face_pos_level(face, 0, 0)
    c, end = start.child_begin(7), start.child_end(7)
    while c != end:
        level7.append(s2.Cell(c).exact_area() * EARTH_RADIUS_KM ** 2)
        c = c.next()
level7 = np.array(level7)
print(f"level 7 enumerated: {len(level7):,} cells, mean {level7.mean():,.1f} km², "
      f"min {level7.min():,.1f}, max {level7.max():,.1f}, max/min {level7.max() / level7.min():.3f}")

cell12 = london.parent(12)
print("London level 12:", cell12.to_token(), "children", [c.to_token() for c in cell12.children()],
      "face", cell12.face())
for label, lat in (("equator", 0.0), ("London", 51.5007), ("70° N", 70.0)):
    c = s2.CellId.from_lat_lng(s2.LatLng.from_degrees(lat, -0.1246)).parent(12)
    print(f"level-12 cell at {label:8}: {s2.Cell(c).exact_area() * EARTH_RADIUS_KM ** 2:.3f} km²")


def centre_km(a, b):
    return a.get_center().angle(b.get_center()) * EARTH_RADIUS_KM


c13 = s2.Cell(london.parent(13))
around = [s2.Cell(n) for n in london.parent(13).get_all_neighbors(13)]
d = [centre_km(c13, n) * 1000 for n in around]
print(f"S2 level 13: {len(around)} neighbours, centre distances {min(d):,.0f}–{max(d):,.0f} m "
      f"(ratio {max(d) / min(d):.3f})")
h8 = h3.latlng_to_cell(51.5007, -0.1246, 8)
hd = [h3.great_circle_distance(h3.cell_to_latlng(h8), h3.cell_to_latlng(n), "m")
      for n in h3.grid_ring(h8, 1)]
print(f"H3 res 8:    {len(hd)} neighbours, centre distances {min(hd):,.0f}–{max(hd):,.0f} m "
      f"(ratio {max(hd) / min(hd):.3f})")

t = time.perf_counter()
big = LeafIndex(million.lat, million.lon)
t_s2 = time.perf_counter() - t
t = time.perf_counter()
[h3.latlng_to_cell(a, b, 9) for a, b in zip(million.lat, million.lon)]
t_h3 = time.perf_counter() - t
c8 = s2.CellId.from_lat_lng(s2.LatLng.from_degrees(51.5007, -0.1246)).parent(8)
t = time.perf_counter()
for _ in range(1000):
    big.count_in(c8)
print(f"1,000,000 points: s2sphere {t_s2:.1f} s, h3 {t_h3:.2f} s; a level-8 count "
      f"{(time.perf_counter() - t) / 1000 * 1e6:.1f} µs; load average {os.getloadavg()[0]:.1f}")

contiguous = 0
parents = sorted({h3.latlng_to_cell(a, b, 5) for a, b in zip(gb.lat, gb.lon)})[:42]
for p5 in parents:
    kids = set(h3.cell_to_children(p5, 9))
    low, high = min(map(h3.str_to_int, kids)), max(map(h3.str_to_int, kids))
    outside = {n for k in kids for n in h3.grid_ring(k, 1)} - kids
    contiguous += not any(low <= h3.str_to_int(n) <= high for n in outside)
print(f"H3: {contiguous} of {len(parents)} resolution-5 parents have no outside neighbour "
      f"inside their children's id range")

500 level-10 cells: children's area differs from the parent's by at most 8.7e-14
level  0: 6 cells, average 8.5011e+07 km²
level  7: 98,304 cells, average 5,188.66 km²
level 12: 100,663,296 cells, average 5.06705 km²
level 30: 6,917,529,027,641,081,856 cells, average 7.37353e-11 km²
level 7 enumerated: 98,304 cells, mean 5,188.7 km², min 3,175.4, max 6,529.1, max/min 2.056
London level 12: 487604d children ['487604c4', '487604cc', '487604d4', '487604dc'] face 2
level-12 cell at equator : 4.316 km²
level-12 cell at London  : 3.804 km²
level-12 cell at 70° N   : 5.168 km²
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)
1,000,000 points: s2sphere 3.6 s, h3 0.61 s; a level-8 count 2.6 µs; load average 1.8
H3: 42 of 42 resolution-5 parents have no outside neighbour inside their children's id range

Squares nest, Hilbert curves make ranges

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 8.7 × 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 3.6 s, against 0.61 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.

A covering never misses part of the region; the budget only decides how much extra it takes in.

Latitude first, leaf cells and signed overflow

LatLng.from_degrees takes 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 level-7 cells range from 3,175 to 6,529 km², a factor of 2.06.

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 on London's meridian, a cell on the equator was 4.32 km², one in London 3.80 km² and one at 70° north 5.17 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 3.6 seconds against 0.61 for h3-py. The C++ library is fast; use a compiled implementation for bulk work, or index once and store the ids.