Пересекающиеся линии, чтобы получить пересечения, используя Python с QGIS?

10

У меня есть набор линий, представляющих автобусные линии. Некоторые линии перекрываются и проходят по одним и тем же дорогам.

введите описание изображения здесь

Я могу извлечь узлы. введите описание изображения здесь

Однако я заинтересован в извлечении только таких переходов: введите описание изображения здесь

Как я могу это сделать? Я ищу способы с QGIS или Python.

Я попробовал метод пересечения из GDAL Python, но это в основном возвращает мне только вершины.

Метод пересечения линий из QGIS возвращает мне пересечения, если две линии пересекаются. Однако в том случае, если две автобусные линии проходят большую часть своего маршрута по одной и той же дороге, это не дает мне точки, где они объединяются.

ustroetz
источник
Пробовали ли вы использовать инструмент пересечения линий в QGIS: Инструмент векторного анализа> Пересечения линий ... Он не даст вам конечные и начальные узлы линии, но все пересечения.
Якоб
Да, я написал об этом в вопросе.
Устроец
Мне непонятно, о чем вы спрашиваете, отчасти потому, что все строки на ваших изображениях одинаково обозначены - я не могу отличить разные маршруты, чтобы понять, на какие узлы вы смотрите или почему их так много в второе изображение. Совпадают ли маршруты на дорогах? Все ли это двухточечные отрезки или непрерывные полилинии? Я отмечаю, что описанное вами поведение аналогично инструменту пересечения ArcGIS: линии / линии с выводом линий дают вам перекрытие, а линии / линии с точечным выводом дают только пересечения.
Крис W
Исходя из этого, чтобы получить то, что я думаю, вы хотите, вам, возможно, придется использовать оба метода. Получите пересечения (line / line = point), затем получите перекрытия (line / line = line) и извлеките начальные / конечные узлы для этих линий перекрытия. Это должны быть все точки / узлы, которые вы ищете.
Крис Вт

Ответы:

20

Узлы:

Вы хотите две вещи, конечные точки полилиний (без промежуточных узлов) и точки пересечения. Есть дополнительная проблема, некоторые конечные точки полилиний также являются точками пересечения:

введите описание изображения здесь

Решением является использование Python и модулей Shapely и Fiona.

1) Прочтите шейп-файл:

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

2) Найти конечные точки линий ( как получить конечные точки ломаной линии ? ):

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) Вычислить пересечения (перебирая пары геометрий в слое с помощью модуля itertools ). Результатом некоторых пересечений являются MultiPoints, и мы хотим получить список точек:

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) Устраните дубликаты между конечными точками и точками пересечения (как вы можете видеть на рисунках)

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

5) Сохраните полученный шейп-файл

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

Конечный результат:

введите описание изображения здесь

Сегменты линии

Если вы хотите также отрезки между узлами, вам необходимо «спланировать» ( планарный график , ребра не пересекаются) ваш шейп-файл. Это можно сделать с помощью функции unary_union в Shapely

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

введите описание изображения здесь

ген
источник
Спасибо @gene за подробный ответ. Я отредактировал ту часть, где описываются разные типы геометрии. В моем случае пересечение также возвращает линии, геометрические коллекции и т. Д. Но это зависит от входных данных. Я не был достаточно ясен в своем вопросе.
Устроец
Отличный ответ. Могу ли я добавить, что нет необходимости делать следующее: result = endpts.extend([pt for pt in inters if pt not in endpts])так как кажется, что .extendметод модифицируется endpt. В моем случае result = Noneпосле этой операции. Это то, endptsчто в итоге содержит результат
урегулирования