Figuring out how to generate an H3 hex grid from a given poly layer, in Qgis python processing environment.

h3-grid-from-poly-test-01.py

# adapted from https://github.com/ThomasG77/30DayMapChallenge/blob/master/day4_hexagons/data/h3-processing.py
import os

from qgis.utils import iface
from qgis.core import (
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsFeature, QgsField, QgsFields,
    QgsGeometry, QgsPointXY, QgsProject,
    QgsProcessingFeedback, QgsMessageLog,
    QgsVectorFileWriter, QgsVectorLayer, QgsWkbTypes)
from qgis.PyQt.QtCore import QVariant
import processing
from h3 import h3

log = QgsMessageLog.logMessage

# https://docs.qgis.org/3.16/en/docs/pyqgis_developer_cookbook/crs.html#crs-transformation
nad83csrs = QgsCoordinateReferenceSystem("EPSG:4617")    # NAD83 (CSRS)
ytalbcsrs = QgsCoordinateReferenceSystem("EPSG:3579")   # YT Albers NAD83 (CSRS)
transformContext = QgsProject.instance().transformContext()
#xform = QgsCoordinateTransform(crsSrc, crsDest, transformContext)

projectPath = os.path.dirname(QgsProject.instance().fileName())
accumulator = {}
r = range(0, 12)
for i in r:
    accumulator[str(i)] = []

mylayer = QgsProject.instance().mapLayersByName('nts-50k')[0]

# h3 assumes all inputs are in geographic lat-long (note reversed X,Y) so convert
#from https://gis.stackexchange.com/questions/316002/pyqgis-reproject-layer
#https://gis.stackexchange.com/questions/350385/pyqgis-script-failing-to-project-layer-to-epsg54002
params = {'INPUT': mylayer, 'TARGET_CRS': nad83csrs,
    'OUTPUT': 'memory:'}
dd_lyr = processing.run('native:reprojectlayer', params)['OUTPUT']
#QgsProject.instance().addMapLayer(dd_lyr)

h3.polyfill([(60,-120),(62,-122)]

for feature in dd_lyr.getFeatures():
    # if point:
        #x = feature.geometry().asPoint()[0]
        #y = feature.geometry().asPoint()[1]
    # if polygon:
    x = feature.geometry().centroid().asPoint()[0]
    y = feature.geometry().centroid().asPoint()[1]
    #log(f"In centroids {x},{y}") # debug
    for i in r:
        # note order is Y,X!
        accumulator[str(i)].append(h3.geo_to_h3(y, x, i)) 

for i in r:
    h_name = f"h_{i:0>2}"
    accumulator[str(i)] = set(accumulator[str(i)])
    fields = QgsFields()
    fields.append(QgsField("id", QVariant.String))
    shpfile = os.path.join(projectPath, f"data/{h_name}.shp")
    writer = QgsVectorFileWriter(shpfile, "UTF8", fields, QgsWkbTypes.Polygon,
        driverName="ESRI Shapefile"
    )
    features = []
    for j in accumulator[str(i)]:
        f = QgsFeature()
        #log("Out hex:" + str(h3.h3_to_geo_boundary(j))) # debug
        f.setGeometry(QgsGeometry.fromPolygonXY([
            # note reversing back to X,Y
            [QgsPointXY(c[1], c[0]) for c in h3.h3_to_geo_boundary(j)] 
        ]))
        f.setAttributes([j])
        features.append(f)
    writer.addFeatures(features)
    processing.run('qgis:definecurrentprojection',{'INPUT': shpfile,
        'CRS': nad83csrs} )

for i in r:
    h_name = f"h_{i:0>2}"
    layer = QgsVectorLayer(
        os.path.join(projectPath, f"data/{h_name}.shp"),
        f"Hex L-{i:0>2}", "ogr"
    )
    QgsProject.instance().addMapLayer(layer)


def count_points_in_hex(h3_layer, in_layer, out_layer):
    ''' Run qgis count points in poly for the specified hex level
        layers: path to data file, e.g. "data/mypoints.shp"
    '''
    processing.run('qgis:countpointsinpolygon', {
        'CLASSFIELD' : None,
        'FIELD' : 'numpoints',
        'POLYGONS': h3_layer, 
        'OUTPUT' : out_layer,
        'POINTS' : in_layer,
        'WEIGHT' : None
    })

def count_and_add_to_map():
    # Count points in polygon and add to map
    for i in r:
        h_name = f"h3_{level:0>2}"
        count_pounts_in_hex(f"data/{h_name}.shp", 
            dd_layer,
            f"data/{h_name}_count.shp")
            
        layer = QgsVectorLayer(f"data/{h_name}_count.shp", 
            f"Hex count {i:0>2}", "ogr")
        QgsProject.instance().addMapLayer(layer)

# Uncomment to run
#count_and_add_to_map()
    


## --- Scrapbook
    
def proj_to_albers(layer, out_name):
        '''Project to Yukon Albers'''
        params = {'INPUT': layer, 'TARGET_CRS': 'EPSG:3579',
            'OUTPUT': f'memory:{out_name}'}
        vlayer = processing.run('native:reprojectlayer', params)['OUTPUT']
        QgsProject.instance().addMapLayer(vlayer)