How to snap lines to points in QGIS using Python

This task is based on my answer to this question poped up on gis.stackexchange.com.

We start from one point vector layer and one line vector layer: the goal for this task is to literally “snap” the line vector layer to the nearest points on the basis of a searching distance (i.e. we want to snap the line to the points which are within a specified distance from it). The result is a memory vector layer which will store the attributes of the input line together with the new geometry.

input

Two problems rise up:

  • How to determine and locate the point from which to start re-drawing the line?
  • How to establish a rule for drawing the line following a non-random path?

The answer to both questions is the using of QgsSpatialIndex.

In more detail, let summarize the workflow:

  1. Create a first spatial index with all the input points;
  2. Create a buffer around the current line using the specified distance;
  3. Find the nearest point to one of the line vertices and call it firstpoint;
  4. Iterate over the features which intersect the buffer and insert them in a second spatial index;
  5. Iterate over the the second spatial index and find the nearest distance from firstpoint. Once the nearest point is set, remove it from the second spatial index;
  6. Create a memory layer and store the attributes of the input line, together with the new geometry.

Here is the code (if you don’t know how to create a new script from Processing, please see here):

##points=vector
##line=vector
##buffer_distance=number 20

from qgis.core import *

layer1 = processing.getObject(points)
layer2 = processing.getObject(line)
crs = layer2.crs().toWkt()

# Create the output layer
outLayer = QgsVectorLayer('Linestring?crs='+ crs, 'snapped' , 'memory')
prov = outLayer.dataProvider()
fields = layer2.pendingFields() # Fields from the input layer
prov.addAttributes(fields) # Add input layer fields to the outLayer
outLayer.updateFields()

all_points = {} # Dictionary containing all the input points
index1 = QgsSpatialIndex() # First spatial index
for feat in layer1.getFeatures():
    index1.insertFeature(feat)
    all_points[feat.id()] = feat

index2 = QgsSpatialIndex() # Second spatial index

first = True
for ft in layer2.getFeatures():
    count = 0
    line_attr = ft.attributes()
    line_geom = ft.geometry()
    bf_geom = line_geom.buffer(buffer_distance, -1) # Buffer using the distance set
    firstpoint = line_geom.interpolate(0) # Start point from the input line
    idsList = index1.intersects(bf_geom.boundingBox())
    for id in idsList:
        inGeom = all_points[id].geometry()
        if bf_geom.intersects(inGeom): # Check if the feature is within the buffer
            index2.insertFeature(all_points[id])
            count += 1
    for num in xrange(0, count):
        nearest = index2.nearestNeighbor(firstpoint.asPoint(), 1)
        index2.deleteFeature(all_points[nearest[0]])
        outGeom = QgsFeature()
        if first:
            seg_start = QgsPoint(firstpoint.asPoint())
            seg_end = QgsPoint(all_points[nearest[0]].geometry().asPoint())
            outGeom.setGeometry(QgsGeometry.fromPolyline([seg_start, seg_end]))
            seg_start = seg_end
            first = False
        else:
            seg_end = QgsPoint(all_points[nearest[0]].geometry().asPoint())
            outGeom.setGeometry(QgsGeometry.fromPolyline([seg_start, seg_end]))
            outGeom.setAttributes(line_attr)
            seg_start = seg_end
            prov.addFeatures([outGeom])

# Add the layer to the Layers panel
QgsMapLayerRegistry.instance().addMapLayer(outLayer)

For example, setting a searching distance of 60 m, the result is the following (I have also inserted the buffer created for more clearness, but it’s not an output):

snapped

Advertisements

2 thoughts on “How to snap lines to points in QGIS using Python

  1. This is was I looking for, it’s very useful. However, for some reason it doesn’t work for me. I’m using QGIS 2.18.11 and I have a point layer and a line layer. I set the threshold to 5 meters, but the lines are snapped to points separated a distance of 20 m. The gaps between some points and lines are sometimes 0.5 m. I am using a custom coordinate system, would that be a problem?

    Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s