From 890c94a730474f057f5237ca07699d6af600ed3f Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 11 Sep 2019 22:47:32 +0200 Subject: createOperations(): make sure sorting function is transitive (a < b and b < c --> a < c), to get consistent results --- src/iso19111/coordinateoperation.cpp | 50 +++++++++++++++++++++++++++--------- 1 file changed, 38 insertions(+), 12 deletions(-) (limited to 'src/iso19111/coordinateoperation.cpp') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 768ef76f..a1612339 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -56,7 +56,9 @@ #include #include -#ifdef DEBUG +// #define DEBUG +// #define DEBUG_SORT +#if defined(DEBUG) || defined(DEBUG_SORT) #include #endif @@ -10309,16 +10311,12 @@ struct SortFunction { return false; } - if (iterA->second.hasGrids_ && iterB->second.hasGrids_) { - // Operations where grids are all available go before other - if (iterA->second.gridsAvailable_ && - !iterB->second.gridsAvailable_) { - return true; - } - if (iterB->second.gridsAvailable_ && - !iterA->second.gridsAvailable_) { - return false; - } + // Operations where grids are all available go before other + if (iterA->second.gridsAvailable_ && !iterB->second.gridsAvailable_) { + return true; + } + if (iterB->second.gridsAvailable_ && !iterA->second.gridsAvailable_) { + return false; } // Operations where grids are all known in our DB go before other @@ -10680,7 +10678,35 @@ struct FilterResults { } // Sort ! - std::sort(res.begin(), res.end(), SortFunction(map)); + SortFunction sortFunc(map); + std::sort(res.begin(), res.end(), sortFunc); + +// Debug code to check consistency of the sort function +#ifdef DEBUG_SORT + constexpr bool debugSort = true; +#elif !defined(NDEBUG) + const bool debugSort = getenv("PROJ_DEBUG_SORT_FUNCT") != nullptr; +#endif +#if defined(DEBUG_SORT) || !defined(NDEBUG) + if (debugSort) { + const bool assertIfIssue = + !(getenv("PROJ_DEBUG_SORT_FUNCT_ASSERT") != nullptr); + for (size_t i = 0; i < res.size(); ++i) { + for (size_t j = i + 1; j < res.size(); ++j) { + if (sortFunc(res[j], res[i])) { +#ifdef DEBUG_SORT + std::cerr << "Sorting issue with entry " << i << "(" + << res[i]->nameStr() << ") and " << j << "(" + << res[j]->nameStr() << ")" << std::endl; +#endif + if (assertIfIssue) { + assert(false); + } + } + } + } + } +#endif } // ---------------------------------------------------------------------- -- cgit v1.2.3 From 63857c92b271bbcd10df0a032304982011acb2a9 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 10 Sep 2019 20:21:36 +0200 Subject: Coordinate transformation: improve transformations from/to WGS84 (Gxxxx) Currently very few transformations from/to WGS84 (Gxxxx) are registered in the EPSG database, and there isn't even transformations between WGS84 EPSG:4326 and those ones. Consequently transformations to those realizations often ended up as no-operation, whereas going through WGS84 EPSG:4326 will bring more meaningful results. So register those EPSG:4326<-->WGS 84 (Gxxx) null transformations, and when having WGS 84 (Gxxx) as source/target, consider EPSG:4326 as an intermediate. This change has no effect on the existing direct transformations from/to WGS 84 (Gxxx). --- src/iso19111/coordinateoperation.cpp | 135 ++++++++++++++++++++++++++++++++++- 1 file changed, 134 insertions(+), 1 deletion(-) (limited to 'src/iso19111/coordinateoperation.cpp') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index a1612339..f7ea385d 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -10171,6 +10171,7 @@ struct CoordinateOperationFactory::Private { const crs::CRSNNPtr &targetCRS; const CoordinateOperationContextNNPtr &context; bool inCreateOperationsWithDatumPivotAntiRecursion = false; + bool inCreateOperationsThroughPreferredHub = false; Context(const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn, @@ -10195,6 +10196,12 @@ struct CoordinateOperationFactory::Private { const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst, Context &context); + static void createOperationsThroughPreferredHub( + std::vector &res, + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst, + Context &context); + static bool hasPerfectAccuracyResult(const std::vector &res, const Context &context); @@ -12016,6 +12023,122 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot( // --------------------------------------------------------------------------- +void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub( + std::vector &res, const crs::CRSNNPtr &sourceCRS, + const crs::CRSNNPtr &targetCRS, const crs::GeodeticCRS *geodSrc, + const crs::GeodeticCRS *geodDst, Private::Context &context) { + + const auto &srcDatum = geodSrc->datum(); + const auto &dstDatum = geodDst->datum(); + + if (!srcDatum || !dstDatum) + return; + const auto &srcDatumIds = srcDatum->identifiers(); + const auto &dstDatumIds = dstDatum->identifiers(); + if (srcDatumIds.empty() || dstDatumIds.empty()) + return; + + const auto &authFactory = context.context->getAuthorityFactory(); + + const auto srcAuthFactory = io::AuthorityFactory::create( + authFactory->databaseContext(), *(srcDatumIds.front()->codeSpace())); + const auto srcPreferredHubs = + srcAuthFactory->getPreferredHubGeodeticReferenceFrames( + srcDatumIds.front()->code()); + + const auto dstAuthFactory = io::AuthorityFactory::create( + authFactory->databaseContext(), *(dstDatumIds.front()->codeSpace())); + const auto dstPreferredHubs = + dstAuthFactory->getPreferredHubGeodeticReferenceFrames( + dstDatumIds.front()->code()); + if (srcPreferredHubs.empty() && dstPreferredHubs.empty()) + return; + + // Currently if we have prefered hubs for both source and target, we + // will use only the one for target, arbitrarily... We could use boths + // but that would complicate things. + if (!srcPreferredHubs.empty() && dstPreferredHubs.empty()) { + std::vector resTmp; + createOperationsThroughPreferredHub(resTmp, targetCRS, sourceCRS, + geodDst, geodSrc, context); + if (!resTmp.empty()) { + resTmp = applyInverse(resTmp); + res.insert(res.end(), resTmp.begin(), resTmp.end()); + } + return; + } + assert(!dstPreferredHubs.empty()); + +#ifdef DEBUG + EnterDebugLevel enterFunction; + debugTrace("createOperationsThroughPreferredHub(" + + objectAsStr(sourceCRS.get()) + "," + + objectAsStr(targetCRS.get()) + ")"); +#endif + + struct AntiRecursionGuard { + Context &context; + + explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) { + assert(!context.inCreateOperationsThroughPreferredHub); + context.inCreateOperationsThroughPreferredHub = true; + } + + ~AntiRecursionGuard() { + context.inCreateOperationsThroughPreferredHub = false; + } + }; + AntiRecursionGuard guard(context); + + std::vector candidatesIntermCRS; + for (const auto &datumHub : dstPreferredHubs) { + auto candidates = + findCandidateGeodCRSForDatum(authFactory, datumHub.get()); + bool addedGeog2D = false; + for (const auto &intermCRS : candidates) { + auto geogCRS = dynamic_cast(intermCRS.get()); + if (geogCRS && + geogCRS->coordinateSystem()->axisList().size() == 2) { + candidatesIntermCRS.emplace_back(intermCRS); + addedGeog2D = true; + break; + } + } + if (!addedGeog2D) { + candidatesIntermCRS.insert(candidatesIntermCRS.end(), + candidates.begin(), candidates.end()); + } + } + + const bool allowEmptyIntersection = true; + for (const auto &intermCRS : candidatesIntermCRS) { +#ifdef DEBUG + EnterDebugLevel loopLevel; + debugTrace("try " + objectAsStr(sourceCRS.get()) + "->" + + objectAsStr(intermCRS.get()) + "->" + + objectAsStr(targetCRS.get()) + ")"); + EnterDebugLevel loopLevel2; +#endif + const auto opsFirst = createOperations(sourceCRS, intermCRS, context); + const auto opsLast = createOperations(intermCRS, targetCRS, context); + for (const auto &opFirst : opsFirst) { + for (const auto &opLast : opsLast) { + if (!opFirst->hasBallparkTransformation() || + !opLast->hasBallparkTransformation()) { + try { + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + {opFirst, opLast}, !allowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { + } + } + } + } + } +} + +// --------------------------------------------------------------------------- + static CoordinateOperationNNPtr createBallparkGeocentricTranslation(const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS) { @@ -12262,7 +12385,17 @@ CoordinateOperationFactory::Private::createOperations( sourceCRS, targetCRS, context.context); res.insert(res.end(), resWithIntermediate.begin(), resWithIntermediate.end()); - doFilterAndCheckPerfectOp = true; + doFilterAndCheckPerfectOp = !res.empty(); + + // If transforming from/to WGS84 (Gxxxx), try through 'neutral' + // WGS84 + if (res.empty() && geodSrc && geodDst && + !context.inCreateOperationsThroughPreferredHub && + !context.inCreateOperationsWithDatumPivotAntiRecursion) { + createOperationsThroughPreferredHub( + res, sourceCRS, targetCRS, geodSrc, geodDst, context); + doFilterAndCheckPerfectOp = !res.empty(); + } } } -- cgit v1.2.3 From eed28e5183579d09e102d1ad72e91fc82005dfe8 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 12 Sep 2019 15:55:04 +0200 Subject: createOperations(): use more candidates when transforming between a geographic and vertical CRS For example when transforming from NAD83+NAVD88 height to WGS84, there is no transformation between NAVD88 height to WGS84. In that case, use all potential transformations from NAVD88 height to another geographic CRS for the vertical part. --- src/iso19111/coordinateoperation.cpp | 301 +++++++++++++++++++++++++++-------- 1 file changed, 238 insertions(+), 63 deletions(-) (limited to 'src/iso19111/coordinateoperation.cpp') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index f7ea385d..aad86410 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -10172,6 +10172,7 @@ struct CoordinateOperationFactory::Private { const CoordinateOperationContextNNPtr &context; bool inCreateOperationsWithDatumPivotAntiRecursion = false; bool inCreateOperationsThroughPreferredHub = false; + bool inCreateOperationsGeogToVertWithIntermediate = false; Context(const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn, @@ -10202,6 +10203,11 @@ struct CoordinateOperationFactory::Private { const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst, Context &context); + static std::vector + createOperationsGeogToVertWithIntermediate(const crs::CRSNNPtr &sourceCRS, + const crs::CRSNNPtr &targetCRS, + Context &context); + static bool hasPerfectAccuracyResult(const std::vector &res, const Context &context); @@ -10950,18 +10956,21 @@ applyInverse(const std::vector &list) { //! @cond Doxygen_Suppress -static void buildSourceAndTargetCRSIds( - const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, - const CoordinateOperationContextNNPtr &context, - std::list> &sourceIds, - std::list> &targetIds) { +static void buildCRSIds(const crs::CRSNNPtr &crs, + const CoordinateOperationContextNNPtr &context, + std::list> &ids) { const auto &authFactory = context->getAuthorityFactory(); assert(authFactory); const auto &authFactoryName = authFactory->getAuthority(); - const auto findObjectInDB = [&authFactory, &authFactoryName]( - const crs::CRSNNPtr &crs, - std::list> &idList) { + for (const auto &id : crs->identifiers()) { + const auto &authName = *(id->codeSpace()); + const auto &code = id->code(); + if (!authName.empty()) { + ids.emplace_back(authName, code); + } + } + if (ids.empty()) { try { const auto tmpAuthFactory = io::AuthorityFactory::create( authFactory->databaseContext(), @@ -10987,35 +10996,37 @@ static void buildSourceAndTargetCRSIds( matches.front().get(), util::IComparable::Criterion::EQUIVALENT) && !matches.front()->identifiers().empty()) { - const auto &ids = matches.front()->identifiers(); - idList.emplace_back(*(ids[0]->codeSpace()), ids[0]->code()); + const auto &tmpIds = matches.front()->identifiers(); + ids.emplace_back(*(tmpIds[0]->codeSpace()), + tmpIds[0]->code()); } } } catch (const std::exception &) { } - }; - - for (const auto &id : sourceCRS->identifiers()) { - const auto &authName = *(id->codeSpace()); - const auto &code = id->code(); - if (!authName.empty()) { - sourceIds.emplace_back(authName, code); - } - } - if (sourceIds.empty()) { - findObjectInDB(sourceCRS, sourceIds); } +} - for (const auto &id : targetCRS->identifiers()) { - const auto &authName = *(id->codeSpace()); - const auto &code = id->code(); - if (!authName.empty()) { - targetIds.emplace_back(authName, code); - } +// --------------------------------------------------------------------------- + +static std::vector +getCandidateAuthorities(const io::AuthorityFactoryPtr &authFactory, + const std::string &srcAuthName, + const std::string &targetAuthName) { + const auto &authFactoryName = authFactory->getAuthority(); + std::vector authorities; + if (authFactoryName == "any") { + authorities.emplace_back(); } - if (targetIds.empty()) { - findObjectInDB(targetCRS, targetIds); + if (authFactoryName.empty()) { + authorities = authFactory->databaseContext()->getAllowedAuthorities( + srcAuthName, targetAuthName); + if (authorities.empty()) { + authorities.emplace_back(); + } + } else { + authorities.emplace_back(authFactoryName); } + return authorities; } // --------------------------------------------------------------------------- @@ -11027,12 +11038,11 @@ findOpsInRegistryDirect(const crs::CRSNNPtr &sourceCRS, const CoordinateOperationContextNNPtr &context) { const auto &authFactory = context->getAuthorityFactory(); assert(authFactory); - const auto &authFactoryName = authFactory->getAuthority(); std::list> sourceIds; std::list> targetIds; - buildSourceAndTargetCRSIds(sourceCRS, targetCRS, context, sourceIds, - targetIds); + buildCRSIds(sourceCRS, context, sourceIds); + buildCRSIds(targetCRS, context, targetIds); for (const auto &idSrc : sourceIds) { const auto &srcAuthName = idSrc.first; @@ -11041,20 +11051,8 @@ findOpsInRegistryDirect(const crs::CRSNNPtr &sourceCRS, const auto &targetAuthName = idTarget.first; const auto &targetCode = idTarget.second; - std::vector authorities; - if (authFactoryName == "any") { - authorities.emplace_back(); - } - if (authFactoryName.empty()) { - authorities = - authFactory->databaseContext()->getAllowedAuthorities( - srcAuthName, targetAuthName); - if (authorities.empty()) { - authorities.emplace_back(); - } - } else { - authorities.emplace_back(authFactoryName); - } + const auto authorities(getCandidateAuthorities( + authFactory, srcAuthName, targetAuthName)); for (const auto &authority : authorities) { const auto tmpAuthFactory = io::AuthorityFactory::create( authFactory->databaseContext(), @@ -11075,6 +11073,81 @@ findOpsInRegistryDirect(const crs::CRSNNPtr &sourceCRS, } return std::vector(); } + +// --------------------------------------------------------------------------- + +// Look in the authority registry for operations from sourceCRS +static std::vector +findOpsInRegistryDirectFrom(const crs::CRSNNPtr &sourceCRS, + const CoordinateOperationContextNNPtr &context) { + const auto &authFactory = context->getAuthorityFactory(); + assert(authFactory); + + std::list> ids; + buildCRSIds(sourceCRS, context, ids); + + for (const auto &id : ids) { + const auto &srcAuthName = id.first; + const auto &srcCode = id.second; + + const auto authorities( + getCandidateAuthorities(authFactory, srcAuthName, srcAuthName)); + for (const auto &authority : authorities) { + const auto tmpAuthFactory = io::AuthorityFactory::create( + authFactory->databaseContext(), + authority == "any" ? std::string() : authority); + auto res = tmpAuthFactory->createFromCoordinateReferenceSystemCodes( + srcAuthName, srcCode, std::string(), std::string(), + context->getUsePROJAlternativeGridNames(), + context->getGridAvailabilityUse() == + CoordinateOperationContext::GridAvailabilityUse:: + DISCARD_OPERATION_IF_MISSING_GRID, + context->getDiscardSuperseded()); + if (!res.empty()) { + return res; + } + } + } + return std::vector(); +} + +// --------------------------------------------------------------------------- + +// Look in the authority registry for operations to targetCRS +static std::vector +findOpsInRegistryDirectTo(const crs::CRSNNPtr &targetCRS, + const CoordinateOperationContextNNPtr &context) { + const auto &authFactory = context->getAuthorityFactory(); + assert(authFactory); + + std::list> ids; + buildCRSIds(targetCRS, context, ids); + + for (const auto &id : ids) { + const auto &targetAuthName = id.first; + const auto &targetCode = id.second; + + const auto authorities(getCandidateAuthorities( + authFactory, targetAuthName, targetAuthName)); + for (const auto &authority : authorities) { + const auto tmpAuthFactory = io::AuthorityFactory::create( + authFactory->databaseContext(), + authority == "any" ? std::string() : authority); + auto res = tmpAuthFactory->createFromCoordinateReferenceSystemCodes( + std::string(), std::string(), targetAuthName, targetCode, + context->getUsePROJAlternativeGridNames(), + context->getGridAvailabilityUse() == + CoordinateOperationContext::GridAvailabilityUse:: + DISCARD_OPERATION_IF_MISSING_GRID, + context->getDiscardSuperseded()); + if (!res.empty()) { + return res; + } + } + } + return std::vector(); +} + //! @endcond // --------------------------------------------------------------------------- @@ -11089,12 +11162,11 @@ static std::vector findsOpsInRegistryWithIntermediate( const auto &authFactory = context->getAuthorityFactory(); assert(authFactory); - const auto &authFactoryName = authFactory->getAuthority(); std::list> sourceIds; std::list> targetIds; - buildSourceAndTargetCRSIds(sourceCRS, targetCRS, context, sourceIds, - targetIds); + buildCRSIds(sourceCRS, context, sourceIds); + buildCRSIds(targetCRS, context, targetIds); for (const auto &idSrc : sourceIds) { const auto &srcAuthName = idSrc.first; @@ -11103,20 +11175,8 @@ static std::vector findsOpsInRegistryWithIntermediate( const auto &targetAuthName = idTarget.first; const auto &targetCode = idTarget.second; - std::vector authorities; - if (authFactoryName == "any") { - authorities.emplace_back(); - } - if (authFactoryName.empty()) { - authorities = - authFactory->databaseContext()->getAllowedAuthorities( - srcAuthName, targetAuthName); - if (authorities.empty()) { - authorities.emplace_back(); - } - } else { - authorities.emplace_back(authFactoryName); - } + const auto authorities(getCandidateAuthorities( + authFactory, srcAuthName, targetAuthName)); for (const auto &authority : authorities) { const auto tmpAuthFactory = io::AuthorityFactory::create( authFactory->databaseContext(), @@ -12139,6 +12199,51 @@ void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub( // --------------------------------------------------------------------------- +std::vector +CoordinateOperationFactory::Private::createOperationsGeogToVertWithIntermediate( + const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS + const crs::CRSNNPtr &targetCRS, // vertical CRS + Private::Context &context) { + + std::vector res; + + struct AntiRecursionGuard { + Context &context; + + explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) { + assert(!context.inCreateOperationsGeogToVertWithIntermediate); + context.inCreateOperationsGeogToVertWithIntermediate = true; + } + + ~AntiRecursionGuard() { + context.inCreateOperationsGeogToVertWithIntermediate = false; + } + }; + AntiRecursionGuard guard(context); + + for (int i = 0; i < 2; i++) { + + // Generally EPSG has operations from GeogCrs to VertCRS + auto ops = + i == 0 ? findOpsInRegistryDirectTo(targetCRS, context.context) + : findOpsInRegistryDirectFrom(targetCRS, context.context); + + for (const auto &op : ops) { + const auto tmpCRS = i == 0 ? op->sourceCRS() : op->targetCRS(); + if (tmpCRS && + dynamic_cast(tmpCRS.get())) { + res.emplace_back(i == 0 ? op : op->inverse()); + } + } + if (!res.empty()) + break; + } + + return res; +} + +// --------------------------------------------------------------------------- + static CoordinateOperationNNPtr createBallparkGeocentricTranslation(const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS) { @@ -12350,6 +12455,21 @@ CoordinateOperationFactory::Private::createOperations( } } + // There's no direct transformation from NAVD88 height to WGS84, + // so try to research all transformations from NAVD88 to another + // intermediate GeographicCRS. + if (res.empty() && + !context.inCreateOperationsGeogToVertWithIntermediate && + geogSrc && vertDst) { + res = createOperationsGeogToVertWithIntermediate( + sourceCRS, targetCRS, context); + } else if (res.empty() && + !context.inCreateOperationsGeogToVertWithIntermediate && + geogDst && vertSrc) { + res = applyInverse(createOperationsGeogToVertWithIntermediate( + targetCRS, sourceCRS, context)); + } + if (res.empty() && !sameGeodeticDatum && !context.inCreateOperationsWithDatumPivotAntiRecursion && geodSrc && geodDst) { @@ -12908,6 +13028,61 @@ CoordinateOperationFactory::Private::createOperations( componentsSrc[1]->extractVerticalCRS()) { verticalTransforms = createOperations(componentsSrc[1], targetCRS, context); + bool foundRegisteredTransformWithAllGridsAvailable = false; + for (const auto &op : verticalTransforms) { + if (!op->identifiers().empty() && authFactory) { + bool missingGrid = false; + const auto gridsNeeded = + op->gridsNeeded(authFactory->databaseContext()); + for (const auto &gridDesc : gridsNeeded) { + if (!gridDesc.available) { + missingGrid = true; + break; + } + } + if (!missingGrid) { + foundRegisteredTransformWithAllGridsAvailable = + true; + break; + } + } + } + if (!foundRegisteredTransformWithAllGridsAvailable && + srcGeogCRS && + !srcGeogCRS->_isEquivalentTo( + geogDst, util::IComparable::Criterion::EQUIVALENT) && + !srcGeogCRS->is2DPartOf3D(NN_NO_CHECK(geogDst))) { + auto verticalTransformsTmp = createOperations( + componentsSrc[1], NN_NO_CHECK(srcGeogCRS), context); + bool foundRegisteredTransform = false; + foundRegisteredTransformWithAllGridsAvailable = false; + for (const auto &op : verticalTransformsTmp) { + if (!op->identifiers().empty() && authFactory) { + bool missingGrid = false; + const auto gridsNeeded = + op->gridsNeeded(authFactory->databaseContext()); + for (const auto &gridDesc : gridsNeeded) { + if (!gridDesc.available) { + missingGrid = true; + break; + } + } + foundRegisteredTransform = true; + if (!missingGrid) { + foundRegisteredTransformWithAllGridsAvailable = + true; + break; + } + } + } + if (foundRegisteredTransformWithAllGridsAvailable) { + verticalTransforms = verticalTransformsTmp; + } else if (foundRegisteredTransform) { + verticalTransforms.insert(verticalTransforms.end(), + verticalTransformsTmp.begin(), + verticalTransformsTmp.end()); + } + } } if (!horizTransforms.empty() && !verticalTransforms.empty()) { for (const auto &horizTransform : horizTransforms) { -- cgit v1.2.3 From e6eae43cf2310c77a466fee257d9974b14ee85fd Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 12 Sep 2019 22:31:07 +0200 Subject: createOperations(): when tranforming from a compoundCRS whose vertical component is a BoundCRS, do not apply the horizontal transformation twice --- src/iso19111/coordinateoperation.cpp | 58 +++++++++++++++++++++++++----------- 1 file changed, 40 insertions(+), 18 deletions(-) (limited to 'src/iso19111/coordinateoperation.cpp') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index aad86410..aea8400c 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -10173,6 +10173,7 @@ struct CoordinateOperationFactory::Private { bool inCreateOperationsWithDatumPivotAntiRecursion = false; bool inCreateOperationsThroughPreferredHub = false; bool inCreateOperationsGeogToVertWithIntermediate = false; + bool skipHorizontalTransformation = false; Context(const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn, @@ -12784,27 +12785,32 @@ CoordinateOperationFactory::Private::createOperations( hubSrcGeog->coordinateSystem()->axisList().size() == 3 && geogDst->coordinateSystem()->axisList().size() == 3) { auto opsFirst = createOperations(sourceCRS, hubSrc, context); - auto opsSecond = createOperations(hubSrc, targetCRS, context); - if (!opsFirst.empty() && !opsSecond.empty()) { - for (const auto &opFirst : opsFirst) { - for (const auto &opLast : opsSecond) { - // Exclude artificial transformations from the hub - // to the target CRS - if (!opLast->hasBallparkTransformation()) { - try { - res.emplace_back( - ConcatenatedOperation:: - createComputeMetadata( - {opFirst, opLast}, - !allowEmptyIntersection)); - } catch ( - const InvalidOperationEmptyIntersection &) { + if (context.skipHorizontalTransformation) { + if (!opsFirst.empty()) + return opsFirst; + } else { + auto opsSecond = createOperations(hubSrc, targetCRS, context); + if (!opsFirst.empty() && !opsSecond.empty()) { + for (const auto &opFirst : opsFirst) { + for (const auto &opLast : opsSecond) { + // Exclude artificial transformations from the hub + // to the target CRS + if (!opLast->hasBallparkTransformation()) { + try { + res.emplace_back( + ConcatenatedOperation:: + createComputeMetadata( + {opFirst, opLast}, + !allowEmptyIntersection)); + } catch ( + const InvalidOperationEmptyIntersection &) { + } } } } - } - if (!res.empty()) { - return res; + if (!res.empty()) { + return res; + } } } } @@ -13026,6 +13032,22 @@ CoordinateOperationFactory::Private::createOperations( std::vector verticalTransforms; if (componentsSrc.size() >= 2 && componentsSrc[1]->extractVerticalCRS()) { + + struct SetSkipHorizontalTransform { + Context &context; + + explicit SetSkipHorizontalTransform(Context &contextIn) + : context(contextIn) { + assert(!context.skipHorizontalTransformation); + context.skipHorizontalTransformation = true; + } + + ~SetSkipHorizontalTransform() { + context.skipHorizontalTransformation = false; + } + }; + SetSkipHorizontalTransform setSkipHorizontalTransform(context); + verticalTransforms = createOperations(componentsSrc[1], targetCRS, context); bool foundRegisteredTransformWithAllGridsAvailable = false; -- cgit v1.2.3