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.

      Note that for PostGIS geometry columns of type GEOGRAPHY, the current
      implementation performs the intersection test on the client side, by dealing
      with geographies as if they were cartesian geometries.

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,23 @@ 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.
+
+      Note that for PostGIS geometry columns of type GEOGRAPHY, the current
+      implementation performs the intersection test on the client side, by dealing
+      with geographies as if they were cartesian geometries.
+
+
 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)