Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
Figuring out how to generate an H3 hex grid from a given poly layer, in Qgis python processing environment.
# adapted from
import os
from qgis.utils import iface
from qgis.core import (
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
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
params = {'INPUT': mylayer, 'TARGET_CRS': nad83csrs,
'OUTPUT': 'memory:'}
dd_lyr ='native:reprojectlayer', params)['OUTPUT']
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
# note reversing back to X,Y
[QgsPointXY(c[1], c[0]) for c in h3.h3_to_geo_boundary(j)]
writer.addFeatures(features)'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"
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"
''''qgis:countpointsinpolygon', {
'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}"
layer = QgsVectorLayer(f"data/{h_name}_count.shp",
f"Hex count {i:0>2}", "ogr")
# Uncomment to run
## --- 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 ='native:reprojectlayer', params)['OUTPUT']
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment