Merge pull request #3248 from rouault/WGS84_to_GDA2020

createOperations(): prefer simpler pipelines / affects WGS 84 to GDA94/GDA2020

include/proj/io.hpp

@@ -479,6 +479,8 @@ class PROJ_GCC_DLL PROJStringFormatter {
 
     PROJ_INTERNAL Convention convention() const;
 
+    PROJ_INTERNAL size_t getStepCount() const;
+
     //! @endcond
 
   protected:

src/iso19111/io.cpp

@@ -8662,6 +8662,13 @@ PROJStringFormatter::Convention PROJStringFormatter::convention() const {
 
 // ---------------------------------------------------------------------------
 
+// Return the number of steps in the pipeline.
+// Note: this value will change after calling toString() that will run
+// optimizations.
+size_t PROJStringFormatter::getStepCount() const { return d->steps_.size(); }
+
+// ---------------------------------------------------------------------------
+
 bool PROJStringFormatter::getUseApproxTMerc() const {
     return d->useApproxTMerc_;
 }

src/iso19111/operation/coordinateoperationfactory.cpp

@@ -704,6 +704,7 @@ struct PrecomputedOpCharacteristics {
     bool gridsAvailable_ = false;
     bool gridsKnown_ = false;
     size_t stepCount_ = 0;
+    size_t projStepCount_ = 0;
     bool isApprox_ = false;
     bool hasBallparkVertical_ = false;
     bool isNullTransformation_ = false;
@@ -712,12 +713,13 @@ struct PrecomputedOpCharacteristics {
     PrecomputedOpCharacteristics(double area, double accuracy,
                                  bool isPROJExportable, bool hasGrids,
                                  bool gridsAvailable, bool gridsKnown,
-                                 size_t stepCount, bool isApprox,
-                                 bool hasBallparkVertical,
+                                 size_t stepCount, size_t projStepCount,
+                                 bool isApprox, bool hasBallparkVertical,
                                  bool isNullTransformation)
         : area_(area), accuracy_(accuracy), isPROJExportable_(isPROJExportable),
           hasGrids_(hasGrids), gridsAvailable_(gridsAvailable),
-          gridsKnown_(gridsKnown), stepCount_(stepCount), isApprox_(isApprox),
+          gridsKnown_(gridsKnown), stepCount_(stepCount),
+          projStepCount_(projStepCount), isApprox_(isApprox),
           hasBallparkVertical_(hasBallparkVertical),
           isNullTransformation_(isNullTransformation) {}
 };
@@ -860,6 +862,18 @@ struct SortFunction {
             return false;
         }
 
+        // Compare number of steps in PROJ pipeline, and prefer the ones
+        // with less operations.
+        if (iterA->second.projStepCount_ != 0 &&
+            iterB->second.projStepCount_ != 0) {
+            if (iterA->second.projStepCount_ < iterB->second.projStepCount_) {
+                return true;
+            }
+            if (iterB->second.projStepCount_ < iterA->second.projStepCount_) {
+                return false;
+            }
+        }
+
         const auto &a_name = a->nameStr();
         const auto &b_name = b->nameStr();
 
@@ -1259,11 +1273,20 @@ struct FilterResults {
 
             bool isPROJExportable = false;
             auto formatter = io::PROJStringFormatter::create();
+            size_t projStepCount = 0;
             try {
-                op->exportToPROJString(formatter.get());
+                const auto str = op->exportToPROJString(formatter.get());
                 // Grids might be missing, but at least this is something
                 // PROJ could potentially process
                 isPROJExportable = true;
+
+                // We exclude pipelines with +proj=xyzgridshift as they
+                // generate more steps, but are more precise.
+                if (str.find("+proj=xyzgridshift") == std::string::npos) {
+                    auto formatter2 = io::PROJStringFormatter::create();
+                    formatter2->ingestPROJString(str);
+                    projStepCount = formatter2->getStepCount();
+                }
             } catch (const std::exception &) {
             }
 
@@ -1282,7 +1305,7 @@ struct FilterResults {
 #endif
             map[op.get()] = PrecomputedOpCharacteristics(
                 area, getAccuracy(op), isPROJExportable, hasGrids,
-                gridsAvailable, gridsKnown, stepCount,
+                gridsAvailable, gridsKnown, stepCount, projStepCount,
                 op->hasBallparkTransformation(),
                 op->nameStr().find(BALLPARK_VERTICAL_TRANSFORMATION) !=
                     std::string::npos,

test/unit/test_operationfactory.cpp

@@ -1732,6 +1732,40 @@ TEST(operation, esri_projectedCRS_to_geogCRS_with_ITRF_intermediate_context) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, geogCRS_to_geogCRS_WGS84_to_GDA2020) {
+    // 2D reduction of use case of /~https://github.com/OSGeo/PROJ/issues/2348
+    auto dbContext = DatabaseContext::create();
+    auto authFactory = AuthorityFactory::create(dbContext, "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    {
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            // GDA2020
+            authFactory->createCoordinateReferenceSystem("7844"),
+            // WGS 84
+            authFactory->createCoordinateReferenceSystem("4326"), ctxt);
+        ASSERT_GE(list.size(), 1U);
+        EXPECT_EQ(list[0]->nameStr(), "GDA2020 to WGS 84 (2)");
+    }
+
+    // Inverse
+    {
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            // WGS 84
+            authFactory->createCoordinateReferenceSystem("4326"),
+            // GDA2020
+            authFactory->createCoordinateReferenceSystem("7844"), ctxt);
+        ASSERT_GE(list.size(), 1U);
+        EXPECT_EQ(list[0]->nameStr(), "Inverse of GDA2020 to WGS 84 (2)");
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 static ProjectedCRSNNPtr createUTM31_WGS84() {
     return ProjectedCRS::create(
         PropertyMap(), GeographicCRS::EPSG_4326,
@@ -5925,6 +5959,42 @@ TEST(operation,
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, compoundCRS_to_geogCRS_3D_WGS84_to_GDA2020_AHD_Height) {
+    // Use case of /~https://github.com/OSGeo/PROJ/issues/2348
+    auto dbContext = DatabaseContext::create();
+    auto authFactory = AuthorityFactory::create(dbContext, "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    {
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            // GDA2020 + AHD height
+            authFactory->createCoordinateReferenceSystem("9463"),
+            // WGS 84 3D
+            authFactory->createCoordinateReferenceSystem("4979"), ctxt);
+        ASSERT_GE(list.size(), 1U);
+        EXPECT_EQ(list[0]->nameStr(), "Inverse of GDA2020 to AHD height (1) + "
+                                      "GDA2020 to WGS 84 (2)");
+    }
+
+    // Inverse
+    {
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            // WGS 84 3D
+            authFactory->createCoordinateReferenceSystem("4979"),
+            // GDA2020 + AHD height
+            authFactory->createCoordinateReferenceSystem("9463"), ctxt);
+        ASSERT_GE(list.size(), 1U);
+        EXPECT_EQ(list[0]->nameStr(), "Inverse of GDA2020 to WGS 84 (2) + "
+                                      "GDA2020 to AHD height (1)");
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, compoundCRS_to_geogCRS_2D_promote_to_3D_context) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");