Bucketing and discretization is important step in any algorithmic training or analysis, today we focus on the ways we can map the continuous spherical surface of Earth into manageable buckets. Why exactly we do this? to index different locations on earth to allow for fast geospatial queries, joins and aggregations. There are three popular geo indexing systems – Geohash, Uber’s H3 and Google’s S2, all let you map (lat, lon) -> a cell id, but they differ in what surface they discretize and therefore what kinds of distortion they introduce

Size and Precision

Grid Cells

To see how these indexing systems behave in practice, let’s look at a specific coordinate in Ottawa (Lat: 45.428, Lon: -75.718) and map it to a “City” resolution across all three systems.

We choose the following resolution / level for each grid system:

"city": {"geohash": 5, "h3": 6, "s2": 11}

"city": {"geohash": 5, "h3": 6, "s2": 11}

which ranges from 16 – 40 km2 (City Ward / Large District)

While Geohash, H3, and S2 all have a scale that conceptually maps to a “City District” their actual geometric footprints vary wildly:

Size

• Geohash (Precision 5)
  • Token: f244k (Parent: f244)
  • Centroid: 45.417480, -75.739746
  • Dimensions: 3433.86m x 4891.99m
  • Estimated Area: ~16.8 million m² (16.8 km²)
• H3 (Resolution 6)
  • Token: 862b83b1fffffff (Parent: 852b83b3fffffff)
  • Centroid: 45.419360, -75.762752
  • Edge Length: 3950.16m
  • Estimated Area: ~39.3 million m² (39.3 km²)
• S2 (Level 11)
  • Token: 4cce04c (Parent: 4cce05)
  • Centroid: 45.443463, -75.710880
  • Edge Length: 3828.80m
  • Estimated Area: ~16.1 million m² (16.1 km²)

Quantization

Quantization error – cell centroid distance – If we are putting a cell around a lat lon – how far is the cell center from the lat lon, we compare this on the smallest scale:

• Geohash (Precision 12)
  • Token: f244kvzv68zp (Parent: f244kvzv68z)
  • Centroid: 45.428000, -75.718000
  • Dimensions: 0.03m x 0.02m
  • Estimated Area: ~0.0005 m²
• H3 (Resolution 15)
  • Token: 8f2b83b1d90b96d (Parent: 8e2b83b1d90b96f)
  • Centroid: 45.427996, -75.718002
  • Edge Length: 0.6174m
  • Estimated Area: ~0.9731 m²
• S2 (Level 30)
  • Token: 4cce04f4d08605cd (Parent: 4cce04f4d08605cc)
  • Centroid: 45.428000, -75.718000
  • Edge Length: 0.0073m
  • Estimated Area: ~0.0001 m²

For the given Lat: 45.428, Lon: -75.718, we’ve computed the distance of the point from the centroid (haversine distance – the details of which you can read in the next section)

the smallest possible quantization error:
Geohash(12) Off-Center Distance: 0.01 meters
H3 (15) Off-Center Distance: 0.52 meters
S2(30) Off-Center Distance: 0.00 meters

city level quantization error (for the cell size displayed in the image):
Geohash(5) Off-Center Distance: 2061.21 meters
H3 (6) Off-Center Distance: 3622.33 meters
S2(11) Off-Center Distance: 1806.92 meters

Notice the massive discrepancy in area. At the city tier, S2 and Geohash yield a similarly sized bounding box (~16 sq km), but H3’s hexagonal resolution jumps to nearly 40 sq km. Furthermore, the centroids of these cells do not perfectly align with the target coordinate. This highlights a crucial rule for geospatial engineering: you cannot blindly swap indexing systems and expect a 1:1 mapping. You must tailor the resolution level to the specific precision required by your dataset.

Distance

Planar

Flat map distance, given latitude and longitude of two points, we employ an equirectangular projection and calculate Euclidean distance.

Point 1: $(\phi_1, \lambda_1)$ , Point 2: $(\phi_2, \lambda_2)$
$\phi$ = latitude (radians), $\lambda$ = longitude (radians), $R$ = Earth radius

Reference point $(\phi_0, \lambda_0)$

Projection Center (reduce distortion):
$\phi_0 = \frac{\phi_1 + \phi_2}{2}$
$\lambda_0 = \frac{\lambda_1 + \lambda_2}{2}$

Equirectangular Projection :

east-west displacement in meters
$x = R \cdot (\lambda – \lambda_0) \cdot \cos(\phi_0)$
$\cos(\phi_0)$ term corrects for the fact that lines of longitude
get closer together toward the poles

north-south displacement in meters
$y = R \cdot (\phi – \phi_0)$

Euclidean Distance in Projected Space (m) :

$d = \sqrt{(x_2 – x_1)^2 + (y_2 – y_1)^2}$

Geodesic

Because the earth is curved, to find distance between two points specified with a latitude or longitude, we employ the haversine formula.

Point 1: $(\phi_1, \lambda_1)$ , Point 2: $(\phi_2, \lambda_2)$
$\phi$ = latitude (radians), $\lambda$ = longitude (radians), $R$ = Earth radius
$\Delta \phi = \phi_2 – \phi_1$, $\Delta \lambda = \lambda_2 – \lambda_1$

Haversine:
$a = \sin^2\left(\frac{\Delta \phi}{2}\right)+ \cos(\phi_1)\cos(\phi_2)\sin^2\left(\frac{\Delta \lambda}{2}\right)$

Angular Distance:
$c = 2 \cdot \arctan2\left(\sqrt{a}, \sqrt{1-a}\right)$

Final Distance:
$d = R \cdot c$

For a collection of points as follows:

test_pairs = {
    "local": [(45.428, -75.718), (45.47, -75.70)], #(Ottawa ~5km)
    "regional": [(45.428, -75.718), (45.5017, -73.5673)], #(Ottawa–Montreal)
    "cross_country": [(45.428, -75.718), (49.2827, -123.1207)], #(Ottawa–Vancouver)
    "equator": [(0.0, 0.0), (0.0, 10.0)], #(0°, 0° - 0°, 10°)
    "poles": [(70.0, 0.0), (70.0, 60.0)], #(70°, 0° - 70°, 60°)
}

levels = {
    # "neighborhood":{"geohash": 6,"h3": 7,"s2": 13},
    "smallest":{"geohash": 12, "h3": 15, "s2": 30},
    "largest":{"geohash": 1, "h3": 0, "s2": 0}, # geohash 0 is the whole world 
    "state": {"geohash": 3, "h3": 2, "s2": 6}, # 100,000+ km2 (State / Province Level)
    "metro": {"geohash": 4, "h3": 4, "s2": 8}, # 500+ km2 (Wide Metro Area / County Level)
    # "city" :{"geohash": 5, "h3": 5, "s2": 10},
    "city": {"geohash": 5, "h3": 6, "s2": 11}, #16 - 40 km2 (City Ward / Large District) 
    "neighborhood": {"geohash": 6, "h3": 7, "s2": 12}, #5 km2 (Neighborhood),
    "sentinel_scale": {"geohash": 8, "h3": 12, "s2": 20} #10m approx
}

test_pairs = {
    "local": [(45.428, -75.718), (45.47, -75.70)], #(Ottawa ~5km)
    "regional": [(45.428, -75.718), (45.5017, -73.5673)], #(Ottawa–Montreal)
    "cross_country": [(45.428, -75.718), (49.2827, -123.1207)], #(Ottawa–Vancouver)
    "equator": [(0.0, 0.0), (0.0, 10.0)], #(0°, 0° - 0°, 10°)
    "poles": [(70.0, 0.0), (70.0, 60.0)], #(70°, 0° - 70°, 60°)
}

levels = {
    # "neighborhood":{"geohash": 6,"h3": 7,"s2": 13},
    "smallest":{"geohash": 12, "h3": 15, "s2": 30},
    "largest":{"geohash": 1, "h3": 0, "s2": 0}, # geohash 0 is the whole world 
    "state": {"geohash": 3, "h3": 2, "s2": 6}, # 100,000+ km2 (State / Province Level)
    "metro": {"geohash": 4, "h3": 4, "s2": 8}, # 500+ km2 (Wide Metro Area / County Level)
    # "city" :{"geohash": 5, "h3": 5, "s2": 10},
    "city": {"geohash": 5, "h3": 6, "s2": 11}, #16 - 40 km2 (City Ward / Large District) 
    "neighborhood": {"geohash": 6, "h3": 7, "s2": 12}, #5 km2 (Neighborhood),
    "sentinel_scale": {"geohash": 8, "h3": 12, "s2": 20} #10m approx
}

we see the geodesic and planar distance are noticeably different on large scales because of earth’s curvature

Now lets compare the error between haversine distance for the actual points vs the cell centroids:

The “Local” Distance Anomaly: If you look at the “Local” test scenario (a distance of roughly 5 km) across all three systems, the error hits exactly 100.0000% at the “State” and “Largest” levels. This happens because the physical distance being measured is vastly smaller than the cell itself. If both points fall into the exact same massive cell, the calculated centroid distance is 0, resulting in a 100% error. This is a vital reminder: never use a spatial index resolution that is physically larger than the phenomena you are trying to measure.

Geohash’s Equator Extreme: At the “Largest” resolution level, Geohash struggles massively at the Equator, showing a staggering 314.0962% error. By comparison, H3 at the Equator for the same level shows a 64.7221% error. Geohash’s rectangular planar projection simply cannot handle extreme longitudinal distances reliably at a low resolution.

S2’s Rigorous Bounding: Look at S2’s “Largest” row: it hits a flat 100.0000% error across all scenarios. Because S2 Level 0 cells cover exactly one-sixth of the Earth’s surface, almost all of our test pairs snapped to the same massive cell centroid. However, as soon as you step down to the “Metro” and “City” levels, S2’s errors drop rapidly and consistently.

Convergence at the Sentinel Scale: By the time we reach the “Sentinel_Scale” (roughly 10-meter resolution, perfect for satellite imagery) and the “Smallest” scale, all three systems converge to near 0.0000% error, regardless of whether you are tracking a cross-country flight or a localized grid outage

The following is how the ‘city’ level cells look like over a regional distance

Intersect

By intersect what we mean is geofencing and point in polygon problem, in most of geospatial problems, we will constantly be checking if telemetry coordinates fall within a specific boundary, and for that we can either do the exact mathematical calculation (which is computationally expensive at scale) or check if the point belongs to a pre computed set of grid cells (which is blazing fast but introduces edge-case errors)

When determining if a point falls inside a complex polygon, we use these test cases (firmly inside the polygon, riding the borderline, and just outside) and compare the exact geometric calculation (using the Shapely library) against native grid set lookups across Geohash, H3, and S2 at the “Neighborhood/City” scale.

geofence_poly = [ #downtown ottawa - bounding box
    [45.430, -75.715],
    [45.432, -75.695],
    [45.418, -75.690],
    [45.415, -75.705],
    [45.430, -75.715] 
]

test_points = {
    "Inside": (45.422, -75.700),
    "Borderline": (45.430, -75.710),
    "Outside": (45.410, -75.720)
}

geofence_poly = [ #downtown ottawa - bounding box
    [45.430, -75.715],
    [45.432, -75.695],
    [45.418, -75.690],
    [45.415, -75.705],
    [45.430, -75.715] 
]

test_points = {
    "Inside": (45.422, -75.700),
    "Borderline": (45.430, -75.710),
    "Outside": (45.410, -75.720)
}

Where In Set is a still a mathematical intersects (using shapely) for the cell polygon and intersects is the the grid systems native intersect function.

The Inside Point: All three systems, as well as the exact calculation, correctly identified that the point was inside the polygon. No surprises here.

The Borderline Point: The exact math confirmed the point was inside. Geohash and S2 successfully flagged it via set lookup. However, H3’s native set lookup returned False. Because hexagons don’t perfectly conform to sharp, straight polygon edges, the H3 polyfill algorithm left a slight gap at the border, resulting in a false negative.

The Outside Point: The exact math correctly flagged the point as outside. Geohash and H3 agreed. However, S2 returned True for both the set lookup and dynamic intersection. Because the S2 cell at this specific resolution level was large enough to overlap the polygon boundary while also swallowing the outside point, it yielded a false positive.

Now you ask if grids introduce false positives and negatives, why use them at all? Speed!

• Exact Geometry (Shapely): Ranged from 0.0053 ms to 0.0500 ms. Exact math is highly variable because the computational cost increases depending on how complex the polygon is and where the point lies relative to its edges. This was a simple bounding box, but geo polygons can be and usually are very detailed
• Geohash Set Lookup: Hovered around ~0.0054 to 0.0065 ms.
• S2 Set Lookup: Ranged from ~0.0142 to 0.0160 ms.
• H3 Set Lookup: Blew the competition away at 0.0018 to 0.0031 ms

(and we can always solve for accuracy by adjusting resolution and/or employing hierarchal checks)

Neighborhood

In geospatial analytics, we rarely look at a single cell in isolation.

Not just optimal drive time solutions but even in Earth Observation problems we apply smoothing filters, calculate local density, or model the spread of phenomena (like a flood or a signal footprint) by radiating outward. This requires pulling a cell’s immediate neighbors.

Based on our analysis at the “Neighborhood” scale (roughly 1 to 5 km²), here is how Geohash, H3, and S2 handle spatial adjacency

Hierarchy & Coverage

we’ve already established how these systems group cells (hierarchy) and how they fill a shape (coverage) are fundamentally different, now we tackle on a problem of indexing a geographic boundary – to retrieve data inside, here we see how at different hierarchies and resolutions, different grid cells cover a geographic area

You can access the code here github