From 0316dc2b42f8354b86e6cc54e2f44d3e29944e57 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') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 9d20b396..74b319ba 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 29a5cae676cf0bcd8226933512ac23a91faa6654 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 | 317 +++++++++++++++++++++++++++-------- src/iso19111/factory.cpp | 69 ++++---- 2 files changed, 283 insertions(+), 103 deletions(-) (limited to 'src') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 74b319ba..d5b766ad 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 inCreateOperationsGeogToVertWithIntermediate = false; Context(const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn, @@ -10195,6 +10196,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); @@ -10943,18 +10949,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(), @@ -10980,35 +10989,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; } // --------------------------------------------------------------------------- @@ -11020,12 +11031,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; @@ -11034,20 +11044,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(), @@ -11068,6 +11066,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 // --------------------------------------------------------------------------- @@ -11082,12 +11155,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; @@ -11096,20 +11168,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(), @@ -12016,6 +12076,51 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot( // --------------------------------------------------------------------------- +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) { @@ -12132,11 +12237,16 @@ CoordinateOperationFactory::Private::createOperations( auto geodSrc = dynamic_cast(sourceCRS.get()); auto geodDst = dynamic_cast(targetCRS.get()); + auto geogSrc = dynamic_cast(sourceCRS.get()); + auto geogDst = dynamic_cast(targetCRS.get()); + auto vertSrc = dynamic_cast(sourceCRS.get()); + auto vertDst = dynamic_cast(targetCRS.get()); // First look-up if the registry provide us with operations. auto derivedSrc = dynamic_cast(sourceCRS.get()); auto derivedDst = dynamic_cast(targetCRS.get()); - if (context.context->getAuthorityFactory() && + auto authFactory = context.context->getAuthorityFactory(); + if (authFactory && (derivedSrc == nullptr || !derivedSrc->baseCRS()->_isEquivalentTo( targetCRS.get(), util::IComparable::Criterion::EQUIVALENT)) && @@ -12170,6 +12280,21 @@ CoordinateOperationFactory::Private::createOperations( util::IComparable::Criterion::EQUIVALENT); } + // 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) { @@ -12228,10 +12353,6 @@ CoordinateOperationFactory::Private::createOperations( "celestial body"); } - auto geogSrc = - dynamic_cast(sourceCRS.get()); - auto geogDst = - dynamic_cast(targetCRS.get()); if (geogSrc && geogDst) { return createOperationsGeogToGeog(res, sourceCRS, targetCRS, geogSrc, geogDst); @@ -12319,7 +12440,6 @@ CoordinateOperationFactory::Private::createOperations( return applyInverse(createOperations(targetCRS, sourceCRS, context)); } - auto geogDst = dynamic_cast(targetCRS.get()); if (boundSrc && geogDst) { const auto &hubSrc = boundSrc->hubCRS(); auto hubSrcGeog = @@ -12508,14 +12628,12 @@ CoordinateOperationFactory::Private::createOperations( } // reverse of previous case - auto geogSrc = dynamic_cast(sourceCRS.get()); if (geogSrc && boundDst) { return applyInverse(createOperations(targetCRS, sourceCRS, context)); } // vertCRS (as boundCRS with transformation to target vertCRS) to // vertCRS - auto vertDst = dynamic_cast(targetCRS.get()); if (boundSrc && vertDst) { auto baseSrcVert = dynamic_cast(boundSrc->baseCRS().get()); @@ -12532,7 +12650,6 @@ CoordinateOperationFactory::Private::createOperations( } // reverse of previous case - auto vertSrc = dynamic_cast(sourceCRS.get()); if (boundDst && vertSrc) { return applyInverse(createOperations(targetCRS, sourceCRS, context)); } @@ -12576,7 +12693,6 @@ CoordinateOperationFactory::Private::createOperations( if (vertSrc && geogDst) { if (vertSrc->identifiers().empty()) { - const auto &authFactory = context.context->getAuthorityFactory(); const auto &vertSrcName = vertSrc->nameStr(); if (authFactory != nullptr && vertSrcName != "unnamed" && vertSrcName != "unknown") { @@ -12727,6 +12843,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) { diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 1f40f1f0..0f6790c0 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -3312,9 +3312,9 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes( bool discardSuperseded) const { auto cacheKey(d->authority()); - cacheKey += sourceCRSAuthName; + cacheKey += sourceCRSAuthName.empty() ? "{empty}" : sourceCRSAuthName; cacheKey += sourceCRSCode; - cacheKey += targetCRSAuthName; + cacheKey += targetCRSAuthName.empty() ? "{empty}" : targetCRSAuthName; cacheKey += targetCRSCode; cacheKey += (usePROJAlternativeGridNames ? '1' : '0'); cacheKey += (discardIfMissingGrid ? '1' : '0'); @@ -3326,27 +3326,30 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes( return list; } - // Look-up first for conversion which is the most precise. - std::string sql("SELECT conversion_auth_name, " - "geodetic_crs_auth_name, geodetic_crs_code FROM " - "projected_crs WHERE auth_name = ? AND code = ?"); - auto params = ListOfParams{targetCRSAuthName, targetCRSCode}; - auto res = d->run(sql, params); - if (!res.empty()) { - const auto &row = res.front(); - bool ok = row[1] == sourceCRSAuthName && row[2] == sourceCRSCode; - if (ok && d->hasAuthorityRestriction()) { - ok = row[0] == d->authority(); - } - if (ok) { - auto targetCRS = d->createFactory(targetCRSAuthName) - ->createProjectedCRS(targetCRSCode); - auto conv = targetCRS->derivingConversion(); - list.emplace_back(conv); - d->context()->d->cache(cacheKey, list); - return list; + if (!targetCRSAuthName.empty()) { + // Look-up first for conversion which is the most precise. + std::string sql("SELECT conversion_auth_name, " + "geodetic_crs_auth_name, geodetic_crs_code FROM " + "projected_crs WHERE auth_name = ? AND code = ?"); + auto params = ListOfParams{targetCRSAuthName, targetCRSCode}; + auto res = d->run(sql, params); + if (!res.empty()) { + const auto &row = res.front(); + bool ok = row[1] == sourceCRSAuthName && row[2] == sourceCRSCode; + if (ok && d->hasAuthorityRestriction()) { + ok = row[0] == d->authority(); + } + if (ok) { + auto targetCRS = d->createFactory(targetCRSAuthName) + ->createProjectedCRS(targetCRSCode); + auto conv = targetCRS->derivingConversion(); + list.emplace_back(conv); + d->context()->d->cache(cacheKey, list); + return list; + } } } + std::string sql; if (discardSuperseded) { sql = "SELECT cov.auth_name, cov.code, cov.table_name, " "ss.replacement_auth_name, ss.replacement_code FROM " @@ -3358,20 +3361,26 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes( "ss.superseded_auth_name = cov.auth_name AND " "ss.superseded_code = cov.code AND " "ss.superseded_table_name = ss.replacement_table_name " - "WHERE source_crs_auth_name = ? AND source_crs_code = ? AND " - "target_crs_auth_name = ? AND target_crs_code = ? AND " - "cov.deprecated = 0"; + "WHERE "; } else { sql = "SELECT cov.auth_name, cov.code, cov.table_name FROM " "coordinate_operation_view cov JOIN area " "ON cov.area_of_use_auth_name = area.auth_name AND " "cov.area_of_use_code = area.code " - "WHERE source_crs_auth_name = ? AND source_crs_code = ? AND " - "target_crs_auth_name = ? AND target_crs_code = ? AND " - "cov.deprecated = 0"; + "WHERE "; + } + ListOfParams params; + if (!sourceCRSAuthName.empty()) { + sql += "source_crs_auth_name = ? AND source_crs_code = ? AND "; + params.emplace_back(sourceCRSAuthName); + params.emplace_back(sourceCRSCode); + } + if (!targetCRSAuthName.empty()) { + sql += "target_crs_auth_name = ? AND target_crs_code = ? AND "; + params.emplace_back(targetCRSAuthName); + params.emplace_back(targetCRSCode); } - params = {sourceCRSAuthName, sourceCRSCode, targetCRSAuthName, - targetCRSCode}; + sql += "cov.deprecated = 0"; if (d->hasAuthorityRestriction()) { sql += " AND cov.auth_name = ?"; params.emplace_back(d->authority()); @@ -3379,7 +3388,7 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes( sql += " ORDER BY pseudo_area_from_swne(south_lat, west_lon, north_lat, " "east_lon) DESC, " "(CASE WHEN accuracy is NULL THEN 1 ELSE 0 END), accuracy"; - res = d->run(sql, params); + auto res = d->run(sql, params); std::set> setTransf; if (discardSuperseded) { for (const auto &row : res) { -- cgit v1.2.3 From f168e2397eb7ac01bb1ec5bd1f08e27ff0a8a8d2 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') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index d5b766ad..f0ebaebd 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 inCreateOperationsGeogToVertWithIntermediate = false; + bool skipHorizontalTransformation = false; Context(const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn, @@ -12599,27 +12600,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; + } } } } @@ -12841,6 +12847,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