I want to create hexagons on my geographic map and want to preserve the digital boundary specified by the shapefile/geojson as well. How do I do it using uber's h3 python library?
I'm new to shapefiles or any other geographic data structures. I'm most comfortable with python.
For convenience, you can use H3-Pandas.
import geopandas as gpd
import h3pandas
# Prepare data
gdf = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
gdf = gdf.loc[gdf['continent'].eq('Africa')]
gdf['gdp_md_per_capita'] = gdf['gdp_md_est'].div(gdf['pop_est'])
resolution = 4
# Resample to H3 cells
gdf.h3.polyfill_resample(resolution)
Use cases like this are precisely why I wrote H3-Pandas, an integration of H3 with Pandas and GeoPandas.
Let's use the naturalearth_lowres dataset included in GeoPandas to demonstrate the usage.
import geopandas as gpd
import h3pandas
path_shapefile = gpd.datasets.get_path('naturalearth_lowres')
gdf = gpd.read_file(path_shapefile)
Let's also create a numeric column for plotting.
gdf = gdf.loc[gdf['continent'].eq('Africa')]
gdf['gdp_md_per_capita'] = gdf['gdp_md_est'].div(gdf['pop_est'])
ax = gdf.plot(figsize=(15, 15), column='gdp_md_per_capita', cmap='RdBu')
ax.axis('off')
We will use H3 resolution 4. See the H3 resolution table for more details.
resolution = 4
We can add H3 hexagons using the function polyfill. This method adds a list of H3 cells whose centroid falls into each polygon.
gdf_h3 = gdf.h3.polyfill(resolution)
print(gdf_h3['h3_polyfill'].head(3))
1 [846aca7ffffffff, 8496b5dffffffff, 847b691ffff...
2 [84551a9ffffffff, 84551cdffffffff, 8455239ffff...
11 [846af51ffffffff, 8496e63ffffffff, 846a803ffff...
Name: h3_polyfill, dtype: object
If we want to explode the values horizontally (ending up with as many rows as there are H3 cells), we can use the parameter explode
gdf_h3 = gdf.h3.polyfill(resolution, explode=True)
print(gdf_h3.head(3))
pop_est continent name iso_a3 gdp_md_est \
1 53950935 Africa Tanzania TZA 150600.0
1 53950935 Africa Tanzania TZA 150600.0
1 53950935 Africa Tanzania TZA 150600.0
geometry gdp_md_per_capita \
1 POLYGON ((33.90371 -0.95000, 34.07262 -1.05982... 0.002791
1 POLYGON ((33.90371 -0.95000, 34.07262 -1.05982... 0.002791
1 POLYGON ((33.90371 -0.95000, 34.07262 -1.05982... 0.002791
h3_polyfill
1 846aca7ffffffff
1 8496b5dffffffff
1 847b691ffffffff
We can then utilize the method h3_to_geo_boundary to obtain the geometries for the H3 cells. It expects that the index already has the H3 cell ids.
gdf_h3 = gdf_h3.set_index('h3_polyfill').h3.h3_to_geo_boundary()
We can now plot the result
ax = gdf_h3.plot(figsize=(15, 15), column='gdp_md_per_capita', cmap='RdBu')
ax.axis('off')
H3-Pandas actually has convenience function that performs all this at once: polyfill_resample
gdf_h3 = gdf.h3.polyfill_resample(resolution)
ax = gdf_h3.plot(figsize=(15, 15), column='gdp_md_per_capita', cmap='RdBu')
ax.axis('off')
h3-py doesn't know how to work with shapefile data directly, but it sounds like you could use a library like https://github.com/GeospatialPython/pyshp to convert you shapefile data to GeoJSON, and then use h3.polyfill() to convert to a collection of H3 hexagons.
There are lots of options for plotting your boundary along with the H3 hexagons. For example, you might use pydeck and its GeoJsonLayer and H3HexagonLayer layers.
If your plotting software can't use H3 hexagons directly, you can convert them to other formats with functions like h3.h3_to_geo_boundary() or h3.h3_set_to_multi_polygon().
Converting a GeoJSON to Uber's h3 is quite simple.
Attaching a sample code snippet and GeoJSON used below:
GeoJSON:
{"type": "FeatureCollection","features": [{"type": "Feature","properties": {},"geometry": {"type": "Polygon", "coordinates": [[[77.5250244140625,13.00857192009273],[77.51266479492188,12.971103764892034],[77.52777099609375,12.94099133483504],[77.57171630859375,12.907528813968185],[77.60604858398438,12.914890953258695],[77.662353515625,12.928276105253065],[77.69874572753906,12.961066692801282],[77.65823364257812,13.00990996390651],[77.58956909179688,13.04469656691112],[77.53944396972656,13.038007215169166],[77.5250244140625,13.00857192009273]]]}}]}
Code:
from h3converter import h3converter
geojson_raw = open("sampleJson.geojson",)
geojson = json.load(geojson_raw)
h3_list = []
h3_resolution = 7
for feature in geojson["features"]:
feature_h3 = h3converter.polyfill(feature["geometry"], h3_resolution)
h3_list.append(feature_h3)
print("H3's created => ", len(feature_h3))
print(h3_list)
Response:
[{'8761892edffffff', '8760145b0ffffff', '87618925effffff', '87618925bffffff', '87618924affffff', '876189256ffffff', '8760145b1ffffff', '8761892eaffffff', '8761892eeffffff', '876189253ffffff', '876189259ffffff', '8760145b2ffffff', '876014586ffffff', '8760145b4ffffff', '8761892e1ffffff', '8760145a2ffffff', '8761892ecffffff', '876189251ffffff', '8760145a4ffffff', '8761892e5ffffff', '87618925affffff', '8761892e9ffffff', '8761892cdffffff', '876189250ffffff', '87618925dffffff', '8760145b6ffffff', '876014595ffffff', '876189252ffffff', '8761892ebffffff', '8760145a3ffffff', '8760145a6ffffff', '876014584ffffff', '876189258ffffff', '8760145b5ffffff', '8760145b3ffffff', '876014594ffffff', '8761892c9ffffff', '87618925cffffff', '8760145a0ffffff', '8761892e8ffffff'}]
Package used: https://pypi.org/project/h3converter/
You could use the above package to convert h3 to GeoJSON as well.
You can extract GeoJSON from Shapley like this:
import geopandas as gpd
import shapely.geometry
# Create Geometry
shapely_polygon = shapely.geometry.Polygon([(0, 0), (0, 1), (1, 0)])
# Extract JSON Feature Collection
gpd.GeoSeries([shapely_polygon]).__geo_interface__
Output:
{"type":"FeatureCollection","features":[{"id":"0","type":"Feature","properties":{},"geometry":{"type":"Polygon","coordinates":[[[0.0,0.0],[0.0,1.0],[1.0,0.0],[0.0,0.0]]]},"bbox":[0.0,0.0,1.0,1.0]}],"bbox":[0.0,0.0,1.0,1.0]}
