ST_Distance returning wrong results near the pole

[matthiasheuser]

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

[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)

[Paul Ramsey <pramsey@…>]

In a3d4f6f3/git:

Adjust geodetic tolerance, references #4835

[Paul Ramsey <pramsey@…>]

In 3afd8a1/git:

Adjust geodetic tolerance, references #4835

[Paul Ramsey <pramsey@…>]

In 356d61e/git:

Adjust geodetic tolerance, references #4835

[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.

[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;

[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

[pramsey]

Milestone renamed

[robe]

In prep for 3.1.3 release

[robe]

Reworking sounds painful. Pushing to 3.2.0. Push back if you feel otherwise.

[pramsey]

Incredibly, I found a fix for this issue that doesn't break the rest of the jenga tower. It's small so I'll back patch it as far as seems reasonable. Basically the point→edge distance routine projects a point onto the plane defined by the edge, and then asks if that point is between the edge-defining points. The code that was being used both tested for between-ness (good) and also for on-the-plane-ness (which since we had already constructed the point to be on-the-plane, we could take as read). Making the test more appropriate for our needs side-steps a precision/tolerance issue with the on-the-plane test.

[Paul Ramsey <pramsey@…>]

In a533bf4/git:

For edge_distance_to_point test, do not check that candidate point is on the target plane: we already know it is because we deliberately constructed it to be there. References #4835

[Paul Ramsey <pramsey@…>]

In 40ca39a/git:

For edge_distance_to_point test, do not check that candidate point is on the target plane: we already know it is because we deliberately constructed it to be there. References #4835

[Paul Ramsey <pramsey@…>]

In 4aad298/git:

For edge_distance_to_point test, do not check that candidate point is on the target plane: we already know it is because we deliberately constructed it to be there. References #4835

[Paul Ramsey <pramsey@…>]

In 1f5ed22/git:

For edge_distance_to_point test, do not check that candidate point is on the target plane: we already know it is because we deliberately constructed it to be there. References #4835