Opened 21 months ago

Closed 3 months ago

#4835 closed defect (fixed)

ST_Distance returning wrong results near the pole

Reported by: matthiasheuser Owned by: pramsey
Priority: high Milestone: PostGIS 3.3.0
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 (18)

comment:1 by pramsey, 20 months ago

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 by Paul Ramsey <pramsey@…>, 20 months ago

In a3d4f6f3/git:

Adjust geodetic tolerance, references #4835

comment:3 by Paul Ramsey <pramsey@…>, 20 months ago

In 3afd8a1/git:

Adjust geodetic tolerance, references #4835

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

In 356d61e/git:

Adjust geodetic tolerance, references #4835

comment:5 by pramsey, 20 months ago

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 by pramsey, 20 months ago

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 by matthiasheuser, 20 months ago

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 by pramsey, 17 months ago

Milestone: PostGIS 3.1.23.1.3

comment:9 by pramsey, 17 months ago

Milestone: 3.1.3PostGIS 3.1.3

Milestone renamed

comment:10 by robe, 15 months ago

Milestone: PostGIS 3.1.3PostGIS 3.1.4

In prep for 3.1.3 release

comment:11 by robe, 13 months ago

Milestone: PostGIS 3.1.4PostGIS 3.2.0

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

comment:12 by pramsey, 12 months ago

Milestone: PostGIS 3.2.0PostGIS 3.3.0

comment:13 by pramsey, 3 months ago

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.

comment:14 by Paul Ramsey <pramsey@…>, 3 months ago

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

comment:15 by Paul Ramsey <pramsey@…>, 3 months ago

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

comment:16 by Paul Ramsey <pramsey@…>, 3 months ago

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

comment:17 by Paul Ramsey <pramsey@…>, 3 months ago

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

comment:18 by pramsey, 3 months ago

Resolution: fixed
Status: newclosed
Note: See TracTickets for help on using tickets.