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

[Claire Palmer]

Thanks for the reply.  Yes, when I run the ST_Distance, _ST_DistanceTree
and the _ST_DistanceUncached I'm getting different results for the adjacent
polygon:

SELECT gid, fld_zone,
ST_Distance(geog_nad83, 'SRID=4269;POINT(-117.876941
33.644253)'::geography) as dist,
_ST_DistanceTree(geog_nad83, 'SRID=4269;POINT(-117.876941
33.644253)'::geography) as dist_tree,
_ST_DistanceUncached(geog_nad83, 'SRID=4269;POINT(-117.876941
33.644253)'::geography) as dist_uncached
FROM nfhl.orange_county_s_fld_haz_ar
WHERE ST_Intersects(geog_nad83, 'SRID=4269;POINT(-117.876941
33.644253)'::geography)

/*
"gid" "fld_zone" "dist" "dist_tree" "dist_uncached"
224245 "X" 0 0 0
224482 "AE" 0 0 190.035004500679
*/

I have tested on PostGIS versions 3.2 + and the result is correct; i.e.
does not return the adjacent polygon (gid 224482) with the ST_Intersects on
geography.

I also see that this issue is documented in other threads.  So apologies
for the duplicate, but thanks for the insight on the internals of the
problem.

Cheers,
Claire

On Sat, Jul 1, 2023 at 4:40 PM Regina Obe <lr at pcorp.us> wrote:

> 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/20230703/ab519a19/attachment.htm>