5

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?

Illustrative scheme of the solution sought

  • 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
6
+100

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[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 20 m, the result is the following (I have also inserted the buffer created for more clearness, but it's not an output):

enter image description here

  • This is a great script! Worked very well for me, going to award you a bounty for this. Hope you don't mind if I edit your script ever so slightly to restrict the vector input layers only to point and line types ;) – Joseph Dec 9 '16 at 10:51
  • Out of curiosity, what was Solution 1? =) – Joseph Dec 9 '16 at 10:53
  • 1
    Ok, I got it. Sincerely, I don't know why this happens because the first point on the next line should be over the buffer distance you initially set (I can't reproduce your problem, so I can't edit the code). A workaround for this problem could be to select the points near a line, then save the selection ("Save selected feature" tool from Processing Toolbox) and the running my code on this selection. If you don't like this workaround, I'm afraid you need to post a new question about it. – mgri Jan 23 '17 at 14:06
  • 1
    Instead, referring to your last point, sorry but it is an evident off-topic request, so please post a new question if it this the case. However, I don't understand what you mean with "lines are break on points": merging or dissolving the lines should be a simple task using the Processing Toolbox. – mgri Jan 23 '17 at 14:08
  • 1
    Thanks a lot, I goto post a new topic about this. But, by the way, I mean merge between each row independently, but I might be able to do this using the Grass toolbox >v.clean> snap. – Franklin Januário Jan 23 '17 at 16:00

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.