How to Turn H3 Cells into a GeoDataFrame of Polygons

On this page

Problem statement

H3 work happens on cell IDs: latlng_to_cell, a groupby, a merge. At some point the cells have to become geometry — to draw a map, to intersect with another layer, to export a GeoPackage, or simply to look at what the numbers mean.

The conversion looks trivial and has four traps, each measured here:

Axis order. h3.cell_to_boundary returns (lat, lng) pairs. Built naively into a shapely polygon, a cell over Bern has its centroid in Ethiopia.
Speed. A Python loop over 500,000 cells took 4.34 s; a vectorised build took 1.72 s.
The antimeridian. 1,547 resolution-5 cells span more than 180° of longitude when drawn naively. One cell near Fiji becomes a sliver 359.9° wide that crosses the entire map.
Area. Measuring the polygons in Web Mercator overstated the area of occupied GeoNames cells by a median 1.35×, and by up to 134× near the poles.

Quick answer

Swap the pairs, and let shapely build all the rings at once:

import geopandas as gpd
import h3
import numpy as np
import shapely


def cells_to_gdf(cells):
    bounds = [h3.cell_to_boundary(c) for c in cells]
    lengths = [len(b) for b in bounds]
    latlng = np.array([pt for b in bounds for pt in b])
    rings = shapely.linearrings(latlng[:, ::-1],
                                indices=np.repeat(np.arange(len(cells)), lengths))
    return gpd.GeoDataFrame({"cell": cells}, geometry=shapely.polygons(rings), crs="EPSG:4326")


cells = h3.grid_disk(h3.latlng_to_cell(46.948, 7.447, 7), 2)
print(cells_to_gdf(cells).head(3))

              cell                                           geometry
0  871f8342cffffff  POLYGON ((7.44113 46.95115, 7.44285 46.93903, ...
1  871f83421ffffff  POLYGON ((7.40954 46.94754, 7.41126 46.93543, ...
2  871f8342effffff  POLYGON ((7.42791 46.93118, 7.42963 46.91906, ...

latlng[:, ::-1] is the axis swap. The indices argument is what lets cells with different vertex counts share one call. For a map that reaches the antimeridian, use the version in Example 1.

The reversal is the whole difference between Bern and Ethiopia.

Reverse the order, expect odd vertex counts, fix the seam

1. Reverse the coordinate order

c = h3.latlng_to_cell(46.948, 7.447, 7)
print(h3.cell_to_boundary(c)[:2])

wrong = shapely.Polygon(h3.cell_to_boundary(c))
right = shapely.Polygon([(lng, lat) for lat, lng in h3.cell_to_boundary(c)])
print(wrong.centroid, right.centroid)

((46.95114682882678, 7.441133457678303), (46.939034899474194, 7.442848983302808))
POINT (46.946889981272605 7.4577903235804675) POINT (7.4577903235804675 46.946889981272605)

H3 follows the (latitude, longitude) convention of its C library; shapely, GeoJSON and GeoPandas use (x, y). The wrong polygon is perfectly valid — a spatial join put its centroid inside Ethiopia — so the mistake only shows on a map or in a join that finds nothing.

2. Do not assume six vertices

A (n, 6, 2) array looks like the obvious vectorisation, and it breaks:

try:
    np.array([h3.cell_to_boundary(x) for x in ["871f8342cffffff", "850069b7fffffff"]])
except ValueError as exc:
    print("ValueError:", str(exc)[:110], "...")

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. T ...

Pentagons have five vertices, but the common case is different: at the odd, "class III" resolutions, cells that straddle an edge of the underlying icosahedron gain extra vertices where the edge bends them. Every one of the 2,016,842 resolution-5 cells, counted:

6 vertices   2,002,160
7 vertices      11,760
8 vertices       2,910
10 vertices         12     (the pentagons)

At resolution 1 a third of all cells have more than six. Flat coordinates plus an indices array, as in the quick answer, handle every case.

3. Attach the values and the CRS

The cell ID is the key, so a merge brings the data across:

import pandas as pd

counts = pd.DataFrame({"cell": cells, "n": range(len(cells))})
gdf = cells_to_gdf(counts["cell"].tolist()).merge(counts, on="cell")
print(gdf.shape, gdf.crs)

(19, 3) EPSG:4326

Set crs="EPSG:4326" when building, not later: the boundary coordinates are degrees on the WGS84 datum, and a GeoDataFrame without a CRS will be silently misread by the next reprojection.

4. Repair cells that cross the antimeridian

fiji = h3.latlng_to_cell(-17.0, 179.99, 5)
naive = shapely.Polygon([(lng, lat) for lat, lng in h3.cell_to_boundary(fiji)])
print(naive.bounds)

(-179.99592759298622, -17.020761499603722, 179.91945817362347, -16.864981747894724)

A cell 0.17° wide has become a polygon from −180° to +180°. Detect these by width (maxx − minx > 180) and split them at 180° into a MultiPolygon — Example 1 does it. On every crossing cell at resolution 5, the split produced 1,545 valid multipolygons out of 1,545 and kept the planar area exactly. The two cells that contain a pole cannot be split, and need a cap built instead.

5. Measure area on the sphere, not in Web Mercator

Calling gdf.area on the degrees raises UserWarning: Geometry is in a geographic CRS. The warning is right, and the tempting fix — to_crs(3857) — is worse than the degrees. On 500 resolution-5 cells over Norway, EPSG:3857 gave a median of 805.7 km² for cells whose real area is 179.0 km², a factor of 4.50. Use h3.cell_area(cell, "km^2"), which is exact on H3's sphere, or a local equal-area projection. Against the WGS84 ellipsoid, cell_area stayed within 0.9912–1.0045 of the geodesic area on 2,000 occupied cells.

6. Dissolve with cell topology, not geometry union

from shapely.geometry import shape

import time

import geopandas as gpd
import shapely

countries = gpd.read_file("ne_10m_admin_0_countries.zip")
switzerland = countries.loc[countries.ADMIN == "Switzerland", "geometry"].iloc[0]
switzerland_cells = h3.geo_to_cells(switzerland, 8)

start = time.perf_counter()
outline = shape(h3.cells_to_geo(switzerland_cells))
topology = time.perf_counter() - start
start = time.perf_counter()
unioned = shapely.union_all(cells_to_gdf(switzerland_cells).geometry.values)
union = time.perf_counter() - start
print(f"{len(switzerland_cells):,} cells: cells_to_geo {topology:.2f} s, union_all {union:.2f} s; "
      f"vertices {len(outline.exterior.coords):,} and {len(unioned.exterior.coords):,}")

58,233 cells: cells_to_geo 0.07 s, union_all 1.05 s; vertices 3,513 and 3,513

cells_to_geo returns a GeoJSON-style dictionary in (lng, lat) order — unlike cell_to_boundary — so shapely reads it directly. On the 58,233 resolution-8 cells of Switzerland it took 0.06 s; union_all on the same polygons took 1.00 s and produced the same 3,513-vertex outline. H3 knows which cell edges are shared; a geometry union has to discover it.

7. Build WKT in DuckDB only when the cells already live there

The DuckDB H3 extension returns the boundary as WKT in (lng, lat) order:

select h3_cell_to_boundary_wkt(cell) as wkt from occupied_cells;

Measured on the same 500,000 cells, including GeoSeries.from_wkt, it took 3.22 s — slower than the vectorised Python build, because the WKT has to be written as text and parsed again. It does not repair the antimeridian either: the Fiji cell came back 359.9° wide. It is worth using when the aggregation already runs in DuckDB and only the output needs geometry.

The antimeridian check adds almost nothing: it tests every polygon's width, and only the few that cross need work.

Cells to polygons with antimeridian and poles handled

The data used on this page

Country outlines from Natural Earth's 1:10m admin-0 layer; every feature in the GeoNames allCountries dump as a point. GeoNames is edited daily, so counts drift a little between downloads.

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")
if not os.path.exists("geonames_all.parquet"):          # every GeoNames feature, about 400 MB zipped
    fetch("https://download.geonames.org/export/dump/allCountries.zip", "allCountries.zip")
    with zipfile.ZipFile("allCountries.zip") as z, z.open("allCountries.txt") as f:
        pd.read_csv(f, sep="\t", header=None, usecols=[4, 5], names=["lat", "lon"],
                    quoting=3).to_parquet("geonames_all.parquet")
print(len(pd.read_parquet("geonames_all.parquet", columns=["lat"])))

13472317

Example 1 — cells to polygons, antimeridian and poles handled

The function, then run on every cell at three resolutions:

import geopandas as gpd
import h3
import numpy as np
import shapely
from shapely.geometry import Polygon, box


def _split_at_180(poly):
    shifted = Polygon([(x + 360 if x < 0 else x, y) for x, y in poly.exterior.coords])
    east = shifted.intersection(box(0, -90, 180, 90))
    west = shapely.affinity.translate(shifted.intersection(box(180, -90, 360, 90)), xoff=-360)
    return shapely.union_all([east, west])


def _polar_cap(cell):
    pts = sorted((lng, lat) for lat, lng in h3.cell_to_boundary(cell))
    pole = 90.0 if pts[0][1] > 0 else -90.0
    (x0, y0), (x1, y1) = pts[-1], (pts[0][0] + 360, pts[0][1])
    edge = y0 + (y1 - y0) * (180 - x0) / (x1 - x0)
    return Polygon([(-180, pole), (-180, edge), *pts, (180, edge), (180, pole)])


def cells_to_gdf(cells, data=None, fix_antimeridian=True):
    """H3 cells (str or int) to a GeoDataFrame of polygons in EPSG:4326."""
    cells = list(cells)
    ids = [h3.int_to_str(c) if isinstance(c, (int, np.integer)) else c for c in cells]
    bounds = [h3.cell_to_boundary(c) for c in ids]
    lengths = np.fromiter((len(b) for b in bounds), dtype=np.int64, count=len(bounds))
    latlng = np.array([pt for b in bounds for pt in b])
    rings = shapely.linearrings(latlng[:, ::-1],
                                indices=np.repeat(np.arange(len(bounds)), lengths))
    geoms = shapely.polygons(rings)

    if fix_antimeridian and len(geoms):
        xmin, _, xmax, _ = shapely.bounds(geoms).T
        for i in np.flatnonzero(xmax - xmin > 180):
            lat, _ = h3.cell_to_latlng(ids[i])
            pole = h3.latlng_to_cell(90 if lat > 0 else -90, 0, h3.get_resolution(ids[i]))
            geoms[i] = _polar_cap(ids[i]) if ids[i] == pole else _split_at_180(geoms[i])

    gdf = gpd.GeoDataFrame({"cell": cells}, geometry=geoms, crs="EPSG:4326")
    if data is not None:
        gdf = gdf.merge(data, on="cell", how="left")
    return gdf


import time

for res in (0, 2, 4):
    every = list(h3.uncompact_cells(h3.get_res0_cells(), res))
    start = time.perf_counter()
    built = cells_to_gdf(every)
    kinds = built.geom_type.value_counts()
    print(f"res {res}: {len(every):>7,} cells  {time.perf_counter() - start:.2f} s  "
          f"invalid {(~built.is_valid).sum()}  Polygon {kinds.get('Polygon', 0):<8,} "
          f"MultiPolygon {kinds.get('MultiPolygon', 0)}")

res 0:     122 cells  0.01 s  invalid 0  Polygon 107      MultiPolygon 15
res 2:   5,882 cells  0.02 s  invalid 0  Polygon 5,796    MultiPolygon 86
res 4: 288,122 cells  1.57 s  invalid 0  Polygon 287,527  MultiPolygon 595

Before the polar cap was added, the south-pole cell at resolution 4 came out as a self-intersecting polygon — one invalid geometry in 288,122, which is exactly the kind that breaks an overlay much later. The caps were valid at every resolution tested: 0, 2, 4, 5 and 7.

Example 2 — areas you can trust, with a built-in cross-check

def add_cell_areas(gdf, cell_col="cell"):
    """Exact spherical cell area from H3, checked against an equal-area projection."""
    out = gdf.copy()
    out["area_km2"] = [h3.cell_area(c if isinstance(c, str) else h3.int_to_str(c), "km^2")
                       for c in out[cell_col]]
    lon, lat = out.geometry.union_all().centroid.coords[0]
    laea = f"+proj=laea +lat_0={lat:.4f} +lon_0={lon:.4f} +units=m"
    projected = out.to_crs(laea).area / 1e6
    worst = (projected / out["area_km2"] - 1).abs().max()
    print(f"{len(out):,} cells; h3 vs local equal-area projection differ by at most {worst:.2%}")
    return out


norway = countries.loc[countries.ADMIN == "Norway", "geometry"].iloc[0]
swiss7 = add_cell_areas(cells_to_gdf(h3.geo_to_cells(switzerland, 7)))
add_cell_areas(cells_to_gdf(h3.geo_to_cells(norway, 6)))
print(f"Swiss resolution-7 cells: {swiss7.area_km2.min():.3f} to {swiss7.area_km2.max():.3f} km²")

8,317 cells; h3 vs local equal-area projection differ by at most 0.29%
14,447 cells; h3 vs local equal-area projection differ by at most 0.86%
Swiss resolution-7 cells: 4.858 to 5.095 km²

The first is Switzerland at resolution 7, the second Norway at resolution 6, where the country's extent stretches a single projection further. Switzerland's resolution-7 cells ranged from 4.858 to 5.095 km² — so dividing a count by a constant area is a small error there, and a much larger one across a continent.

Example 3 — one outline per group, from topology

import pandas as pd
from shapely.geometry import shape


def dissolve_cells(frame, by, cell_col="cell"):
    """One outline per group, built from the cell topology rather than geometry union."""
    rows = []
    for key, group in frame.groupby(by):
        geo = h3.cells_to_geo(group[cell_col].tolist())
        rows.append({by: key, "cells": len(group), "geometry": shape(geo)})
    return gpd.GeoDataFrame(rows, geometry="geometry", crs="EPSG:4326")

On the 26 Swiss cantons filled at resolution 8 (58,233 cells): 0.04 s, returning 19 polygons and 7 multipolygons. GeoDataFrame.dissolve on polygons built from the same cells took 1.08 s. A group whose cells form a ring comes back as a polygon with a hole, and separated groups as a MultiPolygon, so no post-processing is needed.

Example 4 — the measurements quoted on this page

import duckdb
import pandas as pd
from pyproj import Geod

res5 = list(h3.uncompact_cells(h3.get_res0_cells(), 5))
vertices = pd.Series([len(h3.cell_to_boundary(c)) for c in res5]).value_counts().sort_index()
print(f"{len(res5):,} resolution-5 cells by vertex count:", vertices.to_dict())
res1 = list(h3.uncompact_cells(h3.get_res0_cells(), 1))
print(f"resolution 1: {sum(len(h3.cell_to_boundary(c)) > 6 for c in res1) / len(res1):.0%} have more than six")
naive = shapely.polygons(shapely.linearrings(
    np.array([pt for c in res5 for pt in h3.cell_to_boundary(c)])[:, ::-1],
    indices=np.repeat(np.arange(len(res5)), [len(h3.cell_to_boundary(c)) for c in res5])))
wide = np.flatnonzero(np.ptp(shapely.bounds(naive)[:, [0, 2]], axis=1) > 180)
fixed = cells_to_gdf([res5[i] for i in wide])
print(f"{len(wide):,} naive polygons wider than 180°; after the fix "
      f"{fixed.geom_type.value_counts().to_dict()}, all valid: {fixed.is_valid.all()}")

points = pd.read_parquet("geonames_all.parquet").dropna()
occupied = pd.Series([h3.latlng_to_cell(a, b, 7) for a, b in
                      zip(points.lat.tolist()[:3_000_000], points.lon.tolist()[:3_000_000])]).unique()[:500_000]
start = time.perf_counter()
looped = [Polygon([(lng, lat) for lat, lng in h3.cell_to_boundary(c)]) for c in occupied]
loop = time.perf_counter() - start
start = time.perf_counter()
cells_to_gdf(occupied, fix_antimeridian=False)
vector = time.perf_counter() - start
start = time.perf_counter()
cells_to_gdf(occupied)
vector_fix = time.perf_counter() - start
con = duckdb.connect()
con.execute("install h3 from community; load h3")
con.register("occupied_cells", pd.DataFrame({"cell": occupied}))
start = time.perf_counter()
gpd.GeoSeries.from_wkt(con.execute("select h3_cell_to_boundary_wkt(cell) as wkt from occupied_cells").df().wkt)
wkt = time.perf_counter() - start
print(f"{len(occupied):,} occupied resolution-7 cells: Python loop {loop:.2f} s, vectorised {vector:.2f} s, "
      f"with the 180° fix {vector_fix:.2f} s, DuckDB WKT + parse {wkt:.2f} s")

geod = Geod(ellps="WGS84")
res5_occ = pd.Series([h3.latlng_to_cell(a, b, 5) for a, b in zip(points.lat, points.lon)]).unique()
cells5 = cells_to_gdf(res5_occ)
true = np.array([h3.cell_area(c, "km^2") for c in cells5.cell])
mercator = cells5.to_crs(3857).area.to_numpy() / 1e6
centre_lat = np.array([h3.cell_to_latlng(c)[0] for c in cells5.cell])
ok = np.abs(centre_lat) < 85.0511                        # Web Mercator's own latitude limit
print(f"Web Mercator on {ok.sum():,} occupied resolution-5 cells: median {np.median(mercator[ok] / true[ok]):.2f}x, "
      f"max {np.max(mercator[ok] / true[ok]):.1f}x")
over_norway = cells_to_gdf(h3.geo_to_cells(norway, 5)).sample(500, random_state=0)
m = over_norway.to_crs(3857).area / 1e6
t = over_norway.cell.map(lambda c: h3.cell_area(c, "km^2"))
print(f"500 cells over Norway: EPSG:3857 median {m.median():.1f} km², true {t.median():.1f} km² "
      f"({m.median() / t.median():.2f}x)")
check = cells5.sample(2000, random_state=0)
ratio = [h3.cell_area(c, "km^2") / (abs(geod.geometry_area_perimeter(g)[0]) / 1e6)
         for c, g in zip(check.cell, check.geometry)]
print(f"cell_area / WGS84 geodesic area on 2,000 occupied cells: {min(ratio):.4f} .. {max(ratio):.4f}")

2,016,842 resolution-5 cells by vertex count: {6: 2002160, 7: 11760, 8: 2910, 10: 12}
resolution 1: 33% have more than six
1,547 naive polygons wider than 180°; after the fix {'MultiPolygon': 1545, 'Polygon': 2}, all valid: True
500,000 occupied resolution-7 cells: Python loop 5.56 s, vectorised 2.56 s, with the 180° fix 2.52 s, DuckDB WKT + parse 4.23 s
Web Mercator on 450,884 occupied resolution-5 cells: median 1.35x, max 134.5x
500 cells over Norway: EPSG:3857 median 805.7 km², true 179.0 km² (4.50x)
cell_area / WGS84 geodesic area on 2,000 occupied cells: 0.9912 .. 1.0045

Latitude-first boundaries and stretched cells

Why H3 hands back latitude first

H3's C API is written in the geographic tradition — latLngToCell(lat, lng) — and the Python binding keeps that order for every function that takes or returns individual coordinates. The functions that speak GeoJSON (geo_to_cells, cells_to_geo, h3shape.__geo_interface__) use GeoJSON's (lng, lat).

So the rule is by function family, not by library: point-level functions are (lat, lng), GeoJSON-level functions are (lng, lat). Mixing the two in one pipeline is where the Bern-to-Ethiopia bug comes from.

Why some hexagons have seven, eight or ten vertices

H3 projects each of the 20 faces of an icosahedron onto the sphere separately. A cell drawn across a face boundary is bent where the two projections meet, and H3 inserts a vertex at each crossing so the drawn edge follows the real shape.

That only happens at class III resolutions (odd numbers), whose grids are rotated relative to the face edges. At resolution 5, 0.73% of cells had more than six vertices; at even resolutions none did, apart from the twelve five-sided pentagons. Code that indexes boundary[5] or reshapes to six columns works on the test set and fails on the first cell near a face edge.

Why a small cell stretches across the whole map

A polygon in longitude/latitude has no idea that −179.99° and +179.92° are neighbours. Shapely draws the edge the long way round, 359.9° wide, and matplotlib paints a band across the map.

The data itself is not wrong. The geodesic area of the naive Fiji polygon, computed with pyproj's Geod, was 231.77 km² — correct, because geodesic edges take the short route. Only planar operations are wrong: plotting, intersects, clip, and anything else in shapely. Splitting at 180° makes the geometry honest for those operations without changing what it covers.

Why Web Mercator is the wrong place to measure cells

The 4.5× figure is the median measured over Norway, not a worst case.

EPSG:3857 preserves angles and scales area by the square of 1/cos(latitude). A cell at 60° north is drawn four times larger than one of the same size at the equator — close to the 4.5× median measured over Norway.

H3 cells are nearly equal in area already, so any per-area figure should use the cell's own area. h3.cell_area computes it on the sphere, and its disagreement with the WGS84 ellipsoid (under 1%) is far below the variation between cells.

Why topology is faster than union

union_all treats 58,233 hexagons as unrelated polygons: it has to find which edges coincide, snap them and remove them, in floating point. cells_to_geo already knows — two cells either are neighbours or are not — so it only walks the outer edges.

Neighbouring cells do share identical vertex coordinates, which is why the plain union matched the topological outline here. The instinct to make the union "safer" by snapping with grid_size=1e-9 made it worse: 5.01 s, and 3,525 vertices instead of 3,513.

The geodesic area of the naive polygon was already correct; only planar operations are fooled.

Order per function, integer cells and polar cells

cell_to_boundary is (lat, lng); cells_to_geo is (lng, lat). Swap the first, never the second.
Integer cells need h3.int_to_str before the string API, or use h3.api.basic_int throughout.
Two cells per resolution contain a pole. They cannot be split at 180°; build a cap to the pole, as Example 1 does.
The DuckDB boundary function does not repair the antimeridian. Apply the same split to its output.
Do not simplify cell polygons. Simplifying breaks the shared edges, and the dissolve and the map both develop gaps.
GeoParquet keeps the CRS. to_parquet then read_parquet returned EPSG:4326 intact; CSV with WKT does not.
Set the CRS at construction. A GeoDataFrame of degrees with no CRS will be reprojected as if it were anything.

Internal links

H3 explained: how a hexagonal grid indexes the whole Earth — the icosahedron behind the extra vertices
How to fill a polygon with H3 cells — where most cell sets come from
How to aggregate points into H3 hexagons and map them — the map this geometry is for
Fixing H3 hexagons that stretch across the map at the antimeridian — the 359.9° cell in depth
How to use H3 in DuckDB for grid aggregation at scale — the WKT route in context
How to replace a row loop with vectorised GeoPandas operations — the general version of step 2
How to calculate area and distance in GeoPandas correctly — why EPSG:3857 is wrong for area
How to dissolve polygons in Python — the geometry route that topology replaces
WKT, WKB and GeoJSON: how geometry is actually stored — why the text round trip is slow

FAQ

What is the fastest way to build polygons for many H3 cells?

Flatten all boundary coordinates into one array and call shapely.linearrings with an indices array, then shapely.polygons. For 500,000 cells that took 1.72 s against 4.34 s for a Python loop.

Why does one hexagon stretch across my whole world map?

It crosses the antimeridian, so its longitudes jump from about −180° to +180°. Split it at 180° into a multipolygon; 1,547 resolution-5 cells need it.

Why do some H3 cells have more than six vertices?

At odd resolutions, cells crossing an icosahedron face edge get extra vertices where the projection bends them. At resolution 5, 0.73% of cells have seven, eight or ten.

How should I calculate the area of H3 cells?

Use h3.cell_area(cell, "km^2") or a local equal-area projection. Web Mercator overstated cells over Norway by a factor of 4.5.

How do I merge H3 cells into one outline?

Use h3.cells_to_geo and shapely.geometry.shape. It dissolved 58,233 Swiss cells in 0.06 s, against 1.00 s for a geometry union, with an identical result.