QgsCoordinateTransform::transformBoundingBox(): avoid excessive north…

…ings when reprojecting world extent from EPSG:4326 to EPSG:3857 Refs #34518
qgis · Jan 24, 2022 · 04a94d2 · 04a94d2
1 parent 1f1a541
commit 04a94d2
Show file tree

Hide file tree

Showing 4 changed files with 47 additions and 4 deletions.
diff --git a/src/core/proj/qgscoordinatetransform.cpp b/src/core/proj/qgscoordinatetransform.cpp
@@ -514,14 +514,40 @@ QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle &r
     return QgsRectangle( p, p );
   }
 
+  double yMin = rect.yMinimum();
+  double yMax = rect.yMaximum();
+  if ( d->mGeographicToWebMercator &&
+       ( ( direction == Qgis::TransformDirection::Forward && !d->mIsReversed ) ||
+         ( direction == Qgis::TransformDirection::Reverse && d->mIsReversed ) ) )
+  {
+    // Latitudes close to 90 degree project to infinite northing in theory.
+    // We limit to 90 - 1e-1 which reproject to northing of ~ 44e6 m (about twice
+    // the maximum easting of ~20e6 m).
+    // For reference, GoogleMercator tiles are limited to a northing ~85 deg / ~20e6 m
+    // so limiting to 90 - 1e-1 is reasonable.
+    constexpr double EPS = 1e-1;
+    if ( yMin < -90 + EPS )
+    {
+      if ( yMax < -90 + EPS )
+        throw QgsCsException( QObject::tr( "Could not transform bounding box to target CRS" ) );
+      yMin = -90 + EPS;
+    }
+    if ( yMax > 90 - EPS )
+    {
+      if ( yMin > 90 - EPS )
+        throw QgsCsException( QObject::tr( "Could not transform bounding box to target CRS" ) );
+      yMax = 90 - EPS;
+    }
+  }
+
   // 64 points (<=2.12) is not enough, see #13665, for EPSG:4326 -> EPSG:3574 (say that it is a hard one),
   // are decent result from about 500 points and more. This method is called quite often, but
   // even with 1000 points it takes < 1ms.
   // TODO: how to effectively and precisely reproject bounding box?
   const int nPoints = 1000;
-  const double d = std::sqrt( ( rect.width() * rect.height() ) / std::pow( std::sqrt( static_cast< double >( nPoints ) ) - 1, 2.0 ) );
+  const double d = std::sqrt( ( rect.width() * ( yMax - yMin ) ) / std::pow( std::sqrt( static_cast< double >( nPoints ) ) - 1, 2.0 ) );
   const int nXPoints = std::min( static_cast< int >( std::ceil( rect.width() / d ) ) + 1, 1000 );
-  const int nYPoints = std::min( static_cast< int >( std::ceil( rect.height() / d ) ) + 1, 1000 );
+  const int nYPoints = std::min( static_cast< int >( std::ceil( ( yMax - yMin ) / d ) ) + 1, 1000 );
 
   QgsRectangle bb_rect;
   bb_rect.setMinimal();
@@ -537,9 +563,9 @@ QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle &r
   // Populate the vectors
 
   const double dx = rect.width()  / static_cast< double >( nXPoints - 1 );
-  const double dy = rect.height() / static_cast< double >( nYPoints - 1 );
+  const double dy = ( yMax - yMin ) / static_cast< double >( nYPoints - 1 );
 
-  double pointY = rect.yMinimum();
+  double pointY = yMin;
 
   for ( int i = 0; i < nYPoints ; i++ )
   {

diff --git a/src/core/proj/qgscoordinatetransform_p.cpp b/src/core/proj/qgscoordinatetransform_p.cpp
@@ -83,6 +83,7 @@ QgsCoordinateTransformPrivate::QgsCoordinateTransformPrivate( const QgsCoordinat
   , mAvailableOpCount( other.mAvailableOpCount )
   , mIsValid( other.mIsValid )
   , mShortCircuit( other.mShortCircuit )
+  , mGeographicToWebMercator( other.mGeographicToWebMercator )
   , mSourceCRS( other.mSourceCRS )
   , mDestCRS( other.mDestCRS )
   , mSourceDatumTransform( other.mSourceDatumTransform )
@@ -158,6 +159,10 @@ bool QgsCoordinateTransformPrivate::initialize()
     return true;
   }
 
+  mGeographicToWebMercator =
+    mSourceCRS.authid() == QLatin1String( "EPSG:4326" ) &&
+    mDestCRS.authid() == QLatin1String( "EPSG:3857" );
+
   mSourceIsDynamic = mSourceCRS.isDynamic();
   mSourceCoordinateEpoch = mSourceCRS.coordinateEpoch();
   mDestIsDynamic = mDestCRS.isDynamic();

diff --git a/src/core/proj/qgscoordinatetransform_p.h b/src/core/proj/qgscoordinatetransform_p.h
@@ -88,6 +88,9 @@ class QgsCoordinateTransformPrivate : public QSharedData
      */
     bool mShortCircuit = false;
 
+    //! Flag to indicate EPSG:4326 to EPSG:3857 reprojection
+    bool mGeographicToWebMercator = false;
+
     //! QgsCoordinateReferenceSystem of the source (layer) coordinate system
     QgsCoordinateReferenceSystem mSourceCRS;
 

diff --git a/tests/src/python/test_qgscoordinatetransform.py b/tests/src/python/test_qgscoordinatetransform.py
@@ -153,6 +153,15 @@ def testProjectContextProj6(self):
         transform = QgsCoordinateTransform(QgsCoordinateReferenceSystem('EPSG:28356'), QgsCoordinateReferenceSystem('EPSG:3111'), p)
         self.assertEqual(transform.coordinateOperation(), 'proj')
 
+    def testTransformBoundingBoxFullWorldToWebMercator(self):
+        extent = QgsRectangle(-180, -90, 180, 90)
+        transform = QgsCoordinateTransform(QgsCoordinateReferenceSystem('EPSG:4326'), QgsCoordinateReferenceSystem('EPSG:3857'), QgsProject.instance())
+        transformedExtent = transform.transformBoundingBox(extent)
+        self.assertAlmostEqual(transformedExtent.xMinimum(), -20037508.343, delta=1e-3)
+        self.assertAlmostEqual(transformedExtent.yMinimum(), -44927335.427, delta=1e-3)
+        self.assertAlmostEqual(transformedExtent.xMaximum(), 20037508.343, delta=1e-3)
+        self.assertAlmostEqual(transformedExtent.yMaximum(), 44927335.427, delta=1e-3)
+
 
 if __name__ == '__main__':
     unittest.main()