Comment accrocher un réseau routier à une grille hexagonale dans QGIS?

J'essaie d'utiliser QGIS 2.14 pour aligner un réseau routier sur une grille hexagonale, mais j'obtiens d'étranges artefacts.

J'ai créé une grille hexagonale avec MMQGIS , les cellules mesurent environ 20 x 23 m. J'ai tamponné le réseau routier de 1 m et je l'ai densifié pour qu'il y ait un nœud tous les quelques mètres. Vous pouvez voir ce que j'essaie de réaliser ci-dessous. Comme vous pouvez le voir, je peux le faire fonctionner dans certains cas: -

• le bleu est la route densifiée (une ligne tamponnée)
• le rouge est la version «hexifiée» - c'est ce que je veux trouver
• le gris est la grille hexagonale

J'ai ensuite utilisé la nouvelle fonction Snap geometries pour accrocher les nœuds au coin hexagonal le plus proche. Les résultats sont prometteurs, mais il semble y avoir des cas marginaux où la ligne se dilate pour remplir l'hexagone (ou une partie de celui-ci): -

La raison du tampon est que les géométries d'accrochage ne vous permettent pas de vous accrocher à une couche dont la géométrie est différente. Par exemple, vous ne pouvez pas accrocher des nœuds sur une couche LINE à des points sur une couche POINT). Il semble être le plus heureux de claquer POLYGON en POLYGON.

Je soupçonne que les routes se dilatent lorsqu'un côté de la ligne de route tamponnée saute d'un côté de la cellule hexadécimale et que l'autre côté saute de l'autre côté de la cellule hexadécimale. Dans mon exemple, les routes qui traversent l'ouest-est à un angle aigu semblent être les pires.

Des choses que j'ai essayées, sans succès: -

• tampon du réseau routier par une petite quantité, il reste donc un polygone mais est très mince.
• densifier les cellules hexagonales (donc il y a des nœuds le long des bords, pas seulement aux coins)
• faire varier la distance de capture maximale (cela a le plus grand effet, mais je n'arrive pas à trouver une valeur idéale)
• en utilisant des couches LINE, pas des POLYGONs

Je trouve que si je passe à l'utilisation de calques LINE uniquement, cela fonctionne pendant un certain temps, puis se bloque. Il semble enregistrer son travail au fur et à mesure - certaines lignes ont été partiellement traitées.

Quelqu'un connaît-il une autre façon d'aligner les points d'une ligne sur le point le plus proche d'une autre couche de ligne / polygone, idéalement sans avoir besoin d'utiliser postgres / postgis (bien qu'une solution avec postgis soit également la bienvenue)?

ÉDITER

Pour tous ceux qui souhaitent essayer, j'ai mis un projet QGIS de démarrage ici sur Dropbox . Cela inclut les couches Grille hexagonale et Lignes densifiées. (Le réseau routier provient d'OSM, il peut donc être téléchargé à l'aide de QuickOSM, par exemple si vous devez obtenir l'original pour densifier les routes).

Notez que c'est en OSGB (epsg: 27700) qui est une UTM localisée pour le Royaume-Uni, avec des unités en mètres.

Ma solution implique un script PyQGIS qui est plus rapide et plus efficace qu'un workflow impliquant un accrochage (je l'ai également essayé). En utilisant mon algorithme, j'ai obtenu ces résultats:

Vous pouvez exécuter les extraits de code suivants dans l'ordre à partir de QGIS (dans la console QGIS Python). À la fin, vous obtenez une couche mémoire avec les routes cassées chargées dans QGIS.

La seule condition préalable est de créer un Shapefile de route en plusieurs parties (utiliser Processing->Singleparts to multipart, j'ai utilisé le champ fictitiuoscomme Unique ID fieldparamètre). Cela nous donnera un roads_multipart.shpfichier avec une seule fonctionnalité.

Voici l'algorithme expliqué:

1. Obtenez les côtés hexagonaux les plus proches où les routes se croisent. Pour chaque hexagone, nous créons 6 triangles entre chaque paire de sommets voisins et le centroïde correspondant. Si une route croise un triangle, le segment partagé par l'hexagone et le triangle est ajouté à l'itinéraire final capturé. C'est la partie la plus lourde de tout l'algorithme, cela prend 35 secondes de course sur ma machine. Dans les deux premières lignes, il y a 2 chemins Shapefile, vous devez les ajuster pour qu'ils correspondent à vos propres chemins de fichier.

   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. Débarrassez-vous des segments déconnectés (ou «ouverts») en utilisant des listes, des tuples et des dictionnaires Python . À ce stade, il reste des segments déconnectés, c'est-à-dire des segments qui ont un sommet déconnecté mais l'autre connecté à au moins 2 autres segments (voir les segments rouges sur la figure suivante). Nous devons nous en débarrasser.

   

   # 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. Nous pouvons maintenant créer une couche vectorielle à partir de la liste des coordonnées et la charger sur la carte 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)

Une autre partie du résultat:

Si vous avez besoin d'attributs dans les itinéraires capturés, nous pourrions utiliser un index spatial pour évaluer rapidement les intersections (comme dans /gis//a/130440/4972 ), mais c'est une autre histoire.

J'espère que cela t'aides!

Je l'ai fait dans ArcGIS, peut sûrement être implémenté en utilisant QGIS ou simplement python avec un package capable de lire les géométries. Assurez-vous que les routes représentent le réseau, c'est-à-dire qu'elles se croisent aux extrémités uniquement. Vous avez affaire à OSM, je suppose que c'est le cas.

• Convertissez des polygones de proximité en lignes et planarisez-les, afin qu'ils deviennent également un réseau géométrique.
• Placer des points à leurs extrémités - Points Voronoi:
• Placez des points sur la route à intervalles réguliers de 5 m, assurez-vous que les routes du réseau ont un bon nom unique:

• Pour chaque Road Point, trouvez les coordonnées du Voronoi Point le plus proche:
• Créez des «routes» en connectant les points les plus proches dans le même ordre:

Si vous ne voulez pas voir ceci:

N'essayez pas d'utiliser des points de chaînage sur les lignes Voronoi. J'ai peur que cela ne fasse qu'empirer les choses. Ainsi, votre seule option consiste à créer un réseau à partir des lignes de Voronoi et à trouver des itinéraires entre les points d'extrémité de la route, ce n'est pas très grave non plus

Je sais que vous demandez une méthode QGIS, mais soyez indulgent avec moi pour une réponse 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,["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

Remarques:

• Ce script contient de nombreuses boucles dans les boucles et un curseur imbriqué. Il y a certainement place à l'optimisation. J'ai parcouru vos jeux de données en quelques minutes, mais d'autres fonctionnalités aggraveront le problème.

Si vous deviez diviser la ligne de route en segments où chaque segment était entièrement contenu par l'hexagone, votre décision sur les segments de ligne d'hexagone à utiliser serait de savoir si la distance entre le centre de gravité du segment de route divisé et le milieu de chaque côté de l'hexagone était inférieure à la moitié de la diamètre de l'hexagone (ou inférieur au rayon d'un cercle qui s'inscrit à l'intérieur de l'hexagone).

Ainsi, si vous deviez (un segment à la fois) sélectionner des segments de ligne hexagonale (où chaque segment est un côté de l'hexagone) qui sont à une distance du rayon de l'hexagone, vous pourriez copier ces géométries de ligne et les fusionner sur quel que soit l'identifiant unique que vous utilisez pour votre jeu de données routières.

Si vous avez des difficultés à fusionner sur l'identifiant unique, vous pouvez appliquer le tampon et sélectionner par emplacement uniquement sur ces segments pour appliquer les attributs de votre jeu de données de route; De cette façon, vous n'aurez pas à vous soucier de fausses correspondances avec un tampon trop volumineux.

Le problème avec l'outil d'accrochage est qu'il accroche des points sans discernement; il est difficile de trouver cette tolérance parfaite à utiliser. Avec cette méthodologie, vous identifieriez correctement les segments de ligne hexagonaux à utiliser, puis remplaceriez la géométrie de vos données de route (ou insérez les géométries dans un autre ensemble de données).

De plus, si le problème persiste avec les segments de ligne qui sautent d'un côté à l'autre de l'hexagone, vous pouvez diviser la ligne en segments par sommets, calculer la longueur de chaque ligne, puis supprimer tous les segments de ligne plus grands que la longueur moyenne d'un côté de l'hexagone.

Le snapper de géométrie dans qgis 3.0 a été retravaillé et permet désormais d'aligner entre différents types de géométrie. Il a également beaucoup de correctifs. Vous pouvez essayer une version "instantané quotidien" pour accéder au snapper amélioré avant la sortie officielle de la version 3.0.