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.