Berücksichtigen Sie Löcher / Einschränkungen bei der Erstellung von Voronoi-Polygonen in QGIS?

12

Ich versuche, Voronoi-Polygone in QGIS zu erstellen, die "Löcher" im allgemeinen Bereich berücksichtigen. Ein Beispiel wäre:

Geben Sie hier die Bildbeschreibung ein

Ich habe die Voronois in diesem Bild tatsächlich mit QGIS über den Befehl GRASS erstellt und dann mit dem Werkzeug "Differenz" die Löcher erstellt. Als "Differenz" -Ebene wurde ein separates Polygon-Shapefile verwendet, das die Ausdehnung der Löcher enthält. Eine Beispielanwendung wäre das Erstellen von Polygonen um Stichprobenpunkte, die zwischen Strukturen gesammelt wurden, die von der Analyse ausgeschlossen werden sollten.

Hier treten zwei Probleme auf:

  1. Die "Differenz" -Funktion scheint nicht 100% richtig zu funktionieren, da sich einige Polygongrenzen in die "Löcher" erstrecken. Dies kann behoben werden, indem in der Attributtabelle eine Zeile gefunden wird, die keine Polygon-ID-Nummer (oder ID von "0") hat.

  2. Diese Art des nachträglichen "Lochens" kann zu diskontinuierlichen Polygonen führen, wie der rote Pfeil im Bild zeigt.

Meine Frage ist: Gibt es ein Voronoi-Tool oder Plugin, das das Vorhandensein von "Löchern" in der Mitte der Domäne als einen einstufigen Prozess betrachtet und auch die Erzeugung diskontinuierlicher Polygone eliminiert? Ich stelle mir vor, dass ein solches Werkzeug eine Polygongrenze bis zum nächsten Schnittpunkt mit einer anderen Grenze erweitern würde, es sei denn, diese andere Grenze schlägt zuerst gegen eine "Loch" -Grenze.

LeaningCactus
quelle
Dies wäre ähnlich, aber das Gegenteil von (glaube ich) der Verwendung einer Umgebungsmaske in ArcGIS . Auf diese Weise können Sie die erstellten Polygone auf eine bestimmte Grenze beschränken. Mir ist jedoch kein Tool bekannt, das komplexe Grenzen / Löcher verwendet (obwohl die Maske in ArcGIS möglicherweise so komplex sein könnte - ich habe sie nicht getestet und kann sie später ausprobieren, wenn ich Zeit habe).
Chris W
Ich habe die ArcGIS-Theorie getestet und sie funktioniert nicht. Gemäß der verknüpften Frage können Sie die Ergebnisse auf eine äußere Form beschränken. Ein in die Form geschnittenes Loch wird jedoch von den resultierenden Polys ignoriert. Wenn dieses Loch auch einige Punkte enthält, ist das Werkzeug fehlerhaft und kann nicht ausgeführt werden. Ich kann Ihr erstes Problem nicht mit Unterschied erklären, aber das zweite, das zu Splittern führt, ist nicht völlig unerwartet - schließlich würde dieser Bereich immer noch demselben Punkt zugeordnet, selbst wenn das Loch vorhanden ist. Sie können diese Methode verwenden und dann die Splitter mit einer Bereinigungsmethode in ihre Nachbarn integrieren.
Chris W
2
Sie könnten dies möglicherweise lösen, indem Sie zum Raster gehen. Mit einer Rastermaske würde die euklidische Entfernung von Ihren Punkten ausgehen, bis sie entweder Zellen trifft, die von einem anderen Punkt kommen, oder Ihr Maskenraster (Ihre Beschreibung des Boundary Slam). Anschließend führen Sie eine zonale Bereinigung durch und vektorisieren das Ergebnis, um Polygone zu erhalten.
Chris W
1
Ich würde sicherstellen, dass die Voronoi-Geometrie gültig ist, indem ich v.clean ausführe und dann die Geometrie überprüfe. Führen Sie abschließend Difference aus, um die Löcher zu erstellen.
Klewis
Was ist Voronoi an diesen Löchern? Wollen Sie nicht sauber Löcher stanzen? Warum würde keine Polygonebene dies tun?
Mdsumner

Antworten:

3

Dies kann mit Rastern möglich sein. Konvertieren Sie zuerst Ihre Punkte und Grenzpolygone in ein hochauflösendes Raster. Legen Sie mit eine Maske für Ihre Grenzen fest r.mask. Führen Sie dann r.grow.distanceGRASS aus und verwenden Sie die Value= output. Dies gibt Ihnen für jedes Pixel den nächstgelegenen Punkt. Konvertieren Sie dies zurück in Vektorpolygone. Möglicherweise sind zusätzliche Schritte erforderlich, um Splitterpolygone zu entfernen.

Guillaume
quelle
2

Dies ist sicherlich mit Rastern möglich.

Dieser Screenshot zeigt das Problem hoffentlich deutlicher. Der Teil B des Voronoi befindet sich in Luftlinie näher am ursprünglichen Voronoi-Zentrum. Dies berücksichtigt jedoch nicht die Tatsache, dass das Umrunden des Gebäudes länger dauern würde. Mein Verständnis der Frage des OP ist, dass die Voronoi diese zusätzliche Entfernung berücksichtigen müssen, um um das Gebäude herumzugehen.

Geben Sie hier die Bildbeschreibung ein

Ich mag den Vorschlag von @Guillaume. Als ich es versuchte, hatte ich jedoch Probleme, r.grow.distancedie Maske zu ehren (siehe unten. Die Wellen sollten nicht durch die Gebäude gehen).

Mein Graswissen ist nicht so stark wie es sein könnte, also mache ich vielleicht etwas Dummes. Schauen Sie sich diesen Vorschlag auf jeden Fall zuerst an, da er viel weniger Arbeit kostet als meiner ;-)

Geben Sie hier die Bildbeschreibung ein

Schritt 1 - Erstellen Sie eine Kostenfläche

Der erste Schritt besteht darin, eine Kostenfläche zu erstellen. Dies muss nur einmal durchgeführt werden.

  • Erstellen Sie eine bearbeitbare Ebene, Löcher und alles.
  • Fügen Sie ein Feld mit dem Namen "Einheit" hinzu und setzen Sie es auf 1.
  • Verwenden Sie Polygon-zu-Raster auf Ihrer "ausgestanzten" Vektorebene (die Ebene mit den Löchern) und verwenden Sie das Feld "Einheit". Sie haben jetzt eine Ebenen- "Maske", in der 1 freier Speicherplatz ist und 0 erstellt wird.
  • Verwenden Sie den Rasterrechner, um daraus eine Kostenfläche zu machen. Ich werde "draußen" auf 1 und "drinnen" auf 9999 setzen. Dies wird das Bewegen durch Gebäude unerschwinglich erschweren.

    (("Maske @ 1" = 1) * 1) + (("Maske @ 1" = 0) * 9999)

Sie können mehr "organische" Ergebnisse erzielen, indem Sie der Kostenoberfläche ein wenig Rauschen hinzufügen (z. B. verwenden Sie eine Zufallszahl von 1 bis 3 anstelle von nur 1 für Pxiels im Freien.)

Schritt 2. Erstellen Sie kumulative Kostenraster für jedes Voronoi-Zentrum

Jetzt können wir (für jeweils eine Voronoi-Zelle) den GRASS-Algorithmus r.cost.coordinatesgegen unsere Kostenoberflächenschicht ausführen .

Verwenden Sie für die Startkoordinate das Vornoi-Zentrum. Wählen Sie für die Endkoordinate eine der Ecken Ihres Bereichs aus. Ich schlage vor, 'Knights Tour' zu verwenden, da dies zu glatteren Ergebnissen führt.

Das Ergebnis zeigt Linien gleicher Laufzeit von einem Voronoi-Zentrum. Beachten Sie, wie sich die Bänder um die Gebäude wickeln.

Geben Sie hier die Bildbeschreibung ein

Ich bin mir nicht sicher, wie ich das am besten automatisieren kann. Möglicherweise wird der Batch-Modus verarbeitet oder in Pyqgis ausgeführt.

Schritt 3. Führen Sie die Raster zusammen

Dies wird wahrscheinlich Code benötigen. Der Algorithmus wäre

create a raster 'A' to match the size of your cumulative cost images
fill raster 'A' with a suitably high number e.g. 9999
create an array of the same size as the raster.
for each cumulative cost raster number 1..N
    for each cell in image
        if cell < value in raster 'A'
            set value in raster 'A' to cell value
            set corresponding cell in array to cum. cost image number
write out array as a raster

Dieser Ansatz sollte ein Raster ergeben, bei dem jede Zelle unter Berücksichtigung von Hindernissen nach dem Voronoi-Zentrum kategorisiert wird, dem sie am nächsten liegt.

Sie können dann Raster-zu-Polygon verwenden. Sie können dann das Plugin " Verallgemeinern" verwenden , um die Artefakte des "Schritt" -Effekts aus dem Raster zu entfernen.

Entschuldigung für die Unbestimmtheit in den Schritten 2 und 3 ... Ich hoffe, jemand mischt sich mit einer eleganteren Lösung ein :)

Steven Kay
quelle
1
Vielen Dank, Steven, ich habe ein funktionierendes GRASS-Raster-Workaround, aber ich hatte auf eine elegantere Lösung gehofft, wie in der Kopfgeldbeschreibung erwähnt.
Underdark
0

Hinweis 1 : Ich konnte das vorgeschlagene Problem nicht reproduzieren, da das Differenz- Tool in mehreren von mir durchgeführten Tests gut für mich funktioniert hat (möglicherweise aufgrund der einfachen Geometrie des Problems oder weil das Tool seit der Frage verbessert wurde vor 1 Jahr gefragt).

Ich schlage jedoch eine Problemumgehung in PyQGIS vor, um die Verwendung des Differenz- Tools zu vermeiden . Alles basiert auf dem lokalen Schnittpunkt zwischen zwei Eingabeebenen (siehe Abbildung unten):

  1. eine Polygonvektorschicht, die die Voronoi-Polygone darstellt;
  2. eine Polygonvektorebene, die die Löcher / Einschränkungen darstellt, die von der Analyse ausgeschlossen werden müssen.

Geben Sie hier die Bildbeschreibung ein

Hinweis Nr. 2 : Da ich das Differenz- Tool nicht verwenden möchte , kann ich die Erstellung von "Splittern" nicht vermeiden (siehe dann). Daher musste ich das v.cleanTool ausführen , um sie zu entfernen. Darüber hinaus, wie @Chris W sagte,

[...] aber die zweite, die zu Splittern führt, ist nicht ganz unerwartet - schließlich würde dieser Bereich immer noch demselben Punkt zugeordnet, selbst wenn das Loch vorhanden ist. Sie können diese Methode verwenden und dann die Splitter mit einer Bereinigungsmethode in ihre Nachbarn integrieren .

Nach diesen notwendigen Voraussetzungen poste ich meinen Code:

##Voronoi_Polygons=vector polygon
##Constraints=vector polygon
##Voronoi_Cleaned=output vector

from qgis.core import *

voronoi = processing.getObject(Voronoi_Polygons)
crs = voronoi.crs().toWkt()
ex = voronoi.extent()
extent = '%f,%f,%f,%f' % (ex.xMinimum(), ex.xMaximum(), ex.yMinimum(), ex.yMaximum())

constraints = processing.getObject(Constraints)

# Create the output layer
voronoi_mod = QgsVectorLayer('Polygon?crs='+ crs, 'voronoi' , 'memory')
prov = voronoi_mod.dataProvider()
fields = voronoi.pendingFields() # Fields from the input layer
prov.addAttributes(fields) # Add input layer fields to the outLayer
voronoi_mod.updateFields()

# Spatial index containing all the 'constraints'
index_builds = QgsSpatialIndex()
for feat in constraints.getFeatures():
    index_builds.insertFeature(feat)

final_geoms = {}
final_attrs = {}

for feat in voronoi.getFeatures():
    input_geom = feat.geometry()
    input_attrs = feat.attributes()
    final_geom = []
    multi_geom = input_geom.asPolygon()
    input_geoms = [] # edges of the input geometry
    for k in multi_geom:
        input_geoms.extend(k)
    final_geom.append(input_geoms)
    idsList = index_builds.intersects(input_geom.boundingBox())
    mid_geom = [] # edges of the holes/constraints
    if len(idsList) > 0:
        req = QgsFeatureRequest().setFilterFids(idsList)
        for ft in constraints.getFeatures(req):
            geom = ft.geometry()
            hole = []
            res = geom.intersection(input_geom)
            res_geom = res.asPolygon()
            for i in res_geom:
                hole.extend(i)
                mid_geom.append(hole)
        final_geom.extend(mid_geom)
    final_geoms[feat.id()] = final_geom
    final_attrs[feat.id()] = input_attrs

# Add the features to the output layer
outGeom = QgsFeature()
for key, value in final_geoms.iteritems():
    outGeom.setGeometry(QgsGeometry.fromPolygon(value))
    outGeom.setAttributes(final_attrs[key])
    prov.addFeatures([outGeom])

# Add 'voronoi_mod' to the Layers panel
QgsMapLayerRegistry.instance().addMapLayer(voronoi_mod)

# Run 'v.clean'
processing.runalg("grass7:v.clean",voronoi_mod, 2, 0.1, extent, -1, 0.0001, Voronoi_Cleaned, None)

# Remove 'voronoi_mod' to the Layers panel
QgsMapLayerRegistry.instance().removeMapLayer(voronoi_mod)

was zu diesem Ergebnis führt:

Geben Sie hier die Bildbeschreibung ein

Nur aus Gründen der Übersichtlichkeit wäre dies das Ergebnis ohne die Verwendung des v.cleanWerkzeugs:

Geben Sie hier die Bildbeschreibung ein

Der Unterschied zum Ergebnis von @LeaningCactus besteht darin, dass die Geometrien inzwischen nicht gebrochen sind und fehlerfrei "gereinigt" werden können .

mgri
quelle
Machen Sie Löcher länger, z. B. schneiden Sie die gesamte Karte wie einen Fluss durch, und Sie werden das Problem sehen. Durch Hinzufügen von Splittern zu Nachbarn werden Polygone erstellt, die ganz anders aussehen als ein ordnungsgemäß eingeschränktes Voronoi-Diagramm. Das habe ich versucht.
Underdark
Entschuldigung, ich verstehe nicht: Haben Sie einen Fehler in den Ergebnissen gefunden? Ich habe den Code nur für die Fälle getestet, in denen die Polygone den in der Frage vorgeschlagenen ähnlich waren.
Mgri
Sie können den Code derzeit leider nicht testen, aber können Sie die Ergebnisse anzeigen, die mit der Änderung der in i.stack.imgur.com/Jpfra.png skizzierten Löcher erzielt wurden ?
Underdark
Wenn ich die Einschränkung bis zur Funktion auf der rechten Seite erweitere, erhalte ich diese . Wenn ich stattdessen die Einschränkung direkt verschiebe, erhalte ich diese .
Mgri
Das kleine Dreieck, auf das der rote Pfeil in meiner Zeichnung zeigt, ist das Problem. Es sollte nicht da sein, aber es ist auch in Ihren Ergebnissen. Dieser Ansatz scheint das erste Problem der Frage zu lösen, lässt jedoch das zweite Problem ungelöst.
Underdark