geography ST_Intersects results depending on insert order

[gekorob]

We have problems with ST_Intersects on geography feature using postgis 2.1.0 and 2.1.1.

The result of the query

select id, name from images where st_intersects(extent, st_geogfromtext('SRID=4326;POLYGON((-46.625977 81.634149,-46.625977 81.348076,-48.999023 81.348076,-48.999023 81.634149,-46.625977 81.634149))'));

on a table

create table images (
        id integer,
        name varchar,
        extent geography(POLYGON,4326)
);
create index extent_idx on images using gist(extent);

is depending on the insert order. In the following situation

insert into images values(47409, 'TDX-1_2010-10-06T19_44_2375085', st_geogfromtext('SRID=4326;POLYGON((-59.4139571913088 82.9486103943668,-57.3528882462655 83.1123152898828,-50.2302874208478 81.3740574826097,-51.977353304689 81.2431047148532,-59.4139571913088 82.9486103943668))'));
insert into images values(1, 'first_image', st_geogfromtext('SRID=4326;POLYGON((-162.211667 88.046667,-151.190278 87.248889,-44.266389 74.887778,-40.793889 75.043333,-162.211667 88.046667))'));

we can't find 'first_image' (with id=1), the query returns '0 rows'. In the opposite situation, when 'first image' is inserted before, the query, correctly returns the record with id = 1; Adding an additional condition at the end of the query, always returns expected result, for example:

select id, name from images where st_intersects(extent, st_geogfromtext('SRID=4326;POLYGON((-46.625977 81.634149,-46.625977 81.348076,-48.999023 81.348076,-48.999023 81.634149,-46.625977 81.634149))')) and name like 'first%';

We experimented this strange behavior only with version 2.1.x of Postgis. The last one we tried was, on Mac:

PostgreSQL 9.3.1 on x86_64-apple-darwin12.5.0, compiled by Apple LLVM version 5.0 (clang-500.2.75) (based on LLVM 3.3svn), 64-bitPOSTGIS="2.1.1 r12113" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" LIBXML="2.9.1" LIBJSON="UNKNOWN"

and on Ubuntu

PostgreSQL 9.3.1 on x86_64-unknown-linux-gnu, compiled by gcc (Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3, 64-bit POSTGIS="2.1.0 r11822" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" LIBXML="2.7.8" LIBJSON="UNKNOWN"

With Postgis 2.0.x and postgres 9.1.x we can't reproduce the problem (everything seems to work correctly, after one year in production mode).

PostgreSQL 9.1.7 on x86_64-unknown-linux-gnu, compiled by gcc (Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3, 64-bitPOSTGIS="2.0.1 r9979" GEOS="3.3.5-CAPI-1.7.5" PROJ="Rel. 4.8.0, 6 March 2012" LIBXML="2.8.0" LIBJSON="UNKNOWN"
}}}.

[gekorob]

Table ddl and data inserts

[pramsey]

Thank you for the complete and easily replicable ticket. Confirmed, something is amiss. (OSX, pgsql 9.3, postgis 2.1). You note your 2.0 postgis provides correct (though slower) results?

For future reference, when you see order-dependence in results usually you're seeing something wrong in the geometry caching system. The uncached initial results and cached secondary results are disagreeing. In your case, there's something *else* going on, since curiously both the brute force and the tree results agree, but the cached result does not.

ST_Intersects is actually just a synonym for distance == 0, and the algorithms are shared between the two operations, so the first thing I always do when I get a report like this is compare the distance results for various ways of running the query:

select id, name, st_distance(extent, st_geogfromtext('SRID=4326;POLYGON((-46.625977 81.634149,-46.625977 81.348076,-48.999023 81.348076,-48.999023 81.634149,-46.625977 81.634149))')) from images;

select id, name, _st_distancetree(extent, st_geogfromtext('SRID=4326;POLYGON((-46.625977 81.634149,-46.625977 81.348076,-48.999023 81.348076,-48.999023 81.634149,-46.625977 81.634149))')) from images;

select id, name, _ST_DistanceUnCached(extent, st_geogfromtext('SRID=4326;POLYGON((-46.625977 81.634149,-46.625977 81.348076,-48.999023 81.348076,-48.999023 81.634149,-46.625977 81.634149))')) from images;

Usually I find a distance between the two algorithms, the tree and the uncached brute force. In this case, the difference is between the wrapped cache function and the others. This will take some effort to debug, which won't happen for a while as I have other paid work at the moment.

[robe]

Wondering if this ticket is similar to the issue I presented in #2422

[gekorob]

Thank you for your prompt answer, we tested the two native distance functions on our full database (approx 6 milions of images). The answers were correct and the time comparable with the workaround we implemented using st_intersects in union with st_contains (as we found that the problem is where one of the two geometries is completely contained in the other one).

To answer your first question about the behavior in postgis version 2.0.x, the results were correct and the query not so time consuming.

Thank you again for your time bye

Roby

[robe]

another complaint of what looks like same issue: ​http://lists.osgeo.org/pipermail/postgis-users/2014-January/038619.html

[robe]

Very interesting. I took the sample test given in ​http://lists.osgeo.org/pipermail/postgis-users/2014-January/038619.html and simplified it a bit and complicated it. As the original example above its the cached_dist that is screwed up. I think the cache may acutally be right but how the distance is being computed from it. In this case they are all wrong. The first is right of course because on first call, it can't use the cache. I also think this issue might be different from the one I have in #2422 (sinc eit that case the distance tree as I recall doesn't agree with distance uncache)

Here is full example:

CREATE TABLE test_cache (id serial, condition_geo geography);
INSERT INTO test_cache (condition_geo) 
 SELECT  ST_Buffer(ST_MakePoint(20.0,30.0)::geography,10.0) 
  FROM generate_series(1,2) i;

SELECT id , r, ST_Distance(condition_geo,geog) As  cached_dist
 , _ST_DistanceUncached(condition_geo,geog) As uncached_dist, _ST_DistanceTree(condition_geo,geog)
 from test_cache CROSS JOIN 
   (SELECT r, ST_Buffer(ST_GeogFromText('POINT(20.0 30.0)'), r) As geog FROM generate_series(1,14,2) As r) As foo
;

output is:

 id | r  | cached_dist | uncached_dist | _st_distancetree
----+----+-------------+---------------+------------------
  1 |  1 |           0 |             0 |                0
  2 |  1 |           0 |             0 |                0
  1 |  3 |           0 |             0 |                0
  2 |  3 |           0 |             0 |                0
  1 |  5 |           0 |             0 |                0
  2 |  5 |           0 |             0 |                0
  1 |  7 |           0 |             0 |                0
  2 |  7 |           0 |             0 |                0
  1 |  9 |           0 |             0 |                0
  2 |  9 |           0 |             0 |                0
  1 | 11 | 0.995469035 |             0 |                0
  2 | 11 | 0.995469035 |             0 |                0
  1 | 13 |   2.9864071 |             0 |                0
  2 | 13 |   2.9864071 |             0 |                0

Not how both 1 and 2 give the same wrong answer using cache (they are in fact the same geography so the order dependence thing I think in this case is an artifact that the cache is not used for the first record). things go wrong right around 10.4 radius.

SELECT id , r, ST_Distance(condition_geo,geog) As  cached_dist
 , _ST_DistanceUncached(condition_geo,geog) As uncached_dist, _ST_DistanceTree(condition_geo,geog)
 from test_cache CROSS JOIN 
   (SELECT 10*(1 + i*0.01) As r, ST_Buffer(ST_GeogFromText('POINT(20.0 30.0)'), 10*(1 + i*0.01)) As geog FROM generate_series(0,5,1) As i) As foo
   WHERE id = 1;

 id |   r   | cached_dist | uncached_dist | _st_distancetree
----+-------+-------------+---------------+------------------
  1 | 10.00 |           0 |             0 |                0
  1 | 10.10 |           0 |             0 |                0
  1 | 10.20 |           0 |             0 |                0
  1 | 10.30 |           0 |             0 |                0
  1 | 10.40 | 0.398187614 |             0 |                0
  1 | 10.50 | 0.497734518 |             0 |                0
(6 rows)

[robe]

I take back what I said about this being different from the issue I posed in #2422. i think there is something screwy with the circ_tree_distance_tree (or circ_tree_distance_tree_internal ) function.

When comparing for the above case why the _st_distancetree gives the right answer while cached_dist is wrong, I discovered that _st_distancetree does a 2 pronged:

if ( CircTreePIP(circ_tree1, g1, lwgeom2) || CircTreePIP(circ_tree2, g2, lwgeom1) )

CircTreePIP test so since one PIP returns true, it never gets to the circ_tree_distance for this particular case. If I add another prong to the geography_distance_cache function, it too returns 0.

I also tried removing the tree_cache use in geography_distance_cache (just building the tree from scrach in both) similar to what the geography_distance_tree does, and it still returns the wrong same answer. So it does seem to be the tree logic here at fault in both bug tickets and cache is probably fine.

I think that may also explain why for buffer below 10.00 the answer is 0 as expected because the CircTreePIP kicks out before you get the circ_tree_distance_tree test.

[robe]

I suspected as much from the numbers. Not sure where in the code it's doing it, but it seems the cached distance is returning the boundary distance rather than polygonal distance. I recall somewhere in code seeing — this is a polygon use line tree algorithm. Can't remember where that was but I thought it a bit odd and thought there must be a switch somewhere else when computing to do area check but couldn't find it.

SELECT id , r, ST_Distance(condition_geo,geog) As  cached_dist
 , _ST_DistanceUncached(condition_geo,geog) As uncached_dist,-- _ST_DistanceTree(condition_geo,geog), 
  _ST_DistanceUncached(condition_geo,ST_Boundary(geog::geometry)::geography) As bound_dist
 from test_cache CROSS JOIN 
   (SELECT 10*(1 + i*0.0001) As r, ST_Buffer(ST_GeogFromText('POINT(20.0 30.0)'), 10*(1 + i*0.0001)) As geog 
       FROM generate_series(320,324,1) As i) As foo
   WHERE id = 1;

 id |    r    | cached_dist | uncached_dist |    bound_dist
----+---------+-------------+---------------+-------------------
  1 | 10.3200 |           0 |             0 | 0.318550091945215
  1 | 10.3210 |  0.31954556 |             0 |  0.31954556026922
  1 | 10.3220 | 0.320541029 |             0 | 0.320541029428955
  1 | 10.3230 | 0.321536499 |             0 | 0.321536499026094
  1 | 10.3240 | 0.322531967 |             0 |  0.32253196740571

[robe]

I've confirmed same issue with original ticket example so the two are the same issue.

SELECT id, _ST_DistanceUnCached(extent, geog) As uncached_dist, ST_Distance(extent,geog) As dist, _ST_DistanceUncached(ST_Boundary(extent::geometry)::geography, ST_Boundary(geog::Geometry)::geography) As bounddist
from images CROSS JOIN st_geogfromtext('SRID=4326;POLYGON((-46.625977 81.634149,-46.625977 81.348076,-48.999023 81.348076,-48.999023 81.634149,-46.625977 81.634149))') As geog;

 id   |  uncached_dist   |      dist       |    bounddist
------+------------------+-----------------+------------------
47409 | 20623.2731762032 | 20623.273176203 | 20623.2731762032
    1 |                0 | 16361.577605712 | 16361.5776057123

[pramsey]

Can you provide a concise example of where _ST_DistanceTree differs from _ST_DistanceUncached? That is the old standard for "uh oh, we're screwed". I'm currently operating under the assumption that it is the cache logic that is wrong, and hope not to be disproved on that.

[robe]

I think I wrote too much on this thread and it confused you. The _ST_DistanceUncached and _ST_DistanceTree I couldn't get to disagree with these examples (that is because the CircTreePIP kicks out with 0

​http://postgis.net/docs/doxygen/2.2/d8/de3/geography__measurement__trees_8c_aa4c02294d4d0a4bd5822023f9e29cc8e.html#aa4c02294d4d0a4bd5822023f9e29cc8e

so circ_tree_distance_tree is not run.

However with the ST_Distance (cached version), it disagrees with _ST_DistanceTree in these examples (because you only have one CircTreePIP check safety net — so it falls into the circ_tree_distance_tree which returns a boundary distance (so when one is inside another if the PIP you are checking is the right one it will return the right answer and if not it returns the boundary distance.

​http://postgis.net/docs/doxygen/2.2/d8/de3/geography__measurement__trees_8c_a399a5e30071a3979b414ce094ff94de9.html#a399a5e30071a3979b414ce094ff94de9

[pramsey]

Tentative fix in 2.1 branch at r12208 Will need a forward port to trunk of acceptable.

[robe]

r12213 (some missing steps along the way)

[pramsey]

Regression tests added and fixes pulled forward to trunk. Ready to close on confirmation of fix.

[robe]

I'm just going to close this out since no one else has stepped up to confirm and it passes the bug ticket examples