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

13

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

entrez la description de l'image ici

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

entrez la description de l'image ici

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.

entrez la description de l'image ici

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.

Steven Kay
la source
3
Pourriez-vous s'il vous plaît partager un exemple de jeu de données? Je voudrais essayer, mais je ne veux pas revenir sur le processus de création d'échantillons de données à partir de zéro.
Germán Carrillo
@ GermánCarrillo - merci. J'ai ajouté un lien vers un exemple de projet à la question.
Steven Kay

Réponses:

14

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:

entrez la description de l'image ici

entrez la description de l'image ici

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.

    entrez la description de l'image ici

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

entrez la description de l'image ici

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!

Germán Carrillo
la source
1
merci, cela fonctionne parfaitement! J'ai eu des problèmes pour le coller dans la console python ... Je l'ai enregistré en tant que fichier .py dans l'éditeur python qgis, et cela a bien fonctionné à partir de là. L'étape en plusieurs parties supprime les attributs, mais une jointure tampon / spatiale corrigera cela!
Steven Kay
1
Génial! Heureux qu'il ait finalement résolu le problème que vous rencontriez. J'aimerais savoir quel est le cas d'utilisation auquel vous avez affaire. Pensez-vous que nous pourrions en tirer parti pour devenir un plugin QGIS ou peut-être un script inclus dans les scripts de traitement?
Germán Carrillo
1
le cas d'utilisation que j'avais à l'esprit était des cartes de transports en commun comme la Tube Map, où vous devez accrocher des lignes à une grille tesselée ou à un ensemble d'angles restreint. Cela peut être fait manuellement par numérisation, mais j'étais intéressé de voir s'il pouvait être automatisé. J'ai utilisé des hexagones car ils étaient faciles à générer, visuellement intéressants et avaient des angles qui n'étaient pas droits. Je pense que cela vaut la peine d'être examiné plus en détail, surtout si cela pourrait être généralisé pour fonctionner avec d'autres tesselations ...
Steven Kay
1
L'idée derrière le script fonctionnerait sur des grilles de triangles, de carrés, de pentagones, d'hexagones, etc.
Germán Carrillo
6

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: entrez la description de l'image ici
  • 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:

entrez la description de l'image ici entrez la description de l'image ici

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

Si vous ne voulez pas voir ceci: entrez la description de l'image ici

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

FelixIP
la source
c'est super, merci! Vous mentionnez l'utilisation de lignes voronoi, pas trop familières avec cela (Voronois à partir de points, je peux comprendre). Voulez-vous dire que chaque ligne est entourée d'un polygone de tous les points les plus proches de cette ligne? (Je ne connais pas de façon de faire cela dans QGIS). Ou voulez-vous dire les lignes de démarcation d'un maillage voronoï normal, basé sur des points?
Steven Kay
Lignes limites de polygones de proximité. Btw je me suis arrêté trop tôt. Pour terminer la tâche, il suffit de diviser le premier résultat au sommet, d'ajouter un point au milieu et de répéter le processus
FelixIP
4

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

entrez la description de l'image ici

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.
phloème
la source
Merci pour cela, très apprécié. Cela montre exactement l'effet que je visualisais. Les commentaires copieux signifient que je peux obtenir l'essentiel de ce que vous faites même si je ne peux pas exécuter le code. Bien que ce soit arcpy, je suis sûr que ce sera faisable dans pyqgis. Les idées d'algorithmes ici sont intéressantes (en particulier en regardant à la fois dans le sens horaire et antihoraire dans chaque hex, et en choisissant le chemin le plus court)
Steven Kay
2

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.

crld
la source
1

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.

ndawson
la source