[postgis-tickets] [PostGIS] #5025: Inconsistent behavior with ST_Within across versions

[PostGIS]

#5025: Inconsistent behavior with ST_Within across versions
---------------------+--------------------------------
 Reporter:  byrman   |      Owner:  pramsey
     Type:  defect   |     Status:  new
 Priority:  medium   |  Milestone:  PostGIS PostgreSQL
Component:  postgis  |    Version:  3.1.x
 Keywords:           |
---------------------+--------------------------------
 I noticed inconsistent behavior with ST_Within across different
 PostgreSQL-PostGIS versions. To demonstrate the issue, I used
 postgis/postgis Docker images and a minimal example.

 === Steps to reproduce the problem

 {{{
 -- Create table "a" having both a 2D and 3D spatial index on its geometry
 field.
 CREATE TABLE a (id serial, geom geometry(PointZ, 4326) NOT NULL);
 CREATE INDEX a_2d_idx ON a USING GIST(geom);
 CREATE INDEX a_3d_idx ON a USING GIST(geom gist_geometry_ops_nd);

 -- Create table "b" having both a 2D and 3D spatial index on its geometry
 field.
 CREATE TABLE b (id serial, geom geometry(PolygonZ, 4326) NOT NULL);
 CREATE INDEX b_2d_idx ON b USING GIST (geom);
 CREATE INDEX b_3d_idx ON b USING GIST (geom gist_geometry_ops_nd);

 -- Insert a polygon (z=0) into table "b".
 INSERT INTO b (geom) SELECT ST_Force3D(ST_MakeEnvelope(-1, -1, 1, 1,
 4326));

 -- Generate 100 pseudo-random points (z=0) which lie within the polygon
 and insert them into table "a".
 INSERT INTO a (geom) SELECT ST_Force3D((ST_Dump(ST_GeneratePoints((SELECT
 geom FROM b), 100))).geom);

 -- Change the z coordinate of half of the points from 0 into 1.
 UPDATE a SET geom=ST_SetSRID(ST_MakePoint(ST_x(geom), ST_y(geom), 1),
 4326) WHERE MOD(id, 2)=0;

 -- Ensure up-to-date statistics.
 VACUUM ANALYZE;

 -- Finally, run the query below.
 SELECT COUNT(*) FROM a, b WHERE ST_Within(a.geom, b.geom);
 }}}

 === Results
 {{{
 postgis/postgis:11-2.5-alpine: 100
 postgis/postgis:11-3.0-alpine: 100
 postgis/postgis:11-3.1-alpine: 100
 postgis/postgis:12-2.5-alpine: 100
 postgis/postgis:12-3.0-alpine:  50
 postgis/postgis:12-3.1-alpine:  50
 postgis/postgis:13-3.0-alpine:  50
 postgis/postgis:13-3.1-alpine:  50
 postgis/postgis:14-3.1-alpine:  50
 }}}

 1. Clearly, 50 and 100 is the result of a 3D and 2D calculation,
 respectively. For a user, however, it's unclear what calculation to
 expect. Documentation on ST_Within does not help.

 2. Even more confusing is that the result also depends on the indexes
 present. For example, dropping index a_3d_idx in 14-3.1-alpine, changes
 the result from 50 into 100. I expect the presence of an index to affect
 the performance of a query, not its result.

-- 
Ticket URL: <https://trac.osgeo.org/postgis/ticket/5025>
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.