ST_LineInterpolatePoint incorrect output for 3D linestring

[bencaradocdavies]

For a vertical line (a single segment varying only in Z), ST_LineInterpolatePoint incorrectly outputs the the last point for any fraction greater than zero:

$ select st_asewkt(st_lineinterpolatepoint(st_geomfromewkt('LINESTRING(0 0 0, 0 0 1)'), 0.5));
  st_asewkt   
--------------
 POINT(0 0 1)
(1 row)

The output in this case should be:

 POINT(0 0 0.5)

Non-vertical lines appear to be interpolated correctly:

$ select st_asewkt(st_lineinterpolatepoint(st_geomfromewkt('LINESTRING(0 0 0, 1 0 1)'), 0.5));
    st_asewkt     
------------------
 POINT(0.5 0 0.5)
(1 row)
$ select st_asewkt(st_lineinterpolatepoint(st_geomfromewkt('LINESTRING(0 0 0, 0 1 1)'), 0.5));
    st_asewkt     
------------------
 POINT(0 0.5 0.5)
(1 row)

Seen on CentOS 6.5 x86_64:

postgresql93.x86_64 9.3.2-1PGDG.rhel6
postgis2_93.x86_64 2.1.1-1.rhel6

And on Debian sid amd64:

postgresql-9.3 9.3.2-1 amd64
postgis 2.1.1-5 amd64

[bencaradocdavies]

Lines aligned with the X or Y axis also appear to be interpolated correctly:

$ select st_asewkt(st_lineinterpolatepoint(st_geomfromewkt('LINESTRING(0 0 0, 1 0 0)'), 0.5));
   st_asewkt    
----------------
 POINT(0.5 0 0)
(1 row)

$ select st_asewkt(st_lineinterpolatepoint(st_geomfromewkt('LINESTRING(0 0 0, 0 1 0)'), 0.5));
   st_asewkt    
----------------
 POINT(0 0.5 0)
(1 row)

[pramsey]

The interpolation is done in 2D and then the Z and M values of the generated point are created by interpolating between the values of the segment the point is located in. If we do the interpolation in full 3D, you could get a result like this...

ST_LineInterpolatePoint('LINESTRING(0 0 0, 0 0.5 0, 0 1 100)', 0.5) -> 'POINT(0 0.9 50)'

Hard call. I'm not sure we want this to be a real 3d function, or a 2.5d function as it currently is. It's not a clear "bug" though.

[bencaradocdavies]

Paul,

in my view this is a bug. The documentation clearly states "fraction of total linestring length the point has to be located". My understanding of "length" is that for Cartesian coordinate systems the length of a line segment is the Euclidean distance between the end points, and that the length of a LineString is the sum of its constituent segments. If you have a different interpretation of "length", that is very surprising to me. You are implying that this function should interpolate a point corresponding to the fraction of a path in the projection of the linestring into the plane containing its first two dimensions.

Manual calculation in full three dimensions of the example you gives results in:

ST_LineInterpolatePoint('LINESTRING(0 0 0, 0 0.5 0, 0 1 100)', 0.5) -> 'POINT(0 0.748750015624707 49.75000312494141)'

This is expected as the length of the first segment is 0.5 and the length of the second segment is sqrt(0.52 + 1002) which is a little over 100, so the second segment dominates and the interpolated point is close to the result for interpolating on the second segment alone. The current implementation yields 'POINT(0 0.5 0)'.

I also suggest that this operation should give results that commute with rotation or refection. Ignoring the length contribution of the third coordinate when finding the interpolated point will prevent this.

This is a severe bug that prevents the use of this function for 3D geoscience applications. Are other functions affected by similar interpretations?

Kind regards, Ben.

[pramsey]

There's a lot of 2.5D thinking in the code base, in general, unless a function has a "3D" string in the signature. The linear referencing functions are certainly going to think that X/Y are the primary axes and other axes are along for the ride. I'm concerned about touching this in 2.1 branch for fear it creates a "bug" for anyone who is working with 2.5D data expecting X/Y primacy. And generally concerned about changing the behaviour going forward. I don't think your interpretation is a slam dunk for all users.

[robe]

Milestone renamed

[vmo]

adds 3d line interpolate point

[vmo]

Here is the function, it's not a really nice implementation since that's a lot of cpy/paste from the 2.5d function, but works.

[robe]

Vincent, I was going to rewrite this ticket, but given the long comment history and the fact you fixed by implementing a new function, can you instead create a new ticket.

Just call the title ST_3DLineInterpolatePoint and reattach a revised patch.

At a glance one thing that needs changing in your patch is availability.

Since it's a new function, availability should be changed from 2.4.0 to: 2.5.0

[robe]

Also please include some regression tests. A cunit and sql one will be most appreciated.

[komzpa]

In 17097:

ST_3DLineInterpolatePoint.

Patch by Julien Cabieces and Vincent Mora.

Closes #2645
Closes #4171
Closes ​https://github.com/postgis/postgis/pull/344
Closes ​https://github.com/postgis/postgis/pull/346