pyqgissegment
I want to segmentize a MultiLineString layer by a given maximum distance. Result shall be parts of the input linestring with this maximum length. I am using QGIS 3.16 and want a pure PyQGIS solution.
Note: It took a while to research and I could not find a question or answer for PyQGIS, so this is written as self-answer question, therefore no codeattempt. But of course, please add your solution!
It's kind of easy to do with a python script that you can run from the editor of the Qgis console.
First, you need to get your line layer and get the feature inside.
Then, you need to loop on your csvfile to get the distance and create the point with the interpolate method of QgsGeometry.
Finally, add the created point to a new point layer.
The following do the trick, you just need to replace the values in the firsts lines (line_layer, csvfilepath, EPSG) to fit to your data:
import csv
from PyQt4.QtCore import QVariant
#Change as appropriate
layer_name = 'line_truck'
csvfile_path = 'C:/Users/ylecomte/Desktop/test.csv'
EPSG = '29902' # projected in meters units
#get layer and create csv iterator
layer = QgsMapLayerRegistry.instance().mapLayersByName(layer_name)[0]
csvfile = open(csvfile_path, 'rb')
reader = csv.reader(csvfile, delimiter=';')
header = reader.next()
#prepare output layer
mem_layer = QgsVectorLayer("Point?crs=epsg:"+EPSG, 'point', 'memory')
mem_layer.startEditing()
attr = [QgsField(header[0],QVariant.String),QgsField(header[1],QVariant.Double)]
prov =mem_layer.dataProvider()
prov.addAttributes(attr)
mem_layer.updateFields()
#preform the trick by looping on feature and csv
for feat in layer.getFeatures():
for row in reader:
new_feat = QgsFeature()
new_feat.setGeometry(feat.geometry().interpolate(float(row[0])))
new_feat.setAttributes(attr)
new_feat.setAttribute(0, float(row[0]))
new_feat.setAttribute(1,row[1])
mem_layer.addFeatures([new_feat])
#save results and add output to the canvas
mem_layer.commitChanges()
QgsMapLayerRegistry.instance().addMapLayer(mem_layer)
the new layer is added to your map canvas as a memory layer containing the needed points and the csv data in corresponding fields.
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
Best Answer
You can iterate over the parts of the MultiLines and use
QgsLineString.curveSubstring()to create the segments. To know the number of segments as well as start- and endpoints of the segments you only need some simple math: