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_cells returned an empty list.
• A hand-built LatLngPoly with 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 44,940 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).

Degrees, an estimated count, a containment mode

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 as geo_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.

Filling country outlines with checks

The data used on this page

Country outlines from Natural Earth's 1:10m admin-0 layer, and its admin-1 layer of states and provinces for the layer-wide fills.

import os
import zipfile

import pandas as pd
import requests


def fetch(url, path):
    if not os.path.exists(path):
        with requests.get(url, stream=True, timeout=600) as reply:
            reply.raise_for_status()
            with open(path, "wb") as out:
                for chunk in reply.iter_content(1 << 20):
                    out.write(chunk)


fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_0_countries.zip",
      "ne_10m_admin_0_countries.zip")
fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_1_states_provinces.zip",
      "ne_10m_admin_1_states_provinces.zip")
print("ready")

ready

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 cells

On 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

The function, then a run on all the Natural Earth admin-1 polygons at resolution 5:

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, empty


admin1 = gpd.read_file("ne_10m_admin_1_states_provinces.zip")
rows5, empty5 = polygons_to_cells(admin1, 5, "adm1_code")
print(f"{(rows5.adm1_code == 'ATA+00?').sum():,} cells for the retried feature")

4,596 features -> 570,084 rows; 704 empty; retried ['ATA+00?']
45,569 cells for the retried feature

The retry recovered 45,569 cells for Antarctica, 44,940 of them on the mainland (Example 4). 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 = list(set(h3.h3shape_to_cells_experimental(h3shape, res, mode)))   # parts can share a cell
        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%}")


for res in (5, 7):
    print(f">>> coverage_error(switzerland, {res})")
    coverage_error(switzerland, res)
print(">>> coverage_error(ireland, 5)")
coverage_error(countries.loc[countries.ADMIN == "Ireland", "geometry"].iloc[0], 5)

>>> 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%
>>> coverage_error(ireland, 5)
center       289 cells  missed   5.9%  extra   6.2%
full         200 cells  missed  30.5%  extra   0.0%
overlap      377 cells  missed   0.0%  extra  30.9%

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

Example 4 — the other measurements on this page

import time

from shapely.geometry import MultiPolygon, Polygon, box


def km2(geom):
    return abs(_GEOD.geometry_area_perimeter(geom)[0]) / 1e6


def ratio(name, res):
    geom = countries.loc[countries.ADMIN == name, "geometry"].iloc[0]
    expected = km2(geom) / h3.average_hexagon_area(res, "km^2")
    return len(h3.geo_to_cells(geom, res)) / expected - 1


print("cells against area estimate:", {f"{n} r{r}": f"{ratio(n, r):+.1%}" for n, r in
      [("Switzerland", 6), ("Switzerland", 5), ("Chile", 7), ("Iceland", 5), ("Iceland", 7)]})
centre5 = shape(h3.cells_to_geo(h3.geo_to_cells(switzerland, 5)))
print(f"Switzerland r5 centre cells cover {km2(centre5) / km2(switzerland):.1%} of its area")
print(f"compaction at r9 removed {1 - len(compact) / len(cells9):.1%} of the cells")

outer = MultiPolygon([Polygon(p.exterior) for p in za.geoms])
outer_cells = set(h3.geo_to_cells(outer, 6))
print(f"South Africa from exterior rings: {len(outer_cells):,} cells, {len(outer_cells & ls_cells)} of Lesotho's")

# empty fills across the admin-1 layer
a1 = admin1.to_crs(4326)
for res in (4, 5):
    empty, largest = 0, 0
    for fid, geom in zip(a1.adm1_code, a1.geometry):
        try:
            cells = h3.geo_to_cells(geom, res)
        except h3.H3FailedError:
            continue
        if not cells:
            empty += 1
            largest = max(largest, km2(geom))
    print(f"r{res}: {empty:,} of {len(a1):,} empty ({empty / len(a1):.1%}); largest empty {largest:,.0f} km²; "
          f"average cell {h3.average_hexagon_area(res, 'km^2'):,.0f} km²")
sg = countries.loc[countries.ADMIN == "Singapore", "geometry"].iloc[0]
print("Singapore r4 overlap:", len(h3.h3shape_to_cells_experimental(h3.geo_to_h3shape(sg), 4, "overlap")))

# Antarctica: the default fill against the experimental one
ata = countries.loc[countries.ADMIN == "Antarctica", "geometry"].iloc[0]
mainland = max(ata.geoms, key=lambda g: g.area)
print(f"mainland Antarctica: {len(mainland.exterior.coords):,} vertices, x {mainland.bounds[0]:.0f} to {mainland.bounds[2]:.0f}")
for res in (4, 5, 6):
    try:
        print(f"  r{res} geo_to_cells: {len(h3.geo_to_cells(mainland, res)):,} cells")
    except h3.H3FailedError as err:
        print(f"  r{res} geo_to_cells: {type(err).__name__} {str(err)!r}")
print(f"  r5 experimental centre: {len(h3.h3shape_to_cells_experimental(h3.geo_to_h3shape(mainland), 5, 'center')):,} cells")
ata1 = a1.loc[a1.adm1_code == "ATA+00?", "geometry"].iloc[0]
part_by_part = 0
for part in getattr(ata1, "geoms", [ata1]):
    try:
        part_by_part += len(h3.geo_to_cells(part, 5))
    except h3.H3FailedError:
        pass
print(f"part-by-part geo_to_cells retry of ATA+00?: {part_by_part:,} cells")
differ = 0
for fid, geom in zip(a1.adm1_code, a1.geometry):
    if fid == "ATA+00?":
        continue
    shp = h3.geo_to_h3shape(geom)
    differ += set(h3.h3shape_to_cells(shp, 5)) != set(h3.h3shape_to_cells_experimental(shp, 5, "center"))
print(f"experimental centre differs from the default on {differ} of {len(a1) - 1:,} other polygons")

russia = countries.loc[countries.ADMIN == "Russia", "geometry"].iloc[0]
print(f"Russia r4: {len(h3.geo_to_cells(russia, 4)):,} cells, {km2(russia) / h3.average_hexagon_area(4, 'km^2'):,.0f} expected")
total, failed = 0, []
for name, geom in zip(countries.ADMIN, countries.geometry):
    try:
        total += len(h3.geo_to_cells(geom, 5))
    except h3.H3FailedError:
        failed.append(name)
print(f"all {len(countries)} countries at r5: {total:,} cells; failed {failed}")

# shared borders: the Swiss cantons at resolution 8
cantons = a1.loc[a1.admin == "Switzerland"]
for mode in ("center", "overlap"):
    lists = [h3.h3shape_to_cells_experimental(h3.geo_to_h3shape(g), 8, mode) for g in cantons.geometry]
    flat = [c for cells in lists for c in cells]
    print(f"{len(cantons)} cantons r8 {mode}: {len(flat):,} cells, {len(flat) - len(set(flat)):,} duplicates")
print("national r8 fill:", len(h3.geo_to_cells(switzerland, 8)),
      "same set as centre cantons:", set(h3.geo_to_cells(switzerland, 8)) == {c for g in cantons.geometry for c in h3.geo_to_cells(g, 8)})

# a ring that crosses the antimeridian
ring = Polygon([(179, -1), (-179, -1), (-179, 1), (179, 1)])
split = MultiPolygon([box(179, -1, 180, 1), box(-180, -1, -179, 1)])
print(f"179° to -179° ring: {len(h3.geo_to_cells(ring, 5))} cells; split: {len(h3.geo_to_cells(split, 5))}; "
      f"shapely width {ring.bounds[2] - ring.bounds[0]:.0f}°")

cells against area estimate: {'Switzerland r6': '+3.7%', 'Switzerland r5': '+1.9%', 'Chile r7': '+8.1%', 'Iceland r5': '+14.4%', 'Iceland r7': '+14.3%'}
Switzerland r5 centre cells cover 98.4% of its area
compaction at r9 removed 98.2% of the cells
South Africa from exterior rings: 32,749 cells, 779 of Lesotho's
r4: 1,440 of 4,596 empty (31.3%); largest empty 3,798 km²; average cell 1,770 km²
r5: 704 of 4,596 empty (15.3%); largest empty 397 km²; average cell 253 km²
Singapore r4 overlap: 3
mainland Antarctica: 15,955 vertices, x -180 to 180
  r4 geo_to_cells: 6,431 cells
  r5 geo_to_cells: H3FailedError ''
  r6 geo_to_cells: H3FailedError ''
  r5 experimental centre: 44,940 cells
part-by-part geo_to_cells retry of ATA+00?: 629 cells
experimental centre differs from the default on 0 of 4,595 other polygons
Russia r4: 9,245 cells, 9,591 expected
all 258 countries at r5: 524,430 cells; failed ['Antarctica']
26 cantons r8 center: 58,233 cells, 0 duplicates
26 cantons r8 overlap: 62,434 cells, 3,178 duplicates
national r8 fill: 58233 same set as centre cantons: True
179° to -179° ring: 261 cells; split: 261; shapely width 358°

chile = countries.loc[countries.ADMIN == "Chile", "geometry"].iloc[0]
for name, geom in (("Chile", chile), ("Switzerland", switzerland)):
    shp = h3.geo_to_h3shape(geom)
    for mode in ("center", "overlap"):
        started = time.perf_counter()
        n = len(h3.h3shape_to_cells_experimental(shp, 7, mode))
        print(f"{name} r7 {mode}: {n:,} cells in {time.perf_counter() - started:.3f} s")
started = time.perf_counter()
for geom in countries.geometry:
    try:
        h3.geo_to_cells(geom, 5)
    except h3.H3FailedError:
        pass
print(f"all countries at r5: {time.perf_counter() - started:.1f} s")

Chile r7 center: 154,226 cells in 6.757 s
Chile r7 overlap: 163,980 cells in 18.208 s
Switzerland r7 center: 8,317 cells in 0.016 s
Switzerland r7 overlap: 8,709 cells in 0.040 s
all countries at r5: 10.6 s

Centre containment gets totals right and edges wrong

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: 6.8 s for center against 18.2 s for overlap, for 154,226 and 163,980 cells. On a compact polygon like Switzerland the difference was negligible (0.016 s against 0.040 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,955 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 44,940 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 10.6 s for 524,430 cells, and Antarctica was the only failure.

Projected input, coordinate order and shared borders

• A projected polygon returns an empty list. No error, no warning — check the bounds are degrees.
• LatLngPoly is (lat, lng). geo_to_cells and geo_to_h3shape take 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 under overlap and none under center, 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 2°-tall box ring written from 179° to −179° filled with 261 cells at resolution 5, 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.