How to Fill a Polygon with H3 Cells
Problem statement
You have a polygon — a country, a delivery zone, a catchment — and you need the set of H3 cells that represent it. The cell set becomes a join key, a mask for aggregated points, or a region you can store as a list of integers instead of a geometry.
The call itself is one line. The decisions around it are not, and three of the failure modes are silent:
- A polygon in a projected CRS returns zero cells with no error. Measured on Switzerland in EPSG:2056,
geo_to_cellsreturned an empty list. - A hand-built
LatLngPolywith the coordinates in shapefile order returns 1,979 cells instead of 1,189, centred on the Horn of Africa. - The default rule — a cell belongs if its centre is inside — missed 7.1% of Switzerland's area at resolution 5 and added 5.5% outside it. The count looked right because the two errors nearly cancel.
The fourth failure is loud but surprising: the Natural Earth polygon for mainland Antarctica raises H3FailedError at resolution 5, while the experimental fill function in h3 4.5 returns 45,569 cells for the same shape.
Quick answer
Pass the shapely geometry straight to h3.geo_to_cells, in longitude/latitude degrees:
import geopandas as gpd
import h3
countries = gpd.read_file("ne_10m_admin_0_countries.zip")
switzerland = countries.loc[countries.ADMIN == "Switzerland", "geometry"].iloc[0]
cells = h3.geo_to_cells(switzerland, 6)
print(len(cells), cells[:2])1189 ['861f9dd6fffffff', '861f9c2efffffff']
geo_to_cells reads anything with a __geo_interface__, so a shapely Polygon or MultiPolygon works without conversion and the GeoJSON axis order is handled for you. Switzerland's geodesic area divided by the average resolution-6 cell area predicts 1,147 cells; 1,189 is 3.7% more, which is normal (see Explanation).
Step-by-step solution
1. Put the polygon in longitude/latitude degrees
H3 works on the sphere, so the input must be EPSG:4326. A projected polygon is not rejected — its coordinates are simply far outside any valid latitude, and nothing is inside them:
ch = countries.loc[countries.ADMIN == "Switzerland"]
lv95 = ch.to_crs(2056).geometry.iloc[0]
print(lv95.bounds) # metres, not degrees
print(len(h3.geo_to_cells(lv95, 6))) # 0(2485468.010442513, 1075551.7941011817, 2831994.8583251485, 1295112.798875499)
0
Check gdf.crs and call to_crs(4326) before filling. An empty result from a polygon you know is large is almost always this.
2. Estimate the cell count before filling
Geodesic area divided by the average cell area gives the order of magnitude, which tells you whether the resolution is sensible and gives you a number to check the result against:
from pyproj import Geod
geod = Geod(ellps="WGS84")
area_km2 = abs(geod.geometry_area_perimeter(switzerland)[0]) / 1e6
print({res: round(area_km2 / h3.average_hexagon_area(res, "km^2")) for res in range(4, 9)}){4: 23, 5: 164, 6: 1147, 7: 8028, 8: 56197}
Each resolution step multiplies the count by about seven. If the figure is in the millions, decide whether you really need that resolution before generating the list — Switzerland at resolution 9 is 407,632 cells.
3. Choose a containment mode deliberately
geo_to_cells applies the centre rule; h3.geo_to_h3shape followed by h3.h3shape_to_cells is the same operation in two steps, and both gave the identical 1,189-cell set at resolution 6. h3 4.5 exposes four containment rules through h3shape_to_cells_experimental:
shape = h3.geo_to_h3shape(switzerland)
for mode in ("center", "full", "overlap", "bbox_overlap"):
print(mode, len(h3.h3shape_to_cells_experimental(shape, 5, mode)))center 167
full 126
overlap 218
bbox_overlap 268
center— the cell's centre is inside. Same result asgeo_to_cells.full— the whole cell is inside. Nothing outside; the edges are missing.overlap— any part of the cell touches the polygon. Nothing missing; a fringe outside.bbox_overlap— the cell's bounding box touches the polygon. A cheap, loose superset.
The function name says experimental for a reason: the docstring states it has no API stability guarantee across versions. Pin h3 if you depend on it, and keep the mode string in one place. A typo is not a clear error — "intersects" raised KeyError: 'intersects'.
4. Let holes and multipolygons through untouched
Natural Earth's South Africa is a MultiPolygon with a hole where Lesotho sits. H3 honours both:
za = countries.loc[countries.ADMIN == "South Africa", "geometry"].iloc[0]
ls = countries.loc[countries.ADMIN == "Lesotho", "geometry"].iloc[0]
za_cells, ls_cells = set(h3.geo_to_cells(za, 6)), set(h3.geo_to_cells(ls, 6))
print(len(za_cells), len(ls_cells), len(za_cells & ls_cells))31970 779 0
No overlap. Rebuilding the parts from their exterior rings alone — a common side effect of hand-rolled geometry cleaning — gave 32,749 cells, including all 779 of Lesotho's.
5. Build a LatLngPoly yourself only with (lat, lng) order
When coordinates do not come from shapely, h3.LatLngPoly takes (latitude, longitude) pairs — the opposite of shapely, GeoJSON and every shapefile:
ring = list(switzerland.exterior.coords) # (lng, lat) pairs
right = h3.LatLngPoly([(lat, lng) for lng, lat in ring])
wrong = h3.LatLngPoly(ring)
print(len(h3.h3shape_to_cells(right, 6)), len(h3.h3shape_to_cells(wrong, 6)))1189 1979
Both calls succeed. The swapped polygon is Switzerland mirrored across the diagonal, near 9° N 47° E, and it does not even return the same count: a degree of longitude covers more ground near the equator, so the mirrored shape is larger. Nothing short of plotting or checking a centre reveals it.
6. Check the result, then store it compactly
Compare the count with step 2, and look at one centre with h3.cell_to_latlng. Then, if the set is going to be stored or shipped, compact it:
cells9 = h3.geo_to_cells(switzerland, 9)
compact = h3.compact_cells(cells9)
print(len(cells9), len(compact), set(h3.uncompact_cells(compact, 9)) == set(cells9))407632 7252 True
A 98.2% reduction, losslessly: interior blocks of seven children are replaced by their parent, recursively, and only the edge stays fine.
Code examples
Example 1 — a fill with the checks that catch silent failures
import warnings
import h3
from pyproj import Geod
_GEOD = Geod(ellps="WGS84")
def fill_polygon(geom, res, contain="center"):
"""Cells for one lon/lat polygon, with the checks that catch silent failures."""
minx, miny, maxx, maxy = geom.bounds
if not (-180 <= minx <= 180 and -180 <= maxx <= 180
and -90 <= miny <= 90 and -90 <= maxy <= 90):
raise ValueError(f"bounds {geom.bounds} are not degrees — reproject to EPSG:4326 first")
shape = h3.geo_to_h3shape(geom)
if contain == "center":
cells = h3.h3shape_to_cells(shape, res)
else:
cells = h3.h3shape_to_cells_experimental(shape, res, contain)
area_km2 = abs(_GEOD.geometry_area_perimeter(geom)[0]) / 1e6
expected = area_km2 / h3.average_hexagon_area(res, "km^2")
if not cells:
warnings.warn(f"no cells at res {res}: polygon is {expected:.2f} average cells in area")
elif not 0.5 < len(cells) / expected < 2:
warnings.warn(f"{len(cells):,} cells but area suggests ~{expected:,.0f}")
return cellsOn Singapore at resolution 4 it returns an empty list with the warning no cells at res 4: polygon is 0.29 average cells in area; on the EPSG:2056 Switzerland it raises instead of returning nothing.
Example 2 — every feature in a layer, with empties reported
import pandas as pd
def polygons_to_cells(gdf, res, id_col, contain="center"):
"""One row per (feature, cell); reports empty and retried features instead of hiding them."""
gdf = gdf.to_crs(4326)
rows, empty, retried = [], [], []
for fid, geom in zip(gdf[id_col], gdf.geometry):
try:
cells = fill_polygon(geom, res, contain)
except h3.H3FailedError:
# the default fill gives up on some polar rings; the experimental one copes
cells = h3.h3shape_to_cells_experimental(h3.geo_to_h3shape(geom), res, contain)
retried.append(fid)
if not cells:
empty.append(fid)
rows.extend((fid, c) for c in cells)
out = pd.DataFrame(rows, columns=[id_col, "cell"])
print(f"{len(gdf):,} features -> {len(out):,} rows; {len(empty):,} empty; retried {retried}")
return out, emptyOn all 4,596 Natural Earth admin-1 polygons at resolution 5, in 9.8 s:
4,596 features -> 570,084 rows; 704 empty; retried ['ATA+00?']
The retry recovered 45,569 cells for mainland Antarctica. A first version that retried part by part with geo_to_cells recovered only 629, because the failing part was the mainland itself. Across the other 4,595 polygons the experimental center result and the default were identical.
Example 3 — measuring what each mode misses and adds
from shapely.geometry import shape
def coverage_error(geom, res, modes=("center", "full", "overlap")):
"""How much of the polygon each containment mode misses, and how much it adds."""
area = abs(_GEOD.geometry_area_perimeter(geom)[0])
h3shape = h3.geo_to_h3shape(geom)
for mode in modes:
cells = h3.h3shape_to_cells_experimental(h3shape, res, mode)
if not cells:
print(f"{mode:8} {0:>7,} cells missed 100.0% extra 0.0%")
continue
covered = shape(h3.cells_to_geo(cells))
missed = abs(_GEOD.geometry_area_perimeter(geom.difference(covered))[0]) / area
extra = abs(_GEOD.geometry_area_perimeter(covered.difference(geom))[0]) / area
print(f"{mode:8} {len(cells):>7,} cells missed {missed:6.1%} extra {extra:6.1%}")>>> coverage_error(switzerland, 5)
center 167 cells missed 7.1% extra 5.5%
full 126 cells missed 25.8% extra 0.0%
overlap 218 cells missed 0.0% extra 28.5%
>>> coverage_error(switzerland, 7)
center 8,317 cells missed 0.9% extra 0.9%
full 7,927 cells missed 4.7% extra 0.0%
overlap 8,709 cells missed 0.0% extra 4.7%
h3.cells_to_geo builds the outline from cell topology, so the union is instant. Ireland at resolution 5 behaved the same way: full missed 30.5%, overlap added 30.9%.
Explanation
Why the centre rule gets the total right and the edges wrong
A hexagon straddling the boundary is kept when its centre is inside. On a boundary that wiggles at the scale of a cell, roughly as many straddling cells have their centre in as out, so the area kept outside roughly balances the area dropped inside.
That is why Switzerland at resolution 5 came to 98.4% of the true area while being wrong along the whole border — 7.1% missed and 5.5% extra. At resolution 7 both errors fell to 0.9%, because the band of straddling cells shrinks with the cell size.
Use full when a cell must certainly be inside (a guaranteed-interior mask), overlap when nothing inside may be lost (a candidate filter before an exact test), and center for everything statistical.
Why the count is not area divided by average cell area
average_hexagon_area is an average over the whole Earth, and H3 cells are not equal-area. Within one resolution the largest hexagon is about twice the smallest, and the size depends on where the cell sits on the underlying icosahedron, not on latitude.
Measured with the centre rule: Switzerland got 3.7% more cells than the area predicts, Chile 8% more at resolution 7, and Iceland 14% more at resolutions 5 and 7 — Iceland simply sits where resolution-5 and resolution-7 cells are smaller than average. Expect the ratio to vary by ±15%, and treat anything outside a factor of two as a bug.
Why a small polygon gets no cells
With the centre rule, a polygon that contains no cell centre gets nothing. That is not limited to polygons smaller than a cell: an elongated or badly placed polygon can miss every centre.
On the 4,596 admin-1 polygons, 1,440 (31.3%) got no cells at resolution 4 and 704 (15.3%) at resolution 5. The largest empty one at resolution 4 covered 3,798 km², more than twice the 1,770 km² average cell. Singapore at resolution 4 is 0.29 average cells in area and came back empty; overlap gave it 3.
If every feature must be represented, either raise the resolution, use overlap, or fall back to the cell containing a representative point.
Why overlap costs more than centre
The centre rule tests one point per candidate cell. The overlap rule has to intersect the cell's boundary with the polygon's edges, which for a coastline of thousands of vertices is far more work.
Measured on Chile at resolution 7: 7.5 s for center against 19.4 s for overlap, for 154,226 and 163,980 cells. On a compact polygon like Switzerland the difference was negligible (0.017 s against 0.043 s). The cost scales with boundary complexity, not with area.
Why the default fill can fail where the experimental one does not
Natural Earth's mainland Antarctica is one ring of 15,953 vertices that runs from −180° to 180° and encloses the South Pole. geo_to_cells filled it at resolution 4 and raised H3FailedError at resolutions 5 and 6 — the exception carries no message.
h3shape_to_cells_experimental is a separate implementation, and it filled the same ring with 45,569 centre-rule cells at resolution 5. That, and the zero differences on the other 4,595 polygons, is why Example 2 retries with it rather than with a different resolution. The other polygons that touch the antimeridian — Russia and Fiji — are already split at 180° in Natural Earth and filled without complaint: Russia gave 9,245 cells at resolution 4 against 9,591 expected. Across all 258 countries at resolution 5 the default fill took 13.9 s for 524,430 cells, and Antarctica was the only failure.
Edge cases or notes
- A projected polygon returns an empty list. No error, no warning — check the bounds are degrees.
LatLngPolyis (lat, lng).geo_to_cellsandgeo_to_h3shapetake shapely or GeoJSON order and swap it for you.- The experimental function may change. Its docstring says so; pin h3 and wrap the call in one function.
- Shared borders double up under
overlap. The 26 Swiss cantons at resolution 8 gave 3,178 duplicate cells underoverlapand none undercenter, whose 58,233 cells matched the national fill exactly. - Holes are honoured. Stripping interiors during cleaning quietly adds the hole's cells back.
- H3 reads a long edge as crossing 180°. A box ring written from 179° to −179° filled with 712 cells, the same as the split
MultiPolygon— but shapely sees that ring as 358° wide, so split it anyway before any other operation. - Order is not guaranteed. Sort the list if you need a stable file or a reproducible test.
- Integer cells come from
h3.api.basic_int, with the same function names; in pandas they take about a third of the memory of the string form.
Internal links
- H3 explained: how a hexagonal grid indexes the whole Earth — the grid behind every cell in the fill
- Choosing an H3 resolution — cell sizes and counts per level before you pick one
- How to turn H3 cells into a GeoDataFrame of polygons — drawing and checking a fill
- Fixing H3 polygon fill that misses cells or returns nothing — the failure modes in depth
- How to compact an H3 cell set to store a region cheaply — the 98% storage reduction
- H3 hierarchy explained — why compaction works despite imperfect nesting
- How to assign points to H3 cells in Python — the other half of a cell-based join
- How to calculate area and distance in GeoPandas correctly — geodesic area for the count check
FAQ
Which function fills a polygon with H3 cells in h3 v4?
h3.geo_to_cells(geometry, res) for anything with a GeoJSON interface, or h3.polygon_to_cells / h3.h3shape_to_cells for an H3Shape. The v3 name polyfill no longer exists.
Why does filling my polygon return an empty list?
Usually the polygon is in a projected CRS, or it contains no cell centre at that resolution. A Switzerland polygon in EPSG:2056 returned zero cells without an error, and 31.3% of admin-1 polygons were empty at resolution 4.
What is the difference between center, full and overlap?
Whether a cell's centre, its whole area, or any part of it must be inside the polygon. On Switzerland at resolution 5 they gave 167, 126 and 218 cells.
Why do I get more cells than area divided by cell area?
Because cells near your polygon may be smaller than the global average. Iceland got 14% more cells than its area predicts; anything within about ±15% is normal.
Can I pass a shapely polygon directly?
Yes. geo_to_cells and geo_to_h3shape accept any object with __geo_interface__, including shapely polygons and multipolygons, and they handle the axis order.
How do I store a filled polygon compactly?
Call h3.compact_cells. Switzerland at resolution 9 went from 407,632 cells to 7,252, and uncompact_cells restored the identical set.