Always fallback to default coordinate operations when bounding box transfrom fails

Bounding box transforms are inherently approximate, so we can safely
just fallback to the proj default operation if the user-specified
operation fails when transforming a bounding box.

This fixes cases where the user-specified operation involves a
grid shift and the bounds to transform falls outside of the grid.
In this case proj_trans_bounds fails. By falling back to a proj
default operation and re-trying then we get a valid approximate
transformation of the bounding box.

Fixes map rendering issues when users are using a grid shift operation
to transform layers to map CRS

Fixes #60737
Fixes #60753

Author: nyalldawson

src/core/proj/qgscoordinatetransform.cpp

@@ -676,9 +676,27 @@ QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle &r
   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 );
+  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 ) )
+       && ( d->mAvailableOpCount > 1 || d->mAvailableOpCount == -1 ) // only use fallbacks if more than one operation is possible -- otherwise we've already tried it and it failed
+     )
+  {
+    // fail #1 -- try with getting proj to auto-pick an appropriate coordinate operation for the points
+    if ( PJ *transform = d->threadLocalFallbackProjData() )
+    {
+      projResult = proj_trans_bounds( projContext, transform, ( 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 )

tests/src/python/test_qgscoordinatetransform.py

@@ -127,6 +127,51 @@ def testTransformQgsRectangle_Regression17600(self):
             myTransformedExtentReverse.yMinimum(), myExtent.yMinimum()
         )
 
+    def test_transform_bounding_box_grid(self):
+        """
+        This test assumes the ca_nrc_NA83SCRS.tif grid is available on the system!
+        """
+        transform = QgsCoordinateTransform(
+            QgsCoordinateReferenceSystem("EPSG:4269"),
+            QgsCoordinateReferenceSystem("EPSG:3857"),
+            QgsCoordinateTransformContext(),
+        )
+        res = transform.transformBoundingBox(
+            QgsRectangle(
+                -123.65020876249999,
+                45.987175336410544,
+                -101.22289073749998,
+                62.961980263589439,
+            )
+        )
+        self.assertAlmostEqual(res.xMinimum(), -13764678, -2)
+        self.assertAlmostEqual(res.yMinimum(), 5778294, -2)
+        self.assertAlmostEqual(res.xMaximum(), -11268080, -2)
+        self.assertAlmostEqual(res.yMaximum(), 9090934, -2)
+
+        transform = QgsCoordinateTransform(
+            QgsCoordinateReferenceSystem("EPSG:4269"),
+            QgsCoordinateReferenceSystem("EPSG:3857"),
+            QgsCoordinateTransformContext(),
+        )
+        # force use of grid shift operation. This will fail as the source rect is outside of the grid bounds, but we should silently
+        # fall back to the non-grid operation
+        transform.setCoordinateOperation(
+            "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=ca_nrc_NA83SCRS.tif +step +proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84"
+        )
+        res = transform.transformBoundingBox(
+            QgsRectangle(
+                -123.65020876249999,
+                45.987175336410544,
+                -101.22289073749998,
+                62.961980263589439,
+            )
+        )
+        self.assertAlmostEqual(res.xMinimum(), -13764678, -2)
+        self.assertAlmostEqual(res.yMinimum(), 5778294, -2)
+        self.assertAlmostEqual(res.xMaximum(), -11268080, -2)
+        self.assertAlmostEqual(res.yMaximum(), 9090934, -2)
+
     def testContextProj6(self):
         """
         Various tests to ensure that datum transforms are correctly set respecting context