8

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):

Example

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:

Example problem

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:

Result


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.

  • 1
    if its an option: Have you tried to use a multipart polygon as the input polygon? – LaughU Nov 21 '17 at 15:23
  • @LaughU - I just tested it with a multipart polygon but no points are generated. It would be an option to convert the singlepart features into multiparts if the code can be tweaked to work with multipart polygons. – Joseph Nov 21 '17 at 15:34
  • Did you think of dissolving the polygons, followed by multipart-to-singlepart to a temporary layer that you use for the point generation? You can also use a negative buffer for the temporary layer to avoid symbol overlappings at the edges. – Matte Nov 27 '17 at 11:48
  • @Matte - Thanks, I tried dissolving all polygons before. I tried it again and converting it to singlepart (not sure if this does anything as it already was a single feature) but some points near the edges are within points in the other polygons. I would like to avoid using negative buffers as I want to permit points to be near the edges :) – Joseph Nov 27 '17 at 12:12
5
+200

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):

enter image description here

  • 1
    Nice, thanks for noticing! I'm sure this improves efficiency greatly =) – Joseph Nov 24 '17 at 13:24
  • 1
    @Joseph you were right. I've tested it with small areas which worked. I've added some minor changes and it works for me now. Let me know if it still need improvement :) – LaughU Nov 27 '17 at 14:43
  • 1
    Interesting! I will test this tomorrow and report back but the result looks really good ;) – Joseph Nov 27 '17 at 14:45
  • 1
    Yeah, this is good! Your script is also faster and still delivers the results I was looking for. Will award you the bounty eventually, hope you don't mind but I edited your code ever so slightly ;) – Joseph Nov 29 '17 at 10:54
  • 1
    @Joseph No I don't mind ;) I was testing around with the code and forgot to clean up the mess aka whitespaces – LaughU Nov 29 '17 at 12:01
4

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)

    Scenario 1


  • generate_wind_turbines(1000)

    Scenario 2

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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