PostGIS: use full geometry intersection for SetSpatialFilter(), and not only bounding box

Author: rouault

autotest/ogr/ogr_pg.py

@@ -6376,3 +6376,37 @@ def test_ogr_pg_field_truncation(pg_ds):
     lyr.ResetReading()
     f = lyr.GetNextFeature()
     assert f["field"] == b"abcd\xc3\xa9".decode("UTF-8")
+
+
+###############################################################################
+# Test real geometry intersection in spatial filter
+
+
+@gdaltest.enable_exceptions()
+@pytest.mark.parametrize("GEOM_TYPE", ["geometry", "geography"])
+def test_ogr_pg_geometry_intersection_spatial_filter(pg_ds, use_postgis, GEOM_TYPE):
+
+    if use_postgis:
+        srs = osr.SpatialReference(epsg=4326)
+        lyr = pg_ds.CreateLayer("test", srs, options=["GEOM_TYPE=" + GEOM_TYPE])
+    else:
+        lyr = pg_ds.CreateLayer("test")
+    f = ogr.Feature(lyr.GetLayerDefn())
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POINT (0.5 0.5)"))
+    lyr.CreateFeature(f)
+    f = ogr.Feature(lyr.GetLayerDefn())
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POINT (0.9 0.5)"))
+    lyr.CreateFeature(f)
+
+    lyr.ResetReading()
+
+    # Pacman style polygon: almost a square except it doesn't include 0.9, 0.5
+    lyr.SetSpatialFilter(
+        ogr.CreateGeometryFromWkt(
+            "POLYGON ((0 0,0 1,1 1,1 0.6,0.8 0.6,0.8 0.4,1 0.4,1 0,0 0))"
+        )
+    )
+    f = lyr.GetNextFeature()
+    assert f.GetGeometryRef().ExportToWkt() == "POINT (0.5 0.5)"
+    f = lyr.GetNextFeature()
+    assert f is None

ogr/ogrsf_frmts/pg/ogrpgtablelayer.cpp

@@ -1077,10 +1077,24 @@ void OGRPGTableLayer::BuildWhere()
                     sEnvelope.MinY);
         CPLsnprintf(szBox3D_2, sizeof(szBox3D_2), "%.17g %.17g", sEnvelope.MaxX,
                     sEnvelope.MaxY);
-        osWHERE.Printf(
-            "WHERE %s && ST_SetSRID('BOX3D(%s, %s)'::box3d,%d) ",
-            OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef()).c_str(),
-            szBox3D_1, szBox3D_2, poGeomFieldDefn->nSRSId);
+        if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY)
+        {
+            OGRWktOptions sOptions;
+            sOptions.variant = wkbVariantIso;
+            const std::string osWKT = m_poFilterGeom->exportToWkt(sOptions);
+            osWHERE.Printf(
+                "WHERE ST_Intersects(%s, ST_SetSRID(ST_GeomFromText('%s'), "
+                "%d)) ",
+                OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef()).c_str(),
+                osWKT.c_str(), poGeomFieldDefn->nSRSId);
+        }
+        else
+        {
+            osWHERE.Printf(
+                "WHERE %s && ST_SetSRID('BOX3D(%s, %s)'::box3d,%d) ",
+                OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef()).c_str(),
+                szBox3D_1, szBox3D_2, poGeomFieldDefn->nSRSId);
+        }
     }
 
     if (!osQuery.empty())
@@ -1174,7 +1188,6 @@ OGRFeature *OGRPGTableLayer::GetNextFeature()
          * request */
         if (m_poFilterGeom == nullptr || poGeomFieldDefn == nullptr ||
             poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY ||
-            poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY ||
             FilterGeometry(poFeature->GetGeomFieldRef(m_iGeomFieldFilter)))
         {
             if (iFIDAsRegularColumnIndex >= 0)