diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index 035fbab9ce..e13d40711a 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -729,10 +729,12 @@ struct PrecomputedOpCharacteristics { struct SortFunction { const std::map ↦ + const std::string BALLPARK_GEOGRAPHIC_OFFSET_FROM; explicit SortFunction(const std::map &mapIn) - : map(mapIn) {} + : map(mapIn), BALLPARK_GEOGRAPHIC_OFFSET_FROM( + std::string(BALLPARK_GEOGRAPHIC_OFFSET) + " from ") {} // Sorting function // Return true if a < b @@ -860,6 +862,48 @@ struct SortFunction { const auto &a_name = a->nameStr(); const auto &b_name = b->nameStr(); + + // Make sure that + // "Ballpark geographic offset from NAD83(CSRS)v6 to NAD83(CSRS)" + // has more priority than + // "Ballpark geographic offset from ITRF2008 to NAD83(CSRS)" + const auto posA = a_name.find(BALLPARK_GEOGRAPHIC_OFFSET_FROM); + const auto posB = b_name.find(BALLPARK_GEOGRAPHIC_OFFSET_FROM); + if (posA != std::string::npos && posB != std::string::npos) { + const auto pos2A = a_name.find(" to ", posA); + const auto pos2B = b_name.find(" to ", posB); + if (pos2A != std::string::npos && pos2B != std::string::npos) { + const auto pos3A = a_name.find(" + ", pos2A); + const auto pos3B = b_name.find(" + ", pos2B); + const std::string fromA = a_name.substr( + posA + BALLPARK_GEOGRAPHIC_OFFSET_FROM.size(), + pos2A - (posA + BALLPARK_GEOGRAPHIC_OFFSET_FROM.size())); + const std::string toA = + a_name.substr(pos2A + strlen(" to "), + pos3A == std::string::npos + ? pos3A + : pos3A - (pos2A + strlen(" to "))); + const std::string fromB = b_name.substr( + posB + BALLPARK_GEOGRAPHIC_OFFSET_FROM.size(), + pos2B - (posB + BALLPARK_GEOGRAPHIC_OFFSET_FROM.size())); + const std::string toB = + b_name.substr(pos2B + strlen(" to "), + pos3B == std::string::npos + ? pos3B + : pos3B - (pos2B + strlen(" to "))); + const bool similarCRSInA = + (fromA.find(toA) == 0 || toA.find(fromA) == 0); + const bool similarCRSInB = + (fromB.find(toB) == 0 || toB.find(fromB) == 0); + if (similarCRSInA && !similarCRSInB) { + return true; + } + if (!similarCRSInA && similarCRSInB) { + return false; + } + } + } + // The shorter name, the better ? if (a_name.size() < b_name.size()) { return true; @@ -1240,7 +1284,7 @@ struct FilterResults { area, getAccuracy(op), isPROJExportable, hasGrids, gridsAvailable, gridsKnown, stepCount, op->hasBallparkTransformation(), - op->nameStr().find("ballpark vertical transformation") != + op->nameStr().find(BALLPARK_VERTICAL_TRANSFORMATION) != std::string::npos, isNullTransformation(op->nameStr())); } @@ -4482,7 +4526,9 @@ void CoordinateOperationFactory::Private::createOperationsVertToVert( const double factor = convSrc / convDst; if (!equivalentVDatum) { auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()); + name += " ("; name += BALLPARK_VERTICAL_TRANSFORMATION; + name += ')'; auto conv = Transformation::createChangeVerticalUnit( util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name), sourceCRS, targetCRS, @@ -4588,9 +4634,12 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeogBallpark( targetCRSExtent.get(), util::IComparable::Criterion::EQUIVALENT); util::PropertyMap map; - map.set(common::IdentifiedObject::NAME_KEY, - buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) + - BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT) + std::string transfName( + buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr())); + transfName += " ("; + transfName += BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT; + transfName += ')'; + map.set(common::IdentifiedObject::NAME_KEY, transfName) .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY, sameExtent ? NN_NO_CHECK(sourceCRSExtent) : metadata::Extent::WORLD); @@ -5280,45 +5329,84 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( createOperations(sourceCRS, intermGeogSrc, context); const auto opsGeogToTarget = createOperations(intermGeogDst, targetCRS, context); + + // We preferably accept using an intermediate if the operations + // to it do not include ballpark operations, but if they do include + // grids, using such operations result will be better than nothing. + // This will help for example for "NAD83(CSRS) + CGVD28 height" to + // "NAD83(CSRS) + CGVD2013(CGG2013) height" where the transformation + // between "CGVD2013(CGG2013) height" and "NAD83(CSRS)" actually + // involves going through NAD83(CSRS)v6 + const auto hasKnownGrid = + [&dbContext](const CoordinateOperationNNPtr &op) { + const auto grids = op->gridsNeeded(dbContext, true); + if (grids.empty()) { + return false; + } + for (const auto &grid : grids) { + if (grid.available) { + return true; + } + } + return false; + }; + const bool hasNonTrivialSrcTransf = !opsSrcToGeog.empty() && - !opsSrcToGeog.front()->hasBallparkTransformation(); + (!opsSrcToGeog.front()->hasBallparkTransformation() || + hasKnownGrid(opsSrcToGeog.front())); const bool hasNonTrivialTargetTransf = !opsGeogToTarget.empty() && - !opsGeogToTarget.front()->hasBallparkTransformation(); + (!opsGeogToTarget.front()->hasBallparkTransformation() || + hasKnownGrid(opsGeogToTarget.front())); if (hasNonTrivialSrcTransf && hasNonTrivialTargetTransf) { const auto opsGeogSrcToGeogDst = createOperations(intermGeogSrc, intermGeogDst, context); - for (const auto &op1 : opsSrcToGeog) { - if (op1->hasBallparkTransformation()) { - // std::cerr << "excluded " << op1->nameStr() << std::endl; - continue; - } - for (const auto &op2 : opsGeogSrcToGeogDst) { - for (const auto &op3 : opsGeogToTarget) { - if (op3->hasBallparkTransformation()) { - // std::cerr << "excluded " << op3->nameStr() << - // std::endl; - continue; - } - try { - res.emplace_back( - ConcatenatedOperation::createComputeMetadata( - intermGeogSrcIsSameAsIntermGeogDst - ? std::vector< - CoordinateOperationNNPtr>{op1, - op3} - : std::vector< - CoordinateOperationNNPtr>{op1, - op2, - op3}, - disallowEmptyIntersection)); - } catch (const std::exception &) { + // In first pass, exclude (horizontal) ballpark operations, but + // accept them in second pass. + for (int pass = 0; pass <= 1 && res.empty(); pass++) { + for (const auto &op1 : opsSrcToGeog) { + if (pass == 0 && op1->hasBallparkTransformation()) { + // std::cerr << "excluded " << op1->nameStr() << + // std::endl; + continue; + } + if (op1->nameStr().find(BALLPARK_VERTICAL_TRANSFORMATION) != + std::string::npos) { + continue; + } + for (const auto &op2 : opsGeogSrcToGeogDst) { + for (const auto &op3 : opsGeogToTarget) { + if (pass == 0 && op3->hasBallparkTransformation()) { + // std::cerr << "excluded " << op3->nameStr() << + // std::endl; + continue; + } + if (op3->nameStr().find( + BALLPARK_VERTICAL_TRANSFORMATION) != + std::string::npos) { + continue; + } + try { + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + intermGeogSrcIsSameAsIntermGeogDst + ? std::vector< + CoordinateOperationNNPtr>{op1, + op3} + : std::vector< + CoordinateOperationNNPtr>{op1, + op2, + op3}, + disallowEmptyIntersection)); + } catch (const std::exception &) { + } } } } } } + if (!res.empty()) { return; } diff --git a/src/iso19111/operation/oputils.cpp b/src/iso19111/operation/oputils.cpp index b5834edf2c..5efbf07fc5 100644 --- a/src/iso19111/operation/oputils.cpp +++ b/src/iso19111/operation/oputils.cpp @@ -61,10 +61,10 @@ const char *NULL_GEOGRAPHIC_OFFSET = "Null geographic offset"; const char *NULL_GEOCENTRIC_TRANSLATION = "Null geocentric translation"; const char *BALLPARK_GEOGRAPHIC_OFFSET = "Ballpark geographic offset"; const char *BALLPARK_VERTICAL_TRANSFORMATION = - " (ballpark vertical transformation)"; + "ballpark vertical transformation"; const char *BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT = - " (ballpark vertical transformation, without ellipsoid height to vertical " - "height correction)"; + "ballpark vertical transformation, without ellipsoid height to vertical " + "height correction"; // --------------------------------------------------------------------------- diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index a445d8c754..e387b7ab35 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -4518,6 +4518,63 @@ TEST(operation, compoundCRS_to_compoundCRS_issue_2720) { // --------------------------------------------------------------------------- +TEST( + operation, + compoundCRS_to_compoundCRS_concatenated_operation_with_two_vert_transformation_and_ballpark_geog) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + // "NAD83(CSRS) + CGVD28 height" + auto srcObj = createFromUserInput("EPSG:4617+5713", + authFactory->databaseContext(), false); + auto src = nn_dynamic_pointer_cast(srcObj); + ASSERT_TRUE(src != nullptr); + + // "NAD83(CSRS) + CGVD2013(CGG2013) height" + auto dstObj = createFromUserInput("EPSG:4617+6647", + authFactory->databaseContext(), false); + auto dst = nn_dynamic_pointer_cast(dstObj); + ASSERT_TRUE(dst != nullptr); + + // That transformation involves doing CGVD28 height to CGVD2013(CGG2013) + // height by doing: + // - CGVD28 height to NAD83(CSRS): EPSG registered operation + // - NAD83(CSRS) to CGVD2013(CGG2013) height by doing: + // * NAD83(CSRS) to NAD83(CSRS)v6: ballpark + // * NAD83(CSRS)v6 to CGVD2013(CGG2013): EPSG registered operation + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + { + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(src), NN_NO_CHECK(dst), ctxt); + ASSERT_GE(list.size(), 1U); + // Check that we have the transformation using NAD83(CSRS)v6 first + // (as well as the one between NAD83(CSRS) to CGVD28 height) + EXPECT_EQ(list[0]->nameStr(), + "Inverse of NAD83(CSRS) to CGVD28 height (1) + " + "Inverse of Ballpark geographic offset from NAD83(CSRS)v6 to " + "NAD83(CSRS) + " + "NAD83(CSRS)v6 to CGVD2013(CGG2013) height (1) + " + "Inverse of Ballpark geographic offset from NAD83(CSRS) to " + "NAD83(CSRS)v6"); + } + { + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(dst), NN_NO_CHECK(src), ctxt); + ASSERT_GE(list.size(), 1U); + // Check that we have the transformation using NAD83(CSRS)v6 first + // (as well as the one between NAD83(CSRS) to CGVD28 height) + EXPECT_EQ( + list[0]->nameStr(), + "Ballpark geographic offset from NAD83(CSRS) to NAD83(CSRS)v6 + " + "Inverse of NAD83(CSRS)v6 to CGVD2013(CGG2013) height (1) + " + "Ballpark geographic offset from NAD83(CSRS)v6 to NAD83(CSRS) + " + "NAD83(CSRS) to CGVD28 height (1)"); + } +} + +// --------------------------------------------------------------------------- + TEST(operation, vertCRS_to_vertCRS) { auto vertcrs_m_obj = PROJStringParser().createFromPROJString("+vunits=m");