From 48a11f0aefdc91a07cb94014129b06269abf434d Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 1 Feb 2021 16:25:35 +0100 Subject: [PATCH] [Postgres] Use spatial index for tables using geography type (fixes #39453) --- .../postgres/qgspostgresfeatureiterator.cpp | 79 ++++++++++++++++++- tests/src/python/test_provider_postgres.py | 32 ++++++++ tests/testdata/provider/testdata_pg.sh | 1 + .../provider/testdata_pg_geography.sql | 14 ++++ 4 files changed, 124 insertions(+), 2 deletions(-) create mode 100644 tests/testdata/provider/testdata_pg_geography.sql diff --git a/src/providers/postgres/qgspostgresfeatureiterator.cpp b/src/providers/postgres/qgspostgresfeatureiterator.cpp index dac2fce580cd..1c632abc2ad8 100644 --- a/src/providers/postgres/qgspostgresfeatureiterator.cpp +++ b/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 + #include + + 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; } diff --git a/tests/src/python/test_provider_postgres.py b/tests/src/python/test_provider_postgres.py index 17c62791a5a6..e5efc3a8da89 100644 --- a/tests/src/python/test_provider_postgres.py +++ b/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): diff --git a/tests/testdata/provider/testdata_pg.sh b/tests/testdata/provider/testdata_pg.sh index 1ff2cf16efa5..96f59d65eb89 100755 --- a/tests/testdata/provider/testdata_pg.sh +++ b/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=" diff --git a/tests/testdata/provider/testdata_pg_geography.sql b/tests/testdata/provider/testdata_pg_geography.sql new file mode 100644 index 000000000000..93b2f1346689 --- /dev/null +++ b/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 );