diff --git a/docs/source/operations/operations_computation.rst b/docs/source/operations/operations_computation.rst index 86095c14bc..527e967066 100644 --- a/docs/source/operations/operations_computation.rst +++ b/docs/source/operations/operations_computation.rst @@ -566,7 +566,15 @@ Overall logic is similar to the above case. There might be direct operations in the PROJ database, involving grid transformations or simple offsets. The fallback case is to synthesize a ballpark transformation. -This is implemented by the ``createOperationsVertToVert`` method +When no direct operation is found, a search is made for intermediate vertical CRS +entries sharing the same datum as the source or target. If a registered operation +exists between the source (or target) and such an intermediate CRS, it is composed +with the axis/unit conversion to the actual target (or from the actual source). +This is implemented by the ``createOperationsVertToVertWithIntermediateVert`` +method and avoids ballpark fallback for CRS pairs that differ only in axis +direction (height vs depth), units (e.g. metres vs feet), or both. + +The direct case is handled by the ``createOperationsVertToVert`` method. .. code-block:: shell diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index d4c70f0d05..c0c165f1a0 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -560,6 +560,7 @@ struct CoordinateOperationFactory::Private { bool inCreateOperationsWithDatumPivotAntiRecursion = false; bool inCreateOperationsGeogToVertWithAlternativeGeog = false; bool inCreateOperationsGeogToVertWithIntermediateVert = false; + bool inCreateOperationsVertToVertWithIntermediateVert = false; bool skipHorizontalTransformation = false; int nRecLevelCreateOperations = 0; std::map, @@ -653,6 +654,15 @@ struct CoordinateOperationFactory::Private { const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, Context &context); + static std::vector + createOperationsVertToVertWithIntermediateVert( + const crs::CRSNNPtr &sourceCRS, + const util::optional &sourceEpoch, + const crs::CRSNNPtr &targetCRS, + const util::optional &targetEpoch, + const crs::VerticalCRS *vertSrc, const crs::VerticalCRS *vertDst, + Context &context); + static void createOperationsFromDatabaseWithVertCRS( const crs::CRSNNPtr &sourceCRS, const util::optional &sourceEpoch, @@ -4592,15 +4602,17 @@ std::vector CoordinateOperationFactory::Private:: auto candidatesVert = findCandidateVertCRSForDatum( authFactory, vertDst->datumNonNull(dbContext).get()); for (const auto &candidateVert : candidatesVert) { + // Collect all registered source-to-candidate operations, which + // may differ in accuracy or area of use, for later ranking. auto resTmp = createOperations(sourceCRS, sourceEpoch, candidateVert, sourceEpoch, context); if (!resTmp.empty()) { + // The candidate-to-target conversion is expected to be a + // trivial unit or axis change, so we only use the first + // result (opsSecond.front()). const auto opsSecond = createOperations( candidateVert, sourceEpoch, targetCRS, targetEpoch, context); if (!opsSecond.empty()) { - // The transformation from candidateVert to targetCRS should - // be just a unit change typically, so take only the first one, - // which is likely/hopefully the only one. for (const auto &opFirst : resTmp) { if (hasIdentifiers(opFirst)) { if (candidateVert->_isEquivalentTo( @@ -4626,6 +4638,136 @@ std::vector CoordinateOperationFactory::Private:: // --------------------------------------------------------------------------- +// When transforming between two vertical CRS with different datums, check +// if there are registered operations whose target (or source) is a different +// vertical CRS sharing the same datum as our target (or source). If so, +// chain: source to intermediate + intermediate to target (the latter being a +// simple unit/axis conversion handled by createOperationsVertToVert). +// +// Typical example: Baltic 1977 height (EPSG:5705) to Caspian depth (EPSG:5706). +// EPSG has an operation 5705 to 5611 (Caspian height). 5611 and 5706 share the +// same datum but differ only in axis direction. We compose: +// 5705 ->(EPSG:5438)-> 5611 ->(height-to-depth)-> 5706 +std::vector CoordinateOperationFactory::Private:: + createOperationsVertToVertWithIntermediateVert( + const crs::CRSNNPtr &sourceCRS, + const util::optional &sourceEpoch, + const crs::CRSNNPtr &targetCRS, + const util::optional &targetEpoch, + const crs::VerticalCRS *vertSrc, const crs::VerticalCRS *vertDst, + Private::Context &context) { + + ENTER_FUNCTION(); + + std::vector res; + + struct AntiRecursionGuard { + Context &context; + + explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) { + assert(!context.inCreateOperationsVertToVertWithIntermediateVert); + context.inCreateOperationsVertToVertWithIntermediateVert = true; + } + + ~AntiRecursionGuard() { + context.inCreateOperationsVertToVertWithIntermediateVert = false; + } + }; + AntiRecursionGuard guard(context); + + const auto &authFactory = context.context->getAuthorityFactory(); + if (!authFactory) + return res; + const auto dbContext = authFactory->databaseContext().as_nullable(); + + // Try both pivot strategies and collect all candidates. The caller's + // filterAndSort will rank them (by accuracy, area of use, etc.) + + // Strategy 1: pivot through candidates sharing the TARGET's datum. + // source -> candidate (registered CT) -> target (trivial axis/unit change) + { + auto candidatesDst = findCandidateVertCRSForDatum( + authFactory, vertDst->datumNonNull(dbContext).get()); + for (const auto &candidateVert : candidatesDst) { + if (candidateVert->_isEquivalentTo( + targetCRS.get(), + util::IComparable::Criterion::EQUIVALENT)) { + continue; + } + // Collect all registered source-to-candidate operations, which + // may differ in accuracy or area of use, for later ranking. + auto opsFirst = createOperations( + sourceCRS, sourceEpoch, candidateVert, sourceEpoch, context); + if (!opsFirst.empty()) { + // The candidate-to-target conversion is expected to be a + // trivial unit or axis change, so we only use the first + // result (opsSecond.front()). + const auto opsSecond = + createOperations(candidateVert, sourceEpoch, targetCRS, + targetEpoch, context); + if (!opsSecond.empty()) { + for (const auto &opFirst : opsFirst) { + if (hasIdentifiers(opFirst)) { + try { + res.emplace_back( + ConcatenatedOperation:: + createComputeMetadata( + {opFirst, opsSecond.front()}, + context + .disallowEmptyIntersection())); + } catch ( + const InvalidOperationEmptyIntersection &) { + } + } + } + } + } + } + } + + // Strategy 2: pivot through candidates sharing the SOURCE's datum. + // Find vertical CRSs on the same datum as the source, then look for + // operations from each candidate to the target. + if (res.empty()) { + auto candidatesSrc = findCandidateVertCRSForDatum( + authFactory, vertSrc->datumNonNull(dbContext).get()); + for (const auto &candidateVert : candidatesSrc) { + if (candidateVert->_isEquivalentTo( + sourceCRS.get(), + util::IComparable::Criterion::EQUIVALENT)) { + continue; + } + auto opsSecond = createOperations(candidateVert, sourceEpoch, + targetCRS, targetEpoch, context); + if (!opsSecond.empty()) { + const auto opsFirst = + createOperations(sourceCRS, sourceEpoch, candidateVert, + sourceEpoch, context); + if (!opsFirst.empty()) { + for (const auto &opSecond : opsSecond) { + if (hasIdentifiers(opSecond)) { + try { + res.emplace_back( + ConcatenatedOperation:: + createComputeMetadata( + {opsFirst.front(), opSecond}, + context + .disallowEmptyIntersection())); + } catch ( + const InvalidOperationEmptyIntersection &) { + } + } + } + } + } + } + } + + return res; +} + +// --------------------------------------------------------------------------- + std::vector CoordinateOperationFactory::Private:: createOperationsGeogToVertWithAlternativeGeog( const crs::CRSNNPtr &sourceCRS, // geographic CRS @@ -4759,6 +4901,19 @@ void CoordinateOperationFactory::Private:: res = applyInverse(res); } + // Typically to transform between two vertical CRS with different datums + // (e.g. "Baltic 1977 height" to "Caspian depth") by pivoting through + // an intermediate vertical CRS sharing the target's (or source's) datum. + // This handles the case where a registered CT targets a height CRS but + // the user requests the corresponding depth CRS (or vice versa). + if (res.empty() && + !context.inCreateOperationsVertToVertWithIntermediateVert && vertSrc && + vertDst) { + res = createOperationsVertToVertWithIntermediateVert( + sourceCRS, sourceEpoch, targetCRS, targetEpoch, vertSrc, vertDst, + context); + } + // There's no direct transformation from NAVD88 height to WGS84, // so try to research all transformations from NAVD88 to another // intermediate GeographicCRS. diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 35d8141eaf..846416778b 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -7275,6 +7275,82 @@ TEST(operation, vertCRS_to_vertCRS_New_Zealand_context) { // --------------------------------------------------------------------------- +TEST(operation, vertCRS_to_vertCRS_pivot_context) { + // Test that PROJ can chain a registered vertical CT with a + // height-to-depth axis conversion when the target CRS differs from + // the CT's registered target only by axis direction. + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + + auto checkPipeline = [&](const std::string &srcCode, + const std::string &tgtCode, + const std::string &expectedProj) { + auto list = CoordinateOperationFactory::create()->createOperations( + authFactory->createCoordinateReferenceSystem(srcCode), + authFactory->createCoordinateReferenceSystem(tgtCode), ctxt); + ASSERT_GE(list.size(), 1U); + EXPECT_FALSE(list[0]->hasBallparkTransformation()); + EXPECT_EQ( + list[0]->exportToPROJString(PROJStringFormatter::create().get()), + expectedProj); + }; + + // Using Strategy 1 of createOperationsVertToVertWithIntermediateVert() + + // Caspian: EPSG:5705 (Baltic 1977 height) -> EPSG:5706 (Caspian depth) + // via EPSG:5438 (dh=28) + height-to-depth axisswap + checkPipeline("5705", "5706", + "+proj=pipeline +step +proj=geogoffset +dh=28 " + "+step +proj=axisswap +order=1,2,-3"); + + // Caspian: EPSG:5706 (Caspian depth) -> EPSG:5705 (Baltic 1977 height) + // axisswap + inverse of EPSG:5438 (dh=-28) + // Note: yes that pipeline is identical to the above one, since it is its + // own inverse. + checkPipeline("5706", "5705", + "+proj=pipeline +step +proj=geogoffset +dh=28 " + "+step +proj=axisswap +order=1,2,-3"); + + // KOC ft: EPSG:5790 -> EPSG:5614 (KOC WD depth ft) + // via EPSG:7987 (dh=-4.74) + axisswap + unit conversion m->ft + checkPipeline("5790", "5614", + "+proj=pipeline " + "+step +proj=geogoffset +dh=-4.74 " + "+step +proj=axisswap +order=1,2,-3 " + "+step +proj=unitconvert +z_in=m +z_out=ft"); + + // KOC ft: EPSG:5614 (KOC WD depth ft) -> EPSG:5790 (KOC CD height) + // unit conversion ft->m + axisswap + inverse of EPSG:7987 (dh=4.74) + checkPipeline("5614", "5790", + "+proj=pipeline " + "+step +proj=unitconvert +z_in=ft +z_out=m " + "+step +proj=axisswap +order=1,2,-3 " + "+step +proj=geogoffset +dh=4.74"); + + // EPSG:5705 (Baltic 1977 height) to EPSG:5336 (Black Sea depth) + // Strategy 1 composes: EPSG:5447 (5705 to 5735, geogoffset +dh=0.4) + // + height-to-depth (axisswap order=1,2,-3) + checkPipeline("5705", "5336", + "+proj=pipeline " + "+step +proj=geogoffset +dh=0.4 " + "+step +proj=axisswap +order=1,2,-3"); + + // Using Strategy 2 of createOperationsVertToVertWithIntermediateVert() + + // EPSG:5336 (Black Sea depth) to EPSG:5705 (Baltic 1977 height) + // Strategy 2 composes: depth-to-height (5336 to 5735, axisswap) + // + inverse of EPSG:5447 (5735 to 5705, dh=-0.4) + checkPipeline("5336", "5705", + "+proj=pipeline " + "+step +proj=axisswap +order=1,2,-3 " + "+step +proj=geogoffset +dh=-0.4"); +} + +// --------------------------------------------------------------------------- + TEST(operation, projCRS_3D_to_geogCRS_3D) { auto compoundcrs_ft_obj = PROJStringParser().createFromPROJString(