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?
Answer
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:
- Create a first spatial index with all the input points;
- Create a buffer around the current line using the specified distance;
- Find the nearest point to one of the line vertices and call it
firstpoint
; - Iterate over the features which intersect the buffer and insert them in a second spatial index;
- 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; - 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):
No comments:
Post a Comment