Snapping a line over points?

I am looking for a solution that addresses this problem, but I have not been successful in my searches. After vectoring points on a raster, they are minimally aligned, as shown in the appended illustration. Then you get an initial line, the blue line, and you want this line to be adjusted according to solution 1 or solution 2, ignoring points more distant. I do not know if I could do this with PyQGIS or PostGIS. Has anyone had to provide similar solution can you point me a north? • Are you just wanting to snap the line to the closest points? – artwork21 Nov 28 '16 at 13:34
• Exactly, creating a vertex at each points. – Franklin Januário Nov 28 '16 at 14:43
• Have you tried ST_Snap()? You would need to select the points within a given distance of the line, and then create a MultiGeometry to pass to ST_Snap. – amball Nov 30 '16 at 6:03
• Hi, I try this: SELECT ST_Snap(lines.geom, ptos.geom, 0.50) AS linessnapped FROM ptos_simp ptos LEFT JOIN lin_simp lines ON ST_DWithin (ptos.geom, lines.geom, 0.50) WHERE lines.pk IS NOT NULL; but the line is replicate in each point. – Franklin Januário Dec 1 '16 at 11:45
• How many points are within 0.5 units? That depends on your projection. And what do you mean by "the line is replicate in each point?" It would also help to use ST_Collect(ptos.geom) to create a MultiGeometry to pass to ST_Snap. – amball Dec 1 '16 at 22:43

Solution 2

With PyQGIS, my solution is substantially based on the using of QgsSpatialIndex and of a distance (set by the user) used to find the surrounding points.

I summarize my 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.

This is my code:

##points=vector point
##line=vector line
##buffer_distance=number 30

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])
outGeom = QgsFeature()
if first:
seg_start = QgsPoint(firstpoint.asPoint())
seg_end = QgsPoint(all_points[nearest].geometry().asPoint())
outGeom.setGeometry(QgsGeometry.fromPolyline([seg_start, seg_end]))
seg_start = seg_end
first = False
else:
seg_end = QgsPoint(all_points[nearest].geometry().asPoint())
outGeom.setGeometry(QgsGeometry.fromPolyline([seg_start, seg_end]))
outGeom.setAttributes(line_attr)
seg_start = seg_end 