Snapping lines' start_point/end_points to the closest point features to them

Question

I'm trying to make a script/function as part of a QGIS plugin to snap all the lines to their nearest point.

I have a project with 1 LineLayer and 1 PointLayer. Each line feature should be snapped to the nearest point feature from PointLayer. There will be no instance where more than 1 point feature will be in close distance to another.

Here's an example of what it may look like prior to snapping. Using the "Snap geometries to layer" from the processing toolbox I obtain a good result such as this:

Now the main issue here is that is creates a virtual(?) or temporary layer, I need it to edit the actual LineLayer, I'm also trying to compact useful scripts into one major Plugin (or maybe an individual plugin) to make it easier for frequent use for others, which is why I'm trying to write the script myself.

I have searched around and most solutions revolve are made with/for ArcGIS/PostGIS or with Sqlite queries as opposed to basic PyQGIS function (which is what I need). I'm not sure where to look for guides or resources on how to achieve this. I know how to retrieve the concerned layers and edit the needed layer(LineLayer) but not entirely sure how to edit the geometry.

I can retrieve the start and end points of each line(Can vector layer get start point and end point of line using PyQGIS?)

start_point = QgsPoint(geom[0])
end_point = QgsPoint(geom[-1])

I would then need to add a buffer to them and snap the start_point and/or end_point to the point feature closest to it, I assume through getting the coordinates of any point feature overlapping with the buffer. Point features are placed manually and in a way that each line feature should have a point feature snapped to each of its ends.

edit: Following @MrXsquared solution and this thread First and last point of MultiLineString-objects in QGIS 3 I have tried the following code:

with edit(lines):
    # enumerate because we need the feature id to use QgsVectorLayerEditUtils
    for i, feat in enumerate(lines.getFeatures()):
        multiLineString = feat.geometry().constGet()
        first_part = multiLineString.geometryN(0)
        i += 1 # adjust the current feature id
    first_part = multiLineString.geometryN(0)#QgsLineString

    start_point = first_part.pointN(0)#QgsPoint
    end_point = first_part.pointN(first_part.numPoints()-1)#QgsPoint

    # find nearest point of point layer to start and endpoint of line
    nearestneighbor_start = points_idx.nearestNeighbor(QgsPointXY(start_point), neighbors=1, maxDistance=tolerance)
    nearestneighbor_end = points_idx.nearestNeighbor(QgsPointXY(end_point), neighbors=1, maxDistance=tolerance)

    try: # only adjust the startpoint if there is a nearest point within the maxdistance
        nearpoint_startgeom = points.getFeature(nearestneighbor_start[0]).geometry() # get the geometry of the nearest point in pointlayer
        x_start = nearpoint_startgeom.asPoint().x() # extract x value of that point
        y_start = nearpoint_startgeom.asPoint().y() # extract y value of that point
        QgsVectorLayerEditUtils(lines).moveVertex(x_start,y_start,i,0) # use QgsVectorLayerEditUtils to edit the vertex, for more details see https://qgis.org/pyqgis/3.4/core/QgsVectorLayerEditUtils.html#qgis.core.QgsVectorLayerEditUtils
    except: # if there is no nearest point just do nothing
        pass
    try:
        nearpoint_endgeom = points.getFeature(nearestneighbor_end[0]).geometry()
        x_end = nearpoint_endgeom.asPoint().x()
        y_end = nearpoint_endgeom.asPoint().y()
        lastvert = len(feat.geometry().asPolyline())-1 # QgsVectorLayerEditUtils.moveVertex() seems to not accept -1 as index of last vertex, so we need to count for it...
        QgsVectorLayerEditUtils(lines).moveVertex(x_end,y_end,i,lastvert)
    except:
        pass

As start_point and end_point seem to be of type QgsPoint, I was getting an error for the first argument of points_idx.nearestNeighbor() (https://www.qgis.org/api/classQgsSpatialIndex.html#a38804ff7021561070c109abe91342210) so I tried converting them using QgsPointXY(point) and at that point I was getting no error. The code iterated through all the points without but nothing was happening. No nearest point was being detected and x_start/x_end/y_start/y_end were returning nothing if I tried printing them.

Printing nearpoint_startgeom would return something like this:

<QgsGeometry: MultiPoint ((1849531.11318555101752281 1087547.47298211278393865))>
<QgsGeometry: MultiPoint ((1849461.63486740970984101 1087509.80619045789353549))>
<QgsGeometry: MultiPoint ((1849386.41441025864332914 1087465.64150428189896047))>
<QgsGeometry: MultiPoint ((1849337.62416943279094994 1087424.41622101375833154))>
<QgsGeometry: MultiPoint ((1849280.04536620667204261 1087376.88177968817763031))>

As of right now, I am not entirely sure whether it's due to the tolerance variable or something else.

edit2(3?):

I tried converting from multi-type to single type using ConvertToSingleType(), it also execute without returning any error but still does no change to the actual geometry.

with edit(lines):
    # enumerate because we need the feature id to use QgsVectorLayerEditUtils
    for i, feat in enumerate(lines.getFeatures()):
        linegeom = feat.geometry()
        linegeom.convertToSingleType()
        convertedLineGeom = linegeom.asPolyline()
        i += 1 # adjust the current feature id
    # get start and end point of line
    line_startgeom = QgsPointXY(convertedLineGeom[0])
    line_endgeom = QgsPointXY(convertedLineGeom[-1])

    # find nearest point of point layer to start and endpoint of line
    nearestneighbor_start = points_idx.nearestNeighbor(line_startgeom, neighbors=1, maxDistance=tolerance)
    nearestneighbor_end = points_idx.nearestNeighbor(line_endgeom, neighbors=1, maxDistance=tolerance)

Conversion from multi to single part returns true for all elements. Printing multilinestrings .geometry() returns this:

<QgsGeometry: MultiLineString ((1850094.12363056233152747 1085521.28121069446206093, 1850077.50219331961125135 1085563.30667054420337081))>

Priting convertedLineGeom returns

[<QgsPointXY: POINT(1849069.62466584099456668 1086983.79451146279461682)>, <QgsPointXY: POINT(1849106.36500586662441492 1087043.8564368411898613)>]

Which then makes line_startgeom and line_endgeom both receive the respective start and end points of each line. An example of printing line_startgeom:

<QgsPointXY: POINT(1850829.4128553019836545 1082368.52812290168367326)>

I still fail to understand what's at fault.

edit3(4?): The point layer is a Point(MultiPoint) type. Maybe this is what causes an issue?

Ok. Converting the nearpoint_startgeom and nearpoint_endgeom does make the script actually edit geometry so it was indeed due to the point layer geometry type.

I did it this way: (same for nearpoint_endgeom)

nearpoint_startgeom = points.getFeature(nearestneighbor_start[0]).geometry() 
nearpoint_startgeom.convertToSingleType()

The issue here is that the result is kind of messy.

For example, here is a screenshot with the same chunk of the map I used initially, before script execution:

After its execution:

Here's the code looks like:

points = QgsProject.instance().mapLayersByName('PointLayer')[0]
lines = QgsProject.instance().mapLayersByName('LineLayer')[0]
tolerance = 10000 # snapping tolerance in CRS units
create the spatial index of the pointlayer
points_idx = QgsSpatialIndex(points.getFeatures())
with edit(lines):
    # enumerate because we need the feature id to use QgsVectorLayerEditUtils
    for i, feat in enumerate(lines.getFeatures()):
        linegeom = feat.geometry()
        linegeom.convertToSingleType()
        convertedLineGeom = linegeom.asPolyline()
        i += 1 # adjust the current feature id
    # get start and end point of line
    line_startgeom = QgsPointXY(convertedLineGeom[0])
    line_endgeom = QgsPointXY(convertedLineGeom[-1])

    # find nearest point of point layer to start and endpoint of line
    nearestneighbor_start = points_idx.nearestNeighbor(line_startgeom, neighbors=1, maxDistance=tolerance)
    nearestneighbor_end = points_idx.nearestNeighbor(line_endgeom, neighbors=1, maxDistance=tolerance)

    try: # only adjust the startpoint if there is a nearest point within the maxdistance
        nearpoint_startgeom = points.getFeature(nearestneighbor_start[0]).geometry() # get the geometry of the nearest point in pointlayer
        nearpoint_startgeom.convertToSingleType()
        x_start = nearpoint_startgeom.asPoint().x() # extract x value of that point
        y_start = nearpoint_startgeom.asPoint().y() # extract y value of that point
        QgsVectorLayerEditUtils(lines).moveVertex(x_start,y_start,i,0) # use QgsVectorLayerEditUtils to edit the vertex, for more details see https://qgis.org/pyqgis/3.4/core/QgsVectorLayerEditUtils.html#qgis.core.QgsVectorLayerEditUtils
    except: # if there is no nearest point just do nothing
        print("detects nothing for start point")
    try:
        nearpoint_endgeom = points.getFeature(nearestneighbor_end[0]).geometry()
        nearpoint_endgeom.convertToSingleType()
        x_end = nearpoint_endgeom.asPoint().x()
        y_end = nearpoint_endgeom.asPoint().y()
        QgsVectorLayerEditUtils(lines).moveVertex(x_end,y_end,i,-1)
    except:
        print("detects nothing for end point")

I tried meddling with the tolerance value from anywhere between .15 to 20000.

Answer 1

Here is a solution for MultiLineStrings and MultiPoints. It moves the start- and endpoints of every part of the MultiLineString to the closest point.

In case you do not want to move start- and enpoints of every part, but only of every feature, you may want to change these two lines accordingly (simply comment/uncomment in the code below):

• if part.startPoint() == geom.vertexAt(vert): to if feat.geometry().vertexAt(0) == geom.vertexAt(vert):
• elif part.endPoint() == geom.vertexAt(vert): to elif feat.geometry().vertexAt([i for i, f in enumerate(feat.geometry().vertices())][-1]) == geom.vertexAt(vert): (not very convenient, but simply doing feat.geometry().vertexAt(-1) directly doesnt work...)

Here is the full code with further explanations as comments:

points = QgsProject.instance().mapLayersByName('points')[0]
lines = QgsProject.instance().mapLayersByName('lines')[0]
tolerance = 10000 # snapping tolerance in CRS units
if QgsWkbTypes.isMultiType(points.wkbType()): # if point-input is of type multigeometry
    res = processing.run("native:multiparttosingleparts",{'INPUT':points,'OUTPUT':'TEMPORARY_OUTPUT'}) # convert to singlepoint
    singlepoints = res['OUTPUT']
else: # if point-input is of type singlegeometry 
    singlepoints = points # leave the inputlayer as it is
create the spatial index of the pointlayer
points_idx = QgsSpatialIndex(singlepoints.getFeatures())
with edit(lines):
    for feat in lines.getFeatures():
        geom = feat.geometry()
        vert = 0 # reset vertices-count on every new feature
        vert_dict = {} # clear vertices dictionary on every new feature
        for part in geom.parts(): # iterate through all parts of each feature
            for vertex in part.vertices(): # iterate through all vertices of each part
                # if feat.geometry().vertexAt(0) == geom.vertexAt(vert): # if the current vertex is a startpoint of the feature
                if part.startPoint() == geom.vertexAt(vert): # if the current vertex is a startpoint of the features part, then:
                    nearestneighbor_start = points_idx.nearestNeighbor(QgsPointXY(part.startPoint()), neighbors=1, maxDistance=tolerance) # find the closest point to the current startpoint
                    if len(nearestneighbor_start) > 0:# only adjust the startpoint if there is a nearest point within the maxdistance
                        nearpoint_startgeom = singlepoints.getFeature(nearestneighbor_start[0]).geometry() # get the geometry of the nearest point in pointlayer
                        vert_dict[vert] = nearpoint_startgeom.asPoint() # add the index of the current vertex as key and the geometry of the closest point as value to the dictionary
                # elif feat.geometry().vertexAt([i for i, f in enumerate(feat.geometry().vertices())][-1]) == geom.vertexAt(vert): # if the current vertex is an endpoint of the current feature
                elif part.endPoint() == geom.vertexAt(vert): # if the current vertex is an endpoint of the current features part
                    nearestneighbor_end = points_idx.nearestNeighbor(QgsPointXY(part.endPoint()), neighbors=1, maxDistance=tolerance)
                    if len(nearestneighbor_end) > 0:
                        nearpoint_endgeom = singlepoints.getFeature(nearestneighbor_end[0]).geometry()
                        vert_dict[vert] = nearpoint_endgeom.asPoint()
                else:
                    pass # if current vertex is not a start or end point skip it...
                vert += 1 # increase vertices-counter
        for vertindex, newpoint in vert_dict.items(): # for every feature iterate over the just created dictionary (vertindex (=dict key) is the start or endpoint we want to move and newpoint (=dict value) the position we want to move it to)
            QgsVectorLayerEditUtils(lines).moveVertex(newpoint.x(),newpoint.y(),feat.id(),vertindex) # use QgsVectorLayerEditUtils to edit the vertex, for more details see https://qgis.org/pyqgis/3.4/core/QgsVectorLayerEditUtils.html#qgis.core.QgsVectorLayerEditUtils
https://gis.stackexchange.com/questions/410768/snapping-lines-start-point-end-points-to-the-closest-point-features-to-them/410778#410778

For further background you may want to read Modifying specific vertices of MultiLineString using PyQGIS

Example for snapping start- and endpoints of all parts of all feature:

Example for snapping start- and enpoints only for all features:

Answer 2

This solution is for LineStrings and Points, for MultiLineStrings and MultiPoints see my other answer.

You can use QgsSpatialIndex(fi: QgsFeatureIterator, feedback: QgsFeedback = None).nearestNeighbor(self, point: QgsPointXY, neighbors: int) to locate the nearest point to your start- and endpoints. Then get the feature of the pointlayer by using the featureid from the spatial index (pointlayer.getFeature(featureIdIntFromSpatialIndex), read its geometry and use QgsVectorLayerEditUtils(layer: QgsVectorLayer).moveVertex(self, x: float, y: float, atFeatureId: int, atVertex: int) to move the vertex of the linelayers feature.

Here is an example with comments as explanation:

points = QgsProject.instance().mapLayersByName('points')[0]
lines = QgsProject.instance().mapLayersByName('lines')[0]
tolerance = 10000 # snapping tolerance in CRS units
create the spatial index of the pointlayer
points_idx = QgsSpatialIndex(points.getFeatures())
with edit(lines):
    # enumerate because we need the feature id to use QgsVectorLayerEditUtils
    for i, feat in enumerate(lines.getFeatures()):
        linegeom = feat.geometry().asPolyline()
        i += 1 # adjust the current feature id
    # get start and end point of line
    line_startgeom = QgsPointXY(linegeom[0])
    line_endgeom = QgsPointXY(linegeom[-1])

    # find nearest point of point layer to start and endpoint of line
    nearestneighbor_start = points_idx.nearestNeighbor(line_startgeom, neighbors=1, maxDistance=tolerance)
    nearestneighbor_end = points_idx.nearestNeighbor(line_endgeom, neighbors=1, maxDistance=tolerance)

    try: # only adjust the startpoint if there is a nearest point within the maxdistance
        nearpoint_startgeom = points.getFeature(nearestneighbor_start[0]).geometry() # get the geometry of the nearest point in pointlayer
        x_start = nearpoint_startgeom.asPoint().x() # extract x value of that point
        y_start = nearpoint_startgeom.asPoint().y() # extract y value of that point
        QgsVectorLayerEditUtils(lines).moveVertex(x_start,y_start,i,0) # use QgsVectorLayerEditUtils to edit the vertex, for more details see https://qgis.org/pyqgis/3.4/core/QgsVectorLayerEditUtils.html#qgis.core.QgsVectorLayerEditUtils
    except: # if there is no nearest point just do nothing
        pass
    try:
        nearpoint_endgeom = points.getFeature(nearestneighbor_end[0]).geometry()
        x_end = nearpoint_endgeom.asPoint().x()
        y_end = nearpoint_endgeom.asPoint().y()
        lastvert = len(feat.geometry().asPolyline())-1 # QgsVectorLayerEditUtils.moveVertex() seems to not accept -1 as index of last vertex, so we need to count for it...
        QgsVectorLayerEditUtils(lines).moveVertex(x_end,y_end,i,lastvert)
    except:
        pass