from shapely.geometry import MultiPoint
## get reference to the project instance
proj = QgsProject.instance()
#############
## using existing map layers
#bnd = proj.mapLayersByName('boundaries')[0]
#pnt = proj.mapLayersByName('points')[0]
#############
## for a reproducible example
extent = '-0.666666667,-0.166666667,-0.300000000,0.300000000 [EPSG:4326]'
## create an arbitrary grid to represent boundaries
bnd = processing.run("native:creategrid", {'TYPE':2,'EXTENT':extent,'HSPACING':0.1,'VSPACING':0.1,'HOVERLAY':0,'VOVERLAY':0,'CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'OUTPUT':'TEMPORARY_OUTPUT'})['OUTPUT']
## make random points within the grid
pnt = processing.run("native:randompointsinextent", {'EXTENT':extent,'POINTS_NUMBER':300,'MIN_DISTANCE':0,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'MAX_ATTEMPTS':200,'OUTPUT':'TEMPORARY_OUTPUT'})['OUTPUT']
## add grid and random points to canvas
proj.addMapLayer(bnd)
proj.addMapLayer(pnt)
bnd_feats = list(bnd.getFeatures())
pnt_feats = list(pnt.getFeatures())
## initialise empty list to hold cluster centroids
centroids = []
## loop through boundary features, find all points within each grid cell and get the centroid of those points
for b in bnd_feats:
points_within = []
for p in pnt_feats:
if p.geometry().within(b.geometry()):
## append point coordinates to list
g = p.geometry().asPoint()
points_within.append((g.x(), g.y()))
## make a multi point of all points within a cell
points = MultiPoint(points_within)
## get the centroid of the multi point
centroids.append(points.centroid)
## convert shapely objects to QgsPointXY objects
centroid_pnts = [QgsPointXY(point.x, point.y) for point in centroids]
## make empty point layer
centroids_lyr = QgsVectorLayer("Point", "centroids", "memory")
## get reference to data provider
prov = centroids_lyr.dataProvider()
## enter editing mode
centroids_lyr.startEditing()
## add features to empty layer
for point in centroid_pnts:
feat = QgsFeature()
feat.setGeometry( QgsGeometry.fromPointXY(point))
prov.addFeatures( [ feat ] )
## save changes
centroids_lyr.commitChanges()
## name the output layer and add it to the map canvas
centroids_lyr.setName('centroids')
proj.addMapLayer(centroids_lyr)
Bear in mind, for an irregular shaped polygon, the centroid might be outside the boundary. In which case you could use points.representative_point() instead of points.centroid, which would give you the existing point closest to the centroid.
Bear in mind, for an irregular shaped polygon, the centroid might be outside the boundary. In which case you could use points.representative_point() instead of points.centroid, which would give you the existing point closest to the centroid.