Como conectar uma rede viária a uma grade hexagonal no QGIS?

13

Estou tentando usar o QGIS 2.14 para conectar uma rede de estradas a uma grade hexagonal, mas estou obtendo artefatos estranhos.

Eu criei uma grade hexagonal com MMQGIS , as células são aproximadamente 20 x 23 m. Eu tampei a rede rodoviária em 1m e a densifiquei para que haja um nó a cada poucos metros. Você pode ver o que estou tentando alcançar abaixo. Como você pode ver, posso fazê-lo funcionar em alguns casos: -

  • azul é a estrada densificada (uma linha em buffer)
  • vermelho é a versão 'hexificada' - é isso que eu quero encontrar
  • o cinza é a grade hexagonal

insira a descrição da imagem aqui

Em seguida, usei o novo recurso Geometrias de snap para encaixar os nós no canto hexagonal mais próximo. Os resultados são promissores, mas parece haver alguns casos extremos em que a linha se expande para preencher o hexágono (ou parte dele): -

insira a descrição da imagem aqui

O motivo do buffer é que as geometrias de Snap não permitem que você faça snap em uma camada cuja geometria é diferente. Por exemplo, você não pode encaixar nós em uma camada LINE para pontos em uma camada POINT). Parece ser o mais feliz encaixar POLYGON em POLYGON.

Suspeito que as estradas se expandem quando um lado da linha da estrada com buffer salta para um lado da célula hexadecimal e o outro lado pula para o outro lado da célula hexadecimal. No meu exemplo, as estradas que cruzam oeste-leste em um ângulo agudo parecem ser as piores.

Coisas que eu tentei, sem sucesso: -

  • protegendo a rede rodoviária por uma pequena quantidade, permanecendo um polígono, mas é muito fino.
  • densificar as células hexadecimais (para que haja nós ao longo das bordas, não apenas nos cantos)
  • variando a distância máxima de encaixe (isso tem o maior efeito, mas não consigo encontrar um valor ideal)
  • usando camadas LINE, não POLYGONs

Acho que, se eu mudar para o uso de apenas camadas de linha, ele funcionará por um tempo e travará. Parece salvar seu trabalho à medida que avança - algumas linhas foram parcialmente processadas.

insira a descrição da imagem aqui

Alguém sabe de outra maneira de ajustar pontos em uma linha para o ponto mais próximo em outra camada de linha / polígono, idealmente sem precisar usar o postgres / postgis (embora uma solução com o postgis também seja bem-vinda)?

EDITAR

Para quem quiser experimentar, coloquei um projeto inicial do QGIS aqui no Dropbox . Isso inclui as camadas de grade hexadecimal e de densidade. (A rede rodoviária é do OSM, portanto, pode ser baixada usando o QuickOSM, por exemplo, se você precisar fazer com que o original subestime as estradas).

Observe que está no OSGB (epsg: 27700) que é um UTM localizado para o Reino Unido, com unidades em metros.

Steven Kay
fonte
3
Você poderia compartilhar um conjunto de dados de amostra? Gostaria de tentar, mas não quero passar pelo processo de criação de dados de amostra do zero.
Germán Carrillo
@ GermánCarrillo - obrigado. Adicionei um link a um projeto de exemplo para a pergunta.
21416 Steven

Respostas:

14

Minha solução envolve um script PyQGIS que é mais rápido e eficaz do que um fluxo de trabalho que envolve snapping (tentei também). Usando meu algoritmo, obtive estes resultados:

insira a descrição da imagem aqui

insira a descrição da imagem aqui

Você pode executar os seguintes trechos de código em sequência no QGIS (no console do QGIS Python). No final, você obtém uma camada de memória com as rotas encaixadas no QGIS.

O único pré-requisito é criar um Shapefile de estrada com várias partes (use Processing->Singleparts to multipart, usei o campo fictitiuoscomo Unique ID fieldparâmetro). Isso nos dará um roads_multipart.shparquivo com um único recurso.

Aqui está o algoritmo explicado:

  1. Obtenha os lados hexagonais mais próximos onde as rotas se cruzam. Para cada hexágono, criamos 6 triângulos entre cada par de vértices vizinhos e o centróide correspondente. Se alguma estrada cruzar um triângulo, o segmento compartilhado pelo hexágono e pelo triângulo será adicionado à rota final ajustada. Esta é a parte mais pesada de todo o algoritmo, que demora 35 segundos em execução na minha máquina. Nas duas primeiras linhas existem 2 caminhos de Shapefile, você deve ajustá-los para se ajustarem aos seus próprios caminhos de arquivo.

    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. Livre-se dos segmentos desconectados (ou 'abertos') usando listas, tuplas e dicionários do Python . Nesse ponto, existem alguns segmentos desconectados restantes, ou seja, segmentos que têm um vértice desconectado, mas o outro conectado a pelo menos outros 2 segmentos (veja segmentos vermelhos na próxima figura). Precisamos nos livrar deles.

    insira a descrição da imagem aqui

    # 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. Agora podemos criar uma camada vetorial da lista de coordenadas e carregá-la no mapa QGIS :

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

Outra parte do resultado:

insira a descrição da imagem aqui

Se você precisar de atributos nas rotas ajustadas, poderíamos usar um Índice Espacial para avaliar rapidamente interseções (como em /gis//a/130440/4972 ), mas isso é outra história.

Espero que isto ajude!

Germán Carrillo
fonte
1
obrigado, isso funciona perfeitamente! Tive problemas para colá-lo no console python ... salvei-o como um arquivo .py no editor pyg qgis e ele correu bem a partir daí. A etapa multipartes remove os atributos, mas uma junção buffer / espacial corrige isso!
21416 Steven
1
Ótimo! Ainda bem que finalmente resolveu o problema que você estava enfrentando. Estou interessado em saber qual é o caso de uso com o qual você está lidando. Você acha que poderíamos aproveitar isso para se tornar um plug-in QGIS ou talvez um script incluído no Processing scripts?
Germán Carrillo
1
o caso de uso que eu tinha em mente eram mapas de transporte público, como o Mapa do Tubo, em que é necessário ajustar linhas a uma grade em mosaico ou a um conjunto restrito de ângulos. Isso pode ser feito manualmente, digitalizando, mas eu estava interessado em ver se poderia ser automatizado. Usei hexágonos, pois eram fáceis de gerar, visualmente interessantes e tinham ângulos que não eram retos. Eu acho que isso vale a pena olhar com mais detalhes, especialmente se ele poderia ser generalizado para o trabalho com outros Tesselations ...
Steven Kay
1
A idéia por trás do script funcionaria em grades de triângulos, quadrados, pentágonos, hexágonos e assim por diante.
Germán Carrillo
6

Eu fiz isso no ArcGIS, certamente pode ser implementado usando QGIS ou simplesmente python com pacote capaz de ler geometrias. Certifique-se de que as estradas representem a rede, ou seja, se cruzem apenas nas extremidades. Você está lidando com OSM, suponho que seja o caso.

  • Converta polígonos de proximidade em linhas e planarize-os, para que também se tornem uma rede geométrica.
  • Coloque pontos nas suas extremidades - Pontos Voronoi: insira a descrição da imagem aqui
  • Coloque pontos na estrada em intervalos regulares de 5 m, verifique se as estradas da rede têm um bom nome exclusivo:

insira a descrição da imagem aqui insira a descrição da imagem aqui

  • Para cada ponto da estrada, encontre as coordenadas do ponto Voronoi mais próximo: insira a descrição da imagem aqui
  • Crie "Estradas" conectando os pontos mais próximos na mesma ordem: insira a descrição da imagem aqui

Se você não quiser ver isso: insira a descrição da imagem aqui

Não tente usar pontos de corrente nas linhas Voronoi. Eu temia que isso só piorasse. Portanto, sua única opção é criar uma rede a partir das linhas de Voronoi e encontrar rotas entre os pontos finais da estrada, o que também não é grande coisa.

FelixIP
fonte
Isso é ótimo, obrigado! Você menciona o uso de linhas voronoi, que não estão muito familiarizadas com isso (Voronois por pontos, eu posso entender). Você quer dizer que cada linha é cercada por um polígono de todos os pontos mais próximos dessa linha? (Não conheço uma maneira de fazer isso no QGIS). Ou você quer dizer as linhas de contorno de uma malha voronoi normal, com base em pontos?
21416 Steven
Linhas de fronteira dos polígonos de proximidade. Btw eu parei muito cedo. Para tarefa completa é suficiente para dividir 1ª resultado no vértice, adicionar pontos no processo de médio e repita
FelixIP
4

Sei que você está pedindo um método QGIS, mas aceite uma resposta complicada:

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

insira a descrição da imagem aqui

Notas:

  • Este script contém muitos loops dentro de loops e um cursor aninhado. Definitivamente, há espaço para otimização. Examinei seus conjuntos de dados em alguns minutos, mas mais recursos irão compor o problema.
floema
fonte
Obrigado por isso, muito apreciado. Isso mostra exatamente o efeito que eu estava visualizando. Os comentários copiosos significam que eu posso entender o que você está fazendo, mesmo que eu não consiga executar o código. Embora seja um arco, tenho certeza de que isso será possível em pyqgis. As idéias algoritmo aqui são interessantes (especialmente olhando tanto no sentido horário e anti-horário em volta de cada hex, e escolher o mais curto maneira redonda)
Steven Kay
2

Se você dividir a linha da estrada em segmentos em que cada segmento esteja completamente contido pelo hexágono, sua decisão sobre quais segmentos de linha do hexágono usar seria se a distância do centróide do segmento da estrada dividida ao ponto médio de cada lado do hexágono era menos da metade da diâmetro do hexágono (ou menor que o raio de um círculo que se encaixa dentro do hexágono).

Portanto, se você (um segmento de cada vez) seleciona segmentos de linhas hexagonais (em que cada segmento é um lado do hexágono) que estão a uma distância do raio do hexágono, você pode copiar essas geometrias de linhas e mesclá-las em qualquer identificador exclusivo usado para o conjunto de dados da estrada.

Se você tiver problemas para mesclar no identificador exclusivo, poderá aplicar o buffer e selecionar por local somente esses segmentos para aplicar os atributos do seu conjunto de dados da estrada; Dessa forma, você não precisaria se preocupar em fazer correspondências falsas com um buffer muito grande.

O problema com a ferramenta snap é que ela encaixa pontos indiscriminadamente; é difícil encontrar a tolerância perfeita para usar. Com essa metodologia, você identificaria corretamente quais segmentos de linhas hexagonais usar e, em seguida, substituiria a geometria dos dados de sua estrada (ou inseriria as geometrias em um conjunto de dados diferente).

Além disso, se você ainda tiver problemas com os segmentos de linha que saltam de um lado do hexágono para o outro, poderá dividir a linha em segmentos por vértices, calcular o comprimento de cada linha e remover quaisquer segmentos de linha maiores que o comprimento médio de um lado do hexágono.

crld
fonte
1

O snapper de geometria no qgis 3.0 foi reformulado e agora permite o snap entre diferentes tipos de geometria. Ele também tem muitas correções. Você pode tentar uma versão "instantâneo diário" para obter acesso ao snapper aprimorado antes que o 3.0 seja lançado oficialmente.

ndawson
fonte