Skip to content
Permalink
Browse files

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

…39453)

(cherry picked from commit b936829)
  • Loading branch information
rouault authored and nyalldawson committed Feb 19, 2021
1 parent 1805d16 commit 32867532dd477c207207332d6b0f59ae84d50cfa
@@ -437,11 +437,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
{
@@ -450,7 +451,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 ||
@@ -461,6 +462,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
@@ -488,6 +562,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;
}

@@ -2528,6 +2528,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):

@@ -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="
@@ -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 );

0 comments on commit 3286753

Please sign in to comment.