From aaf0d2c2be3c4c04e072c704bcb6eb0802853362 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 5 Oct 2022 19:17:21 +0200 Subject: [PATCH] createOperations(): emulate PROJ < 6 behavior when doing geocentric <--> geographic transformation between datum with unknown transformation (fixes #3360) --- .../operation/coordinateoperationfactory.cpp | 39 +++++++++++++++++-- test/unit/test_io.cpp | 21 +++++----- test/unit/test_operationfactory.cpp | 27 ++++++------- 3 files changed, 61 insertions(+), 26 deletions(-) diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index 7eaac8dcc3..c27a9dcff3 100644 --- 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::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 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 index ac88a81670..1fd56e6fad 100644 --- 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 index cdfe851663..cc423e2041 100644 --- 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"); } // ---------------------------------------------------------------------------