ST_Intersects() returns false while ST_Intersection returns a geography

[mhiper3pg]

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?

[mhiper3pg]

Iowa shapefile to reproduce the problem

[pramsey]

There are many things going on here.

First, probably ST_Intersects is actually *right*, the polygon does intersect Iowa, see ​http://blog.opengeo.org/2012/04/30/the-earth-is-not-flat-volume-2/ 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).

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

[pramsey]

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

[pramsey]

View of the inputs

[pramsey]

View of the inputs

[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(
  geog, 
  ST_GeographyFromText('SRID=4326;POINT(-95.9 42.9)')
) 
from iowa 
 where name10 = 'Iowa';

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

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

select 
generate_series(0,10),
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

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

[pramsey]

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.