Jak podłączyć sieć drogową do heksagonalnej siatki w QGIS?

Próbuję użyć QGIS 2.14 do przyciągnięcia sieci drogowej do sześciokątnej siatki, ale otrzymuję dziwne artefakty.

Utworzyłem siatkę heksadecymalną z MMQGIS , komórki mają około 20 x 23 m. Buforowałem sieć dróg o 1 m i zagęszcziłem ją, aby węzeł pojawiał się co kilka metrów. Poniżej możesz zobaczyć, co próbuję osiągnąć. Jak widać, w niektórych przypadkach mogę go uruchomić: -

• niebieski to zagęszczona droga (linia buforowana)
• czerwony to wersja „szesnastkowa” - to właśnie chcę znaleźć
• szary to siatka heksadecymalna

Następnie użyłem nowej funkcji Snap geometrii, aby przyciągnąć węzły do ​​najbliższego rogu sześciokąta. Wyniki są obiecujące, ale wydaje się, że istnieją przypadki, w których linia rozszerza się, wypełniając sześciokąt (lub jego część):

Powodem tego bufora jest to, że przyciąganie geometrii nie pozwala na przyciąganie do warstwy, której geometria jest inna. Na przykład nie można przyciągać węzłów na warstwie LINE do punktów na warstwie POINT). Najszczęśliwsze wydaje się przyciąganie POLYGONA do POLYGONA.

Podejrzewam, że drogi rozszerzają się, gdy jedna strona buforowanej linii drogi przeskakuje na jedną stronę komórki heksadecymalnej, a druga strona przeskakuje na drugą stronę komórki heksadecymalnej. W moim przykładzie najgorsze wydają się drogi, które przecinają zachód-wschód pod ostrym kątem.

Rzeczy, których próbowałem, bez powodzenia: -

• buforując sieć dróg o niewielką ilość, więc pozostaje wielokątem, ale jest bardzo cienki.
• zagęszczanie komórek szesnastkowych (więc na krawędziach znajdują się węzły, nie tylko w rogach)
• zmienianie maksymalnej odległości przyciągania (ma to największy efekt, ale nie wydaje mi się, aby znaleźć idealną wartość)
• używając warstw LINE, a nie POLYGONÓW

Uważam, że jeśli przejdę na używanie tylko warstw LINE, to działa przez chwilę, a następnie ulega awarii. Wydaje się, że zapisuje swoją pracę na bieżąco - niektóre linie zostały częściowo przetworzone.

Czy ktoś zna inny sposób przyciągania punktów na linii do najbliższego punktu na innej warstwie linii / wielokąta, idealnie bez potrzeby korzystania z postgres / postgis (chociaż rozwiązanie z postgis również byłoby mile widziane)?

EDYTOWAĆ

Dla każdego, kto chciałby spróbować, umieściłem tutaj na Dropbox startowy projekt QGIS . Obejmuje to warstwę siatki heksadecymalnej i warstwy linii zagęszczonych. (Sieć drogowa pochodzi z OSM, więc można ją pobrać za pomocą QuickOSM, np. Jeśli potrzebujesz uzyskać oryginał, aby nie zagęścić dróg).

Pamiętaj, że jest to OSGB (epsg: 27700), który jest zlokalizowanym UTM dla Wielkiej Brytanii, z jednostkami w metrach.

Moje rozwiązanie dotyczy skryptu PyQGIS, który jest szybszy i bardziej efektywny niż przepływ pracy obejmujący przyciąganie (próbowałem też). Korzystając z mojego algorytmu uzyskałem następujące wyniki:

Możesz uruchomić następujące fragmenty kodu w sekwencji z poziomu QGIS (w konsoli QGIS Python). Na koniec dostajesz warstwę pamięci z przyciąganymi trasami załadowanymi do QGIS.

Jedynym warunkiem jest utworzenie wieloczęściowego pliku Shapefile (użyj Processing->Singleparts to multipart, użyłem pola fictitiuosjako Unique ID fieldparametru). To da nam roads_multipart.shpplik z jedną funkcją.

Oto wyjaśniony algorytm:

1. Zdobądź najbliższe boki sześciokąta, w którym krzyżują się trasy Dla każdego sześciokąta tworzymy 6 trójkątów między każdą parą sąsiednich wierzchołków i odpowiadającą centroidem. Jeśli jakakolwiek droga przecina trójkąt, fragment dzielony przez sześciokąt i trójkąt jest dodawany do końcowej zatrzasniętej trasy. Jest to cięższa część całego algorytmu, zajmuje 35 sekund na mojej maszynie. W pierwszych dwóch wierszach znajdują się 2 ścieżki Shapefile, należy je dopasować, aby pasowały do ​​własnych ścieżek plików.

   hexgrid = QgsVectorLayer("/docs/borrar/hex_grid_question/layers/normal-hexgrid.shp", "hexgrid", "ogr")
   roads = QgsVectorLayer("/docs/borrar/hex_grid_question/layers/roads_multipart.shp", "roads", "ogr")  # Must be multipart!
   
   roadFeat = roads.getFeatures().next() # We just have 1 geometry
   road = roadFeat.geometry() 
   indicesHexSides = ((0,1), (1,2), (2,3), (3,4), (4,5), (5,0))
   
   epsilon = 0.01
   # Function to compare whether 2 segments are equal (even if inverted)
   def isSegmentAlreadySaved(v1, v2):
       for segment in listSegments:        
           p1 = QgsPoint(segment[0][0], segment[0][1])
           p2 = QgsPoint(segment[1][0], segment[1][1])
           if v1.compare(p1, epsilon) and v2.compare(p2, epsilon) \
               or v1.compare(p2, epsilon) and v2.compare(p1, epsilon):
               return True
       return False
   
   # Let's find the nearest sides of hexagons where routes cross
   listSegments = []
   for hexFeat in hexgrid.getFeatures():
       hex = hexFeat.geometry()
       if hex.intersects( road ):
           for side in indicesHexSides:
               triangle = QgsGeometry.fromPolyline([hex.centroid().asPoint(), hex.vertexAt(side[0]), hex.vertexAt(side[1])])
               if triangle.intersects( road ):
                   # Only append new lines, we don't want duplicates!!!
                   if not isSegmentAlreadySaved(hex.vertexAt(side[0]), hex.vertexAt(side[1])): 
                       listSegments.append( [[hex.vertexAt(side[0]).x(), hex.vertexAt(side[0]).y()], [hex.vertexAt(side[1]).x(),hex.vertexAt(side[1]).y()]] )  

2. Pozbądź się odłączonych (lub „otwartych”) segmentów, używając list, krotek i słowników w języku Python . W tym momencie pozostały niektóre odłączone segmenty, tj. Segmenty, które mają jeden odłączony wierzchołek, ale drugi jest połączony z co najmniej innymi 2 segmentami (patrz czerwone segmenty na następnym rysunku). Musimy się ich pozbyć.

   

   # Let's remove disconnected/open segments
   lstVertices = [tuple(point) for segment in listSegments for point in segment]
   dictConnectionsPerVertex = dict((tuple(x),lstVertices.count(x)-1) for x in set(lstVertices))
   
   # A vertex is not connected and the other one is connected to 2 segments
   def segmentIsOpen(segment):
       return dictConnectionsPerVertex[tuple(segment[0])] == 0 and dictConnectionsPerVertex[tuple(segment[1])] >= 2 \
           or dictConnectionsPerVertex[tuple(segment[1])] == 0 and dictConnectionsPerVertex[tuple(segment[0])] >= 2
   
   # Remove open segments
   segmentsToDelete = [segment for segment in listSegments if segmentIsOpen(segment)]        
   for toBeDeleted in segmentsToDelete:
       listSegments.remove( toBeDeleted )

3. Teraz możemy utworzyć warstwę wektorową z listy współrzędnych i załadować ją na mapę QGIS :

   # Create a memory layer and load it to QGIS map canvas
   vl = QgsVectorLayer("LineString", "Snapped Routes", "memory")
   pr = vl.dataProvider()
   features = []
   for segment in listSegments:
       fet = QgsFeature()
       fet.setGeometry( QgsGeometry.fromPolyline( [QgsPoint(segment[0][0], segment[0][1]), QgsPoint(segment[1][0], segment[1][1])] ) )
       features.append(fet)
   
   pr.addFeatures( features )
   vl.updateExtents()
   QgsMapLayerRegistry.instance().addMapLayer(vl)

Kolejna część wyniku:

Jeśli potrzebujesz atrybutów w przyciąganych trasach, moglibyśmy skorzystać z indeksu przestrzennego, aby szybko ocenić skrzyżowania (takie jak na /gis//a/130440/4972 ), ale to już inna historia.

Mam nadzieję że to pomoże!

Zrobiłem to w ArcGIS, na pewno można go zaimplementować za pomocą QGIS lub po prostu python z pakietem zdolnym do odczytu geometrii. Upewnij się, że drogi reprezentują sieć, tzn. Przecinają się tylko na końcach. Masz do czynienia z OSM, tak przypuszczam.

• Konwertuj wielokąty zbliżeniowe na linie i planuj je płasko, aby stały się również siecią geometryczną.
• Umieść punkty na ich końcach - Punkty Voronoi:
• Umieść punkty na drodze w regularnych odstępach 5 m, upewnij się, że drogi sieciowe mają dobrą unikalną nazwę:

• Dla każdego punktu drogowego znajdź współrzędne najbliższego punktu Voronoi:
• Utwórz „Drogi”, łącząc najbliższe punkty w tej samej kolejności:

Jeśli nie chcesz tego widzieć:

Nie próbuj używać punktów pikowania na liniach Voronoi. Obawiam się, że tylko pogorszy to sytuację. Zatem twoją jedyną opcją jest tworzenie sieci z linii Voronoi i znajdowanie tras między punktami końcowymi dróg, to też nie jest wielka sprawa

Zdaję sobie sprawę, że prosisz o metodę QGIS, ale proszę o zwięzłą odpowiedź:

roads = 'clipped roads' # roads layer
hexgrid = 'normal-hexgrid' # hex grid layer
sr = arcpy.Describe('roads').spatialReference # spatial reference
outlines = [] # final output lines
points = [] # participating grid vertices
vert_dict = {} # vertex dictionary
hex_dict = {} # grid dictionary
with arcpy.da.SearchCursor(roads,["SHAPE@","OID@"], spatial_reference=sr) as r_cursor: # loop through roads
    for r_row in r_cursor:
        with arcpy.da.SearchCursor(hexgrid,["SHAPE@","OID@"], spatial_reference=sr) as h_cursor: # loop through hex grid
            for h_row in h_cursor:
                if not r_row[0].disjoint(h_row[0]): # check if the shapes overlap
                    hex_verts = []
                    for part in h_row[0]:
                        for pnt in part:
                            hex_verts.append(pnt) # add grid vertices to list
                    int_pts = r_row[0].intersect(h_row[0],1) # find all intersection points between road and grid
                    hex_bnd = h_row[0].boundary() # convert grid to line
                    hex_dict[h_row[1]] = hex_bnd # add grid geometry to dictionary
                    for int_pt in int_pts: # loop through intersection points
                        near_dist = 1000 # arbitrary large number
                        int_pt = arcpy.PointGeometry(int_pt,sr)
                        for hex_vert in hex_verts: # loop through hex vertices
                            if int_pt.distanceTo(hex_vert) < near_dist: # find shortest distance between intersection point and grid vertex
                                near_vert = hex_vert # remember geometry
                                near_dist = int_pt.distanceTo(hex_vert) # remember distance
                        vert_dict.setdefault(h_row[1],[]).append(arcpy.PointGeometry(near_vert,sr)) # store geometry in dictionary
                        points.append(arcpy.PointGeometry(near_vert,sr)) # add to points list
for k,v in vert_dict.iteritems(): # loop through participating vertices
    if len(v) < 2: # skip if there was only one vertex
        continue
    hex = hex_dict[k] # get hex grid geometry
    best_path = hex # longest line possible is hex grid boundary
    for part in hex:
        for int_vert in v: # loop through participating vertices
            for i,pnt in enumerate(part): # loop through hex grid vertices
                if pnt.equals(int_vert): # find vertex index on hex grid corresponding to current point
                    start_i = i
                    if start_i == 6:
                        start_i = 0
                    for dir in [[0,6,1],[5,-1,-1]]: # going to loop once clockwise, once counter-clockwise
                        past_pts = 0 # keep track of number of passed participating vertices
                        cur_line_arr = arcpy.Array() # polyline coordinate holder
                        cur_line_arr.add(part[start_i]) # add starting vertex to growing polyline
                        for j in range(dir[0],dir[1],dir[2]): # loop through hex grid vertices
                            if past_pts < len(v): # only make polyline until all participating vertices have been visited
                                if dir[2] == 1: # hex grid vertex index bookkeeping
                                    if start_i + j < 6:
                                        index = start_i + j
                                    else:
                                        index = (start_i - 6) + j
                                else:
                                    index = j - (5 - start_i)
                                    if index < 0:
                                        index += 6
                                cur_line_arr.add(part[index]) # add current vertex to growing polyline
                                for cur_pnt in v:
                                    if part[index].equals(cur_pnt): # check if the current vertex is a participating vertex
                                        past_pts += 1 # add to counter
                        if cur_line_arr.count > 1:
                            cur_line = arcpy.Polyline(cur_line_arr,sr)
                            if cur_line.length < best_path.length: # see if current polyline is shorter than any previous candidate
                                best_path = cur_line # if so, store polyline
    outlines.append(best_path) # add best polyline to list
arcpy.CopyFeatures_management(outlines, r'in_memory\outlines') # write list
arcpy.CopyFeatures_management(points, r'in_memory\mypoints') # write points, if you want

Uwagi:

• Ten skrypt zawiera wiele pętli w pętlach i zagnieżdżony kursor. Zdecydowanie jest miejsce na optymalizację. Przejrzałem twoje zestawy danych w ciągu kilku minut, ale więcej funkcji skomplikuje problem.

Jeśli podzielisz linię drogi na segmenty, w których każdy segment był całkowicie zawarty w sześciokącie, twoją decyzją, które segmenty linii sześciokąta zastosować, będzie to, czy odległość od środka ciężkości podzielonego odcinka drogi do punktu środkowego każdej strony sześciokąta będzie mniejsza niż połowa średnica sześciokąta (lub mniejsza niż promień koła pasującego do sześciokąta).

Tak więc, jeśli (jeden segment na raz) wybierzesz segmenty linii sześciokąta (gdzie każdy segment jest bokiem sześciokąta), które znajdują się w odległości od promienia sześciokąta, możesz skopiować te geometrie linii i scalić je niezależnie od unikalnego identyfikatora używanego w zbiorze danych o drogach.

Jeśli masz problemy z scaleniem unikalnego identyfikatora, możesz zastosować bufor i wybrać według lokalizacji tylko w tych segmentach, aby zastosować atrybuty zestawu danych o drodze; w ten sposób nie będziesz musiał się martwić o tworzenie fałszywych dopasowań przy zbyt dużym buforze.

Problem z narzędziem przyciągania polega na tym, że przyciąga punkty bez rozróżnienia; trudno jest znaleźć idealną tolerancję w użyciu. Dzięki tej metodologii poprawnie identyfikujesz, które segmenty linii sześciokąta mają być używane, a następnie zastępujesz geometrię danych o drodze (lub wstawiasz geometrie do innego zestawu danych).

Ponadto, jeśli nadal masz problem z segmentami linii, które przeskakują z jednej strony sześciokąta na drugą, możesz podzielić linię na segmenty według wierzchołków, obliczyć długość każdej linii, a następnie usunąć wszelkie segmenty linii, które są większe niż średnia długość jednego boku sześciokąta.

Snapper geometrii w qgis 3.0 został przerobiony i umożliwia teraz przyciąganie między różnymi typami geometrii. Ma również wiele poprawek. Możesz wypróbować wersję „codziennej migawki”, aby uzyskać dostęp do ulepszonej wersji migawki przed oficjalnym wydaniem wersji 3.0.