Linien schneiden, um mit Python mit QGIS Kreuzungen zu erhalten?

10

Ich habe eine Reihe von Linien, die Buslinien darstellen. Einige der Linien überlappen sich und nehmen die gleichen Straßen.

Geben Sie hier die Bildbeschreibung ein

Ich kann die Knoten extrahieren. Geben Sie hier die Bildbeschreibung ein

Ich bin jedoch daran interessiert, nur Kreuzungen wie diese zu extrahieren: Geben Sie hier die Bildbeschreibung ein

Wie kann ich das machen? Ich suche nach Wegen mit QGIS oder Python.

Ich habe die Schnittmethode von GDAL Python ausprobiert, aber dies gibt mir im Grunde nur die Eckpunkte zurück.

Die Line Intersections-Methode von QGIS gibt mir die Kreuzungen zurück, wenn sich zwei Linien kreuzen. Für den Fall, dass zwei Buslinien einen weiten Teil ihrer Route auf derselben Straße fahren, gibt es mir nicht den Punkt, an dem sie zusammenführen.

ustroetz
quelle
Haben Sie das Linienkreuzungswerkzeug in QGIS ausprobiert: Vektoranalyse-Werkzeug> Linienschnittpunkte ... Es gibt Ihnen nicht die End- und Startknoten einer Linie, sondern alle Schnittpunkte.
Jakob
Ja, ich habe darüber in der Frage geschrieben.
Ustroetz
Ich bin mir nicht sicher, was Sie fragen, zum Teil, weil alle Linien in Ihren Bildern auf die gleiche Weise symbolisiert sind. Ich kann nicht verschiedene Routen voneinander unterscheiden, um zu verstehen, welche Knoten Sie betrachten oder warum es so viele in der gibt zweites Bild. Stimmen die Routen auf Straßen überein? Sind sie alle Zweipunkt-Liniensegmente oder durchgehende Polylinien? Ich stelle fest, dass das von Ihnen beschriebene Verhalten das gleiche ist wie das Schnittpunkt- Werkzeug von ArcGIS. Linien / Linien mit Linienausgabe ergeben Überlappungen, Linien / Linien mit Punktausgabe ergeben jedoch nur Kreuzungen.
Chris W
Auf dieser Grundlage müssen Sie möglicherweise beide Methoden anwenden, um das zu erreichen, was Sie meiner Meinung nach wollen. Holen Sie sich die Kreuzungen (Linie / Linie = Punkt) und dann die Überlappungen (Linie / Linie = Linie) und extrahieren Sie die Start- / Endknoten für diese Überlappungslinien. Dies sollten alle Punkte / Knoten sein, nach denen Sie suchen.
Chris W

Antworten:

19

Die Knoten:

Sie möchten zwei Dinge, die Endpunkte der Polylinien (ohne Zwischenknoten) und die Schnittpunkte. Es gibt ein zusätzliches Problem, einige Polylinienendpunkte sind auch Schnittpunkte:

Geben Sie hier die Bildbeschreibung ein

Eine Lösung besteht darin, Python und die Module Shapely und Fiona zu verwenden

1) Lesen Sie das Shapefile:

from shapely.geometry import Point, shape
import fiona
lines = [shape(line['geometry']) for line in fiona.open("your_shapefile.shp")]

2) Finden Sie die Endpunkte der Linien ( wie würde man die Endpunkte einer Polylinie erhalten? ):

endpts = [(Point(list(line.coords)[0]), Point(list(line.coords)[-1])) for line  in lines]
# flatten the resulting list to a simple list of points
endpts= [pt for sublist in endpts  for pt in sublist] 

Geben Sie hier die Bildbeschreibung ein

3) Berechnen Sie die Schnittpunkte (Iteration durch Geometriepaare in der Ebene mit dem itertools- Modul). Das Ergebnis einiger Kreuzungen sind MultiPoints und wir möchten eine Liste von Punkten:

import itertools
inters = []
for line1,line2 in  itertools.combinations(lines, 2):
  if  line1.intersects(line2):
    inter = line1.intersection(line2)
    if "Point" == inter.type:
        inters.append(inter)
    elif "MultiPoint" == inter.type:
        inters.extend([pt for pt in inter])
    elif "MultiLineString" == inter.type:
        multiLine = [line for line in inter]
        first_coords = multiLine[0].coords[0]
        last_coords = multiLine[len(multiLine)-1].coords[1]
        inters.append(Point(first_coords[0], first_coords[1]))
        inters.append(Point(last_coords[0], last_coords[1]))
    elif "GeometryCollection" == inter.type:
        for geom in inter:
            if "Point" == geom.type:
                inters.append(geom)
            elif "MultiPoint" == geom.type:
                inters.extend([pt for pt in geom])
            elif "MultiLineString" == geom.type:
                multiLine = [line for line in geom]
                first_coords = multiLine[0].coords[0]
                last_coords = multiLine[len(multiLine)-1].coords[1]
                inters.append(Point(first_coords[0], first_coords[1]))
                inters.append(Point(last_coords[0], last_coords[1]))

Geben Sie hier die Bildbeschreibung ein

4) Beseitigen Sie Duplikate zwischen Endpunkten und Schnittpunkten (wie Sie in den Abbildungen sehen können).

result = endpts.extend([pt for pt in inters  if pt not in endpts])

5) Speichern Sie das resultierende Shapefile

from shapely.geometry import mapping
# schema of the shapefile
schema = {'geometry': 'Point','properties': {'test': 'int'}}
# creation of the shapefile
with fiona.open('result.shp','w','ESRI Shapefile', schema) as output:
    for i, pt in enumerate(result):
        output.write({'geometry':mapping(pt), 'properties':{'test':i}})

Endergebnis:

Geben Sie hier die Bildbeschreibung ein

Die Liniensegmente

Wenn Sie auch die Segmente zwischen den Knoten möchten, müssen Sie Ihr Shapefile "planarisieren" ( planares Diagramm , keine Kanten kreuzen sich). Dies kann durch die unary_union- Funktion von Shapely erfolgen

from shapely.ops import unary_union
graph = unary_union(lines)

Geben Sie hier die Bildbeschreibung ein

Gen
quelle
Danke @gene für die ausführliche Antwort. Ich habe den Teil bearbeitet, in dem die verschiedenen Geometrietypen behandelt werden. In meinem Fall gibt der Schnittpunkt auch Linien, Geometriesammlungen usw. zurück. Dies hängt jedoch von den Eingabedaten ab. Ich war in meiner Frage nicht klar genug.
Ustroetz
Gute Antwort. Darf ich hinzufügen, dass Folgendes nicht erforderlich ist: result = endpts.extend([pt for pt in inters if pt not in endpts])Da sich die .extendMethode anscheinend ändert endpt. In meinem Fall result = Nonenach dieser Operation. Es ist das, endptswas am Ende die Ergebniseinstellung enthält
user32882