[postgis-tickets] [PostGIS] #4835: ST_Distance returning wrong results near the pole

[PostGIS]

#4835: ST_Distance returning wrong results near the pole
----------------------------+---------------------------
 Reporter:  matthiasheuser  |      Owner:  pramsey
     Type:  defect          |     Status:  new
 Priority:  high            |  Milestone:  PostGIS 3.1.2
Component:  postgis         |    Version:  3.1.x
 Keywords:                  |
----------------------------+---------------------------
 Hi,

 I experienced an issue in the calculation of the distance between the
 north pole and a line that passes nearby.

 When running the following query:
 {{{
 select st_distance(st_geogfromtext('POINT (0
 90)'),st_geogfromtext('LINESTRING (-166.11 68.875,15.55 78.216667)'));
 }}}

 I get the value "`1315941.39494886`", where it should be around
 "`25115.89021614`".

 Further invstigation showed that this issue happens in PostGIS versions
 2.4.6 and later (including 3.1.0). The versions before (i.e. 2.4.5) are
 returning the correct result.

 I tried to identify more details about this issue and found out that using
 ST_Segmentize with the line as geography object works, where using
 ST_Segmentize with a geometry object results is the same value as returned
 by the original calculation.

 {{{
 SELECT  ST_Distance(p1, l1) as dist_p1_l1, ST_Distance(p1, l2) as
 dist_p1_l2, ST_Distance(p1, l3) as dist_p1_l3
 FROM (SELECT
         ST_GeogFromText('POINT(0 90)') as p1,
         ST_GeogFromText('LINESTRING(15.55 78.216667, -166.11 68.875)') as
 l1,
         ST_Segmentize(ST_GeomFromText('LINESTRING(15.55 78.216667, -166.11
 68.875)'),10) as l2,
         ST_Segmentize(ST_GeogFromText('LINESTRING(15.55 78.216667, -166.11
 68.875)'),10) as l3
         ) As foo;
 }}}

 Results in:
 {{{
     dist_p1_l1    |    dist_p1_l2    |   dist_p1_l3
 ------------------+------------------+----------------
  1315941.39494886 | 1315941.39494886 | 25115.89026901
 }}}

 The slight difference between the expected value and the value in
 dist_p1_l3 are likely a cause of the segmentation.

 Further investigation using the following SQL showed that this issue does
 not happen only if the line is very close to the pole, it seems to happen
 sporadically at some coordinates, where others that are nearer to the pole
 work.

 {{{
 do $$
 declare
    coord decimal := -170.0;
    end_coord decimal := -159.0;
    inc decimal := 0.01 ;
    dist decimal;
    sqlquery text;
 begin
         loop
                 exit when coord > end_coord ;
                 sqlquery=format('select
 trunc(ST_Distance(ST_GeogFromText(''LINESTRING(15.55 78.216667, %s
 68.875)''), ST_GeogFromText(''POINT(0 90)'')))', coord);
                 EXECUTE sqlquery into dist;
                 raise notice 'Coord: %, Dist: %', coord, dist;
                 coord = coord + inc;
         end loop;
 end; $$
 }}}

 Gives the following output (I only added a small section here):
 {{{
 NOTICE:  Coord: -166.23, Dist: 1315941
 NOTICE:  Coord: -166.22, Dist: 26780
 NOTICE:  Coord: -166.21, Dist: 26628
 NOTICE:  Coord: -166.20, Dist: 1315941
 NOTICE:  Coord: -166.19, Dist: 1315941
 NOTICE:  Coord: -166.18, Dist: 1315941
 NOTICE:  Coord: -166.17, Dist: 26023
 NOTICE:  Coord: -166.16, Dist: 25872
 NOTICE:  Coord: -166.15, Dist: 1315941
 NOTICE:  Coord: -166.14, Dist: 25569
 NOTICE:  Coord: -166.13, Dist: 1315941
 NOTICE:  Coord: -166.12, Dist: 1315941
 NOTICE:  Coord: -166.11, Dist: 1315941
 NOTICE:  Coord: -166.10, Dist: 1315941
 NOTICE:  Coord: -166.09, Dist: 1315941
 NOTICE:  Coord: -166.08, Dist: 1315941
 NOTICE:  Coord: -166.07, Dist: 1315941
 NOTICE:  Coord: -166.06, Dist: 24359
 NOTICE:  Coord: -166.05, Dist: 1315941
 NOTICE:  Coord: -166.04, Dist: 24056
 NOTICE:  Coord: -166.03, Dist: 1315941
 }}}

 The same issue seems to be present in version 2.4.5 as well, but happens
 not as often as in the newer versions.
 The same SQL executed on 2.4.5 gives the following result (again only a
 part of the output is included here):

 {{{
 NOTICE:  Coord: -164.51, Dist: 907
 NOTICE:  Coord: -164.50, Dist: 756
 NOTICE:  Coord: -164.49, Dist: 605
 NOTICE:  Coord: -164.48, Dist: 1315941
 NOTICE:  Coord: -164.47, Dist: 302
 NOTICE:  Coord: -164.46, Dist: 151
 NOTICE:  Coord: -164.45, Dist: 0
 NOTICE:  Coord: -164.44, Dist: 151
 NOTICE:  Coord: -164.43, Dist: 302
 NOTICE:  Coord: -164.42, Dist: 1315941
 NOTICE:  Coord: -164.41, Dist: 605
 NOTICE:  Coord: -164.40, Dist: 756
 NOTICE:  Coord: -164.39, Dist: 907
 }}}

 The occurences of this issue shown above are the only ones I could
 identify for version 2.4.5.

 I did some variations of the above SQL and varied the coordinates of the
 point and experienced this issue not only for latitude 90 but for anything
 above around 89.7.

 The latest PostGIS version I used to test the issue is:

 {{{
 POSTGIS="3.1.0 5e2af69" [EXTENSION] PGSQL="130" GEOS="3.9.0-CAPI-1.16.2"
 PROJ="7.2.1" LIBXML="2.9.1" LIBJSON="0.11"
 }}}

 Best regards,

 Matthias

-- 
Ticket URL: <https://trac.osgeo.org/postgis/ticket/4835>
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.