[Backport 7.6] PostGIS: fix ST_Intersects() with bounding box that is a point

[rouault]

(follow-up of fixes #6181, fixes #6230) (fixes #6347 (comment))

Backport of PR #6354

[sdlime]

Just a FYI. Ran into another regression that seems like it might be related. Will create a ticket when I can isolate the problem more...

Ok, a little more. In my case the "shape" that ultimately gets intersected with the bounding box is created using ST_Collect(), so there's a "GROUP BY" and then an aggregation using ST_Collect().

The new bounding box check ST_Intersects("shape", ST_GeomFromText('POLYGON((125000 4785000,125000 5489000,789000 5489000,789000 4785000,125000 4785000))',26915)) fails against that geometry with an error like:

ERROR:  Relate Operation called with a LWGEOMCOLLECTION type.  This is unsupported.
HINT:  Change argument 1: 'GEOMETRYCOLLECTION(MULTILINESTRING((515124.6335 5277837.6223,515130.0913 5277...'
SQL state: XX000

The old bounding box check "shape" && ST_GeomFromText('POLYGON((125000 4785000,125000 5489000,789000 5489000,789000 4785000,125000 4785000))',26915) works as expected.

This might not be something that can be addressed in the code but I suspect that other users might be using aggregation functions that might break as a result of the change.

--Steve

[sdlime]

Note that if I switch ST_Collect() to ST_Union() the error goes away - but at a performance cost.

[jmckenna]

Should we wait to release 7.6.4, for this fix?

[sdlime]

Should we wait to release 7.6.4, for this fix?

I don't know, but it is a regression. The issue, more generally, is that a GeometryCollection fails. In my case I'm just doing a query to generate a light JSON response so I don't use the resulting geometry for anything (e.g. drawing or outputting geometries) so the actual type (GeometryCollection vs MultiLinestring) doesn't matter. We're using PostGIS 2.4 and it looks like this is no longer an issue in version 2.5.

Perhaps a note to -dev to see what folks think.

--Steve

[rouault]

@sdlime Does the following fix your issue ? (it is likely not sufficient if the geometry collection contains heterogenous single geometry types)

diff --git a/mappostgis.cpp b/mappostgis.cpp
index 0ce749a04..342bcfd4e 100644
--- a/mappostgis.cpp
+++ b/mappostgis.cpp
@@ -1902,9 +1902,9 @@ static std::string msPostGISBuildSQLWhere(layerObj *layer, const rectObj *rect,
         // otherwise if find_srid() would return 0, ST_Intersects() would not
         // work at all, which breaks the msautotest/query/query_postgis.map
         // tests, releated to bdry_counpy2 layer that has no SRID
-        strRect = "ST_Intersects(\"";
+        strRect = "ST_Intersects(ST_CollectionHomogenize(\"";
         strRect += layerinfo->geomcolumn;
-        strRect += "\", ";
+        strRect += "\"), ";
         strRect += strBox;
         strRect += ')';
     }

[rouault]

That said I'm concerned it might prevent PostGIS from using its spatial index... We should probably check if the geometry column is a true geometry column, and if so use ST_Intersects(), and if not probably add a && filter in addition

[rouault]

Could @pramsey potentially advise what would be a robust & efficient strategy for calling ST_Intersects( some_geometry_that_might_be_a_column_or_the_result_of_a_query, ST_GeomFromText(...)) as a spatial filter, if such strategy exists :-) ?

[pramsey]

Your intuition that wrapping the column in any ST_Function() is going to invalidate the index is correct, so don't do that.

The lack of support of geometrycollection in ST_Intersects in some versions is the key issue, and it might be the simplest thing is to just go with && instead of ST_Intersects. You'll get an over-determined result set, however, and that might be another regression.

This has other knock-on performance impacts, but ST_DWithin(a, b, 0) is logically equivalent to ST_Interects(a, b) and doesn't have the GEOS-related "no collections" limitation, and also is index assisted. Even better might be a && b AND ST_Distance(a, b) = 0 as that results in a simpler plan than what happens under the covers of ST_DWithin()

[sdlime]

What about && if < 2.5 and ST_Intersects if >= 2.5?

[pramsey]

Then I'd have to be definitive about what version the behavior changed in :) maybe people with older things running can quickly crowdsource if this runs

select st_intersects('POINT(1 1)', 'GEOMETRYCOLLECTION(LINESTRING(2 2, 3 3), POINT(1 1))');

[pramsey]

Crowd confirms that collections work >= 2.5.