H3-Python: How to generate hex bins over shapefile?

Question

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
APERTURE_SIZE = 6

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

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

# 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
dfh.plot(x='lng',y='lat',style='.',figsize=(48,48));

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.