PG: add a SPATIAL_FILTER_INTERSECTION=LOCAL/SERVER open option

-  .. oo:: SPATIAL_FILTER_INTERSECTION
      :choices: LOCAL, SERVER
      :default: LOCAL
      :since: 3.13

      Since GDAL 3.13, spatial filters are evaluated using full geometry intersection
      rather than only bounding box intersection, as in earlier versions.
      The ``SPATIAL_FILTER_INTERSECTION`` open option controls where this
      intersection operation is performed. Even when the intersection is evaluated
      locally (which is the default), the server still applies a bounding
      box-based spatial filter to efficiently retrieve candidate features.

Author: rouault

autotest/ogr/ogr_pg.py

@@ -6384,13 +6384,19 @@ def test_ogr_pg_field_truncation(pg_ds):
 
 @gdaltest.enable_exceptions()
 @pytest.mark.parametrize("GEOM_TYPE", ["geometry", "geography"])
-def test_ogr_pg_geometry_intersection_spatial_filter(pg_ds, use_postgis, GEOM_TYPE):
+@pytest.mark.parametrize("open_options", [None, ["SPATIAL_FILTER_INTERSECTION=SERVER"]])
+@pytest.mark.require_geos
+def test_ogr_pg_geometry_intersection_spatial_filter(
+    pg_ds, use_postgis, GEOM_TYPE, open_options
+):
+
+    ds = reconnect(pg_ds, open_options=open_options)
 
     if use_postgis:
         srs = osr.SpatialReference(epsg=4326)
-        lyr = pg_ds.CreateLayer("test", srs, options=["GEOM_TYPE=" + GEOM_TYPE])
+        lyr = ds.CreateLayer("test", srs, options=["GEOM_TYPE=" + GEOM_TYPE])
     else:
-        lyr = pg_ds.CreateLayer("test")
+        lyr = ds.CreateLayer("test")
     f = ogr.Feature(lyr.GetLayerDefn())
     f.SetGeometry(ogr.CreateGeometryFromWkt("POINT (0.5 0.5)"))
     lyr.CreateFeature(f)
@@ -6413,7 +6419,7 @@ def test_ogr_pg_geometry_intersection_spatial_filter(pg_ds, use_postgis, GEOM_TY
 
     lyr.SetSpatialFilter(None)
 
-    with pg_ds.ExecuteSQL("SELECT * FROM test") as sql_lyr:
+    with ds.ExecuteSQL("SELECT * FROM test") as sql_lyr:
         sql_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))"

doc/source/drivers/vector/pg.rst

@@ -243,6 +243,19 @@ The following open options are supported:
       :oo::`PRELUDE_STATEMENTS`, the appropriate CLOSING_STATEMENTS would be
       "COMMIT".
 
+-  .. oo:: SPATIAL_FILTER_INTERSECTION
+      :choices: LOCAL, SERVER
+      :default: LOCAL
+      :since: 3.13
+
+      Since GDAL 3.13, spatial filters are evaluated using full geometry intersection
+      rather than only bounding box intersection, as in earlier versions.
+      The ``SPATIAL_FILTER_INTERSECTION`` open option controls where this
+      intersection operation is performed. Even when the intersection is evaluated
+      locally (which is the default), the server still applies a bounding
+      box-based spatial filter to efficiently retrieve candidate features.
+
+
 Dataset Creation Options
 ~~~~~~~~~~~~~~~~~~~~~~~~
 

ogr/ogrsf_frmts/pg/ogr_pg.h

@@ -604,6 +604,8 @@ class OGRPGDataSource final : public GDALDataset
     bool m_bHasWritePermissionsOnMetadataTableRun = false;
     bool m_bHasWritePermissionsOnMetadataTableSuccess = false;
 
+    bool m_bSpatialFilterIntersectionIsLocal = true;
+
     void LoadTables();
 
     CPLString osDebugLastTransactionCommand{};
@@ -714,6 +716,11 @@ class OGRPGDataSource final : public GDALDataset
     bool CreateMetadataTableIfNeeded();
     bool HasOgrSystemTablesMetadataTable();
     bool HasWritePermissionsOnMetadataTable();
+
+    bool IsSpatialFilterIntersectionLocal() const
+    {
+        return m_bSpatialFilterIntersectionIsLocal;
+    }
 };
 
 #endif /* ndef OGR_PG_H_INCLUDED */

ogr/ogrsf_frmts/pg/ogrpgdatasource.cpp

@@ -261,6 +261,11 @@ int OGRPGDataSource::Open(const char *pszNewName, int bUpdate, int bTestOpen,
     CPLAssert(nLayers == 0);
     papszOpenOptions = CSLDuplicate(papszOpenOptionsIn);
 
+    m_bSpatialFilterIntersectionIsLocal =
+        EQUAL(CSLFetchNameValueDef(papszOpenOptions,
+                                   "SPATIAL_FILTER_INTERSECTION", "LOCAL"),
+              "LOCAL");
+
     const char *pszPreludeStatements =
         CSLFetchNameValue(papszOpenOptions, "PRELUDE_STATEMENTS");
     if (pszPreludeStatements)

ogr/ogrsf_frmts/pg/ogrpgdrivercore.cpp

@@ -79,6 +79,12 @@ void OGRPGDriverSetCommonMetadata(GDALDriver *poDriver)
         "  <Option name='SKIP_VIEWS' type='boolean' description='Whether "
         "views should be omitted from the list' "
         "default='NO'/>"
+        "  <Option name='SPATIAL_FILTER_INTERSECTION' type='string-select' "
+        "description='Whether geometry intersection of spatial filter is "
+        "evaluated locally or on the server' default='LOCAL'>"
+        "    <Value>LOCAL</Value>"
+        "    <Value>DATABASE</Value>"
+        "  </Option>"
         "</OpenOptionList>");
 
     poDriver->SetMetadataItem(GDAL_DMD_CREATIONOPTIONLIST,

ogr/ogrsf_frmts/pg/ogrpgresultlayer.cpp

@@ -275,7 +275,8 @@ OGRFeature *OGRPGResultLayer::GetNextFeature()
 
         if ((m_poFilterGeom == nullptr || poGeomFieldDefn == nullptr ||
              (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY &&
-              poGeomFieldDefn->nSRSId > 0) ||
+              poGeomFieldDefn->nSRSId > 0 &&
+              !poDS->IsSpatialFilterIntersectionLocal()) ||
              FilterGeometry(poFeature->GetGeomFieldRef(m_iGeomFieldFilter))) &&
             (m_poAttrQuery == nullptr || m_poAttrQuery->Evaluate(poFeature)))
             return poFeature;
@@ -305,7 +306,8 @@ OGRErr OGRPGResultLayer::ISetSpatialFilter(int iGeomField,
             if (m_poFilterGeom != nullptr && poDS->sPostGISVersion.nMajor >= 0)
             {
                 if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY &&
-                    poGeomFieldDefn->nSRSId > 0)
+                    poGeomFieldDefn->nSRSId > 0 &&
+                    !poDS->IsSpatialFilterIntersectionLocal())
                 {
                     char *pszHexEWKB = OGRGeometryToHexEWKB(
                         m_poFilterGeom, poGeomFieldDefn->nSRSId,

ogr/ogrsf_frmts/pg/ogrpgtablelayer.cpp

@@ -1060,7 +1060,8 @@ void OGRPGTableLayer::BuildWhere()
          poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY))
     {
         poGeomFieldDefn->GetSpatialRef();  // make sure nSRSId is resolved
-        if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY)
+        if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY &&
+            !poDS->IsSpatialFilterIntersectionLocal())
         {
             char *pszHexEWKB = OGRGeometryToHexEWKB(
                 m_poFilterGeom, poGeomFieldDefn->nSRSId,
@@ -1188,12 +1189,11 @@ OGRFeature *OGRPGTableLayer::GetNextFeature()
             return nullptr;
 
         /* We just have to look if there is a geometry filter */
-        /* If there's a PostGIS geometry column, the spatial filter */
-        /* is already taken into account in the select request */
         /* The attribute filter is always taken into account by the select
          * request */
         if (m_poFilterGeom == nullptr || poGeomFieldDefn == nullptr ||
-            poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY ||
+            (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY &&
+             !poDS->IsSpatialFilterIntersectionLocal()) ||
             FilterGeometry(poFeature->GetGeomFieldRef(m_iGeomFieldFilter)))
         {
             if (iFIDAsRegularColumnIndex >= 0)