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:
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()]] )
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 )
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!
ST_LineLocatePoint( geometry, ST_EndPoint( geometry)) don't have much sense, because is searching for the closest point of a line to the endpoint of the same line, and will always return 1.
And ST_LineInterpolatePoint( geometry, 1) will always return the endpoint.
Since you want to snap geometry vertices, use ST_Snap instead. Just find the distance threshold to search a vertex to snap.
In this example, I am snapping the geometries of the lines table to the vertices of the geometries in the smoothed_lines table, that are at most to 0.5 meters away.
Both tables have an id column, and a geometry column with the geometry( LINESTRING, 54004) type. I am casting the ST_Snap() output to the same type.
The inner join is returning only the candidate pairs, and from that pairs ST_Snap() is returning all lines geometries, modified or not.
DROP TABLE IF EXISTS snapped_lines;
CREATE TABLE IF NOT EXISTS snapped_lines AS (
SELECT
l.id AS id_l,
s.id AS id_s,
ST_Snap(
l.geometry,
s.geometry,
0.5
)::geometry( LineString, 54004) AS geometry
FROM
lines AS l
INNER JOIN
smoothed_lines AS s
ON
ST_DWithin( l.geometry, s.geometry, 0.5)
);
UPDATE:
If you need to snap lines within the same layer, you can perform the inner join to quickly exclude the pairs of lines that are far from each other, snap between each pair of selected lines, but filter the results for those lines that are more than zero apart from their pair (to avoid the lines that are generated by snapping to themselves or their adjacents).
One criterion to choose which lines must be snapped to which other, is the distance (<= 0.5) between them, adjust the value of 0.5 (reference system of the geometry column units) to the one that fits your needs.
When the lines from and to snapping, belong to the same layer, the same criterion that is applied to snapping line 1 to line 2, is applied to snapping line 2 to line 1. So let's add one more criterion than be deterministic to decide which line should snap to the other, for example, that the shortest line must snaps to the longest one.
DROP TABLE IF EXISTS snapped_lines;
CREATE TABLE IF NOT EXISTS snapped_lines AS (
SELECT
l.id AS id_l,
s.id AS id_s,
ST_Snap(
l.geometry,
s.geometry,
0.5
)::geometry( LineString, 54004) AS geometry
FROM
smoothed_lines AS l
INNER JOIN
smoothed_lines AS s
ON
ST_DWithin( l.geometry, s.geometry, 0.5)
WHERE
ST_Distance( l.geometry, s.geometry) > 0
AND
ST_Length( l.geometry) < ST_Length( s.geometry)
);
Best Answer
The behaviour of this tool is a bit tricky, indeed. So instead of snapping, better create the closest point on the line layer using QGIS expressions with Geometry generator or Geometry by Expression (see here for details):