PyQGIS – Modifying Specific Vertices of MultiLineString in PyQGIS

multipartpyqgis

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:

MultiLineString Layer
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.

Best Answer

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

Related Solutions

QGIS PyQGIS: Snapping Lines’ Start/End Points to Nearest Features

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:

PyQGIS – Modifying Attributes with Map Selection Using PyQGIS

Try the script below. Note that the 'A' and 'B' fields must be of numeric type (not string).

class customMapTool(QgsMapToolIdentifyFeature):
    
    def __init__(self, iface):
        self.iface = iface
        self.canvas = self.iface.mapCanvas()
        self.layer = self.iface.activeLayer()
        QgsMapToolIdentifyFeature.__init__(self, self.canvas, self.layer)
        self.iface.currentLayerChanged.connect(self.active_changed)
        
        self.first_click = True
        
        self.A_Val = None
        self.B_Val = None
        
        
    def active_changed(self, layer):
        if layer.type() == QgsMapLayerType.VectorLayer and layer.isSpatial():
            self.layer = layer
            self.setLayer(self.layer)
            
    def canvasReleaseEvent(self, event):
        if event.button() == Qt.LeftButton:
            if self.layer.geometryType() != QgsWkbTypes.PointGeometry:
                self.iface.messageBar().pushMessage('Please select a point layer', Qgis.Info, 2)
                return
            found_features = self.identify(event.x(), event.y(), [self.layer], QgsMapToolIdentify.ActiveLayer)
            if found_features:
                feat = found_features[0].mFeature
                if self.first_click:
                    if feat['A'] is not NULL and feat['B'] is not NULL:
                        self.A_Val = feat['A']
                        self.B_Val = feat['B']
                        self.first_click = False
                    return
                elif not self.first_click:
                    if self.A_Val and self.B_Val:
                        self.A_Val += 1
                        self.B_Val += 1
                        self.update_attributes(feat)
                    return
        elif event.button() == Qt.RightButton:
            self.iface.actionPan().trigger()
    
    def update_attributes(self, feat):
        a_idx = self.layer.fields().lookupField('A')
        b_idx = self.layer.fields().lookupField('B')
        if a_idx != -1 and b_idx != -1:
            if self.A_Val and self.B_Val:
                self.layer.dataProvider().changeAttributeValues({feat.id(): {a_idx: self.A_Val,
                                                                        b_idx: self.B_Val}})
            self.layer.triggerRepaint()
        
    def deactivate(self):
        self.iface.currentLayerChanged.disconnect(self.active_changed)
            
        
t = customMapTool(iface)
iface.mapCanvas().setMapTool(t)

See screencast below for expected results. I labelled the points with the expression: "A" ||', '|| "B" and added a call to self.layer.triggerRepaint() to illustrate the feature attributes being updated on each mouse click.

Best Answer

Related Solutions

QGIS PyQGIS: Snapping Lines’ Start/End Points to Nearest Features

PyQGIS – Modifying Attributes with Map Selection Using PyQGIS

Related Question