Interseção de linhas para obter cruzamentos usando Python com QGIS?

10

Eu tenho um conjunto de linhas representando linhas de ônibus. Algumas das linhas estão sobrepostas e seguem as mesmas estradas.

insira a descrição da imagem aqui

Eu sou capaz de extrair os nós. insira a descrição da imagem aqui

No entanto, estou interessado em extrair apenas cruzamentos como este: insira a descrição da imagem aqui

Como posso fazer isso? Estou procurando maneiras com QGIS ou Python.

Eu tentei o método de interseção do GDAL Python, mas isso basicamente me retorna apenas os vértices.

O método Intersecções de linha do QGIS retorna as passagens se duas linhas se cruzarem. No entanto, no caso de duas linhas de ônibus percorrerem parte da rota na mesma estrada, isso não indica que eles apontam para onde se fundem.

ustroetz
fonte
Você já experimentou a ferramenta de interseção de linhas no QGIS: Ferramenta de Análise de Vetor> Interseções de Linhas ... Ela não fornecerá os nós finais e iniciais de uma linha, mas todas as interseções.
Jakob
Sim, eu escrevi sobre isso na pergunta.
ustroetz
Não sei ao certo o que você está perguntando, em parte porque todas as linhas são simbolizadas da mesma maneira em suas imagens - não consigo distinguir rotas diferentes para entender quais nós você está vendo ou por que existem tantos no segunda imagem. As rotas são coincidentes nas estradas? São todos segmentos de linha de dois pontos ou polilinhas contínuas? Observo que o comportamento que você descreve é ​​o mesmo que a ferramenta Intersect do ArcGIS - linhas / linhas com saída de linhas dão sobreposição, mas linhas / linhas com saída de ponto apenas fornecem cruzamentos.
22415 Chris N,
Com base nisso, para obter o que acho que você deseja, pode ser necessário usar os dois métodos. Obtenha os cruzamentos (linha / linha = ponto) e, em seguida, obtenha as sobreposições (linha / linha = linha) e extraia os nós de início / fim dessas linhas de sobreposição. Esses devem ser todos os pontos / nós que você está procurando.
27515 Chris

Respostas:

20

Os nós:

Você deseja duas coisas: os pontos finais das polilinhas (sem nós intermediários) e os pontos de interseção. Há um problema adicional, alguns pontos finais das polilinhas também são pontos de interseção:

insira a descrição da imagem aqui

Uma solução é usar Python e os módulos Shapely e Fiona

1) Leia o shapefile:

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

2) Encontre os pontos finais das linhas ( como obter os pontos finais de uma polilinha? ):

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] 

insira a descrição da imagem aqui

3) Calcule as interseções (iterando através de pares de geometrias na camada com o módulo itertools ). O resultado de algumas interseções são MultiPoints e queremos uma lista de pontos:

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

insira a descrição da imagem aqui

4) Elimine duplicatas entre os pontos finais e os pontos de interseção (como você pode ver nas figuras)

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

5) Salve o shapefile resultante

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

Resultado final:

insira a descrição da imagem aqui

Os segmentos de linha

Se você quiser também os segmentos entre os nós, precisará "planarizar" ( Gráfico Planar , nenhuma aresta se cruza) no seu arquivo de forma. Isso pode ser feito pela função unary_union do Shapely

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

insira a descrição da imagem aqui

gene
fonte
Obrigado @gene pela resposta detalhada. Eu editei a parte em que ela aborda os diferentes tipos de geometria. No meu caso, a interseção também retorna linhas, geometrias, etc. Mas isso depende dos dados de entrada. Eu não estava claro o suficiente na minha pergunta.
ustroetz
Ótima resposta. Posso acrescentar que não é necessário fazer o seguinte: result = endpts.extend([pt for pt in inters if pt not in endpts])pois parece que o .extendmétodo é modificado endpt. No meu caso result = Nonedepois dessa operação. É endptsque acaba contendo o resultado sett
user32882