ST_Distance yields inconsistent results when used in a statement involving multiple Point and Multipolygon Geography pairs

[ewcz]

The function ST_Distance returns inconsistent results when calculating the mutual geographical distance of several points with respect to given (multi)polygon.

Using the tables as in the attached minimal example, following SQL statement works as expected:

sample=# SELECT N.node_id, P.iid, ST_Distance(P.outline::geography, N.location::geography)
FROM N, P WHERE (N.node_id IN (2) AND P.iid=1);
 node_id | iid |   st_distance    
---------+-----+------------------
       2 |   1 | 1730626.41376862
(1 row)

While a similar statement taking into account all three points gives wrong distance for the second point, i.e.,

sample=# SELECT N.node_id, P.iid, ST_Distance(P.outline::geography, N.location::geography)
FROM N, P WHERE (N.node_id IN (1,2,3) AND P.iid=1);
 node_id | iid |   st_distance    
---------+-----+------------------
       1 |   1 | 1951140.06556506
       2 |   1 |                0
       3 |   1 | 1727662.98350914
(3 rows)

However, using _ST_DistanceUnCached in the latter statement recovers the expected result:

sample=# SELECT N.node_id, P.iid, _ST_DistanceUnCached(P.outline::geography, N.location::geography)
FROM N, P WHERE (N.node_id IN (1,2,3) AND P.iid=1);
 node_id | iid | _st_distanceuncached 
---------+-----+----------------------
       1 |   1 |     1951140.06556506
       2 |   1 |     1730626.41376862
       3 |   1 |     1727662.98350914
(3 rows)

Statements above were tested on:

PostgreSQL 9.5.0 on x86_64-pc-linux-gnu, compiled by gcc (GCC) 4.8.5 20150623
(Red Hat 4.8.5-4), 64-bit

POSTGIS="2.2.1 r14555" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.8.0, 6 March 2012"
GDAL="GDAL 1.11.2, released 2015/02/10" LIBXML="2.9.1"
LIBJSON="0.11" RASTER

[ewcz]

sample SQL data

[pramsey]

Confirmed on OSX and against trunk,

select 
  _st_distancetree(n.location, p.outline),
  _st_distanceuncached(n.location, p.outline), 
  n.node_id, p.iid 
from n, p;

 _st_distancetree | _st_distanceuncached | node_id | iid 
------------------+----------------------+---------+-----
 1951140.06557644 |     1951140.06557644 |       1 |   1
                0 |      1730626.4137779 |       2 |   1
  1727662.9835187 |      1727662.9835187 |       3 |   1
(3 rows)

[pramsey]

The point-in-polygon test is short-circuiting the calculation early with an erroneous result, hence the zero distance return. Why? Dunno, probably a robustness failure in the spherical geometry core algorithms, gah.

[dbaston]

Interestingly, the issue seems to disappear if we extract the polygon from p, even though p is a single-member multipolygon:

 select 
  _st_distancetree(n.location, ST_GeometryN(p.outline, 1)),
  _st_distanceuncached(n.location, p.outline), 
  n.node_id, p.iid 
from n, p;
 _st_distancetree | _st_distanceuncached | node_id | iid 
------------------+----------------------+---------+-----
 1951140.06556506 |     1951140.06556506 |       1 |   1
 1730626.41376862 |     1730626.41376862 |       2 |   1
 1727662.98350914 |     1727662.98350914 |       3 |   1

[pramsey]

Oooh, that's a nice clue...

[robe]

A clue you will not do anything with for a while.

[dbaston]

This seems important enough to fix back to 2.2, since it causes PostGIS to report erroneous results.

[robe]

Dan,

I'm planning to release 2.2.3 sometime next week. I feel we should wait for 2.2.4 for this, since I want to be able to say that 2.2.3 and 2.3.0 are on parallel tracks.

[robe]

Sadly we still have not fixed this. My 2.3.1dev shows the same issue.

Pushing this to 2.3.2. If pramsey or dbaston have the energy to fix this by Monday or Tuesday fill free to push back to 2.2.4 and let me know, otherwise going to have to wait till next release.