[Postgres] Use spatial index for tables using geography type (fixes #39453)

src/providers/postgres/qgspostgresfeatureiterator.cpp

@@ -430,11 +430,12 @@ QString QgsPostgresFeatureIterator::whereClauseRect()
   }
 
   QString qBox;
+  const QString bboxSrid = mSource->mRequestedSrid.isEmpty() ? mSource->mDetectedSrid : mSource->mRequestedSrid;
   if ( mConn->majorVersion() < 2 )
   {
     qBox = QStringLiteral( "setsrid('BOX3D(%1)'::box3d,%2)" )
            .arg( rect.asWktCoordinates(),
-                 mSource->mRequestedSrid.isEmpty() ? mSource->mDetectedSrid : mSource->mRequestedSrid );
+                 bboxSrid );
   }
   else
   {
@@ -443,7 +444,7 @@ QString QgsPostgresFeatureIterator::whereClauseRect()
                  qgsDoubleToString( rect.yMinimum() ),
                  qgsDoubleToString( rect.xMaximum() ),
                  qgsDoubleToString( rect.yMaximum() ),
-                 mSource->mRequestedSrid.isEmpty() ? mSource->mDetectedSrid : mSource->mRequestedSrid );
+                 bboxSrid );
   }
 
   bool castToGeometry = mSource->mSpatialColType == SctGeography ||
@@ -454,6 +455,79 @@ QString QgsPostgresFeatureIterator::whereClauseRect()
                               castToGeometry ? "::geometry" : "",
                               qBox );
 
+  // For geography type, using a && filter with the geography column cast as
+  // geometry prevents the use of a spatial index. So for "small" filtering
+  // bounding boxes, use the filtering in the geography space. But a bbox in
+  // geography uses geodesic arcs, and not rhumb lines, so we must expand a bit the QGIS
+  // bounding box (which assumes rhumb lines) to the south in the northern hemisphere
+  // See https://trac.osgeo.org/postgis/ticket/2495 for some background
+  if ( mConn->majorVersion() >= 2 &&
+       mSource->mSpatialColType == SctGeography &&
+       bboxSrid == QLatin1String( "4326" ) &&
+       std::fabs( rect.yMaximum()  - rect.yMinimum() ) <= 10 &&
+       std::fabs( rect.xMaximum()  - rect.xMinimum() ) <= 10 &&
+       std::fabs( rect.yMaximum() ) <= 70 )
+  {
+    /* The following magic constant has been obtained by running :
+        #include "geodesic.h"
+        #include <math.h>
+        #include <stdio.h>
+
+        int main()
+        {
+            struct geod_geodesic g;
+            struct geod_geodesicline l;
+            int i;
+            geod_init(&g, 6378137, 1/298.257223563);
+            double lat1 = 60;
+            double lon1 = 0;
+            double lat2 = 70;
+            double lon2 = 10;
+            geod_inverseline(&l, &g, lat1, lon1, lat2, lon2, 0);
+            double maxdlat = 0;
+            for (i = 0; i <= 100; ++i)
+            {
+                double lat, lon;
+                geod_position(&l, i * l.s13 * 0.01, &lat, &lon, 0);
+                double alpha = (lon - lon1) / (lon2 - lon1);
+                double lat_rhumb = lat1 + (lat2 - lat1) * alpha;
+                double dlat = lat - lat_rhumb;
+                if( fabs(dlat) > fabs(maxdlat) )
+                    maxdlat = dlat;
+                //printf("%f: %f %f delta=%f\n", lon, lat, lat_rhumb, dlat);
+            }
+            printf("maxdlat = %f\n", maxdlat);
+            return 0;
+        }
+    */
+    // And noticing that the difference between the rhumb line and the geodesics
+    // increases with higher latitude maximum, and differences of longitude and latitude.
+    // We could perhaps get a formula that would give dlat as a function of
+    // delta_lon, delta_lat and max(lat), but those maximum values should be good
+    // enough for now.
+    double dlat = 1.04;
+
+    // For smaller filtering bounding box, use smaller bbox expansion
+    if ( std::fabs( rect.yMaximum()  - rect.yMinimum() ) <= 1 &&
+         std::fabs( rect.xMaximum()  - rect.xMinimum() ) <= 1 )
+    {
+      // Value got by changing lat1 to 69 and lon2 to 1 in the above code snippet
+      dlat = 0.013;
+    }
+    // In the northern hemisphere, extends the geog bbox to the south
+    const double yminGeog = rect.yMinimum() >= 0 ? std::max( 0.0, rect.yMinimum() - dlat ) : rect.yMinimum();
+    const double ymaxGeog = rect.yMaximum() >= 0 ? rect.yMaximum() : std::min( 0.0, rect.yMaximum() + dlat );
+    const QString qBoxGeog = QStringLiteral( "st_makeenvelope(%1,%2,%3,%4,%5)" )
+                             .arg( qgsDoubleToString( rect.xMinimum() ),
+                                   qgsDoubleToString( yminGeog ),
+                                   qgsDoubleToString( rect.xMaximum() ),
+                                   qgsDoubleToString( ymaxGeog ),
+                                   bboxSrid );
+    whereClause += QStringLiteral( " AND %1 && %2" )
+                   .arg( QgsPostgresConn::quotedIdentifier( mSource->mBoundingBoxColumn ),
+                         qBoxGeog );
+  }
+
   if ( mRequest.flags() & QgsFeatureRequest::ExactIntersect )
   {
     QString curveToLineFn; // in PostGIS < 1.5 the st_curvetoline function does not exist
@@ -481,6 +555,7 @@ QString QgsPostgresFeatureIterator::whereClauseRect()
     whereClause += QStringLiteral( " AND %1" ).arg( QgsPostgresConn::postgisTypeFilter( mSource->mGeometryColumn, ( QgsWkbTypes::Type )mSource->mRequestedGeomType, castToGeometry ) );
   }
 
+  QgsDebugMsgLevel( QStringLiteral( "whereClause = %1" ).arg( whereClause ), 4 );
   return whereClause;
 }
 

tests/src/python/test_provider_postgres.py

@@ -2530,6 +2530,38 @@ def testHasSpatialIndex(self):
                 self.assertTrue(vl.isValid())
                 self.assertEqual(vl.hasSpatialIndex(), spatial_index)
 
+    def testBBoxFilterOnGeographyType(self):
+        """Test bounding box filter on geography type"""
+
+        vl = QgsVectorLayer(
+            self.dbconn +
+            ' sslmode=disable key=\'pk\' srid=4326 type=POINT table="qgis_test"."testgeog" (geog) sql=',
+            'test', 'postgres')
+
+        self.assertTrue(vl.isValid())
+
+        def _test(vl, extent, ids):
+            request = QgsFeatureRequest().setFilterRect(extent)
+            values = {feat['pk']: 'x' for feat in vl.getFeatures(request)}
+            expected = {x: 'x' for x in ids}
+            self.assertEqual(values, expected)
+
+        _test(vl, QgsRectangle(40 - 0.01, -0.01, 40 + 0.01, 0.01), [1])
+        _test(vl, QgsRectangle(40 - 5, -5, 40 + 5, 5), [1])
+        _test(vl, QgsRectangle(40 - 5, 0, 40 + 5, 5), [1])
+        _test(vl, QgsRectangle(40 - 10, -10, 40 + 10, 10), [1])  # no use of spatial index currently
+        _test(vl, QgsRectangle(40 - 5, 0.01, 40 + 5, 5), [])  # no match
+
+        _test(vl, QgsRectangle(40 - 0.01, 60 - 0.01, 40 + 0.01, 60 + 0.01), [2])
+        _test(vl, QgsRectangle(40 - 5, 60 - 5, 40 + 5, 60 + 5), [2])
+        _test(vl, QgsRectangle(40 - 5, 60 - 0.01, 40 + 5, 60 + 9.99), [2])
+
+        _test(vl, QgsRectangle(40 - 0.01, -60 - 0.01, 40 + 0.01, -60 + 0.01), [3])
+        _test(vl, QgsRectangle(40 - 5, -60 - 5, 40 + 5, -60 + 5), [3])
+        _test(vl, QgsRectangle(40 - 5, -60 - 9.99, 40 + 5, -60 + 0.01), [3])
+
+        _test(vl, QgsRectangle(-181, -90, 181, 90), [1, 2, 3])  # no use of spatial index currently
+
 
 class TestPyQgsPostgresProviderCompoundKey(unittest.TestCase, ProviderTestCase):
 

tests/testdata/provider/testdata_pg.sh

@@ -17,6 +17,7 @@ SCRIPTS="
   tests/testdata/provider/testdata_pg_pointcloud.sql
   tests/testdata/provider/testdata_pg_bigint_pk.sql
   tests/testdata/provider/testdata_pg_hasspatialindex.sql
+  tests/testdata/provider/testdata_pg_geography.sql
 "
 
 SCRIPTS12="

tests/testdata/provider/testdata_pg_geography.sql

@@ -0,0 +1,14 @@
+
+DROP TABLE IF EXISTS qgis_test.testgeog CASCADE;
+
+CREATE TABLE qgis_test.testgeog
+(
+  pk SERIAL NOT NULL PRIMARY KEY,
+  geog GEOGRAPHY
+);
+
+INSERT INTO qgis_test.testgeog(geog) VALUES ('POINT(40 0)'::geography);
+INSERT INTO qgis_test.testgeog(geog) VALUES ('POINT(40 60)'::geography);
+INSERT INTO qgis_test.testgeog(geog) VALUES ('POINT(40 -60)'::geography);
+
+CREATE INDEX testgeog_geog_idx ON qgis_test.testgeog USING GIST ( geog );