Generieren Sie Punkte in mehreren Features unter Berücksichtigung des Mindestabstands

Ich habe eine Funktion, die Windturbinen erzeugt, die als Punkte dargestellt werden. Im Wesentlichen wird der Code aus dem Werkzeug Zufällige Punkte in Polygonen (fest) verwendet, allerdings mit einigen geringfügigen Änderungen.

Ziel ist es, zufällige Punkte innerhalb von Polygonen unter Berücksichtigung des angegebenen Mindestabstands zu erstellen. Dies funktioniert sehr gut, insbesondere bei Polygonen, die nicht nahe beieinander liegen (z. B. ein einzelnes Polygon):

Wenn sich das Polygon jedoch in der Nähe eines anderen Polygons befindet oder an dieses angrenzt (z. B. wie unten gezeigt), können sich die Punkte von jedem Polygon innerhalb des Mindestabstands befinden, wie in rot dargestellt:

Wie kann ich den Code so ändern, dass diese Punkte in Rot nicht in der Nähe eines anderen aus einem nahe gelegenen Polygon liegen?

Idealerweise möchte ich, dass mehrere Punkte durch einen einzelnen Punkt ersetzt werden:

Hier ist der Code, der in der Python-Konsole reproduziert werden kann. Vor dem Ausführen der Funktion muss eine Polygonebene mit einem relevanten CRS ausgewählt werden:

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)

Bearbeiten:

Das Auflösen und / oder Konvertieren der Polygone in einzelne Teile scheint nicht viel zu helfen, da die erzeugten Punkte immer noch innerhalb des Mindestabstands zu liegen scheinen.

Getestet mit QGIS 2.18.3 .

Sie müssen zwei Dinge ändern, damit dies funktioniert. Sie erhalten jedoch nicht das Maximum an Windkraftanlagen pro Gebiet. Dazu müssten Sie für jeden Wert einige Iterationen ausführen und die maximale Punktzahl ermitteln.

Ich habe die index = QgsSpatialIndex()und points = dict()außerhalb der for-Schleife verschoben . Der Code würde folgendermaßen aussehen:

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

BEARBEITEN:

Joseph hatte recht. Meine Änderungen funktionierten nur für einen wirklich kleinen Bereich. Ich habe herum getestet und eine neue Lösung gefunden, indem ich zwei Variablen aus der for-Schleife verschoben und die pointCountVariable geändert habe .

Ich habe es mit 500m getestet und dies ist das Ergebnis (zwei verschiedene Versuche):

Eine Methode könnte darin bestehen, eine andere Funktion zu erstellen, die Folgendes ausführt:

1. Generieren Sie Puffer mit demselben Radius wie der Abstand um alle Punkte und speichern Sie diese in einer Polygonebene.
2. Erstellen Sie eine Liste mit der ID aller Puffer, die sich überschneiden.
3. Löscht die Puffer mithilfe dieser ID-Liste.
4. Generieren Sie einen Schwerpunkt der verbleibenden Puffer und speichern Sie diese in einer Punktebene.

Dies ist die Funktion, die ich verwendet habe:

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)

Die Funktion kann sofort am Ende der generate_wind_turbines()Funktion ausgeführt werden, indem Folgendes verwendet wird:

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

Dies ergibt die Ergebnisse wie im Bild in der Frage gezeigt. Wahrscheinlich nicht die effizienteste Lösung, aber es scheint zu funktionieren.

Einige Beispiele, bei denen die roten Punkte die anfänglich erzeugten sind und die Punkte als Turbinen mit einer Grenze angezeigt werden, sind das Endergebnis:

• generate_wind_turbines(500)

  

• generate_wind_turbines(1000)