Opened 5 months ago

Last modified 4 weeks ago

#4835 new defect

ST_Distance returning wrong results near the pole

Reported by: matthiasheuser Owned by: pramsey
Priority: high Milestone: PostGIS 3.1.3
Component: postgis Version: 3.1.x
Keywords: Cc:

Description

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

Change History (9)

comment:1 Changed 4 months ago by pramsey

I fear this is in general going to require a re-write of the core geographic routines, as the "solution" is jimmying with the tolerances to find particular very very very small dot-products to be zero or not. (Changed FP_TOLERANCE to 5e-14 in lwgeodetic.c)

comment:2 Changed 4 months ago by Paul Ramsey <pramsey@…>

In a3d4f6f3/git:

Adjust geodetic tolerance, references #4835

comment:3 Changed 4 months ago by Paul Ramsey <pramsey@…>

In 3afd8a1/git:

Adjust geodetic tolerance, references #4835

comment:4 Changed 4 months ago by Paul Ramsey <pramsey@…>

In 356d61e/git:

Adjust geodetic tolerance, references #4835

comment:5 Changed 4 months ago by pramsey

If you're capable of testing I'd be interested to know just how much better it is. My guess is "a little better but not perfect". I think for perfect I have to re-work this stuff with a higher-precision calculation.

comment:6 Changed 4 months ago by pramsey

Just going to store this here for future use

WITH foo AS (
	SELECT
	'POINT(0 90)'::geography as p1,
	'LINESTRING(15.55 78.216667, -166.11 68.875)'::geography as l1,
	ST_Segmentize('LINESTRING(15.55 78.216667, -166.11 68.875)'::geography, 5) as l2,
	ST_Segmentize('LINESTRING(15.55 78.216667, -166.11 68.875)'::geography, 10) as l3
	)
SELECT 
_ST_DistanceUnCached(p1, l1) as dist_p1_l1, 
_ST_DistanceUnCached(p1, l2) as dist_p1_l2, 
_ST_DistanceUnCached(p1, l3) as dist_p1_l3
FROM foo;

comment:7 Changed 4 months ago by matthiasheuser

I compiled the new version today and did some test runs of the script mentioned in the issue. As you already said, it is better but not perfect.

Here are the results (showing only the first occurrences of the issue):

Version 3.1.0 (at about 74 km distance):

NOTICE:  Coord: -169.35, Dist: 1315941
NOTICE:  Coord: -159.55, Dist: 1315941

Version 3.1.2dev (at about 14 km distance):

NOTICE:  Coord: -165.37, Dist: 1315941
NOTICE:  Coord: -163.53, Dist: 1315941

Version 2.4.5 (at about 0.5 km distance):

NOTICE:  Coord: -164.48, Dist: 1315941
NOTICE:  Coord: -164.42, Dist: 1315941

comment:8 Changed 4 weeks ago by pramsey

Milestone: PostGIS 3.1.23.1.3

comment:9 Changed 4 weeks ago by pramsey

Milestone: 3.1.3PostGIS 3.1.3

Milestone renamed

Note: See TracTickets for help on using tickets.