Opened 6 years ago

Closed 6 years ago

#1976 closed defect (fixed)

ST_Intersects() returns false while ST_Intersection returns a geography

Reported by: mhiper3pg Owned by: pramsey
Priority: medium Milestone: PostGIS 1.5.6
Component: postgis Version: 1.5.X
Keywords: Cc:


When running ST_Intersects() with multipolygons, sometimes ST_Intersects() returns false for the intersection, while ST_Intersection() returns a geography (not empty).

How to reproduce: loading the geography for the state of Iowa (attached) in a table (with no indexes), and querying for the intersection with a polygon fully contained in the state:

select st_intersects(geog, ST_GeographyFromText('SRID=4326;POLYGON((-95.9 42.9,-95.9 43,-95.8 43,-95.8 42.9,-95.9 42.9))')), st_astext(st_intersection(geog,ST_GeographyFromText('SRID=4326;POLYGON((-95.9 42.9,-95.9 43,-95.8 43,-95.8 42.9,-95.9 42.9))'))) from states where name10='Iowa'

yields this result:

 st_intersects |                                                           st_astext                                                           
 f             | POLYGON((-95.9 42.8999999999999,-95.9 42.9999999999999,-95.8 42.9999999999999,-95.8 42.8999999999999,-95.9 42.8999999999999))
(1 row)

So, the query with ST_Intersects() returns false (which is wrong), but using ST_Intersection() provides the correct result.

Is this behavior normal?

Attachments (3) (239.0 KB) - added by mhiper3pg 6 years ago.
Iowa shapefile to reproduce the problem
screenshot_01.jpg (148.8 KB) - added by pramsey 6 years ago.
View of the inputs
screenshot_01.2.jpg (148.8 KB) - added by pramsey 6 years ago.
View of the inputs

Download all attachments as: .zip

Change History (11)

Changed 6 years ago by mhiper3pg

Attachment: added

Iowa shapefile to reproduce the problem

comment:1 Changed 6 years ago by pramsey

Resolution: invalid
Status: newclosed

There are many things going on here.

First, probably ST_Intersects is actually *right*, the polygon does intersect Iowa, see for an explanation.

Second, because ST_Intersection is not actually running in spherical coordinates (it's a cheater wrap to a planar projection) the results are not consistent with ST_Intersects, which *is* running in spherical.

In all, I think I need to either make the planar hack more sophisticated (rathole) or start thinking about how to do constructive geometry on the sphere (ratmine).

comment:2 Changed 6 years ago by mhiper3pg

I see how ST_Intersects() and ST_Intesection() may yield different restuls. However, I am not so sure that ST_Intersects() is returning the correct result here (sorry if I am being obtuse).

I see the issue with the flat Earth, as I also have several of those cases. However, what I thought weird in this case is that I had expected the great circles to actually "expand" the shape of the north border, which is the closest one, and the one that is defined only by a small set of points (and therefore, more likely to have long straight lines in the linear interpolation). In that case, the small polygon should still intersect with the state.

Another issue is that, when I find this situation, I can usually extend the small polygon in the direction of the closest border of the large polygon, and ST_Intersects() will still return false (meaning that the small polygon is "outside" the large one). However, I cannot do so in this case. (I realize this "test" may fail depending on the shape of the large polygon, but I don't see how it would in this case).

comment:3 Changed 6 years ago by pramsey

Resolution: invalid
Status: closedreopened

I'll have to look closer and see if ST_Intersects is doing the right thing.

Changed 6 years ago by pramsey

Attachment: screenshot_01.jpg added

View of the inputs

Changed 6 years ago by pramsey

Attachment: screenshot_01.2.jpg added

View of the inputs

comment:4 Changed 6 years ago by pramsey

Confirmed there is something fishy in st_intersects in 1.5/2.0. It's in the point-in-polygon test, as seen here, the first point of the polygon returns false

select st_intersects(
  ST_GeographyFromText('SRID=4326;POINT(-95.9 42.9)')
from iowa 
 where name10 = 'Iowa';

comment:5 Changed 6 years ago by pramsey

As an added piece of WTF, 1.5 returns a wrong answer in 120ms, while 2.0/2.1 return wrong answers in 500ms.

comment:6 Changed 6 years ago by pramsey

And another piece: the new tree-based point-in-polygon code returns the *right* answer for that case (first return is uncached, non-tree answer, the rest are tree answers)

st_intersects(geog, ST_GeographyFromText('SRID=4326;POINT(-95.9 42.9)')) 
from iowa where name10='Iowa';

 generate_series | st_intersects 
               0 | f
               1 | t
               2 | t

comment:7 Changed 6 years ago by pramsey

OK, so the speed difference is because between 1.5 and 2.0 we switched from passing boxes into the distance function as part of the signature to letting them ride in implicitly on the lwgeom objects. The trouble with the latter is that the distance function recursively unwraps the geometries, so that when collections get unwrapped, the lwgeom of their sub-geometries that get passed into the next call DO NOT HAVE boxes defined. So they get calculated on the fly, expensively. We need an lwgeom_add_box_deep to add boxes all the way down into lwgeoms for each element.

comment:8 Changed 6 years ago by pramsey

Resolution: fixed
Status: reopenedclosed

Committed patches to trunk (r10527), 2.0 branch (r10528) and 1.5 branch (r10529). The problem was the point in ring algorithm flagging a small edge as being a crossing edge even though it wasn't, because of relatively loose tolerances and error introduced in the edge intersection calculations. The "solution" was a completely new edge intersection routine that hopefully introduces less error and is more robust for more cases.

Note: See TracTickets for help on using tickets.