[PostGIS] #6133: Geography line interpolation: NaN heading along a meridian, wrong meridian past a pole, uninitialised Z and M

[PostGIS]

#6133: Geography line interpolation: NaN heading along a meridian, wrong meridian
past a pole, uninitialised Z and M
-------------------------------------------------+-------------------------
 Reporter:  ezimanyi                             |      Owner:  pramsey
     Type:  defect                               |     Status:  new
 Priority:  medium                               |  Milestone:  PostGIS
                                                 |  3.6.5
Component:  postgis                              |    Version:  3.6.x
 Keywords:  geography ST_LineInterpolatePoint    |
  ST_LineInterpolatePoints sphere_direction      |
  sphere_project                                 |
-------------------------------------------------+-------------------------
 The interpolation of geography lines, which came with #5460, goes
 through `interpolate_point4d_spheroid()` in
 `liblwgeom/lwgeodetic_measures.c` and, on the sphere, through
 `sphere_direction()` and `sphere_project()` in `liblwgeom/lwgeodetic.c`.
 Each of the three has a defect of its own. All three reproduce on
 PostGIS 3.6.3, and master 2bb71185a8 carries the same code.

 == 1. sphere_direction() returns NaN for two points on one meridian ==

 `sphere_direction()` computes the initial heading from its cosine,

 {{{
 f = (sin(lat2) - sin(lat1) * cos(d)) / (sin(d) * cos(lat1))
 heading = acos(f)
 }}}

 For two points on one meridian the exact value of `f` is 1 or -1, and
 rounding leaves it outside [-1, 1] by more than `FP_TOLERANCE` (5e-14
 for the geodetic code), so `acos()` returns NaN; the `fabs(f) > 1.0`
 branch logs the value and still calls `acos()`. Over 43800 segments
 along the prime meridian between 80S and 80N, up to 0.02 degrees long,
 these expressions give NaN for 17791 of them. `sphere_project()` then
 fails, and (see 3) the result is the caller's uninitialised point:

 {{{#!sql
 SELECT ST_AsText(ST_LineInterpolatePoint(
   'LINESTRING(0 -79.269, 0 -79.239)'::geography, 0.5, false));
 -- POINT(6.95314402316196e-310 6.9229486544248e-310)

 SELECT ST_AsText(ST_LineInterpolatePoint(
   'LINESTRING(0 -79.269, 0 -79.239)'::geography, 0.5, true));
 -- POINT(0 -79.25400000727056)
 }}}

 The circular tree already guards against the NaN:
 `circ_center_spherical()` in `lwgeodetic_tree.c` tests the heading
 with `isnan()` ("Catch sphere_direction when it barfs") and falls back
 to a Cartesian centre.

 == 2. sphere_project() keeps the start meridian past a pole ==

 For an azimuth of 0 or pi, `sphere_project()` keeps the start
 longitude, which holds up to the pole and not beyond it; the source
 marks it with "TODO: this isn't quite true, what if we're going over
 the pole?". Past the pole the path runs down the opposite meridian:

 {{{#!sql
 SELECT ST_AsText(ST_LineInterpolatePoint(
   'LINESTRING(0 88, 180 88)'::geography, 0.75, false));
 -- POINT(0 89.0000000000001)

 SELECT ST_AsText(ST_LineInterpolatePoint(
   'LINESTRING(0 88, 180 88)'::geography, 0.75, true));
 -- POINT(180 89.00000307849999)
 }}}

 == 3. interpolate_point4d_spheroid() leaves Z, M, and a failed x and y
 uninitialised ==

 The function sets only x and y, and only when the projection succeeds:

 {{{
 /* If success, use newly computed lat and lon,
 * otherwise return precomputed cartesian result */
 if (success == LW_SUCCESS)
 {
       p->x = rad2deg(longitude_radians_normalize(g.lon));
       p->y = rad2deg(latitude_radians_normalize(g.lat));
 }
 }}}

 Nothing precomputes the Cartesian result the comment names, and both
 callers, `geography_substring()` and `geography_interpolate_points()`,
 pass an uninitialised `POINT4D pt`. So every interpolated Z and M is
 whatever the stack held, on the sphere and on the spheroid alike:

 {{{#!sql
 SELECT ST_AsText(ST_LineInterpolatePoint(
   'LINESTRING Z(0 0 0, 10 0 10)'::geography, 0.5, true));
 -- POINT Z (5.000000000201743 0 6.95291892163796e-310)

 SELECT ST_AsText(ST_LineInterpolatePoint(
   'LINESTRING Z(0 0 0, 10 0 10)'::geography, 0.5, false));
 -- POINT Z (5.000000000202194 3.057732552914371e-16 5.0170810976182e-310)
 }}}

 The function comes from MobilityDB (#5460), whose version, now in
 `meos/src/geo/tpoint_spatialfuncs.c`, calls `interpolate_point4d(p1,
 p2, p, f)` before the success test; that call is absent from the
 PostGIS version.

 == Fix ==

  1. `sphere_direction()` reads the heading with `atan2()` from both
     components of the initial bearing,
     {{{
 heading = atan2(sin(dlon) * cos(lat2),
                 cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dlon))
     }}}
     The sine and cosine of the heading are these two components
     divided by sin(d) > 0, so this is the heading wherever the two
     points are distinct, with no cosine to leave its range and no
     clamp. It answers exactly 0 or pi along a meridian and keeps the
     sign convention. The distance argument is no longer read and stays
     in the signature.
  2. `sphere_project()` keeps its shortcut and answers the opposite
     meridian past the pole,
     {{{
 lon2 = ( cos(d) - sin(lat1) * sin(lat2) < 0.0 ) ? lon1 + M_PI : lon1;
     }}}
     `cos(d) - sin(lat1) * sin(lat2)` is the second component of the
     general formula below the shortcut and equals `cos(lat1) *
     cos(lat2) * cos(dlon)`, so along a meridian its sign is that of
     `cos(dlon)`: positive on the start meridian, negative on the
     opposite one. The TODO goes.
  3. `interpolate_point4d_spheroid()` computes the Cartesian
     interpolation first, which also interpolates z and m, and a
     successful projection overwrites x and y.

 == Tests ==

 The patch adds to `liblwgeom/cunit/cu_geodetic.c` the meridian pair
 (0 -79.269) and (0 -79.239) in both directions to
 `test_sphere_direction()`, expecting 0 and pi, and a projection due
 north from 88N over 3 degrees to `test_sphere_project()`, expecting 89N
 on the opposite meridian. It adds to `regress/core/geography.sql` the
 cases `lrs_lip_pole_sphere`, `lrs_lip_meridian_sphere`, `lrs_lip_z` and
 `lrs_lip_zm`.

 Measured on master 2bb71185a8 with PostgreSQL 18.3 and CUnit 2.1-3:
  * with the patch, the cunit geodetic suite passes 25 of 25 tests and
    536 of 536 assertions, and every suite 375 tests and 5901
    assertions, none failing;
  * with the new tests but without the fixes in `lwgeodetic.c` and
    `lwgeodetic_measures.c`, the geodetic suite fails 2 of 25 tests and
    3 of 536 assertions, exactly the new ones: `cu_geodetic.c:96` and
    `:98` (the heading along the meridian, NaN) and `:185` (the
    longitude past the pole, 0 instead of pi);
  * with the patch, `regress/core/geography` passes under
    `run_test.pl --extension`, the four new cases answering
    `POINT(180 89)`, `POINT(0 -79.254)`, `POINT Z (5 0 5)` and
    `MULTIPOINT ZM ((5 0 5 10),(10 0 10 20))`.

 These were found while fixing the same functions in the copy of
 liblwgeom that MobilityDB vendors:
 [https://github.com/MobilityDB/MobilityDB/pull/2722 MobilityDB#2722]
 (the heading) and
 [https://github.com/MobilityDB/MobilityDB/pull/2720 MobilityDB#2720]
 (the pole).
-- 
Ticket URL: <https://trac.osgeo.org/postgis/ticket/6133>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.