4

First of, this question is related to this one: Snapping lines' start_point/end_points to the closest point features to them. I am trying to modify my answer to make it work with MultiLineStrings.

Basically I have two inputs:

  1. MultiLineString Layer
  2. MultiPoint Layer

I am now trying to snap the start and end points of each part of the MultiLines to the closest Point. Therefore the MultiPoints are converted to SinglePoints and the nearest points geometry is beeing read. Now the issue is, that it seems like I am not able to modify the Parts of the MultiLineString by using

  • part.setXAt()
  • part.setYAt()
  • part.moveVertex()

Link to the QGIS-Docs. Debugging all this using print() statements actually return the correct results. Start and endpoints seem to be modified and the correctly modified LineString is printed.

However, there are no changes made in the actual geometries. The input lines just stay as they were before running the script, no error, no nothing. So my question is: why doesnt my script commit the changes? And in case setXAt() and setYAt() are not meant to change things, what are they for and what could I use instead?

Here is the entire code:

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:
    singlepoints = points


creat the spatial index of the pointlayer


points_idx = QgsSpatialIndex(singlepoints.getFeatures())


with edit(lines):
    for i, feat in enumerate(lines.getFeatures()):
        for j, part in enumerate(geom.parts()):
            # get start and end point of line
            line_startgeom = QgsPointXY(part.startPoint())
            line_endgeom = QgsPointXY(part.endPoint())
            print('\n')
            print('Feat: ' + str(i) + ' Part: ' + str(j) + ' OriginalLine: ' + str(part))


        # 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)

        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
            part.setXAt(0, nearpoint_startgeom.asPoint().x())
            part.setYAt(0, nearpoint_startgeom.asPoint().y())

        if len(nearestneighbor_end) > 0:
            nearpoint_endgeom = singlepoints.getFeature(nearestneighbor_end[0]).geometry()
            part.setXAt(-1, nearpoint_endgeom.asPoint().x())
            part.setYAt(-1, nearpoint_endgeom.asPoint().y())
            #part.moveVertex(QgsVertexId(-1,-1,-1),QgsPoint(nearpoint_endgeom.asPoint())) # honestly, QgsVertexId is not well documented, so I have no idea how this is supposed to work...
            #print(str(part.xAt(-1))+ '\n' + str(nearpoint_endgeom.asPoint())+'\n')

        print('Feat: ' + str(i) + ' Part: ' + str(j) + ' ModifiedLine: ' + str(part))
        lines.updateFeature(feat)



lines.commitChanges()

And an example output from the print statement:

Feat: 2 Part: 3 OriginalLine: <QgsLineString: LineString (689753.69139794236980379 5338632.41991709917783737, 697525.0708597571356222 5339841.7065342990681529, 695706.26743028790224344 5337553.34959724545478821)>
Feat: 2 Part: 3 ModifiedLine: <QgsLineString: LineString (688633.31491136807017028 5338762.2247987687587738, 697525.0708597571356222 5339841.7065342990681529, 695385.03714443475473672 5337152.7669952092692256)>

As you can see it prints different line strings, but does not change geometries.

MrXsquared
  • 34,292
  • 21
  • 67
  • 117

1 Answers1

6

About why it does not modify the actual vertices, here a quote from the docs --> its a read only copy:

vertices(self) → QgsVertexIterator

Returns a read-only, Java-style iterator for traversal of vertices of all the geometry, including all geometry parts and rings.

The iterator returns a copy of individual vertices, and accordingly geometries cannot be modified using the iterator. See transformVertices() for a safe method to modify vertices “in-place”.

I guess, even thoug I do not use vertices(), part.setXAt(), part.setYAt() and part.moveVertex() apply to it.

What you can do instead:

I guess there is a more intuitive way to do this, but for now I am stepping back and make use of QgsVectorLayerEditUtils(layer).moveVertex() again. Luckily Python has dictionarys which can be used for a lot of stuff, just like this:

So basically I am iterating over each vertex of each part of each feature and check wheter the current vertex is a start- or an endpoint of the current part. If so, I'll search for the nearest point, read its geometry and store this geometry as value in the dictionary and the index of the current vertex as key. So I get a dictionary of vertexindizes I want to move (keys) and points where these vertices should be moved to (values). When iterating over this dictionary, I can go back using QgsVectorLayerEditUtils(layer).moveVertex(value.x(),value.y(),feat.id(),key) instead of the not working part.moveVertex(), part.setXAt() or part.setYAt() methods. Simple as that; here is the complete code for reference with comments as further explanations:

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 part.startPoint() == geom.vertexAt(vert): # if the current vertex is a startpoint, 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 part.endPoint() == geom.vertexAt(vert):
                    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
MrXsquared
  • 34,292
  • 21
  • 67
  • 117