Sunday, 13 November 2016

qgis - 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?


Illustrative scheme of the solution sought



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:



  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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...