QGIS – How to Snap a Road Network to a Hexagonal Grid in QGIS

geometryhexagonal gridpyqgisqgissnapping

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

enter image description here

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

enter image description here

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.

enter image description here

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.

Best Answer

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:

enter image description here

enter image description here

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.

    enter image description here

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

enter image description here

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!

Related Question