# Moving points onto lines (~neighborhood)

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

• There would be a PostGIS solution: PostGIS : nearest point on a linestring to a given point Jan 3 '13 at 23:23
• Are you looking for a totally automated function to do this, or would some sort of snapping tool to do it by hand be OK? Aug 25 '13 at 8:01
• I asked a similar question, I was trying to snap line to points but never found an easy solution. gis.stackexchange.com/questions/52232/… Sep 18 '13 at 13:42
• What about triangulation and distance matching? Jan 29 '14 at 21:14
• I found this question about a method that works in ArcGIS using Near. Went searching for QGIS Near equivalent and found this forum post where someone suggested GRASS v.distance. That lead me to this tutorial that may identify a method. Perhaps somewhere in there someone has written a plugin by now? Apr 26 '14 at 23:39

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
nearestIterator = lineLayer.getFeatures(QgsFeatureRequest().setFilterFid(featureId))
nearFeature = QgsFeature()
nearestIterator.nextFeature(nearFeature)
#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
snapGeometry = QgsGeometry.fromPoint(QgsPoint(closePoint,closePoint))

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

dist = math.hypot(p2.x() - p1.x(), p2.y() - p1.y())
print dist
``````
• running into nightmare trying to format this python code..argh!! Apr 29 '14 at 9:04

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') #iface.mapCanvas().layer(0)
lineLayer = QgsProject.instance().mapLayersByName('lines') # 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:
closePoint = closeSegResult
mindist = closeSegResult
side = closeSegResult
feat = QgsFeature()
feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(closePoint,closePoint)))
feat.setAttributes([i,mindist,side])
feats.append(feat)