Merge pull request #3361 from rouault/fix_3360

createOperations(): emulate PROJ < 6 behavior when doing geocentric <--> geographic transformation between datum with unknown transformation (fixes #3360)
OSGeo · Oct 6, 2022 · 632eb41 · 632eb41
2 parents a75ae93 + aaf0d2c
commit 632eb41
Show file tree

Hide file tree

Showing 3 changed files with 61 additions and 26 deletions.
diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -4037,6 +4037,10 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod(
             res.emplace_back(
                 Conversion::createGeographicGeocentric(sourceCRS, targetCRS));
         } else if (isSrcGeocentric && geogDst) {
+#if 0
+            // The below logic was used between PROJ >= 6.0 and < 9.2
+            // It assumed that the geocentric origin of the 2 datums
+            // matched.
             std::string interm_crs_name(geogDst->nameStr());
             interm_crs_name += " (geocentric)";
             auto interm_crs =
@@ -4053,15 +4057,42 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod(
                 createBallparkGeocentricTranslation(sourceCRS, interm_crs);
             auto opSecond =
                 Conversion::createGeographicGeocentric(interm_crs, targetCRS);
-            res.emplace_back(ConcatenatedOperation::createComputeMetadata(
-                {opFirst, opSecond}, disallowEmptyIntersection));
+#else
+            // The below logic is used since PROJ >= 9.2. It emulates the
+            // behavior of PROJ < 6 by converting from the source geocentric CRS
+            // to its corresponding geographic CRS, and then doing a null
+            // geographic offset between that CRS and the target geographic CRS
+            std::string interm_crs_name(geodSrc->nameStr());
+            interm_crs_name += " (geographic)";
+            auto interm_crs = util::nn_static_pointer_cast<crs::CRS>(
+                crs::GeographicCRS::create(
+                    addDomains(util::PropertyMap().set(
+                                   common::IdentifiedObject::NAME_KEY,
+                                   interm_crs_name),
+                               geodSrc),
+                    geodSrc->datum(), geodSrc->datumEnsemble(),
+                    cs::EllipsoidalCS::createLongitudeLatitudeEllipsoidalHeight(
+                        common::UnitOfMeasure::DEGREE,
+                        common::UnitOfMeasure::METRE)));
+            auto opFirst =
+                Conversion::createGeographicGeocentric(sourceCRS, interm_crs);
+            auto opsSecond = createOperations(interm_crs, targetCRS, context);
+            for (const auto &opSecond : opsSecond) {
+                try {
+                    res.emplace_back(
+                        ConcatenatedOperation::createComputeMetadata(
+                            {opFirst, opSecond}, disallowEmptyIntersection));
+                } catch (const InvalidOperationEmptyIntersection &) {
+                }
+            }
+#endif
         } else {
             // Apply previous case in reverse way
             std::vector<CoordinateOperationNNPtr> resTmp;
             createOperationsGeodToGeod(targetCRS, sourceCRS, context, geodDst,
                                        geodSrc, resTmp);
-            assert(resTmp.size() == 1);
-            res.emplace_back(resTmp.front()->inverse());
+            resTmp = applyInverse(resTmp);
+            res.insert(res.end(), resTmp.begin(), resTmp.end());
         }
 
         return;

diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
@@ -11105,10 +11105,11 @@ TEST(io, projparse_cart_unit) {
         GeographicCRS::EPSG_4326, NN_NO_CHECK(crs));
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-              "+ellps=WGS84 +step +proj=unitconvert +xy_in=m +z_in=m "
-              "+xy_out=km +z_out=km");
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m "
+              "+step +proj=cart +ellps=WGS84 "
+              "+step +proj=unitconvert +xy_in=m +z_in=m +xy_out=km +z_out=km");
 }
 
 // ---------------------------------------------------------------------------
@@ -11123,11 +11124,13 @@ TEST(io, projparse_cart_unit_numeric) {
     auto op = CoordinateOperationFactory::create()->createOperation(
         GeographicCRS::EPSG_4326, NN_NO_CHECK(crs));
     ASSERT_TRUE(op != nullptr);
-    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-              "+ellps=WGS84 +step +proj=unitconvert +xy_in=m +z_in=m "
-              "+xy_out=500 +z_out=500");
+    EXPECT_EQ(
+        op->exportToPROJString(PROJStringFormatter::create().get()),
+        "+proj=pipeline "
+        "+step +proj=axisswap +order=2,1 "
+        "+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m "
+        "+step +proj=cart +ellps=WGS84 "
+        "+step +proj=unitconvert +xy_in=m +z_in=m +xy_out=500 +z_out=500");
 }
 
 // ---------------------------------------------------------------------------

diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
@@ -550,7 +550,7 @@ TEST(operation, geogCRS_to_geogCRS_context_helmert_geog3D_to_geocentirc) {
         authFactory->createCoordinateReferenceSystem("4939"),
         // GDA2020 geocentric
         authFactory->createCoordinateReferenceSystem("7842"), ctxt);
-    ASSERT_EQ(list.size(), 1U);
+    ASSERT_GE(list.size(), 1U);
 
     // Check there is no push / pop of v_3
     EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
@@ -1331,13 +1331,13 @@ TEST(operation, geocentricCRS_to_geogCRS_different_datum) {
         createGeocentricDatumWGS84(), GeographicCRS::EPSG_4269);
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->nameStr(),
-              "Ballpark geocentric translation from WGS 84 to NAD83 "
-              "(geocentric) + Conversion from NAD83 "
-              "(geocentric) to NAD83");
+              "Conversion from WGS 84 to WGS 84 (geographic) + "
+              "Ballpark geographic offset from WGS 84 (geographic) to NAD83");
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +inv +proj=cart +ellps=GRS80 +step "
-              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
-              "+order=2,1");
+              "+proj=pipeline "
+              "+step +inv +proj=cart +ellps=WGS84 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
 }
 
 // ---------------------------------------------------------------------------
@@ -1347,13 +1347,14 @@ TEST(operation, geogCRS_to_geocentricCRS_different_datum) {
     auto op = CoordinateOperationFactory::create()->createOperation(
         GeographicCRS::EPSG_4269, createGeocentricDatumWGS84());
     ASSERT_TRUE(op != nullptr);
-    EXPECT_EQ(op->nameStr(), "Conversion from NAD83 to NAD83 (geocentric) + "
-                             "Ballpark geocentric translation from NAD83 "
-                             "(geocentric) to WGS 84");
+    EXPECT_EQ(op->nameStr(),
+              "Ballpark geographic offset from NAD83 to WGS 84 (geographic) + "
+              "Conversion from WGS 84 (geographic) to WGS 84");
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-              "+ellps=GRS80");
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m "
+              "+step +proj=cart +ellps=WGS84");
 }
 
 // ---------------------------------------------------------------------------