Interpolation of Z values in GEOSInterpolate and GEOSInterpolateNormalized

[theroggy]

For both GEOSInterpolate and GEOSInterpolateNormalized, it seems that the interpolation distance along the line is treated as a 2D distance rather than a 3D distance.

Possibly this is by design, but in that case it might be useful to document it. However, especially for a perfectly vertical line this gives some odd results...

In addition, GEOSInterpolateNormalized on a vertical line always seems to return the starting point, even if the relative interpolation distance is 1, which doesn't seem logical even if the distance is treated as 2D.

edit: I added the "length" as calculated by GEOS for each line in the output as well... and it seems that the Z dimension is being ignored in the length calculation.

Test script:

from shapely import LineString

for line in [
    LineString([(0, 0, 0), (0, 0, 1), (0, 0, 2)]),
    LineString([(0, 0, 0), (0, 1, 1), (0, 2, 2)]),
]:
    for normalized in [True, False]:
        for dist in [0.5, 1]:
            result = line.interpolate(dist, normalized=normalized)
            print(f"interpolate {line=} ({line.length=}), {dist=}, {normalized=}: {result}")

Results:

interpolate line=<LINESTRING Z (0 0 0, 0 0 1, 0 0 2)> (line.length=0.0), dist=0.5, normalized=True: POINT Z (0 0 0)
interpolate line=<LINESTRING Z (0 0 0, 0 0 1, 0 0 2)> (line.length=0.0), dist=1, normalized=True: POINT Z (0 0 0)
interpolate line=<LINESTRING Z (0 0 0, 0 0 1, 0 0 2)> (line.length=0.0), dist=0.5, normalized=False: POINT Z (0 0 2)
interpolate line=<LINESTRING Z (0 0 0, 0 0 1, 0 0 2)> (line.length=0.0), dist=1, normalized=False: POINT Z (0 0 2)
interpolate line=<LINESTRING Z (0 0 0, 0 1 1, 0 2 2)> (line.length=2.0), dist=0.5, normalized=True: POINT Z (0 1 1)
interpolate line=<LINESTRING Z (0 0 0, 0 1 1, 0 2 2)> (line.length=2.0), dist=1, normalized=True: POINT Z (0 2 2)
interpolate line=<LINESTRING Z (0 0 0, 0 1 1, 0 2 2)> (line.length=2.0), dist=0.5, normalized=False: POINT Z (0 0.5 0.5)
interpolate line=<LINESTRING Z (0 0 0, 0 1 1, 0 2 2)> (line.length=2.0), dist=1, normalized=False: POINT Z (0 1 1)

Originally reported here: shapely/shapely#1195

[pramsey]

What should/do length and distance return? The trouble here is answering the question of what a 2.5D library should do. I would tend to argue the answer is "never calculate with the extra dimensions, only carry them". We may have already wandered into the soup of having different answers for different questions, but this madness eventually drove PostGIS to say "all calculations are 2D except those where the function signature has a '3D' in it".

[theroggy]

I added to the script above that the length of each linestring according to GEOS is also plotted... so its seems that at the moment, the Z dimension is ignored in the calculation.

What length and distance should return I leave up to you guys to decide...