Opened 8 years ago

Closed 8 years ago

#884 closed defect (fixed)

Unstable results with ST_Within, ST_Intersects

Reported by: strk Owned by: chodgson
Priority: high Milestone: PostGIS 1.5.3
Component: postgis Version: 1.5.X
Keywords: history Cc: lreeder

Description (last modified by robe)

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

Attachments (1)

foo_reload.sql (106.2 KB) - added by strk 8 years ago.

Download all attachments as: .zip

Change History (26)

comment:1 Changed 8 years ago by 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

comment:2 Changed 8 years ago by 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

comment:3 Changed 8 years ago by strk

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

comment:4 Changed 8 years ago by strk

Milestone: PostGIS 2.0.0PostGIS 1.5.3
Version: trunk1.5.X

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)

comment:5 Changed 8 years ago by robe

Description: modified (diff)
Summary: Unstable results from ST_WithinUnstable results with Prepared Geometry (ST_Within, ST_Intersects)

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

comment:6 Changed 8 years ago by strk

Summary: Unstable results with Prepared Geometry (ST_Within, ST_Intersects)Unstable results with ST_Within, ST_Intersects

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

comment:7 Changed 8 years ago by strk

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

comment:8 Changed 8 years ago by robe

what about ST_Intersects?

comment:9 Changed 8 years ago by strk

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

comment:10 Changed 8 years ago by strk

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

comment:11 Changed 8 years ago by strk

ST_CoveredBy is another such case. As well pointing at point_in_multipolygon_rtree

comment:12 Changed 8 years ago by 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

Changed 8 years ago by strk

Attachment: foo_reload.sql added

comment:13 Changed 8 years ago by strk

dateset attached.

comment:14 Changed 8 years ago by strk

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

comment:15 Changed 8 years ago by lreeder

Cc: lreeder added

comment:16 Changed 8 years ago by 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.

comment:17 Changed 8 years ago by chodgson

Owner: changed from pramsey to chodgson
Status: newassigned

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.

comment:18 Changed 8 years ago by strk

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

comment:19 Changed 8 years ago by chodgson

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

comment:20 Changed 8 years ago by chodgson

Fix merged into trunk and committed in r7137.

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

comment:21 Changed 8 years ago by chodgson

Resolution: fixed
Status: assignedclosed

comment:22 Changed 8 years ago by chodgson

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

comment:23 Changed 8 years ago by robe

Keywords: history added

comment:24 Changed 8 years ago by chodgson

Resolution: fixed
Status: closedreopened

lreeder found a new case which breaks my earlier fix.

comment:25 Changed 8 years ago by chodgson

Resolution: fixed
Status: reopenedclosed

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

Note: See TracTickets for help on using tickets.