13

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.

  • 3
    Could you please share a sample dataset? I'd like to give it a try but don't want to go over the process of creating sample data from scratch. – Germán Carrillo Apr 13 '16 at 23:20
  • @GermánCarrillo - thanks. I've added a link to a sample project to the question. – Steven Kay Apr 14 '16 at 18:48
14
+100

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!

  • 1
    thank you, that works perfectly! Had problems pasting it into python console... I saved it as a .py file in the qgis python editor, and it ran fine from there. The multipart step removes the attributes, but a buffer/spatial join will fix that! – Steven Kay Apr 16 '16 at 19:08
  • 1
    Great! Glad it finally solved the problem you were facing. I'm interested in knowing what's the use case you're dealing with. Do you think we could leverage this to become a QGIS plugin or perhaps a script that is included into Processing scripts? – Germán Carrillo Apr 16 '16 at 21:16
  • 1
    the use case I had in mind was public transport maps like the Tube Map, where you need to snap lines to a tesselated grid, or to a restricted set of angles. This can be done manually by digitizing, but I was interested to see if it could be automated. I used hexes as they were easy to generate, visually interesting and had angles which weren't right angles. I think this is worth looking at in more detail, especially if it could be generalised to work with other tesselations... – Steven Kay Apr 17 '16 at 8:29
  • 1
    The idea behind the script would work on grids of triangles, squares, pentagons, hexagons, and so on. – Germán Carrillo Apr 19 '16 at 20:58
6

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: enter image description here
  • Place points on the road at regular interval of 5 m, make sure network roads have good unique name:

enter image description here enter image description here

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

If you don't want to see this: enter image description here

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

  • this is great, thank you! You mention using voronoi lines, not too familiar with that (Voronois from points, I can understand). Do you mean each line is surrounded by a polygon of all points closest to that line? (I'm not aware of a way of doing that in QGIS). Or do you mean the boundary lines from a normal voronoi mesh, based on points? – Steven Kay Apr 15 '16 at 20:21
  • Boundary lines of proximity polygons. Btw I stopped too early. To complete task it's enough to split 1st result at vertex, add point at the middle and repeat process – FelixIP Apr 15 '16 at 20:39
4

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

enter image description here

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.
  • Thank you for this, much appreciated. This shows exactly the effect I was visualising. The copious comments mean I can get the gist of what you're doing even if I can't run the code. Although it's arcpy, I'm sure this will be doable in pyqgis. The algorithm ideas in here are interesting (especially looking both clockwise and anticlockwise round each hex, and choosing the shortest way round) – Steven Kay Apr 15 '16 at 20:36
2

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.

1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.