Split 3D polyline at point using QGIS

Question

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.

Answer 1

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 QgsPointSequencein 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.

Answer 2

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.