Biorąc pod uwagę dziury / ograniczenia w tworzeniu wielokątów Voronoi w QGIS?

Usiłuję utworzyć wielokąty voronoi w QGIS, które uwzględniają „dziury” w domenie ogólnej. Przykładem może być:

W rzeczywistości stworzyłem Voronoi na tym obrazie za pomocą QGIS za pomocą polecenia GRASS, a następnie za pomocą narzędzia „Różnica”, aby utworzyć dziury. Oddzielny wielokątny plik kształtu, który zawiera zakres otworów, został użyty jako warstwa „Różnica”. Przykładową aplikacją byłoby tworzenie wielokątów wokół punktów próbkowania, które zostały zebrane między strukturami, które należy wykluczyć z analizy.

Pojawiają się tutaj dwa problemy:

1. Funkcja „różnicy” nie działa w 100% poprawnie, a niektóre granice wielokątów rozciągają się na „dziury”. Można to naprawić, znajdując wiersz w tabeli atrybutów, który nie ma numeru identyfikacyjnego wielokąta (lub identyfikatora „0”).

2. Ten rodzaj „dziurkowania” po fakcie może skutkować nieciągłymi wielokątami, jak pokazuje czerwona strzałka na obrazie.

Moje pytanie brzmi: czy istnieje narzędzie lub wtyczka Voronoi, które mogą uznać obecność „dziur” w środku domeny za proces jednoetapowy, a także wyeliminować generowanie nieciągłych wielokątów? Wyobrażam sobie, że takie narzędzie przedłużyłoby granicę wielokąta do najbliższego przecięcia z inną granicą, chyba że ta inna granica najpierw uderzy o granicę „dziury”.

Może to być możliwe przy użyciu rastrów. Najpierw przekonwertuj punkty i wielokąty graniczne na raster wysokiej rozdzielczości. Ustaw maskę dla swoich granic za pomocą r.mask. Następnie uruchom r.grow.distanceGRASS i użyj Value= output. To da ci za każdy piksel, który jest najbliższym punktem. Przekształć to z powrotem w wielokąty wektorowe. Konieczne mogą być dodatkowe kroki, aby pozbyć się wielokątów taśmy.

Jest to z pewnością możliwe w przypadku rastrów.

Mam nadzieję, że ten zrzut ekranu lepiej pokazuje problem. Część B voronoi jest bliżej „w linii prostej” od oryginalnego centrum voronoi, ale nie bierze to pod uwagę faktu, że obejście budynku zajęłoby więcej czasu. Rozumiem pytanie PO, że voronoi muszą wziąć pod uwagę tę dodatkową odległość, aby obejść budynek.

Podoba mi się sugestia @ Guillaume. Jednak gdy spróbowałem, miałem problemy r.grow.distancez uhonorowaniem maski (patrz poniżej. Fale nie powinny przechodzić przez budynki).

Moja wiedza na temat Trawy nie jest tak silna, jak mogłaby być, więc może robię coś głupiego. Zdecydowanie sprawdź najpierw tę sugestię, ponieważ będzie to o wiele mniej pracy niż moja ;-)

Krok 1 - Utwórz powierzchnię kosztów

Pierwszym krokiem jest stworzenie powierzchni kosztowej. Tę czynność należy wykonać tylko raz.

• utwórz edytowalną warstwę, dziury i wszystko.
• dodaj pole o nazwie „jednostka”, ustaw je na 1.
• używając wielokąta do rastra na „wykrawanej” warstwie wektorowej (tej, która ma dziury), używając pola „jednostki”. Masz teraz „maskę” warstwy, gdzie 1 oznacza wolne miejsce, a 0 buduje.

• użyj kalkulatora rastrowego, aby przekształcić go w powierzchnię kosztów. Ustawię „na zewnątrz” na 1, a „wewnątrz” na 9999. To sprawi, że poruszanie się po budynkach będzie wyjątkowo trudne.

  ((„mask @ 1” = 1) * 1) + ((„mask @ 1” = 0) * 9999)

Możesz uzyskać więcej „organicznych” wyników, dodając trochę szumu do powierzchni kosztów (np. Użyj losowej liczby od 1 do 3, a nie tylko 1 dla pikseli zewnętrznych).

Krok 2. Utwórz skumulowane rastry kosztów dla każdego centrum Voronoi

Teraz możemy uruchomić (dla jednej komórki voronoi na raz) algorytm GRASS r.cost.coordinateswzględem naszej warstwy powierzchni kosztowej.

Jako współrzędną początkową użyj centrum Vornoi. Jako współrzędną końcową wybierz jeden z rogów swojego obszaru. Sugeruję użycie „Knights Tour”, ponieważ daje to płynniejsze wyniki.

Wynik pokazuje linie równego czasu podróży z jednego centrum voronoi. Zwróć uwagę, jak opaski owijają się wokół budynków.

Nie wiem, jak najlepiej to zautomatyzować. Może przetwarzanie w trybie wsadowym lub w pyqgis.

Krok 3. Scal rastry

Prawdopodobnie będzie to wymagać kodu. Algorytm byłby

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

Takie podejście powinno dać raster, w którym każda komórka jest podzielona na kategorie według najbliższego centrum voronoi, biorąc pod uwagę przeszkody.

Następnie możesz użyć rastra do wielokąta. Następnie możesz użyć wtyczki Uogólnij, aby usunąć artefakty efektu „kroku” z rastra.

Przepraszamy za niejasność za kroki 2 i 3 ... Mam nadzieję, że ktoś wpadnie na bardziej eleganckie rozwiązanie :)

Uwaga 1 : Nie byłem w stanie odtworzyć proponowanego problemu, ponieważ narzędzie Różnica działało dla mnie dobrze w kilku testach, które przeprowadziłem (być może było to spowodowane prostą geometrią problemu lub dlatego, że narzędzie zostało ulepszone od czasu pytania zapytał 1 rok temu).

Proponuję jednak obejście w PyQGIS, aby uniknąć użycia narzędzia Różnica . Wszystko opiera się na lokalnym przecięciu dwóch warstw wejściowych (patrz rysunek poniżej):

1. warstwa wektora wielokąta reprezentująca wielokąty Voronoi;
2. warstwa wektora wielokąta reprezentująca otwory / ograniczenia, które należy wykluczyć z analizy.

Uwaga 2 : Ponieważ nie chcę używać narzędzia Różnica , nie jestem w stanie uniknąć tworzenia „taśm” (patrz wtedy), więc musiałem uruchomić to v.cleannarzędzie, aby je wyeliminować. Ponadto, jak powiedział @Chris W,

[...] ale drugi wynik w postaci taśm nie jest całkowicie nieoczekiwany - w końcu obszar ten nadal byłby przydzielony do tego samego punktu, nawet gdyby dziura była obecna. Możesz użyć tej metody, a następnie włączyć taśmy do sąsiadów za pomocą metody czyszczenia .

Po tych niezbędnych przesłaniach umieszczam kod:

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

co prowadzi do tego wyniku:

Dla jasności byłby to wynik bez użycia v.cleannarzędzia:

Różnica w stosunku do wyniku @LeaningCactus polega na tym, że do tej pory geometrie nie są zepsute i można je „wyczyścić” bez błędów .