[postgis-users] ST_Intersects on geography in WHERE clause returns wrong results

[Regina Obe]

It’s hard to tell from below if there is an issue or not because something intersects in planar space or spheroidal space doesn’t guarantee the other

Especially when you are talking about fairly large areas.

 

Can you provide a sample of the data set you are working with and I can test with a newer version of PostGIS.

 

We have run into issues with geography in the past and usually the culprit is tolerance or some issue with distancetree vs. uncached answers.

 

So one way to check is to compare what ST_Distance, _ST_DistanceTree  and _ST_DistanceUncached return as answers.  They should all return the same answer.

 

Geography ST_Intersects just uses Distance with a very small tolerance as discussed in these tickets, https://trac.osgeo.org/postgis/ticket/4223 , https://trac.osgeo.org/postgis/ticket/4932

 

All which have been fixed since your release version.

 

Geometry on the other hand uses DE9I matrix - https://postgis.net/docs/using_postgis_query.html#DE-9IM  but that doesn’t mean the answer geometry gives you is right either since it is assuming straight lines.  One other option to confirm is to transform your data to a suitable space preserving California planar (meter or foot) spatial ref sys that covers your fload zone and confirm the answers are still the same.

 

From: postgis-users [mailto:postgis-users-bounces at lists.osgeo.org] On Behalf Of Claire Palmer
Sent: Saturday, July 1, 2023 7:18 PM
To: postgis-users at lists.osgeo.org
Subject: [postgis-users] ST_Intersects on geography in WHERE clause returns wrong results

 

Dear PostGIS users,

 

I am running a simple point-in-polygon analysis on Orange County flood zone polygons.  I noticed that I'm getting some false positives when running ST_Intersects on geography in the where clause.  

 

My postgis_full_version is:

POSTGIS-"2.4.3 R16312" PGSQL=100 GEOS="3.6.2-CAPI-1.10.2 4d2925d" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.3, released 2017/11/20 GDAL_DATA not found" LIBXML="2.7.8" LIBJSON="0.12" LIBROTOBUF="1.2.1" RASTER

 

The following query should return one flood zone, but instead it's returning two:

SELECT gid, fld_zone, ST_Intersects(geog_nad83, 'SRID=4269;POINT(-117.876941 33.644253)'::geography)
FROM nfhl.orange_county_s_fld_haz_ar
WHERE ST_Intersects(geog_nad83, 'SRID=4269;POINT(-117.876941 33.644253)'::geography);

 

It's selecting the polygon the point intersects with, and also the adjacent polygon.

Here's a visual of the two polygons together with the test point:

 

 

 

If I remove the where clause and query ST_Intersects for the adjacent polygon, the result is correctly FALSE:

SELECT gid, fld_zone, ST_Intersects(geog_nad83, 'SRID=4269;POINT(-117.876941 33.644253)'::geography)
FROM nfhl.orange_county_s_fld_haz_ar
WHERE gid = 224482;

 

To test further, I generated 1000 points on a 200 meter buffer around the test point.  The results can be seen in the following visual.  Some points return only the polygon they intersect with (grey, green) but others return an additional adjacent polygon (blue, red):

 

 

 

Also, I should note that running the same query on geometry does not introduce any errors.  We'd like to stick with using geography if possible and were a bit surprised with these false positives.

 

Does anyone know what's going on?  Does it have to do with the earlier PostGIS version?  We are planning an upgrade to the latest version, but our API is stuck for now with the version I posted above.  

 

Any insight as to why the false positives are occurring would be appreciated!

 

Thanks,

Claire

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20230701/3c382f94/attachment.htm>