[PostGIS] #6128: ST_intersects Gives false Positives related to meridian / antimeridian

[PostGIS]

#6128: ST_intersects Gives false Positives related to meridian / antimeridian
----------------------------------------------+---------------------------
 Reporter:  BoulderPika                       |      Owner:  pramsey
     Type:  defect                            |     Status:  new
 Priority:  medium                            |  Milestone:  PostGIS 3.6.5
Component:  postgis                           |    Version:  3.6.x
 Keywords:  geospatial false positive lookup  |
----------------------------------------------+---------------------------
 I am seeing some squirrely behavior with ST_intersects when determining if
 a point intersects a polygon when a polygon touches (but does not span)
 the antimeridian and the lookup point lies on the meridian. This is easily
 reproducible in the latest docker image of postgis I could find.

 **Versions:**\\
 SELECT version();

 {{{
                                                       version
 --------------------------------------------------------------------------------------------------------------------
  PostgreSQL 18.6 (Debian 18.6-1.pgdg13+2) on x86_64-pc-linux-gnu, compiled
 by gcc (Debian 14.2.0-19) 14.2.0, 64-bit
 (1 row)

 }}}

 SELECT postgis_full_version();

 {{{
 postgis_full_version
 ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  POSTGIS="3.6.4 94d984b" [EXTENSION] PGSQL="180" GEOS="3.14.1-CAPI-1.20.5"
 (compiled against GEOS 3.13.1) PROJ="9.8.1 NETWORK_ENABLED=OFF
 URL_ENDPOINT=https://cdn.proj.org
 USER_WRITABLE_DIRECTORY=/var/lib/postgresql/.local/share/proj
 DATABASE_PATH=/usr/share/proj/proj.db" (compiled against PROJ 9.6.0)
 LIBXML="2.9.14" LIBJSON="0.18" LIBPROTOBUF="1.5.1" WAGYU="0.5.0
 (Internal)" TOPOLOGY
 (1 row)
 }}}



 **Below are the steps to setup the geography:**\\
 docker pull postgis/postgis:18-3.6\\
 docker run --name postgis-container -e POSTGRES_PASSWORD=password -d
 postgis/postgis:18-3.6\\
 docker exec -it postgis-container psql -U postgres\\

 -- (1) Initialize\\

 {{{
 CREATE EXTENSION IF NOT EXISTS postgis;
 }}}

 -- (2) Versions before 3.12 had a bug:
 https://trac.osgeo.org/postgis/ticket/5838\\

 {{{
 SELECT PostGIS_GEOS_Version();
 }}}


 -- (3) Make a table with (id, name, geography)\\

 {{{
 CREATE TABLE geojson_features ( id SERIAL PRIMARY KEY, name TEXT, geom
 geography(POLYGON, 4326) NOT NULL );
 }}}


 -- (4) Add an index\\

 {{{
 CREATE INDEX idx_geojson_features_geom ON geojson_features USING GIST
 (geom);
 }}}


 -- (5) Add a polygon in the northern hemisphere, which contains the point
 (lon, lat) = (0, 25)\\

 {{{
 INSERT INTO geojson_features (name, geom) VALUES('poly1',
 ST_GeomFromGeojson(' { "type": "Polygon", "coordinates": [ [ [ -10, 10 ],
 [ -5, 10 ], [ 5, 20 ], [ 5, 30 ], [ -5, 40 ], [ -10, 40 ], [ -10, 10 ] ] ]
 } ' )::geography);
 }}}


 -- (6) The below southern hemisphere polygon is the western portion of
 "poly1" when translated to span the antemeridian and then sliced there\\

 {{{
 INSERT INTO geojson_features (name, geom) VALUES('poly2_west',
 ST_GeomFromGeojson(' { "type": "Polygon", "coordinates": [ [ [ -180,
 -35.4101725284 ], [ -175, -30 ], [ -175, -20 ], [ -180, -15.1725999559 ],
 [ -180, -35.4101725284 ] ] ] } ' )::geography);
 }}}

 -- (7) The below southern hemisphere polygon is the eastern portion of
 "poly1" when translated to span the antemeridian and then sliced there\\

 {{{
 INSERT INTO geojson_features (name, geom) VALUES('poly2_east',
 ST_GeomFromGeojson(' { "type": "Polygon", "coordinates": [ [ [ 180,
 -15.1725999559 ], [ 175, -10 ], [ 170, -10 ], [ 170, -40 ], [ 175, -40 ],
 [ 180, -35.4101725284 ], [ 180, -15.1725999559 ] ] ] } ' )::geography);
 }}}


 =====================================================

 At this point, we have three polygons setup. I now wish to perform some
 point intersections. The below two execute as expected with a point very
 close to zero lon:

 {{{
 SELECT name, ST_intersects(geom, ST_makepoint(-0.00001, 25)::geography)
 FROM geojson_features;
 }}}
 name | st_intersects\\
 -------+---------------\\
 poly1 | t\\
 poly2 | f\\
 poly3 | f\\
 (3 rows)\\


 {{{
 SELECT name, ST_intersects(geom, ST_makepoint(0.00001, 25)::geography)
 FROM geojson_features;
 }}}
 name | st_intersects\\
 -------+---------------\\
 poly1 | t\\
 poly2 | f\\
 poly3 | f\\
 (3 rows)\\

 But on the next query, I get false positive results:

 {{{
 SELECT name, ST_intersects(geom, ST_makepoint(0, 25)::geography) FROM
 geojson_features;
 }}}
 name | st_intersects\\
 -------+---------------\\
 poly1 | t\\
 poly2 | t\\
 poly3 | t\\
 (3 rows)\\


 Now the polygons in the southern hemisphere are included when the
 longitude is exactly 0. But they were excluded when close to zero.

 **Additional Details:**
 If I were to setup the two southern hemisphere polygons such that they
 were very close to the antimeridian (e.g. edges at 179.9999 instead of 180
 and -179.9999 instead of -180) then the intersections behave as expected.

 This also appears to be a symptom of not using the index. When I have 5000
 more polygons in the table, then I confirmed the index is being used in
 the optimizer and I get the results I expect. With just the three
 polygons, a sequetial scan is performed and I get the false positives.
-- 
Ticket URL: <https://trac.osgeo.org/postgis/ticket/6128>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.