How to convert H3 cell boundaries to Shapely polygons in Python

Uber H3 grid over Boulder, CO

I frequently use Uber’s H3 grid system for my data analysis tasks because it makes aggregation fast and easy. In this post, I will demonstrate how to convert an H3 cell into a Shapely geometry in Python.

Gotcha #1: Python version

I use conda to manage my Python environments, and at time of writing, the latest version of h3-py on conda-forge is 3.7.4. However, version 4.0.0 is in beta, and several breaking changes are slated.

Introducing H3 Version 4.0.0

Be sure to check your version and keep these changes in mind before proceeding.

Gotcha #2: Latitude = y, Longitude = x

Ahh, yes… The all time gotcha; number two in this article, but number one in the hearts of Geographers. You spend your whole life learning about Euclidean space in terms of x then y, and then you’re expected to make the mental flip to lat then long. *le sigh*

With those two caveats out of the way, we can begin!

Sample Data

I first ran into this issue while analyzing emergency response times for the fire department at the City of Boulder, CO. Data downloads can be done through their Open Data Portal, but here’s a randomly generated sample of data to get you going.

from shapely.geometry import Point
import numpy as np
import random

# Define the range for ELAPSETIME
start = 0
stop = 500

# Define the bounding box for Boulder, CO
minx, miny = -105.301758, 39.964069
maxx, maxy = -105.178925, 40.094555

# Generate a list of random lat-long points within the bounding box
points = [Point(random.uniform(minx, maxx), random.uniform(miny, maxy)) for _ in np.arange(10000)]

# Create a geodataframe with the specified columns and attributes
gdf = gpd.GeoDataFrame({'ELAPSETIME': [random.randint(start, stop) for _ in np.arange(10000)],
                        'geometry': points}, crs='EPSG:4326')

H3 Cell Identifier

Let’s enrich the dataframe by labelling each record with the unique identifier of the hex in which it resides. To do this, we use the geo_to_h3 function (renamed latlng_to_cell in v4.0.0):

import h3

# x = longitude, y = latitude!!!
res = 10
col = f"H3_{res}"
gdf[col] = gdf.apply(lambda row: str(h3.geo_to_h3(row.geometry.y, row.geometry.x, res)), axis=1)

Aggregate by H3 cell

With labelling out of the way, let’s aggregate the data based on those labels!

h3_df = gdf.groupby(col)['ELAPSETIME'].describe().reset_index()

Convert to Shapely Polygon

We need the boundary of a hex in order to convert it to a Shapely geometry object. Luckily, we can do this with the H3 library’s h3_to_geo_boundary function (renamed cell_to_boundary in v4.0.0). This returns a tuple-of-tuples representing a ring of (lat, long) coordinate pairs for a specific cell identifier, with identical coordinates in the first and last positions.

Yay, we have a representation of a polygon as coordinate pairs! Just shove that into a shapely Polygon and we’re good, right? WRONG!! Remember Gotcha #2! Lat = y, Long =x, and shapely Polygon expects pairs as xy. Let’s define a simple function to flip each coordinate pair so they are (long, lat) pairs:

from shapely.geometry import Polygon

def cell_to_shapely(cell):
    coords = h3.h3_to_geo_boundary(cell)
    flipped = tuple(coord[::-1] for coord in coords)
    return Polygon(flipped)

Then we create a brand new geodataframe from the aggregated information:

h3_geoms = h3_df[col].apply(lambda x: cell_to_shapely(x))
h3_gdf = gpd.GeoDataFrame(data=h3_df, geometry=h3_geoms, crs=4326)

There you have it! From here, you can keep manipulating the data in Python, or export and play around in your favorite GIS software. Here’s what actual median response times look like after following this process and importing into QGIS:

Median Response Times for the Boulder Fire Department.

Thanks for reading!