2

I have multiple sets of around 200-300 points and I have the following requirements:

  1. Cluster the points in each set into a group of points with a minimum & maximum number (as defined by me).
  2. In the process the clustering needs to take care not to cross a set of lines, which is in another layer.

How to go about it? I was able to get the clustering part using K Means / DB Scan, but it does not take into consideration of my second requirement

I am using QGIS to get the desired result

Points are in red & lines in blue in sample photo

enter image description here

Expected outcome is given below

enter image description here

0

1 Answer 1

3

If your boundaries were polygons, you could find all points within each polygon and get the centroid of each demarcated group.

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)

There is probably a cleaner way of dealing with the geometries, but this works.

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.

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.