I am trying to generate hexbins over my shapefile to eventually cluster other geospatial events to them using H3. Here is my process, but I am wondering if there is a better way.

# packages
from h3 import h3
import pandas as pd
import geopandas as gpd
from shapely import geometry, ops
from shapely.geometry.polygon import Polygon
from shapely.geometry.multipolygon import MultiPolygon
from collections import Counter
from shapely.ops import cascaded_union

# transform multipolygon to polygon
def explode(indata):
    indf = gpd.GeoDataFrame.from_file(indata)
    outdf = gpd.GeoDataFrame(columns=indf.columns)
    for idx, row in indf.iterrows():
        if type(row.geometry) == Polygon:
            outdf = outdf.append(row,ignore_index=True)
        if type(row.geometry) == MultiPolygon:
            multdf = gpd.GeoDataFrame(columns=indf.columns)
            recs = len(row.geometry)
            multdf = multdf.append([row]*recs,ignore_index=True)
            for geom in range(recs):
                multdf.loc[geom,'geometry'] = row.geometry[geom]
            outdf = outdf.append(multdf,ignore_index=True)
    return outdf

# transform multipolygon to polygon
explode_vals = explode('C:\\Users\\Andrew\Desktop\\af_shape\\9c21fbef-c5ee-4eb6-88f3-ae4b938014bf202049-1-xml0al.1eimn.shp')

# set hex size

# Unify the boundries
union_poly = cascaded_union(explode_vals.geometry)

# Find the hexs within the shape boundary using PolyFill
for n,g in enumerate(union_poly):
    temp  = geometry.mapping(g)
    temp['coordinates']=[[[j[1],j[0]] for j in i] for i in temp['coordinates']]

# create hex dataframe
hex_col = 'hex{}'.format(APERTURE_SIZE)
dfh = pd.DataFrame(hex_list,columns=[hex_col])
print('Sanity Check\nnumber of hexes:', len(hex_list))
print('number of duplicates:', len(hex_list) - len(dfh.drop_duplicates()))

# add lat & lng of center of hex
dfh['lat']=dfh[hex_col].apply(lambda x: h3.h3_to_geo(x)[0])
dfh['lng']=dfh[hex_col].apply(lambda x: h3.h3_to_geo(x)[1])

# create Point object based on hex latlng
dfh['geometry'] = dfh.apply(lambda x: geometry.Point(x.lng,x.lat),1)
dfh.crs = {"init": "epsg:4326"}

# plot hex latlng

The second problem I have, is, if I want to filter my country level shape file to just a province, I get an error saying that my polygon is not iterable.

I am filtering like so, and following the same process.

explode_vals = explode_vals.loc[explode_vals['PROV_34_NA'] == 'Hilmand'].reset_index(drop=True)

I have also tried to use QGIS to export just the province of interest to geoJSON and a shapefile, but something about cutting it off from the rest of the country makes this whole process not work.

  • 1
    Could you share your shapefile (link etc.)? So that we perhaps may solve the problem. When I use my shapefile I don't get that error. Dec 4, 2020 at 16:53
  • Could GDAL with SpatiaLite dialect be any option for you? Usage example with a geometry that is created with SQL as well ogrinfo -dialect SQLite -sql "SELECT ST_HexagonalGrid(geom, 10) FROM (SELECT ST_Buffer(ST_GeomFromText('LINESTRING(50 50,150 150,150 50)'),10) AS geom) AS s" foo.jml. Ogrinfo and ogr2ogr require that some valid data source exists. Use your own instead of foo.jml for testing.
    – user30184
    Dec 7, 2020 at 8:41
  • I'll look for a reproducible shapefile or perhaps I will just host somewhere since they are all hosted inside zips. I am less familiar with using SQL to accomplish this.
    – John Stud
    Dec 7, 2020 at 15:08
  • Have you tried QGIS Create Grid tool? docs.qgis.org/3.16/en/docs/user_manual/processing_algs/qgis/… Dec 8, 2020 at 5:41
  • I'll give that a shot to see how QGIS does it!
    – John Stud
    Dec 8, 2020 at 15:22


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.