[Backport release-3_40] Delegate bounding box transformation to proj on >= 8.2

python/PyQt6/core/auto_generated/proj/qgscoordinatetransform.sip.in

@@ -280,6 +280,10 @@ the returned rectangle.
 
 :return: rectangle in destination CRS
 
+.. warning::
+
+   Do not call this method if the transformation involves geocentric CRS -- in this situation transformation of a 2D bounding box is meaningless! Calling this method with a geocentric CRS will result in a :py:class:`QgsCsException` being thrown.
+
 :raises QgsCsException: if the transformation fails
 %End
 

python/core/auto_generated/proj/qgscoordinatetransform.sip.in

@@ -280,6 +280,10 @@ the returned rectangle.
 
 :return: rectangle in destination CRS
 
+.. warning::
+
+   Do not call this method if the transformation involves geocentric CRS -- in this situation transformation of a 2D bounding box is meaningless! Calling this method with a geocentric CRS will result in a :py:class:`QgsCsException` being thrown.
+
 :raises QgsCsException: if the transformation fails
 %End
 

python/plugins/processing/tests/GdalAlgorithmsRasterTest.py

@@ -3498,7 +3498,7 @@ def testRasterize(self):
                 ),
                 [
                     "gdal_rasterize",
-                    "-l polys2 -a id -ts 0.0 0.0 -te -1.000000001857055 -2.9999999963940835 10.000000000604244 5.99999999960471 -ot Float32 -of JPEG "
+                    "-l polys2 -a id -ts 0.0 0.0 -te -1.000000001857055 -2.9999999963940835 10.000000000604246 5.999999999604708 -ot Float32 -of JPEG "
                     + source
                     + " "
                     + outdir

src/core/maprenderer/qgsmaprendererjob.cpp

@@ -409,9 +409,9 @@ bool QgsMapRendererJob::reprojectToLayerExtent( const QgsMapLayer *ml, const Qgs
         extent = approxTransform.transformBoundingBox( extent, Qgis::TransformDirection::Reverse );
     }
   }
-  catch ( QgsCsException & )
+  catch ( QgsCsException &e )
   {
-    QgsDebugError( QStringLiteral( "Transform error caught" ) );
+    QgsDebugError( QStringLiteral( "Transform error caught: %1" ).arg( e.what() ) );
     extent = QgsRectangle( std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max() );
     r2 = QgsRectangle( std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max() );
     res = false;

src/core/proj/qgscoordinatetransform.cpp

@@ -594,6 +594,26 @@ QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle &r
     return QgsRectangle( p, p );
   }
 
+#ifdef QGISDEBUG
+  if ( !mHasContext )
+  {
+    QgsDebugMsgLevel( QStringLiteral( "No QgsCoordinateTransformContext context set for transform" ), 4 );
+  }
+#endif
+
+  // we can't calculate if transform involves a geocentric CRS. This is silly anyway,
+  // as transformation of a 2d bounding box makes no sense when a geocentric CRS is involved!
+  if ( d->mSourceCRS.type() == Qgis::CrsType::Geocentric )
+  {
+    throw QgsCsException( QObject::tr( "Could not transform bounding box for geocentric CRS %1" ).arg( d->mSourceCRS.authid() ) );
+  }
+  if ( d->mDestCRS.type() == Qgis::CrsType::Geocentric )
+  {
+    throw QgsCsException( QObject::tr( "Could not transform bounding box for geocentric CRS %1" ).arg( d->mDestCRS.authid() ) );
+  }
+
+  const double xMin = rect.xMinimum();
+  const double xMax = rect.xMaximum();
   double yMin = rect.yMinimum();
   double yMax = rect.yMaximum();
   if ( d->mGeographicToWebMercator &&
@@ -620,6 +640,109 @@ QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle &r
     }
   }
 
+  // delegate logic to proj if version >= 8.2 available
+#if PROJ_VERSION_MAJOR>8 || (PROJ_VERSION_MAJOR==8 && PROJ_VERSION_MINOR>=2)
+
+  QgsScopedProjSilentLogger errorLogger;
+
+  QgsDebugMsgLevel( QStringLiteral( "Entering transformBoundingBox..." ), 4 );
+
+  ProjData projData = d->threadLocalProjData();
+  PJ_CONTEXT *projContext = QgsProjContext::get();
+
+#if PROJ_VERSION_MAJOR< 9 || (PROJ_VERSION_MAJOR==9 && PROJ_VERSION_MINOR<6)
+  // if source or destination crs include vertical components, we need to demote them to
+  // 2d crs first, otherwise proj_trans_bounds fails on proj < 9.6 (see https://github.com/OSGeo/PROJ/pull/4333)
+
+  QgsProjUtils::proj_pj_unique_ptr srcCrs( proj_get_source_crs( projContext, projData ) );
+  QgsProjUtils::proj_pj_unique_ptr destCrs( proj_get_target_crs( projContext, projData ) );
+
+  QgsProjUtils::proj_pj_unique_ptr srcCrsHorizontal;
+  QgsProjUtils::proj_pj_unique_ptr destCrsHorizontal;
+  QgsProjUtils::proj_pj_unique_ptr transform2D;
+  if ( QgsProjUtils::hasVerticalAxis( srcCrs.get() ) ||
+       QgsProjUtils::hasVerticalAxis( destCrs.get() ) )
+  {
+    srcCrsHorizontal = QgsProjUtils::crsToHorizontalCrs( srcCrs.get() );
+    destCrsHorizontal = QgsProjUtils::crsToHorizontalCrs( destCrs.get() );
+    transform2D.reset( proj_create_crs_to_crs_from_pj( projContext, srcCrsHorizontal.get(), destCrsHorizontal.get(), nullptr, nullptr ) );
+    projData = transform2D.get();
+  }
+#endif
+
+  double transXMin = 0;
+  double transYMin = 0;
+  double transXMax = 0;
+  double transYMax = 0;
+
+  proj_errno_reset( projData );
+  // proj documentation recommends 21 points for densification
+  constexpr int DENSIFY_POINTS = 21;
+  const int projResult = proj_trans_bounds( projContext, projData, ( direction == Qgis::TransformDirection::Forward && !d->mIsReversed ) || ( direction == Qgis::TransformDirection::Reverse && d->mIsReversed ) ? PJ_FWD : PJ_INV,
+                         xMin, yMin, xMax, yMax,
+                         &transXMin, &transYMin, &transXMax, &transYMax, DENSIFY_POINTS );
+  if ( projResult != 1
+       || !std::isfinite( transXMin )
+       || !std::isfinite( transXMax )
+       || !std::isfinite( transYMin )
+       || !std::isfinite( transYMax ) )
+  {
+    const QString projErr = QString::fromUtf8( proj_context_errno_string( projContext, proj_errno( projData ) ) );
+    const QString dir = ( direction == Qgis::TransformDirection::Forward ) ? QObject::tr( "Forward transform" ) : QObject::tr( "Inverse transform" );
+    const QString msg = QObject::tr( "%1 (%2 to %3) of bounding box failed: %4" )
+                        .arg( dir,
+                              ( direction == Qgis::TransformDirection::Forward ) ? d->mSourceCRS.authid() : d->mDestCRS.authid(),
+                              ( direction == Qgis::TransformDirection::Forward ) ? d->mDestCRS.authid() : d->mSourceCRS.authid(),
+                              projErr );
+    QgsDebugError( msg );
+
+    throw QgsCsException( msg );
+  }
+
+  // check if result bbox is geographic and is crossing 180/-180 line: ie. min X is before the 180° and max X is after the -180°
+  bool doHandle180Crossover = false;
+
+  if ( handle180Crossover
+       && ( ( direction == Qgis::TransformDirection::Forward && d->mDestCRS.isGeographic() ) ||
+            ( direction == Qgis::TransformDirection::Reverse && d->mSourceCRS.isGeographic() ) )
+       && ( transXMax < transXMin ) )
+  {
+    //if crossing the date line, temporarily add 360 degrees to -ve longitudes
+    std::swap( transXMax, transXMin );
+    if ( transXMin < 0 )
+      transXMin += 360;
+    if ( transXMax < 0 )
+      transXMax += 360;
+    doHandle180Crossover = true;
+  }
+
+  QgsRectangle boundingBoxRect{ transXMin, transYMin, transXMax, transYMax };
+  if ( boundingBoxRect.isNull() )
+  {
+    // something bad happened when reprojecting the filter rect... no finite points were left!
+    throw QgsCsException( QObject::tr( "Could not transform bounding box to target CRS" ) );
+  }
+
+  if ( doHandle180Crossover )
+  {
+    //subtract temporary addition of 360 degrees from longitudes
+    if ( boundingBoxRect.xMinimum() > 180.0 )
+      boundingBoxRect.setXMinimum( boundingBoxRect.xMinimum() - 360.0 );
+    if ( boundingBoxRect.xMaximum() > 180.0 )
+      boundingBoxRect.setXMaximum( boundingBoxRect.xMaximum() - 360.0 );
+  }
+
+  QgsDebugMsgLevel( "Projected extent: " + boundingBoxRect.toString(), 4 );
+
+  if ( boundingBoxRect.isEmpty() )
+  {
+    QgsDebugMsgLevel( "Original extent: " + rect.toString(), 4 );
+  }
+
+  return boundingBoxRect;
+#else
+  // this logic is buggy! See https://github.com/qgis/QGIS/issues/59821
+
   // 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.
@@ -649,7 +772,7 @@ QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle &r
   {
 
     // Start at right edge
-    double pointX = rect.xMinimum();
+    double pointX = xMin;
 
     for ( int j = 0; j < nXPoints; j++ )
     {
@@ -733,6 +856,7 @@ QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle &r
   }
 
   return bb_rect;
+#endif
 }
 
 void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *y, double *z, Qgis::TransformDirection direction ) const

src/core/proj/qgscoordinatetransform.h

@@ -257,6 +257,7 @@ class CORE_EXPORT QgsCoordinateTransform
      * \param handle180Crossover set to TRUE if destination CRS is geographic and handling of extents
      * crossing the 180 degree longitude line is required
      * \returns rectangle in destination CRS
+     * \warning Do not call this method if the transformation involves geocentric CRS -- in this situation transformation of a 2D bounding box is meaningless! Calling this method with a geocentric CRS will result in a QgsCsException being thrown.
      * \throws QgsCsException if the transformation fails
      */
     QgsRectangle transformBoundingBox( const QgsRectangle &rectangle, Qgis::TransformDirection direction = Qgis::TransformDirection::Forward, bool handle180Crossover = false ) const SIP_THROW( QgsCsException );

src/core/proj/qgsprojutils.cpp

@@ -279,6 +279,11 @@ bool QgsProjUtils::hasVerticalAxis( const PJ *crs )
       return false;
     }
 
+    case PJ_TYPE_BOUND_CRS:
+    {
+      return hasVerticalAxis( proj_get_source_crs( context, crs ) );
+    }
+
     // maybe other types to handle like this??
 
     default:
@@ -363,6 +368,10 @@ void QgsProjUtils::proj_collecting_logger( void *user_data, int /*level*/, const
   dest->append( messageString );
 }
 
+void QgsProjUtils::proj_silent_logger( void * /*user_data*/, int /*level*/, const char * /*message*/ )
+{
+}
+
 void QgsProjUtils::proj_logger( void *, int level, const char *message )
 {
 #ifdef QGISDEBUG
@@ -623,3 +632,18 @@ QgsScopedProjCollectingLogger::~QgsScopedProjCollectingLogger()
   // reset logger back to terminal output
   proj_log_func( QgsProjContext::get(), nullptr, QgsProjUtils::proj_logger );
 }
+
+//
+// QgsScopedProjSilentLogger
+//
+
+QgsScopedProjSilentLogger::QgsScopedProjSilentLogger()
+{
+  proj_log_func( QgsProjContext::get(), nullptr, QgsProjUtils::proj_silent_logger );
+}
+
+QgsScopedProjSilentLogger::~QgsScopedProjSilentLogger()
+{
+  // reset logger back to terminal output
+  proj_log_func( QgsProjContext::get(), nullptr, QgsProjUtils::proj_logger );
+}

src/core/proj/qgsprojutils.h

@@ -257,6 +257,11 @@ class CORE_EXPORT QgsProjUtils
      */
     static void proj_collecting_logger( void *user_data, int level, const char *message );
 
+    /**
+     * QGIS proj log function which ignores errors.
+     */
+    static void proj_silent_logger( void *user_data, int level, const char *message );
+
 #if 0 // not possible in current Proj 6 API
 
     /**
@@ -311,6 +316,38 @@ class CORE_EXPORT QgsProjContext
 #endif
 };
 
+/**
+ * \ingroup core
+ *
+ * \brief Scoped object for temporary suppression of PROJ logging output.
+ *
+ * Temporarily sets the PROJ log function to one which suppresses errors for the lifetime of the object,
+ * before returning it to the default QGIS proj logging function on destruction.
+ *
+ * \note The collecting logger ONLY applies to the current thread.
+ *
+ * \note Not available in Python bindings
+ * \since QGIS 3.42
+ */
+class CORE_EXPORT QgsScopedProjSilentLogger
+{
+  public:
+
+    /**
+     * Constructor for QgsScopedProjSilentLogger.
+     *
+     * PROJ errors will be ignored.
+     */
+    QgsScopedProjSilentLogger();
+
+    /**
+     * Returns the PROJ logger back to the default QGIS PROJ logger.
+     */
+    ~QgsScopedProjSilentLogger();
+
+};
+
+
 
 /**
  * \ingroup core

tests/src/core/testqgscoordinatetransform.cpp

@@ -611,21 +611,11 @@ void TestQgsCoordinateTransform::transformBoundingBox()
   QgsCoordinateTransform tr( sourceSrs, destSrs, QgsProject::instance() );
   const QgsRectangle crossingRect( 6374985, -3626584, 7021195, -3272435 );
   QgsRectangle resultRect = tr.transformBoundingBox( crossingRect, Qgis::TransformDirection::Forward, true );
-  QgsRectangle expectedRect;
-  expectedRect.setXMinimum( 175.771 );
-  expectedRect.setYMinimum( -39.7222 );
-  expectedRect.setXMaximum( -176.549 );
-  expectedRect.setYMaximum( -36.3951 );
-
-  qDebug( "BBox transform x min: %.17f", resultRect.xMinimum() );
-  qDebug( "BBox transform x max: %.17f", resultRect.xMaximum() );
-  qDebug( "BBox transform y min: %.17f", resultRect.yMinimum() );
-  qDebug( "BBox transform y max: %.17f", resultRect.yMaximum() );
-
-  QGSCOMPARENEAR( resultRect.xMinimum(), expectedRect.xMinimum(), 0.001 );
-  QGSCOMPARENEAR( resultRect.yMinimum(), expectedRect.yMinimum(), 0.001 );
-  QGSCOMPARENEAR( resultRect.xMaximum(), expectedRect.xMaximum(), 0.001 );
-  QGSCOMPARENEAR( resultRect.yMaximum(), expectedRect.yMaximum(), 0.001 );
+
+  QGSCOMPARENEAR( resultRect.xMinimum(), 175.771, 0.001 );
+  QGSCOMPARENEAR( resultRect.yMinimum(), -39.7222, 0.001 );
+  QGSCOMPARENEAR( resultRect.xMaximum(), -176.549, 0.001 );
+  QGSCOMPARENEAR( resultRect.yMaximum(), -36.3951, 0.001 );
 
   // test transforming a bounding box, resulting in an invalid transform - exception must be thrown
   tr = QgsCoordinateTransform( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:28356" ) ), QgsProject::instance() );