Wie kann ein Straßennetz in QGIS an einem hexagonalen Raster ausgerichtet werden?

Ich versuche, mit QGIS 2.14 ein Straßennetz an einem hexagonalen Raster auszurichten, aber es treten merkwürdige Artefakte auf.

Ich habe mit MMQGIS ein Hex-Gitter erstellt , die Zellen sind ca. 20 x 23 m groß . Ich habe das Straßennetz um 1 m gepuffert und verdichtet , sodass alle paar Meter ein Knoten vorhanden ist. Sie können unten sehen, was ich zu erreichen versuche. Wie Sie sehen, kann ich es in einigen Fällen zum Laufen bringen:

• blau ist die verdichtete Straße (eine gepufferte Linie)
• rot ist die "hexifizierte" Version - das ist, was ich finden möchte
• Das Grau ist das Hex-Gitter

Anschließend habe ich die neue Funktion " Geometrien fangen" verwendet, um die Knoten an der nächstgelegenen Sechseckecke zu fangen. Die Ergebnisse sind vielversprechend, aber es scheint einige Randfälle zu geben, bei denen sich die Linie ausdehnt, um das Sechseck (oder einen Teil davon) auszufüllen:

Der Grund für den Puffer ist , dass Snap - Geometrien nicht , dass Sie auf eine Schicht schnappen sich lassen , deren Geometrie ist anders. Sie können beispielsweise keine Knoten auf einer LINE-Ebene an Punkten auf einer POINT-Ebene ausrichten. Es scheint am glücklichsten zu sein, POLYGON an POLYGON zu knipsen.

Ich vermute, die Straßen dehnen sich aus, wenn eine Seite der gepufferten Straßenlinie zur einen Seite der Hex-Zelle und die andere Seite zur anderen Seite der Hex-Zelle springt. In meinem Beispiel scheinen die Straßen, die sich in einem spitzen Winkel nach West-Ost kreuzen, die schlechtesten zu sein.

Dinge, die ich erfolglos versucht habe:

• Puffern des Straßennetzes um einen winzigen Betrag, so dass es ein Polygon bleibt, aber sehr dünn ist.
• Verdichtung der hexadezimalen Zellen (also gibt es Knoten entlang der Kanten, nicht nur an den Ecken)
• Variieren der maximalen Fangentfernung (dies hat den größten Effekt, aber ich kann scheinbar keinen idealen Wert finden)
• unter Verwendung von LINE-Ebenen, nicht von POLYGONen

Ich finde, wenn ich nur LINE-Ebenen verwende, funktioniert es eine Weile und stürzt dann ab. Es scheint seine Arbeit zu retten, da einige Zeilen teilweise verarbeitet wurden.

Kennt jemand eine andere Möglichkeit, um Punkte auf einer Linie zum nächsten Punkt auf einer anderen Linie / Polygonebene zu fangen, idealerweise ohne die Verwendung von Postgres / Postgis (obwohl eine Lösung mit Postgis auch wünschenswert wäre)?

BEARBEITEN

Für alle, die es ausprobieren möchten, habe ich hier auf Dropbox ein QGIS- Startprojekt eingestellt . Dies schließt die Ebenen Hex Grid und Densified lines ein. (Das Straßennetz stammt von OSM und kann mit QuickOSM heruntergeladen werden, z. B. wenn Sie das Original benötigen, um die Straßen zu verdichten.)

Beachten Sie, dass es sich bei OSGB (epsg: 27700) um eine lokalisierte UTM für Großbritannien mit Einheiten in Metern handelt.

Meine Lösung beinhaltet ein PyQGIS-Skript, das schneller und effektiver ist als ein Workflow, bei dem ein Snap ausgeführt wird (ich habe es auch ausprobiert). Mit meinem Algorithmus habe ich folgende Ergebnisse erhalten:

Sie können die folgenden Code-Snippets nacheinander in QGIS (in der QGIS Python-Konsole) ausführen. Am Ende erhalten Sie eine Speicherebene mit den in QGIS geladenen Routen.

Die einzige Voraussetzung ist die Erstellung eines mehrteiligen Straßen-Shapefiles (verwenden Sie Processing->Singleparts to multipart, ich habe das Feld fictitiuosals Unique ID fieldParameter verwendet). Dadurch erhalten wir eine roads_multipart.shpDatei mit einer einzigen Funktion.

Hier ist der Algorithmus erklärt:

1. Holen Sie sich die nächsten Sechseckseiten, wo sich die Routen kreuzen. Für jedes Sechseck erstellen wir 6 Dreiecke zwischen jedem Paar von Nachbarscheitelpunkten und dem entsprechenden Schwerpunkt. Wenn eine Straße ein Dreieck schneidet, wird das Segment, das das Sechseck und das Dreieck gemeinsam haben, zur endgültigen Route hinzugefügt. Dies ist der schwerere Teil des gesamten Algorithmus. Auf meinem Computer dauert es 35 Sekunden. In den ersten beiden Zeilen befinden sich 2 Shapefile-Pfade, die Sie an Ihre eigenen Dateipfade anpassen sollten.

   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. Befreien Sie sich mit Python-Listen, -Tupeln und -Wörterbüchern von nicht verbundenen (oder "offenen") Segmenten . Zu diesem Zeitpunkt sind noch einige getrennte Segmente vorhanden, dh Segmente, bei denen ein Scheitelpunkt getrennt, der andere jedoch mit mindestens zwei anderen Segmenten verbunden ist (siehe rote Segmente in der nächsten Abbildung). Wir müssen sie loswerden.

   

   # 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. Jetzt können wir einen Vektorlayer aus der Liste der Koordinaten erstellen und in die QGIS-Karte laden :

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

Ein weiterer Teil des Ergebnisses:

Sollten Sie Attribute für die gefangenen Routen benötigen, können Sie mit einem räumlichen Index schnell Kreuzungen auswerten (z. B. in /gis//a/130440/4972 ), aber das ist eine andere Geschichte.

Hoffe das hilft!

Ich habe es in ArcGIS gemacht, es kann sicherlich mit QGIS oder einfach mit Python mit einem Paket implementiert werden, das Geometrien lesen kann. Stellen Sie sicher, dass Straßen ein Netzwerk darstellen, dh sich nur an den Enden kreuzen. Sie beschäftigen sich mit OSM, ich nehme an, dass es der Fall ist.

• Wandeln Sie Proximity-Polygone in Linien um und planarisieren Sie sie, sodass sie auch zu einem geometrischen Netzwerk werden.
• Platziere Punkte an ihren Enden - Voronoi Punkte:
• Platzieren Sie die Punkte in regelmäßigen Abständen von 5 m auf der Straße und vergewissern Sie sich, dass die Netzwerkstraßen einen eindeutigen Namen haben:

• Finden Sie für jeden Straßenpunkt die Koordinaten des nächsten Voronoi-Punkts:
• Erstellen Sie „Straßen“, indem Sie die nächstgelegenen Punkte in derselben Reihenfolge verbinden:

Wenn Sie das nicht sehen wollen:

Versuchen Sie nicht, auf Voronoi-Linien Verkettungspunkte zu verwenden. Ich fürchte, es wird es nur noch schlimmer machen. Sie können also nur ein Netzwerk aus Voronoi-Linien erstellen und Routen zwischen Straßenendpunkten finden, auch das ist keine große Sache

Mir ist klar, dass Sie nach einer QGIS-Methode fragen, aber ich bitte Sie um eine bogenförmige Antwort:

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

Anmerkungen:

• Dieses Skript enthält viele Schleifen in Schleifen und einen verschachtelten Cursor. Es gibt definitiv Raum für Optimierungen. Ich habe Ihre Datensätze in ein paar Minuten durchgearbeitet, aber weitere Funktionen werden das Problem verschlimmern.

Wenn Sie die Straßenlinie in Segmente aufteilen würden, in denen jedes Segment vollständig vom Sechseck umschlossen ist, wäre Ihre Entscheidung, welche Sechseckliniensegmente verwendet werden sollen, ob der Abstand zwischen dem Schwerpunkt des aufgeteilten Straßensegments und dem Mittelpunkt der einzelnen Sechseckseiten weniger als die Hälfte beträgt Durchmesser des Sechsecks (oder kleiner als der Radius eines Kreises, der in das Sechseck passt).

Wenn Sie also (jeweils ein Segment) Sechskantliniensegmente auswählen (wobei jedes Segment eine Seite des Sechskants ist), die sich in einem Abstand zum Radius des Sechskants befinden, können Sie diese Liniengeometrien kopieren und zusammenführen Welche eindeutige Kennung Sie auch für Ihren Straßendatensatz verwenden.

Wenn beim Zusammenführen des eindeutigen Bezeichners Probleme auftreten, können Sie den Puffer anwenden und nach Position nur in diesen Segmenten auswählen, um die Attribute Ihres Straßendatensatzes anzuwenden. Auf diese Weise müssen Sie sich keine Gedanken darüber machen, ob mit einem zu großen Puffer falsche Übereinstimmungen hergestellt werden.

Das Problem mit dem Ausrichtungswerkzeug ist, dass es wahllos Punkte ausrichtet. Es ist schwierig, diese perfekte Toleranz zu finden. Mit dieser Methode können Sie die zu verwendenden Sechskantliniensegmente korrekt identifizieren und dann die Geometrie Ihrer Straßendaten ersetzen (oder die Geometrien in einen anderen Datensatz einfügen).

Wenn Sie weiterhin Probleme mit den Liniensegmenten haben, die von einer Seite des Sechsecks zur anderen springen, können Sie die Linie durch Scheitelpunkte in Segmente aufteilen, die Länge jeder Linie berechnen und dann alle Liniensegmente entfernen, die größer als sind die durchschnittliche Länge einer Seite des Sechsecks.

Der Geometrie-Snapper in qgis 3.0 wurde überarbeitet und ermöglicht nun das Einrasten zwischen verschiedenen Geometrietypen. Es hat auch viele Korrekturen. Sie können eine "Daily Snapshot" -Version ausprobieren, um Zugriff auf den verbesserten Snapper zu erhalten, bevor 3.0 offiziell veröffentlicht wird.