Cómo complemento de una red de carreteras a una rejilla hexagonal en QGIS?

Question

Estoy tratando de usar QGIS 2.14 a la presión de una red de carreteras a una rejilla hexagonal, pero me estoy poniendo de artefactos extraños.

He creado un hex de cuadrícula con MMQGIS, las células son aprox 20 x 23 m. He tamponada de la red de carreteras por 1m y densificado para que haya un nodo cada pocos metros. Usted puede ver lo que estoy tratando de conseguir a continuación. Como se puede ver, puedo conseguir que funcione en algunos casos:-

• el azul es el densificado de la carretera (en el buffer de línea)
• el rojo es el 'hexified' versión - esto es lo que quiero encontrar
• el gris es el hex de la cuadrícula

A continuación, utiliza el nuevo Complemento de geometrías en función a la presión de los nodos a los más cercanos hexágono de la esquina. Los resultados son prometedores, pero parece ser que hay algunos casos extremos en que la línea se expande para llenar el hexágono (o parte de ella):-

La razón para el búfer es que el Complemento de las geometrías no se puede ajustar a una capa cuya geometría es diferente. Por ejemplo, no puedes tomar los nodos en una capa de LÍNEA a los puntos de una capa de puntos). Parece ser más felices ajuste POLÍGONO a POLÍGONO.

Sospecho que los caminos se expanda cuando uno de los lados del búfer de la carretera de la línea de saltos a un lado de la hexagonal de la célula, y por el otro lado, salta hacia el otro lado de la celda hexagonal. En mi ejemplo, las carreteras que cruzan de oeste a este en un ángulo agudo, parecen ser los peores.

Cosas que he probado, sin éxito:-

• el búfer de la red de carreteras por parte de una pequeña cantidad, por lo que sigue siendo un polígono, pero es muy delgada.
• densificar el hexagonal de las celdas (por lo que hay nodos a lo largo de los bordes, no sólo en las esquinas)
• la variación de la máxima distancia de ajuste (este tiene el mayor efecto, pero me parece que no puede encontrar un valor ideal)
• usando la LÍNEA de las capas, no de Polígonos

Me parece que si me cambio a utilizar sólo la LÍNEA de capas, funciona por un tiempo, luego se bloquea. Parece guardar su trabajo como se va - algunas de las líneas que han sido parcialmente procesados.

¿Alguien sabe de alguna otra forma a los puntos de anclaje en una línea al punto más cercano de otro línea/polígono de la capa, idealmente, sin necesidad de utilizar postgres/postgis (aunque una solución con postgis sería bienvenido demasiado)?

EDITAR

Para cualquier persona que desea ir, le he puesto un motor de arranque de QGIS proyecto aquí en Dropbox. Esto incluye el Hexagonal de la Cuadrícula y Densificado líneas de capas. (La red de carreteras es de OSM, así que puede ser descargado a través de QuickOSM por ejemplo, si usted necesita para obtener el original undensify los caminos).

Tenga en cuenta que en OSGB (epsg:27700) que un localizada UTM para el reino unido, con unidades en metros.

Answer 1

Mi solución implica una PyQGIS secuencia de comandos que es más rápido y más eficaz que un flujo de trabajo que implican ajuste (me dio una oportunidad demasiado). Usando mi algoritmo he obtenido estos resultados:

Puede ejecutar los siguientes fragmentos de código en secuencia desde dentro de QGIS (en el QGIS consola de Python). Al final, usted obtiene una capa de memoria con el acoplada rutas de carga en QGIS.

El único requisito es crear un multipart carretera Shapefile (uso Processing->Singleparts to multipart, he utilizado el campo fictitiuos como Unique ID field del parámetro). Esto nos dará un roads_multipart.shp archivo con una sola característica.

Aquí es el algoritmo explicado:

1. Obtener el más cercano de los lados del hexágono, donde las rutas de la cruz. Para cada hexágono creamos 6 triángulos entre cada par de vecino vértices y el correspondiente centro de gravedad. Si cualquier camino se cruza con un triángulo, el segmento compartido por el hexágono y el triángulo se agrega al final rompió la ruta. Esta es la parte más pesada de todo el algoritmo, se tarda 35 segundos corriendo en mi máquina. En las dos primeras líneas hay 2 Shapefile rutas de acceso, usted debe ajustar para adaptarse a sus propias rutas de acceso de archivo.

   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. Deshacerse de desconectado (o 'abrir') segmentos mediante el uso de Python listas, tuplas y diccionarios. En este punto, hay algunos desconectado segmentos de la izquierda, es decir, los segmentos que tienen un vértice desconectado, pero el otro conectado a al menos otros 2 segmentos (ver el rojo de los segmentos en la siguiente figura). Tenemos que 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 de la lista de coordenadas y de carga para 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:

Si usted necesita atributos en el acoplada rutas, podríamos utilizar un Índice Espacial para evaluar rápidamente las intersecciones (como en http://gis.stackexchange.com/a/130440/4972 ), pero esa es otra historia.

Espero que esto ayude!

Answer 2

Yo lo hice en ArcGIS, seguramente puede ser implementado utilizando QGIS o simplemente python con el paquete capaz de leer geometrías. Asegúrese de que las carreteras representan la red, es decir, se intersecan en los extremos sólo. Usted está tratando con OSM, supongo que es el caso.

• Convertir la proximidad a las líneas y los polígonos planarise ellos, por lo que convertirse en una red geométrica así.
• Colocar puntos en sus extremos – de Voronoi Puntos:
• Colocar puntos en el camino en intervalos regulares de 5 m, asegúrese de que la red de carreteras tiene buen nombre único:

• Para cada Punto de la Carretera encontrar las coordenadas de la más cercana de Voronoi Punto:
• Crear "Vías" conectando los puntos más próximos en el mismo orden:

Si no quieres ver esto:

No intente usar pk puntos de Voronoi Líneas. Me temo que esto sólo va a empeorar las cosas. Por lo tanto, su única opción es la creación de la red de Voronoi líneas y encontrar las rutas entre el final de la carretera puntos, que no es la gran cosa tampoco

Answer 3

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

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

Notas:

• Este script contiene muchos bucles dentro de los bucles y un cursor anidado. Definitivamente hay espacio para la optimización. Corrí a través de sus bases de datos en un par de minutos, pero más características agravará el problema.

Answer 4

Si se divide la carretera de la línea en segmentos, donde cada segmento fue completamente contenida por el hexágono, en su decisión en la que hexágono segmentos de línea a utilizar sería si la distancia desde la carretera dividida segmento del centroide para cada hexágono lado del punto medio fue de menos de la mitad del diámetro del hexágono (o menor que el radio de un círculo dentro del que se inscribe el hexágono).

Por lo tanto, si usted fuera a (un segmento en una hora) seleccione hexágono segmentos de línea (donde cada segmento es un lado del hexágono) que se encuentra a una distancia de la radio de la hexagonal, puede copiar los de la línea de las geometrías y la combinación de ellos en cualquiera identificador único que se utiliza para el trazado de la carretera conjunto de datos.

Si tiene problemas para combinar en el identificador único, se podría aplicar el búfer y seleccionar por ubicación en los segmentos de aplicar los atributos de la carretera conjunto de datos, de esta manera usted no tiene que preocuparse de hacer falsas coincide con un colchón que es demasiado grande.

El problema con el complemento de la herramienta es que se asiente puntos indiscriminadamente; es difícil encontrar la perfecta tolerancia a utilizar. Con esta metodología sería la correcta identificación de los cuales hexágono segmentos de línea a utilizar, a continuación, en sustitución de la geometría de la carretera de datos (o la inserción de geometrías en un conjunto de datos diferente).

También, si usted todavía tiene el problema con los segmentos de línea que saltar de un lado del hexágono para el otro, se podría dividir la línea en segmentos por vértices, calcular la longitud de cada línea, a continuación, quitar cualquiera de los segmentos de línea que son más grandes que el promedio de la longitud de uno de los lados del hexágono.