# Split 3D polyline at point using QGIS

I want to split a `QgsPolyline` geometry into two separate lines at a given `QgsPoint` (by projection onto the line), adding the point to the line as end point of the incoming line and start point of the outgoing one.

I got the following code that works for 2D geometries (creating a `QgsPointXY` point and a `QgsPolylineXY` poly from the input `QgsFeature`):

``````def breaklineatpoint(pointfeat, polyfeat):
point = pointfeat.geometry().asPoint()
poly = polyfeat.geometry().asPolyline()
geo = QgsGeometry().fromPolylineXY(poly)

(sqdist, projpoint, aftervertex, leftOf) = geo.closestSegmentWithContext(point)
newline1 = poly[:aftervertex]  # from start to break point
newline1.append(projpoint)  # append point
newline2 = poly[aftervertex:]  # from break point to end
newline2.insert(0, projpoint)  # insert at position 0

return newline1, newline2
``````

However, it drops the feature's Z and M values. How do I translate this function so that it can handle 3D geometries (`QgsPolyline` instead of `QgsPolylineXY` and `QgsPoint` instead of `QgsPointXY`)?

I already fail creating a `QgsPolyline` from a `QgsGeometry`. Using a `QgsLineString` by `QgsGeometry().constGet()` is no alternative as I then can't apply the Python list cutting and slicing functions.

After some further research, I found a solution:

The 3D equivalent to `QgsGeometry().asPolyline` (which returns a list of `QgsPointXY`) is `QgsGeometry().constGet().points()`, providing a list of `QgsPoint`. Such list is called a `QgsPointSequence`in `QgsAbstractGeometry` class, but called a `QgsPolyline` in `QgsGeometry` class. That's definetely confusing.

However, my code is now:

``````def breaklineatpoint(pointfeat, polyfeat):
point = pointfeat.geometry().asPoint()
poly = polyfeat.geometry().constGet().points()
geo = QgsGeometry().fromPolyline(poly)

(sqdist, projpoint, aftervertex, leftOf) = geo.closestSegmentWithContext(QgsPointXY(point))
pp = QgsPoint(projpoint)
pp.setZ(interpolateZ(projpoint, poly, aftervertex))
newline1 = poly[:aftervertex]  # from start to break point
newline1.append(pp)  # append point
newline2 = poly[aftervertex:]  # from break point to end
newline2.insert(0, pp)  # insert at position 0
``````

I just need an additional function to interpolate a Z value for my projected point (as this is returned as a 2D QsqPointXY). I use the neighbouring points p1 and p2:

``````def interpolateZ(p, polyline, aftervertex):
linegeo = QgsGeometry().fromPolyline(polyline)
d = linegeo.lineLocatePoint(QgsGeometry().fromPointXY(p))
p1 = polyline[aftervertex-1]
p2 = polyline[aftervertex]
z1 = p1.z()
z2 = p2.z()
sl = p1.distance(p2)
return - d/sl * z1 + z1 + d/sl * z2
``````

Please comment if you know a more simple solution. And I'm very interested in any explaination for the confusing wording in the 3d geometry part of the API.

Here is a different method to split a `MultiLineZM` by several points at once. Unfortunately I was not able to find a method that also works fine in not-projected crs'. So make sure to run this on layers with units meters. I guess `lineLocatePoint()` as well as every other method of `QgsGeometry` to find nearby points just gets confused with degrees.

Basically its searching for the nearest points by using a spatial index and interpolates a point along the line, which is next to the nearest point. Since we dont know the order of the nearest points along the line yet, we store it in a dictionary and sort it. Then we iterate over the sorted dictionary and densify the line with the interpolated points. We need to do this, because afterwards we will iterate over the new parts we want to keep and delete everything else of the line. Why? Because this way its a lot easier to keep gaps in multi line strings. The only thing I am not absolutely sure about is, whether the spatial index always returns the correct nearest neighbors when we feed it with PointZ's and not PointXY's. Seems like the answewr is no.

Here is the example code with an attempt of explanations in the comments:

``````lyr_line = QgsProject.instance().mapLayersByName('linez')[0]
lyr_point = QgsProject.instance().mapLayersByName('pointz')[0]
max_dist = 2500

idx_point = QgsSpatialIndex(lyr_point.getFeatures(), flags=QgsSpatialIndex.FlagStoreFeatureGeometries)

for feat_line in lyr_line.getFeatures():
vertices_dict = {}
line_geom = feat_line.geometry()
# get the nearest points
nearest_point_ids = idx_point.nearestNeighbor(feat_line.geometry(),-1,max_dist) # int
for nearest_point_id in nearest_point_ids:
nearest_point_geom = idx_point.geometry(nearest_point_id)
# get the distance along the line to the point
dist_along_line = line_geom.lineLocatePoint(nearest_point_geom)
# interpolate a point at that distance
point_on_line = line_geom.interpolate(dist_along_line)
# get the vertex id after the interpolated point
vertex_after_id = line_geom.constGet().closestSegment(point_on_line.vertices().next(),10)[2]
# get the vertex index from the vertex id
vertex_after_nr_new = line_geom.vertexNrFromVertexId(vertex_after_id)
# add the distance from the line part as key to a dictionary. We will use the distance from start as sorting criteria afterwards
# add the interpolated point and the vertex index as values
# use try/except in case there is more than one point at the same distance
try:
vertices_dict[dist_along_line] = [point_on_line,vertex_after_nr_new]
except:
pass

# sort the dictionary by distance from line start, so we dont get confused where to insert the new vertices
vertices_dict = dict(sorted(vertices_dict.items()))
from_to_list = []
# insert vertices at the distances to densify the line
densified_geom = feat_line.geometry()
for i, (k, v) in enumerate(vertices_dict.items()):
v[1] += i # count up the vertex index, where we want to insert the new vertex
densified_geom.insertVertex(v[0].vertices().next(),v[1])
from_to_list.append(v[1])
print(densified_geom.asWkt())
# get the index of the last vertex and add indices 0 and last to a list
max_vert = densified_geom.constGet().nCoordinates() - 1
from_to_list.append(0)
from_to_list.append(max_vert)
from_to_list.sort()
# from_to_list contains the vertex indices we build our new lines from and to
print(from_to_list)
for from_to_index, from_to_value in enumerate(from_to_list):
if from_to_value == max_vert:
break
from_vert = from_to_value
try:
to_vert = from_to_list[from_to_index+1]
except IndexError:
to_vert = densified_geom.constGet().nCoordinates() - 1
print('Build line from ',str(from_vert),' to ',str(to_vert))
# create a deepcopy of the geometry we can safely modify
new_geom = QgsGeometry(densified_geom.constGet().clone())
# create a list of the vertices that shall not get deleted
keep_vert_index = list(range(from_vert,to_vert+1))
vertices_deleted = 0
del_vert = 0
# delete index 0 as long as we dont get to the part we do not want to delete
# then delete the first index after the wart we want to keeep until the rest of the line is deleted
for vert_index in range(0,densified_geom.constGet().nCoordinates()):
if vert_index in keep_vert_index:
del_vert += 1 # count up the delete index
continue
# if the second last vertex of a part gets deleted, two vertices are deleted at the same time
# because there is no response about this, we need to somehow figure out when that happened
del_diff = densified_geom.constGet().nCoordinates() - (new_geom.constGet().nCoordinates() + vertices_deleted)
if del_diff > 0:
vertices_deleted += 1
continue
# finally delete the vertex at index
new_geom.deleteVertex(del_vert)
vertices_deleted += 1
# debug
print('vert_index: ',str(vert_index),
' del_vert: ',str(del_vert),
' diff: ',str(del_diff),
' vertices_deleted: ',str(vertices_deleted)
)
# clean up the geometry we created
new_geom.removeDuplicateNodes(4,True)
remove_parts = []
for k, part in enumerate(new_geom.parts()):
if part.length() == 0:
remove_parts.append(k)
# iterate in reverse so we dont mess with part indices
for k in reversed(remove_parts):
new_geom.deletePart(k)
# do whatever with the line part; add it to a layer or just print it
print(new_geom.asWkt())
``````

You can find a full implementation of this method in this algorithm. And why not using `QgsGeometry().splitGeometry()`? Thats because this method takes a list of points as input and not a geometry, so we can only have one huge split line, because we have multiple split points. Means we have no control over possible gaps in the split-line. This may result in unwanted split locations of the line due to unwanted crossings of the split line. It can also not be done one by one because we want to use multiple split points at once to cut the line into pieces.