Unstable results with ST_Within, ST_Intersects

[strk]

This is pretty puzzling. Load attached dataset, then: (I suspect it affects all functions that used prepared geometry - does same with ST_Intersects)

-- This query says none of the 3 polygons contain the point
=# select gid, ST_Within(
ST_GeomFromText('POINT (-115.04252 36.05137)', -1), the_geom
) from foo_reload;

-- This query says polygon 3 does contain the point
=# select gid, ST_Within(
ST_GeomFromText('POINT (-115.04252 36.05137)', -1), the_geom
) from foo_reload where gid = 3; 

Storing the point geometry into a table makes no difference. Geometry number 3 really contains the point. The other two don't.

I've tested this with PostgreSQL 8.4.3, PostGIS 2.0.0SVN and GEOS 3.3.0dev

[robe]

Just to add to strk — me too, and see ST_DWithin works fine and since ST_Within works fine with just one record, and I don't think ST_DWithin uses prepared geometry, that rules out Geos and any weirdness in PostgreSQL. So I'm guessing prepared geometry hashing. Check this out:

select gid, ST_IsValid(the_geom), ST_Within(ST_GeomFromText('POINT (-115.04252 36.05137)',  -1), the_geom) from foo_reload where gid !=2;


 gid | st_within
-----+-----------
   1 | f
   3 | t

select gid, ST_NPoints(the_geom), ST_Distance(the_geom,ST_GeomFromText('POINT (-115.04252 36.05137)',  -1)) As dist , ST_DWithin(ST_GeomFromText('POINT (-115.04252 36.05137)',  -1), the_geom,0.0000001) from foo_reload

 gid | st_npoints |        dist         | st_dwithin
-----+------------+---------------------+------------
   1 |        111 |   0.149589748004329 | f
   2 |        342 | 0.00350800000001072 | f
   3 |       2916 |                   0 | t

[robe]

Oops wrong output - but still the same result

select gid, ST_IsValid(the_geom), 
ST_Within(ST_GeomFromText('POINT (-115.04252 36.05137)',  -1), the_geom) 
from foo_reload where gid !=2;


 gid | st_isvalid | st_within
-----+------------+-----------
   1 | t          | f
   3 | t          | t

[strk]

I think it has to do with prepared geometries too. You get true where gid > 2 and two falses where gid > 1.

[strk]

Reporter says this bug is in 1.5.2, calling for an 1.5.3 on bugfix here (but not planning to fix this myself, not in the short term at least)

[robe]

I think this highly points the finger at prepared geometry

select gid, ST_IsValid(the_geom), ST_Intersects(ST_GeomFromText('POINT (-115.04252 36.05137)',  -1), the_geom) from foo_reload;

 gid | st_isvalid | st_intersects
-----+------------+---------------
   1 | t          | f
   2 | t          | f
   3 | t          | f

[strk]

Lookin at the code the ST_Within case isn't strictly using geos prepared geometries, but a custom caching thing, takin to point_in_multipolygon_rtree

[strk]

Actually, also ST_Interscets has a non-GEOS implementation, based on cached rings.

[robe]

what about ST_Intersects?

[strk]

As I commented earlier, ST_Intersects is hitting the cache code (out of GEOS). We're probably talking about lwgeom_rtree

[strk]

both ST_Intersects and ST_Within, ONLY in the failing case, are callin point_in_multipolygon_rtree.

[strk]

ST_CoveredBy is another such case. As well pointing at point_in_multipolygon_rtree

[nicklas]

Hallo

I don't see the attached dataset.

It seems to be very close in the code to #852, but I cannot see how that would give this behavior.

But what I found in #852 is that the tolerance doesn't work on vertical and horizontal boundary lines and that the tolerance is not a tolerance of 1e-12 map units as expected but depends on the length of the boundary segment that it is checked against.

Regards Nicklas

[strk]

dateset attached.

[strk]

regress-testing this becomes of upmost importance. Easier if outside of postgresql, with unit testing.

[lreeder]

Have there been any updates on this bug? I was able to workaround this problem by creating thin bridges from the island sub-polygon to the mainland sub-polygon, thus converting the multipolygon to a polygon. However, we've come across this problem again in our production system with a very complicated multi-polygon that will be difficult to fix with the workaround.

[chodgson]

The logic used in point_in_multipolygon_rtree() was different that the logic in point_in_multipolygon(), and gave the wrong answer for the specific case of a point inside a polygon that is inside a hole of a larger polygon (island in a lake inside a larger area).

I modified the multipolygon rtree cache structure to use the standard polygon order for multipolygons instead of the goofy list of all outer rings first followed by a list of all inner rings, which then make it possible to use the right logic in the rtree case.

Committed fix for 1.5 branch in r7136.

[strk]

Chris: can you add a regression test for this in regress/tickets.sql ?

[chodgson]

Yah I'm working on it right now - that's why the ticket's not closed

[chodgson]

Fix merged into trunk and committed in r7137.

Regression tests for 1.5 branch and trunk added in r7138 and 7139, respectively.

[chodgson]

A small error in the earlier patch, fixed in r7140 and r7141.

[chodgson]

lreeder found a new case which breaks my earlier fix.

[chodgson]

My previous bug fix failed to test polygons after the first one due to using continue and skipping over the index incrementing. This is now corrected in r7460 (1.5) and r7461 (trunk/2.0).