Moving points onto lines (~neighborhood)

Question

I do have two vector layer, of which one is a point layer based on "events" by remote sensing and the second one is a line layer from local research.

In my case these are earthquakes and tectonic faults, but I guess one could simply choose "car-accidents and roads" as a general example.

So what I'd like to do is move/copy the points onto the closest point of the lines, as long as its within a tolerance distance (say 1-2km or 0.0xx°), with the new point layer (+attr moved y/n).

Any ideas ?

Linux, QGIS 1.8

Answer 1

Posted a code snippet(tested in python console) that doest the below

1. Use QgsSpatialIndex to find the nearest line feature to a point
2. Find the nearest point on this line to the point. I used shapely package as a shortcut for this. I found the QGis methods for this as insufficient(or most probably i do not understand them properly)
3. Added rubberbands to the snap locations

from shapely.wkt import *
from shapely.geometry import *
from qgis.gui import *
from PyQt4.QtCore import Qt
lineLayer = iface.mapCanvas().layer(0)
pointLayer =  iface.mapCanvas().layer(1)
canvas =  iface.mapCanvas()
spIndex = QgsSpatialIndex() #create spatial index object
lineIter =  lineLayer.getFeatures()
for lineFeature in lineIter:
    spIndex.insertFeature(lineFeature)        
pointIter =  pointLayer.getFeatures()
for feature in pointIter:
    ptGeom = feature.geometry()
    pt = feature.geometry().asPoint()
    nearestIds = spIndex.nearestNeighbor(pt,1) # we need only one neighbour
    featureId = nearestIds[0]
    nearestIterator = lineLayer.getFeatures(QgsFeatureRequest().setFilterFid(featureId))
    nearFeature = QgsFeature()
    nearestIterator.nextFeature(nearFeature)
    shplyLineString = shapely.wkt.loads(nearFeature.geometry().exportToWkt())
    shplyPoint = shapely.wkt.loads(ptGeom.exportToWkt())
    #nearest distance from point to line
    dist = shplyLineString.distance(shplyPoint)
    print dist
    #the point on the road where the point should snap
    shplySnapPoint = shplyLineString.interpolate(shplyLineString.project(shplyPoint))
    #add rubber bands to the new points
    snapGeometry = QgsGeometry.fromWkt(shapely.wkt.dumps(shplySnapPoint))
    r = QgsRubberBand(canvas,QGis.Point)
    r.setColor(Qt.red)
    r.setToGeometry(snapGeometry,pointLayer)

Edit: Just now found that @radouxju method using closestSegmentWithContext gives the same results in less lines of code. I wonder why they came up with this weird method name? should have been something like closestPointOnGeometry.

So we can avoid shapely and do like,

nearFeature = QgsFeature()
nearestIterator.nextFeature(nearFeature)   

closeSegResult = nearFeature.geometry().closestSegmentWithContext(ptGeom.asPoint())
closePoint = closeSegResult[1]
snapGeometry = QgsGeometry.fromPoint(QgsPoint(closePoint[0],closePoint[1])) 

p1 = ptGeom.asPoint()
p2 = snapGeometry.asPoint()

dist = math.hypot(p2.x() - p1.x(), p2.y() - p1.y())
print dist

Answer 2

here is a pseudo-code to start with. I hope this helps and that someone will have time to provide full code (I don't have at the moment)

first thing to do is to loop on the point and select the lines that are located within the threshold distance to each point. Thi can be done with QgsSpatialIndex

Within the first loop, the second thing to do is to loop on the selected lines and find the closest point on the line. This can be done directly based on QgsGeometry::closestSegmentWithContext

double QgsGeometry::closestSegmentWithContext ( const QgsPoint & point, QgsPoint & minDistPoint, int & afterVertex, double * leftOf = 0, double epsilon = DEFAULT_SEGMENT_EPSILON )

Searches for the closest segment of geometry to the given point.

Parameters point Specifies the point for search

minDistPoint  Receives the nearest point on the segment

afterVertex   Receives index of the vertex after the closest segment. The vertex before the closest segment is always afterVertex -

1 leftOf Out: Returns if the point lies on the left of right side of the segment ( < 0 means left, > 0 means right ) epsilon epsilon for segment snapping (added in 1.8)

the third step (within the first loop) would consist in updating the geometry of the point with the geometry of the minDistPoint with the smallest distance

update with some code (on QGIS3)

pointlayer = QgsProject.instance().mapLayersByName('point')[0] #iface.mapCanvas().layer(0)
lineLayer = QgsProject.instance().mapLayersByName('lines')[0] # iface.mapCanvas().layer(1)

epsg = pointlayer.crs().postgisSrid()
uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer&field=distance:double(20,2)&field=left:integer&index=yes"
snapped = QgsVectorLayer(uri,'snapped', 'memory')

prov = snapped.dataProvider()

testIndex = QgsSpatialIndex(lineLayer)
i=0

feats=[]

for p in pointlayer.getFeatures():
    i+=1
    mindist = 10000.
    near_ids = testIndex.nearestNeighbor(p.geometry().asPoint(),4) #nearest neighbor works with bounding boxes, so I need to take more than one closest results and further check all of them. 
    features = lineLayer.getFeatures(QgsFeatureRequest().setFilterFids(near_ids))
    for tline in features:
        closeSegResult = tline.geometry().closestSegmentWithContext(p.geometry().asPoint())
        if mindist > closeSegResult[0]:
            closePoint = closeSegResult[1]
            mindist = closeSegResult[0]
            side = closeSegResult[3]
    feat = QgsFeature()
    feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(closePoint[0],closePoint[1])))
    feat.setAttributes([i,mindist,side])
    feats.append(feat)

prov.addFeatures(feats)
QgsProject.instance().addMapLayer(snapped)