Przecinające się linie, aby uzyskać skrzyżowania za pomocą Pythona z QGIS?

Mam zestaw linii reprezentujących linie autobusowe. Niektóre linie nachodzą na siebie i biegną tymi samymi drogami.

Jestem w stanie wyodrębnić węzły.

Jednak interesuje mnie wyodrębnianie tylko takich przejść:

W jaki sposób mogę to zrobić? Szukam sposobów na QGIS lub Python.

Próbowałem metody przecięcia z GDAL Python, ale to w zasadzie zwraca mi tylko wierzchołki.

Metoda przecięcia linii z QGIS zwraca mi skrzyżowania, jeśli przecinają się dwie linie. Jednak w przypadku, gdy dwie linie autobusowe idą daleko na części swojej trasy na tej samej drodze, nie daje mi to punktu, w którym się łączą.

Węzły:

Potrzebujesz dwóch rzeczy, punktów końcowych polilinii (bez węzłów pośrednich) i punktów przecięcia. Istnieje dodatkowy problem, niektóre punkty końcowe polilinii są również punktami przecięcia:

Rozwiązaniem jest użycie Pythona i modułów Shapely i Fiona

1) Przeczytaj plik kształtu:

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

2) Znajdź punkty końcowe linii ( jak uzyskać punkty końcowe polilinii? ):

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] 

3) Oblicz przecięcia (iteracja przez pary geometrii w warstwie z modułem itertools ). Wynikiem niektórych skrzyżowań są MultiPoints i chcemy listę punktów:

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

4) Wyeliminuj duplikaty między punktami końcowymi i punktami przecięcia (jak widać na rysunkach)

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

5) Zapisz wynikowy plik kształtu

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

Ostateczny wynik:

Segmenty linii

Jeśli chcesz także segmentów między węzłami, musisz „planować” ( wykres planarny , żadne krawędzie nie krzyżują się) swój plik kształtu. Można tego dokonać za pomocą funkcji unary_union programu Shapely

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