How to snap a road network to a hexagonal grid in QGIS?

Question

I'm trying to use QGIS 2.14 to snap a road network to a hexagonal grid, but I'm getting strange artifacts.

I've created a hex grid with MMQGIS, cells are approx 20 x 23 m. I've buffered the road network by 1m and densified it so there is a node every few meters. You can see what I'm trying to achieve below. As you can see, I can get it to work in some cases:-

• blue is the densified road (a buffered line)
• red is the 'hexified' version - this is what I want to find
• the grey is the hex grid

I then used the new Snap geometries feature to snap the nodes to the closest hexagon corner. The results are promising, but there seem to be some edge cases where the line expands out to fill the hexagon (or part of it):-

The reason for the buffer is that Snap geometries doesn't let you snap to a layer whose geometry is different. For example, you can't snap nodes on a LINE layer to points on a POINT layer). It seems to be happiest snapping POLYGON to POLYGON.

I suspect the roads expand out when one side of the buffered road line jumps to one side of the hex cell, and the other side jumps to the other side of the hex cell. In my example, the roads which cross west-east at an acute angle seem to be the worst.

Things I've tried, without success:-

• buffering the road network by a tiny amount, so it remains a polygon but is very thin.
• densifying the hex cells (so there are nodes along the edges, not just at the corners)
• varying the maximum snapping distance (this has the largest effect, but I can't seem to find an ideal value)
• using LINE layers, not POLYGONs

I find that if I change to using just LINE layers, It works for a while, then crashes. It seems to save its work as it goes - some lines have been partially processed.

Does anyone know of any other way to snap points on a line to the nearest point on another line/polygon layer, ideally without needing to use postgres/postgis (although a solution with postgis would be welcome too)?

EDIT

For anyone who'd like to have a go, I've put a starter QGIS project here on Dropbox. This includes the Hex Grid and Densified lines layers. (The road network is from OSM, so can be downloaded using QuickOSM e.g. if you need to get the original to undensify the roads).

Note that it's in OSGB (epsg:27700) which a localised UTM for the UK, with units in meters.

Answer 1

My solution involves a PyQGIS script that is faster and more effective than a workflow involving snapping (I gave it a try too). Using my algorithm I've obtained these results:

You can run the following code snippets in sequence from within QGIS (in the QGIS Python console). At the end you get a memory layer with the snapped routes loaded into QGIS.

The only prerequisite is to create a multipart road Shapefile (use Processing->Singleparts to multipart, I used the field fictitiuos as Unique ID field parameter). This will give us a roads_multipart.shp file with a single feature.

Here is the algorithm explained:

1. Get the nearest hexagon sides where routes cross. For each hexagon we create 6 triangles between each pair of neighbour vertices and the corresponding centroid. If any road intersects a triangle, the segment shared by the hexagon and the triangle is added to the final snapped route. This is the heavier part of the whole algorithm, it takes 35 seconds running on my machine. In the first two lines there are 2 Shapefile paths, you should adjust them to fit your own file paths.

   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. Get rid of disconnected (or 'open') segments by using Python lists, tuples, and dictionaries. At this point, there are some disconnected segments left, i.e., segments that have one vertex disconnected but the other one connected to at least other 2 segments (see red segments in the next figure). We need to get rid of them.

   

   # 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. Now we can create a vector layer from the list of coordinates and load it to the QGIS map:

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

Another part of the result:

Should you need attributes in the snapped routes, we could use a Spatial Index to rapidly evaluate intersections (such as in https://gis.stackexchange.com/a/130440/4972 ), but that's another story.

Hope this helps!

Answer 2

I did it in ArcGIS, surely can be implemented using QGIS or simply python with package capable of reading geometries. Make sure that roads represent network, i.e. intersect each other at the ends only. You are dealing with OSM, I suppose it is the case.

• Convert proximity polygons to lines and planarise them, so they become a geometric network as well.
• Place points at their ends – Voronoi Points:
• Place points on the road at regular interval of 5 m, make sure network roads have good unique name:

• For every Road Point find coordinates of nearest Voronoi Point:
• Create “Roads” by connecting nearest points in the same order:

If you don't want to see this:

Don't try to use chainage points on Voronoi Lines. I afraid it will only make it worse. Thus your only option is creating network from Voronoi lines and find routes between road end points, that is no big deal either

Answer 3

I realize you're asking for a QGIS method, but bear with me for an arcpy answer:

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

Notes:

• This script contains many loops within loops and a nested cursor. There is definitely room for optimization. I ran through your datasets in a couple minutes, but more features will compound the issue.

Answer 4

If you were to split the road line into segments where each segment was completely contained by the hexagon, your decision on which hexagon line segments to use would be whether the distance from the split road segment's centroid to each hexagon side's midpoint was less than half the diameter of the hexagon (or less than the radius of a circle that fits inside the hexagon).

Thus, if you were to (one segment at a time) select hexagon line segments (where each segment is a side of the hexagon) that are within a distance of the radius of the hexagon, you could copy those line geometries and merge them on whatever unique identifier you use for your road dataset.

If you have trouble merging on the unique identifier, you could apply the buffer and select by location on only those segments to apply the attributes of your road dataset; that way you wouldn't have to worry about making false matches with a buffer that is too large.

The problem with the snap tool is that it snaps points indiscriminately; it's difficult to find that perfect tolerance to use. With this methodology you would be correctly identifying which hexagon line segments to use, then replacing the geometry of your road data (or inserting the geometries into a different dataset).

Also, if you still have the issue with the line segments that jump from one side of the hexagon to the other, you could split the line into segments by vertices, calculate the length of each line, then remove any line segments that are larger than the average length of one side of the hexagon.

Answer 5

The geometry snapper in qgis 3.0 has been reworked and now allows snapping between different geometry types. It also has lots of fixes. You could try a "daily snapshot" version to get access to the improved snapper before 3.0 is officially released.