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

12

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

wprowadź opis zdjęcia tutaj

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”.

LeaningCactus
źródło
Byłoby to podobne, ale przeciwne (myślę) za pomocą maski środowiska w ArcGIS . To pozwoli ci ograniczyć utworzone wielokąty do określonej granicy. Jednak nie znam żadnego narzędzia, które wykorzystywałoby złożone granice / dziury (chociaż może w ArcGIS maska ​​może być tak złożona - nie przetestowałem jej i mógłbym spróbować później, jeśli będę miał czas).
Chris W
Testowałem teorię ArcGIS i to nie zadziała. Zgodnie z połączonym pytaniem możesz ograniczyć wyniki do zewnętrznego kształtu. Jednak powstałe w wyniku wycięcia dziury w kształcie są ignorowane. Ponadto, jeśli w tym otworze znajdują się również niektóre punkty, narzędzie jest błędne i nie działa. Nie potrafię wyjaśnić z różnicą twojego pierwszego problemu, ale drugi, w wyniku którego powstają taśmy, nie jest zupełnie nieoczekiwany - w końcu obszar ten byłby nadal przydzielony do tego samego punktu, nawet jeśli dziura jest obecna. Możesz użyć tej metody, a następnie włączyć taśmy do sąsiadów za pomocą metody czyszczenia.
Chris W
2
Możesz potencjalnie rozwiązać ten problem, przechodząc do rastra. W przypadku maski rastrowej odległość euklidesowa wychodziłaby z twoich punktów, dopóki nie uderzyłaby w komórki wychodzące z innego punktu lub w twoją maskę rastrową (twój opis uderzenia granicy). Następnie wykonujesz oczyszczanie stref i wektoryzujesz wynik, aby uzyskać wielokąty.
Chris W
1
Upewnij się, że geometria voronoi jest poprawna, uruchamiając v.clean, a następnie sprawdź geometrię. Na koniec wykonaj Różnicę, aby utworzyć dziury.
klewis
Co jest Voronoi w tych dziurach? Nie chcesz czysto dziurawić otworów? Dlaczego żadna warstwa wielokąta nie miałaby tego robić?
mdsumner,

Odpowiedzi:

3

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.

Guillaume
źródło
2

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.

wprowadź opis zdjęcia tutaj

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

wprowadź opis zdjęcia tutaj

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.

wprowadź opis zdjęcia tutaj

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

Steven Kay
źródło
1
Dzięki Steven, mam działającego rastra GRASS, ale miałem nadzieję na bardziej eleganckie rozwiązanie, jak wspomniano w opisie nagrody.
podmrok
0

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.

wprowadź opis zdjęcia tutaj

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:

wprowadź opis zdjęcia tutaj

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

wprowadź opis zdjęcia tutaj

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 .

mgri
źródło
Wykonuj dziury dłużej, np. Przecinając całą mapę, np. Rzekę, a zobaczysz problem. Dodanie taśm do sąsiadów tworzy wielokąty, które wyglądają zupełnie inaczej niż odpowiednio ograniczony diagram Voronoi. Próbowałem tego.
podmrok
Przepraszam, nie rozumiem: czy znalazłeś błąd w wynikach? Testowałem tylko kod dla przypadków, w których wielokąty były podobne do tych zaproponowanych w pytaniu.
mgri
Niestety, nie możesz teraz przetestować kodu, ale czy możesz pokazać wyniki uzyskane ze zmianą otworów opisaną w i.stack.imgur.com/Jpfra.png ?
podmrok
Jeśli rozszerzę ograniczenie do funkcji po prawej stronie, uzyskam to . Zamiast tego, jeśli bezpośrednio przesunę ograniczenie, uzyskam to .
mgri
Problemem jest mały trójkąt, na który wskazuje czerwona strzałka na moim rysunku. Nie powinno tam być, ale jest również w twoich wynikach. Wydaje się, że takie podejście rozwiązuje problem nr 1 pytania, ale pozostawia problem nr 2 nierozwiązany.
podmrok