Skip to content

Instantly share code, notes, and snippets.

Embed
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 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment