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

8

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

Beispiel

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:

Beispiel Problem

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:

Ergebnis


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 .

Joseph
quelle
1
Wenn dies eine Option ist: Haben Sie versucht, ein mehrteiliges Polygon als Eingabepolygon zu verwenden?
LaughU
@LaughU - Ich habe es gerade mit einem mehrteiligen Polygon getestet, aber es werden keine Punkte generiert. Es wäre eine Option, die Einzelteil-Features in Mehrteiler umzuwandeln, wenn der Code so angepasst werden kann, dass er mit Mehrteiler-Polygonen funktioniert.
Joseph
Haben Sie daran gedacht, die Polygone aufzulösen, gefolgt von Multipart-to-Singlepart zu einer temporären Ebene, die Sie für die Punktgenerierung verwenden? Sie können auch einen negativen Puffer für die temporäre Ebene verwenden, um Symbolüberlappungen an den Rändern zu vermeiden.
Matte
@Matte - Danke, ich habe zuvor versucht, alle Polygone aufzulösen. Ich habe es erneut versucht und es in ein einzelnes Teil konvertiert (nicht sicher, ob dies etwas bewirkt, da es bereits ein einzelnes Feature war), aber einige Punkte in der Nähe der Kanten befinden sich innerhalb von Punkten in den anderen Polygonen. Ich möchte die Verwendung negativer Puffer vermeiden, da ich zulassen möchte, dass sich Punkte in der Nähe der Ränder befinden :)
Joseph

Antworten:

5

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

Geben Sie hier die Bildbeschreibung ein

LaughU
quelle
1
Schön, danke, dass du es bemerkt hast! Ich bin sicher, dies verbessert die Effizienz erheblich =)
Joseph
1
@ Joseph du hattest recht. Ich habe es mit kleinen Bereichen getestet, die funktionierten. Ich habe einige kleinere Änderungen hinzugefügt und es funktioniert jetzt für mich. Lassen Sie mich wissen, ob es noch verbessert werden muss :)
LaughU
1
Interessant! Ich werde dies morgen testen und zurückmelden, aber das Ergebnis sieht wirklich gut aus;)
Joseph
1
Ja, das ist gut! Ihr Skript ist auch schneller und liefert immer noch die Ergebnisse, die ich gesucht habe. Werde dir irgendwann das Kopfgeld geben, hoffe es macht dir nichts aus, aber ich habe deinen Code ganz leicht bearbeitet;)
Joseph
1
@ Joseph Nein, es macht mir nichts aus;) Ich habe mit dem Code getestet und vergessen, das Chaos, auch bekannt als Whitespaces, zu
beseitigen
4

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)

    Szenario 1


  • generate_wind_turbines(1000)

    Szenario 2

Joseph
quelle