2

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.

enter image description here

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:

enter image description here

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:

enter image description here

After its execution:

enter image description here

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.

8
  • 1
    You can use a spatial index to efficiently locate the nearest point, then modify the end/startgeometry of the line to that nearest point.
    – MrXsquared
    Sep 2 at 12:44
  • What prevents you from saving the output of the snap-tool to a file and working with this data?
    – Erik
    Sep 2 at 12:48
  • 1
    @Erik I would much rather save the changes to the exact same .shp file I'm initially working with. This is supposed to be used by various users, the less meddling with files the better Sep 2 at 12:53
  • @MrXsquared Could you please explain to me how to use spatial indices ? I'm fairly now to Python and PyQGIS in general Sep 2 at 12:57
  • 1
    @Erik As of right now, we get projects like that from another software which doesn't allow snapping so the tool/script has to be ran until things change. Until then, I do have to make a script that takes cares of it and I would like to make it part of a bigger plugin so its done from a simple button click without extra work from other users. Sep 2 at 13:17
2

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:

enter image description here

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

enter image description here

4
  • It does work just about perfect, thanks a lot with a thoroughly commented/explained script. I'm a tad curious though, doesn't the processing script create a new temporary layer with updated type (since I did use it with my multi-type layers) ? Also, as of right now since all of my layers have their name end with the project's code, I'm using title = QgsProject.instance().title() points = QgsProject.instance().mapLayersByName('Points '+title)[0] But the title variable prints empty while it did work for other projects. I'm wondering what causes the issue and what's a workaround Sep 13 at 12:02
  • I have asked a little while ago how to retrieve a layer's name through a part of its name, you have provided me with a solution but as I try the same method in a different project, I get an error that 'layername' is referenced before assignment, exact same code that works for another project. I'm genuinely confused. I cannot reply to that other thread to update my question as I just created a legit account and can't comment on the other one. Should I create another thread over the same question ? Sep 13 at 12:49
  • @IKindaNeedAHand yes, please create two new separate questions for these topics. If you need to provide context, just include a link to the previous questions. Background is, that we use a "focused Q/A model" on GIS SE, where a Question should only contain one single question and only answers to this specific question. So in contrast to other forums we encourage to open multiple threads. About your second account you can read up here: meta.stackexchange.com/questions/18232/…
    – MrXsquared
    Sep 13 at 13:08
  • Eventually, all I had to do was retrieve the file name using .baseName() using the project's file path. Everything works, thanks a lot for the help as well as the explanations. Sep 13 at 13:34
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

enter image description here

11
  • Sorry for a late response and thanks. As I ran the script, I got the intitial error: TypeError: MultiLineString geometry cannot be converted to a polyline. Only single line or curve types are permitted. So I tried replacing asPolyline with asMultiPolyline but then I got another error as follows: TypeError: QgsPointXY(): arguments did not match any overloaded call: overload 1: too many arguments overload 2: argument 1 has unexpected type 'list' (same until overload 6) line_startgeom = QgsPointXY(linegeom[0]) (line of error) Sep 6 at 8:57
  • Ok yes, MultiLineStrings need to be handled a little different. How should they be treated in terms of Start- and Endpoints? Should each part of that MultiLine have a Start- and Endpoint or only each MultiLine overall one Start- and Endpoint? Can take a look into this when I have time again. Meanwhile, if thats possible / makes sense for your data, you can try to convert your MultiLineStrings into SingleLines.
    – MrXsquared
    Sep 6 at 9:07
  • 1
    I'm currently looking at this thread that seems to explain the process to retrieve the start and end points out of a MultiLineString gis.stackexchange.com/questions/304755/… I will try to make something out of it and if I find a way, I'll add a comment with the edited script. Sep 6 at 9:13
  • Sir, I have edited my initial post with more details over what I have attempted. It seems that x_start/x_end/y_start/y_end return nothing Sep 6 at 10:16
  • Second attempt added to my initial post with ConvertToSingleType() (qgis.org/pyqgis/3.6/core/…) still executes but no result Sep 6 at 13:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.