¿Cómo ajustar una red de carreteras a una cuadrícula hexagonal en QGIS?

13

Estoy tratando de usar QGIS 2.14 para ajustar una red de carreteras a una cuadrícula hexagonal, pero estoy recibiendo artefactos extraños.

He creado una cuadrícula hexagonal con MMQGIS , las celdas miden aproximadamente 20 x 23 m. He amortiguado la red de carreteras en 1 m y la ha densificado para que haya un nodo cada pocos metros. Puedes ver lo que estoy tratando de lograr a continuación. Como puede ver, puedo hacer que funcione en algunos casos: -

  • azul es la carretera densificada (una línea amortiguada)
  • rojo es la versión 'hexificada', esto es lo que quiero encontrar
  • el gris es la cuadrícula hexadecimal

LuegoutilicélanuevafunciónAjustargeometríasparaajustarlosnodosalaesquinahexagonalmáscercana.Losresultadossonprometedores,peroparecehaberalgunoscasosdebordeenlosquelalíneaseexpandeparallenarelhexágono(opartedeél):-

ElmotivodelbúferesqueGeometríasinstantáneasnotepermiteajustaraunacapacuyageometríaesdiferente.Porejemplo,nopuedeajustarlosnodosenunacapaLINEalospuntosenunacapaPOINT).ParecesermásfelizcambiarPOLYGONaPOLYGON.

Sospechoquelascarreterasseexpandencuandounladodelalíneadelacarreteraamortiguadasaltaaunladodelaceldahexagonal,yelotroladosaltaalotroladodelaceldahexagonal.Enmiejemplo,lascarreterasquecruzandeoesteaesteenunánguloagudoparecenserlaspeores.

Cosasqueheprobado,sinéxito:-

  • amortiguarlareddecarreterasenunapequeñacantidad,porloquesiguesiendounpolígonoperoesmuydelgado.
  • densificarlasceldashexadecimales(porloquehaynodosalolargodelosbordes,nosoloenlasesquinas)
  • variandoladistanciamáximadeajuste(estotieneelmayorefecto,peroparecequenopuedoencontrarunvalorideal)
  • utilizandocapasLINE,noPOLYGONs

EncuentroquesicambioausarsolocapasdeLÍNEA,funcionaporuntiempoyluegofalla.Parecequeguardasutrabajoamedidaqueavanza,algunaslíneassehanprocesadoparcialmente.

  

¿Alguiensabedealgunaotraformadeajustarpuntosenunalíneaalpuntomáscercanoenotra?  capadelínea/polígono,idealmentesinnecesidaddeusarpostgres/postgis(aunqueunasoluciónconpostgistambiénseríabienvenida)?

EDIT

Paracualquierpersonaquequieratenerunaoportunidad,hepuestounproyectoQGISdeinicio aquí en Dropbox . Esto incluye las capas Hex Grid y Densified lines. (La red de carreteras es de OSM, por lo que puede descargarse utilizando QuickOSM, por ejemplo, si necesita obtener el original para no compensar las carreteras).

Tenga en cuenta que está en OSGB (epsg: 27700), que es un UTM localizado para el Reino Unido, con unidades en metros.

    
pregunta Steven Kay 09.04.2016 - 21:47

5 respuestas

14

Mi solución involucra un script PyQGIS que es más rápido y más efectivo que un flujo de trabajo que involucra ajustes (también lo intenté). Usando mi algoritmo he obtenido estos resultados:

PuedeejecutarlossiguientesfragmentosdecódigoensecuenciadesdeQGIS(enlaconsolaPythondeQGIS).Alfinal,obtienesunacapadememoriaconlasrutasajustadascargadasenQGIS.

ElúnicorequisitoprevioescrearunShapefiledecarreteradevariaspartes(useProcessing->Singlepartstomultipart,uséelcampofictitiuoscomoparámetroUniqueIDfield).Estonosdaráunarchivoroads_multipart.shpconunasolafunción.

Aquíestáelalgoritmoexplicado:

  1. Obtengalosladosdelhexágonomáscercanosdondesecruzanlasrutas.Paracadahexágonocreamos6triángulosentrecadapardevérticesvecinosyelcentroidecorrespondiente.Sialgunacarreteracruzauntriángulo,elsegmentocompartidoporelhexágonoyeltriánguloseagregaalarutafinalajustada.Estaeslapartemáspesadadetodoelalgoritmo,tarda35segundosenejecutarseenmimáquina.Enlasdosprimeraslíneashay2rutasdeShapefile,debesajustarlasparaqueseajustenatuspropiasrutasdearchivo.

    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. Deshazte de los segmentos desconectados (o 'abiertos') utilizando listas, tuplas y diccionarios de Python . En este punto, quedan algunos segmentos desconectados, es decir, segmentos que tienen un vértice desconectado pero el otro conectado a al menos otros 2 segmentos (vea los segmentos rojos en la siguiente figura). Necesitamos deshacernos de ellos.

    # 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. Ahora podemos crear una capa vectorial desde la lista de coordenadas y cargarla en el mapa de 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)
    

Otra parte del resultado:

Sinecesitaatributosenlasrutasajustadas,podríamosusarunÍndiceespacialparaevaluarrápidamentelasintersecciones(comoen enlace ), pero esa es otra historia.

Espero que esto ayude!

    
respondido por el Germán Carrillo 16.04.2016 - 01:02
6

Lo hice en ArcGIS, seguramente se puede implementar utilizando QGIS o simplemente python con un paquete capaz de leer geometrías. Asegúrese de que las carreteras representan la red, es decir, se intersectan entre sí solo en los extremos. Usted está tratando con OSM, supongo que es el caso.

  • Convierta los polígonos de proximidad en líneas y planévelos, para que Conviértete también en una red geométrica.
  • Coloque puntos en sus extremos - Puntos Voronoi:
  • Coloquepuntosenlacarreteraaintervalosregularesde5m,asegúresedeLoscaminosdelaredtienenunbuennombreúnico:

  • Paracadapuntodecarretera,encuentrelascoordenadasdelpuntoVoronoimáscercano:
  • Crea"Carreteras" conectando los puntos más cercanos en el mismo orden:

Sinoquieresveresto:

No intentes utilizar puntos de PK en líneas Voronoi. Me temo que solo lo empeorará. Por lo tanto, su única opción es crear una red a partir de las líneas de Voronoi y encontrar rutas entre puntos finales de carreteras, eso tampoco es un gran problema.

    
respondido por el FelixIP 14.04.2016 - 23:13
4

Me doy cuenta de que estás pidiendo un método QGIS, pero ten paciencia conmigo para una respuesta arcpy:

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,["[email protected]","[email protected]"], spatial_reference=sr) as r_cursor: # loop through roads
    for r_row in r_cursor:
        with arcpy.da.SearchCursor(hexgrid,["[email protected]","[email protected]"], 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

Notas:

  • Este script contiene muchos bucles dentro de bucles y un cursor anidado. Definitivamente hay espacio para la optimización. Revisé tus conjuntos de datos en un par de minutos, pero más características agravarán el problema.
respondido por el phloem 15.04.2016 - 20:54
2

Si tuviera que dividir la línea de la carretera en segmentos donde cada segmento estuviera completamente contenido por el hexágono, su decisión sobre qué segmentos de la línea del hexágono usar sería si la distancia desde el centroide del segmento de la carretera dividida al punto medio de cada lado del hexágono era menor menos de la mitad del diámetro del hexágono (o menos del radio de un círculo que encaja dentro del hexágono).

Por lo tanto, si tuviera que (un segmento a la vez) seleccionar segmentos de línea hexagonales (donde cada segmento es un lado del hexágono) que están dentro de una distancia del radio del hexágono, podría copiar esas geometrías de línea y fusionarlos en cualquier identificador único que use para su conjunto de datos de carreteras.

Si tiene problemas para fusionar el identificador único, puede aplicar el búfer y seleccionar solo por ubicación en esos segmentos para aplicar los atributos de su conjunto de datos de carreteras; de esa manera, no tendría que preocuparse por hacer coincidencias falsas con un búfer demasiado grande.

El problema con la herramienta de ajuste es que ajusta puntos indiscriminadamente; Es difícil encontrar esa tolerancia perfecta para usar. Con esta metodología, estaría identificando correctamente qué segmentos de línea hexagonal utilizar y luego reemplazaría la geometría de sus datos de carretera (o insertaría las geometrías en un conjunto de datos diferente).

Además, si aún tiene el problema con los segmentos de línea que saltan de un lado del hexágono al otro, puede dividir la línea en segmentos por vértices, calcular la longitud de cada línea y luego eliminar cualquier segmento de línea que son más grandes que la longitud promedio de un lado del hexágono.

    
respondido por el crld 15.04.2016 - 00:36
1

El snapper de geometría en qgis 3.0 se ha reelaborado y ahora permite el ajuste entre diferentes tipos de geometría. También tiene muchas correcciones. Puede probar una versión de "instantánea diaria" para obtener acceso al snapper mejorado antes del lanzamiento oficial de la versión 3.0.

    
respondido por el ndawson 13.12.2016 - 19:49

Lea otras preguntas en las etiquetas