Generate points in multiple features taking minimum distance into account

Question

I have a function which creates Wind Turbines represented as points. Essentially, it uses the code from the Random points inside polygons (fixed) tool albeit with some slight changes.

The goal is to create random points inside polygons taking into account the minimum distance specified. This works very well especially with polygons which are not close to another (e.g. a single polygon):

However, if the polygon is close to or adjoined to another polygon (e.g. as shown below), the points from each polygon can be within the minimum distance as shown in red:

How could I change the code so that those points in red are not close to another from a nearby polygon?

Ideally, I would like multiple points to be replaced by a single point:

---

Here is the code which can be reproduced in the Python Console, a polygon layer must be selected with a relevant CRS before running the function:

import random
from PyQt4.QtCore import QVariant

def checkMinDistance(point, index, distance, points):
    if distance == 0:
        return True
    neighbors = index.nearestNeighbor(point, 1)
    if len(neighbors) == 0:
        return True
    if neighbors[0] in points:
        np = points[neighbors[0]]
        if np.sqrDist(point) < (distance * distance):
            return False
    return True

def generate_wind_turbines(spacing):
    layer = iface.activeLayer()
    crs = layer.crs()
    # Memory layer
    memory_lyr = QgsVectorLayer("Point?crs=epsg:" + unicode(crs.postgisSrid()) + "&index=yes", "Wind turbines for " + str(layer.name()), "memory")
    QgsMapLayerRegistry.instance().addMapLayer(memory_lyr)
    memory_lyr.startEditing()
    provider = memory_lyr.dataProvider()
    provider.addAttributes([QgsField("ID", QVariant.Int)])
    # Variables
    point_density = 0.0001
    fid = 1
    distance_area = QgsDistanceArea()
    # List of features
    fts = []
    # Create points
    for f in layer.getFeatures():
        fGeom = QgsGeometry(f.geometry())
        bbox = fGeom.boundingBox()
        pointCount = int(round(point_density * distance_area.measure(fGeom)))
        index = QgsSpatialIndex()
        points = dict()
        nPoints = 0
        fid += 1
        nIterations = 0
        maxIterations = pointCount * 200
        random.seed()
        while nIterations < maxIterations and nPoints < pointCount:
            rx = bbox.xMinimum() + bbox.width() * random.random()
            ry = bbox.yMinimum() + bbox.height() * random.random()
            pnt = QgsPoint(rx, ry)
            geom = QgsGeometry.fromPoint(pnt)
            if geom.within(fGeom) and checkMinDistance(pnt, index, spacing, points):
                f = QgsFeature(nPoints)
                f.setAttributes([fid])
                f.setGeometry(geom)
                fts.append(f)
                index.insertFeature(f)
                points[nPoints] = pnt
                nPoints += 1
            nIterations += 1
    provider.addFeatures(fts)
    memory_lyr.updateFields()
    memory_lyr.commitChanges()

generate_wind_turbines(500)

---

Edit:

Dissolving and/or converting the polygons to singleparts does not seem to help much as the generated points still seem to fall within the minimum distance.

Tested on QGIS 2.18.3.

Answer 1

You have to change two things in order to get this to work. However, you won't get the maximum of wind turbines per area. For this, you would need to run some iterations for each value and get the max point count.

I've moved the index = QgsSpatialIndex() and points = dict() outside the for loop. The code would look like this:

import random
def generate_wind_turbines(spacing):
    layer = self.iface.activeLayer()
    crs = layer.crs()
    # Memory layer
    memory_lyr = QgsVectorLayer("Point?crs=epsg:" + unicode(crs.postgisSrid()) + "&index=yes", "Wind turbines for " + str(layer.name()), "memory")
    QgsMapLayerRegistry.instance().addMapLayer(memory_lyr)
    memory_lyr.startEditing()
    provider = memory_lyr.dataProvider()
    provider.addAttributes([QgsField("ID", QVariant.Int)])
    # Variables
    point_density = 0.0001
    fid = 1
    distance_area = QgsDistanceArea()
    # List of features
    fts = []
    # Create points
    points = dict() # changed from me 
    index = QgsSpatialIndex()# changend from me 
    nPoints = 0 # changed in the edit 
    pointCount = 0 # changed in the edit 

    for f in layer.getFeatures():
        fGeom = QgsGeometry(f.geometry())
        bbox = fGeom.boundingBox()
        # changed here as well 
        pointCount = int(round(point_density * distance_area.measure(fGeom))) + int(pointCount)
        fid += 1
        nIterations = 0
        maxIterations = pointCount * 200
        random.seed()
        while nIterations < maxIterations and nPoints < pointCount:
            rx = bbox.xMinimum() + bbox.width() * random.random()
            ry = bbox.yMinimum() + bbox.height() * random.random()
            pnt = QgsPoint(rx, ry)
            geom = QgsGeometry.fromPoint(pnt)
            if geom.within(fGeom) and checkMinDistance(pnt, index, spacing, points):
                f = QgsFeature(nPoints)
                f.setAttributes([fid])
                f.setGeometry(geom)
                fts.append(f)
                index.insertFeature(f)
                points[nPoints] = pnt
                nPoints += 1
            nIterations += 1
    provider.addFeatures(fts)
    memory_lyr.updateFields()
    memory_lyr.commitChanges()

def checkMinDistance( point, index, distance, points):
    if distance == 0:
        return True
    neighbors = index.nearestNeighbor(point, 1)
    if len(neighbors) == 0:
        return True
    if neighbors[0] in points:
        np = points[neighbors[0]]
        if np.sqrDist(point) < (distance * distance):
            return False
    return True

EDIT:

Joseph was right. my changes only worked for a realy small area. I've tested around and found a new solution by moving two variable out of the for loop and changed the pointCount variable.

I've tested it with 500m and this is the result ( two different tries):

Answer 2

One method could be to create another function which does the following:

1. Generate buffers with the same radius as the spacing around all points and stores these in a polygon layer.
2. Create a list containing the id of all buffers which intersect one another.
3. Deletes the buffers using this id list.
4. Generate centroid of the remaining buffers and stores these in a point layer.

This is the function I used:

def wind_turbine_spacing_checker(layer, spacing):
    # Create buffers for points
    poly_layer =  QgsVectorLayer("Polygon?crs=epsg:27700", 'Buffers' , "memory")
    pr = poly_layer.dataProvider() 
    pr.addAttributes([QgsField("ID", QVariant.Int)])
    feat_list = []
    for f in layer.getFeatures():
        poly = QgsFeature()
        f_buffer = f.geometry().buffer((spacing / 2), 99)
        f_poly = poly.setGeometry(QgsGeometry.fromPolygon(f_buffer.asPolygon()))
        poly.setAttributes([1])
        poly.setGeometry(f_buffer)
        feat_list.append(poly)

    pr.addFeatures(feat_list)
    poly_layer.updateExtents()
    poly_layer.updateFields()
    QgsMapLayerRegistry.instance().addMapLayers([poly_layer])

    # Get pairs of intersecting buffer features
    features = [feat for feat in poly_layer.getFeatures()]
    ids = []
    for feat in poly_layer.getFeatures():
        for geat in features:
            if feat.id() != geat.id():
                if geat.geometry().intersects(feat.geometry()):
                    ids.append([feat.id(), geat.id()])

    # Set/sort list and get id of intersecting feature(s)
    for x in ids:
        x.sort()

    ids_sort = set(tuple(x) for x in ids)
    ids_list = [list(x) for x in ids_sort]
    ids_firstItem = [item[0] for item in ids_list]
    final_list = list(set(ids_firstItem))

    # Use ids from final_list to remove intersecting buffer features
    with edit(poly_layer):
        poly_layer.deleteFeatures(final_list)

    # Create new point layer and get centroids from buffers
    # (using final_list to delete the original features may not delete those points where the buffer interesects
    # so best to obtain the centroid of the buffers and output this as a new file)
    result_layer = QgsVectorLayer('Point?crs=epsg:27700&field=id:string', 'Result' , 'memory')
    result_layer.startEditing()
    for feat in poly_layer.getFeatures():
        centroid = feat.geometry().centroid()
        name = feat.attribute("ID")
        centroid_feature = QgsFeature(poly_layer.fields())
        centroid_feature.setGeometry(centroid)
        centroid_feature['ID'] = name
        result_layer.addFeature(centroid_feature)

    result_layer.commitChanges()
    QgsMapLayerRegistry.instance().addMapLayer(result_layer)

The function can be executed immediately at the end of the generate_wind_turbines() function by using:

...
memory_lyr.commitChanges()
wind_turbine_spacing_checker(memory_lyr, spacing)

This gives the results as shown in the image in the question. Probably not the most efficient of solutions but it seems to work.

---

Some examples where the red points are those generated initially and the points displayed as turbines with a boundary are the final result:

• generate_wind_turbines(500)

  

---

• generate_wind_turbines(1000)