Skip to content

Commit 869ca44

Browse files
phaarnesphaarnes
committed
Enhance vertical CRS operations chaining and add more tests
1 parent d03b464 commit 869ca44

2 files changed

Lines changed: 200 additions & 78 deletions

File tree

‎src/iso19111/operation/coordinateoperationfactory.cpp‎

Lines changed: 66 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -4584,74 +4584,86 @@ std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
45844584
return res;
45854585
const auto dbContext = authFactory->databaseContext().as_nullable();
45864586

4587+
// Try both pivot strategies and collect all candidates. The caller's
4588+
// filterAndSort will rank them (by accuracy, area of use, etc.)
4589+
45874590
// Strategy 1: pivot through candidates sharing the TARGET's datum.
4588-
// Find vertical CRSs on the same datum as the target, then look for
4589-
// operations from the source to each candidate (using the full
4590-
// createOperations machinery so multi-step chains are also found).
4591-
auto candidatesDst = findCandidateVertCRSForDatum(
4592-
authFactory, vertDst->datumNonNull(dbContext).get());
4593-
for (const auto &candidateVert : candidatesDst) {
4594-
if (candidateVert->_isEquivalentTo(
4595-
targetCRS.get(), util::IComparable::Criterion::EQUIVALENT)) {
4596-
continue; // Skip the target itself (already tried by the caller)
4597-
}
4598-
// Collect all registered source-to-candidate operations, which
4599-
// may differ in accuracy or area of use, for later ranking.
4600-
auto opsFirst = createOperations(sourceCRS, sourceEpoch, candidateVert,
4601-
sourceEpoch, context);
4602-
if (!opsFirst.empty()) {
4603-
// The candidate-to-target conversion is expected to be a
4604-
// trivial unit or axis change, so we only use the first
4605-
// result (opsSecond.front()).
4606-
const auto opsSecond = createOperations(
4607-
candidateVert, sourceEpoch, targetCRS, targetEpoch, context);
4608-
if (!opsSecond.empty()) {
4609-
for (const auto &opFirst : opsFirst) {
4610-
if (hasIdentifiers(opFirst)) {
4611-
try {
4612-
res.emplace_back(
4613-
ConcatenatedOperation::createComputeMetadata(
4614-
{opFirst, opsSecond.front()},
4615-
context.disallowEmptyIntersection()));
4616-
} catch (const InvalidOperationEmptyIntersection &) {
4591+
// source -> candidate (registered CT) -> target (trivial axis/unit change)
4592+
{
4593+
auto candidatesDst = findCandidateVertCRSForDatum(
4594+
authFactory, vertDst->datumNonNull(dbContext).get());
4595+
for (const auto &candidateVert : candidatesDst) {
4596+
if (candidateVert->_isEquivalentTo(
4597+
targetCRS.get(),
4598+
util::IComparable::Criterion::EQUIVALENT)) {
4599+
continue;
4600+
}
4601+
// Collect all registered source-to-candidate operations, which
4602+
// may differ in accuracy or area of use, for later ranking.
4603+
auto opsFirst = createOperations(sourceCRS, sourceEpoch,
4604+
candidateVert, sourceEpoch,
4605+
context);
4606+
if (!opsFirst.empty()) {
4607+
// The candidate-to-target conversion is expected to be a
4608+
// trivial unit or axis change, so we only use the first
4609+
// result (opsSecond.front()).
4610+
const auto opsSecond = createOperations(
4611+
candidateVert, sourceEpoch, targetCRS, targetEpoch,
4612+
context);
4613+
if (!opsSecond.empty()) {
4614+
for (const auto &opFirst : opsFirst) {
4615+
if (hasIdentifiers(opFirst)) {
4616+
try {
4617+
res.emplace_back(
4618+
ConcatenatedOperation::
4619+
createComputeMetadata(
4620+
{opFirst, opsSecond.front()},
4621+
context
4622+
.disallowEmptyIntersection()));
4623+
} catch (
4624+
const InvalidOperationEmptyIntersection &) {
4625+
}
46174626
}
46184627
}
46194628
}
4620-
if (!res.empty())
4621-
return res;
46224629
}
46234630
}
46244631
}
46254632

46264633
// Strategy 2: pivot through candidates sharing the SOURCE's datum.
46274634
// Find vertical CRSs on the same datum as the source, then look for
46284635
// operations from each candidate to the target.
4629-
auto candidatesSrc = findCandidateVertCRSForDatum(
4630-
authFactory, vertSrc->datumNonNull(dbContext).get());
4631-
for (const auto &candidateVert : candidatesSrc) {
4632-
if (candidateVert->_isEquivalentTo(
4633-
sourceCRS.get(), util::IComparable::Criterion::EQUIVALENT)) {
4634-
continue;
4635-
}
4636-
auto opsSecond = createOperations(candidateVert, sourceEpoch, targetCRS,
4637-
targetEpoch, context);
4638-
if (!opsSecond.empty()) {
4639-
const auto opsFirst = createOperations(
4640-
sourceCRS, sourceEpoch, candidateVert, sourceEpoch, context);
4641-
if (!opsFirst.empty()) {
4642-
for (const auto &opSecond : opsSecond) {
4643-
if (hasIdentifiers(opSecond)) {
4644-
try {
4645-
res.emplace_back(
4646-
ConcatenatedOperation::createComputeMetadata(
4647-
{opsFirst.front(), opSecond},
4648-
context.disallowEmptyIntersection()));
4649-
} catch (const InvalidOperationEmptyIntersection &) {
4636+
{
4637+
auto candidatesSrc = findCandidateVertCRSForDatum(
4638+
authFactory, vertSrc->datumNonNull(dbContext).get());
4639+
for (const auto &candidateVert : candidatesSrc) {
4640+
if (candidateVert->_isEquivalentTo(
4641+
sourceCRS.get(),
4642+
util::IComparable::Criterion::EQUIVALENT)) {
4643+
continue;
4644+
}
4645+
auto opsSecond = createOperations(candidateVert, sourceEpoch,
4646+
targetCRS, targetEpoch, context);
4647+
if (!opsSecond.empty()) {
4648+
const auto opsFirst = createOperations(
4649+
sourceCRS, sourceEpoch, candidateVert, sourceEpoch,
4650+
context);
4651+
if (!opsFirst.empty()) {
4652+
for (const auto &opSecond : opsSecond) {
4653+
if (hasIdentifiers(opSecond)) {
4654+
try {
4655+
res.emplace_back(
4656+
ConcatenatedOperation::
4657+
createComputeMetadata(
4658+
{opsFirst.front(), opSecond},
4659+
context
4660+
.disallowEmptyIntersection()));
4661+
} catch (
4662+
const InvalidOperationEmptyIntersection &) {
4663+
}
46504664
}
46514665
}
46524666
}
4653-
if (!res.empty())
4654-
return res;
46554667
}
46564668
}
46574669
}

‎test/unit/test_operationfactory.cpp‎

Lines changed: 134 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -7279,57 +7279,167 @@ TEST(operation, vertCRS_to_vertCRS_height_depth_pivot_context) {
72797279
// Test that PROJ can chain a registered vertical CT with a
72807280
// height-to-depth axis conversion when the target CRS differs from
72817281
// the CT's registered target only by axis direction.
7282+
auto authFactory =
7283+
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
7284+
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
7285+
7286+
auto checkPipeline = [&](const std::string &srcCode,
7287+
const std::string &tgtCode,
7288+
const std::string &expectedProj) {
7289+
auto list = CoordinateOperationFactory::create()->createOperations(
7290+
authFactory->createCoordinateReferenceSystem(srcCode),
7291+
authFactory->createCoordinateReferenceSystem(tgtCode), ctxt);
7292+
ASSERT_GE(list.size(), 1U);
7293+
EXPECT_FALSE(list[0]->hasBallparkTransformation());
7294+
EXPECT_EQ(
7295+
list[0]->exportToPROJString(PROJStringFormatter::create().get()),
7296+
expectedProj);
7297+
};
7298+
7299+
// Caspian: EPSG:5705 (Baltic 1977 height) -> EPSG:5706 (Caspian depth)
7300+
// via EPSG:5438 (dh=28) + height-to-depth axisswap
7301+
checkPipeline("5705", "5706",
7302+
"+proj=pipeline +step +proj=geogoffset +dh=28 "
7303+
"+step +proj=axisswap +order=1,2,-3");
7304+
7305+
// AIOC95: EPSG:5705 -> EPSG:5734 (AIOC95 depth)
7306+
// via EPSG:5443 (dh=26.3) + height-to-depth axisswap
7307+
checkPipeline("5705", "5734",
7308+
"+proj=pipeline +step +proj=geogoffset +dh=26.3 "
7309+
"+step +proj=axisswap +order=1,2,-3");
7310+
7311+
// KOC: EPSG:5790 (KOC CD height) -> EPSG:5789 (KOC WD depth)
7312+
// via EPSG:7986 (dh=-4.74) + height-to-depth axisswap
7313+
checkPipeline("5790", "5789",
7314+
"+proj=pipeline +step +proj=geogoffset +dh=-4.74 "
7315+
"+step +proj=axisswap +order=1,2,-3");
7316+
7317+
// Kuwait PWD: EPSG:5788 -> EPSG:5789 (KOC WD depth)
7318+
// via EPSG:7981 (dh=-4.25) + height-to-depth axisswap
7319+
checkPipeline("5788", "5789",
7320+
"+proj=pipeline +step +proj=geogoffset +dh=-4.25 "
7321+
"+step +proj=axisswap +order=1,2,-3");
7322+
7323+
// KOC ft: EPSG:5790 -> EPSG:5614 (KOC WD depth ft)
7324+
// via EPSG:7987 (dh=-4.74) + axisswap + unit conversion m->ft
7325+
checkPipeline("5790", "5614",
7326+
"+proj=pipeline "
7327+
"+step +proj=geogoffset +dh=-4.74 "
7328+
"+step +proj=axisswap +order=1,2,-3 "
7329+
"+step +proj=unitconvert +z_in=m +z_out=ft");
7330+
}
7331+
7332+
// ---------------------------------------------------------------------------
7333+
7334+
TEST(operation, vertCRS_to_vertCRS_depth_height_pivot_context) {
7335+
// Test the reverse direction: depth to height.
7336+
// When the source is depth (axis DOWN) and the target is height (axis UP),
7337+
// the code prefers pivoting through the SOURCE's datum, composing:
7338+
// depth-to-height axisswap + inverse of the registered forward op.
7339+
auto authFactory =
7340+
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
7341+
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
7342+
7343+
auto checkPipeline = [&](const std::string &srcCode,
7344+
const std::string &tgtCode,
7345+
const std::string &expectedProj) {
7346+
auto list = CoordinateOperationFactory::create()->createOperations(
7347+
authFactory->createCoordinateReferenceSystem(srcCode),
7348+
authFactory->createCoordinateReferenceSystem(tgtCode), ctxt);
7349+
ASSERT_GE(list.size(), 1U);
7350+
EXPECT_FALSE(list[0]->hasBallparkTransformation());
7351+
EXPECT_EQ(
7352+
list[0]->exportToPROJString(PROJStringFormatter::create().get()),
7353+
expectedProj);
7354+
};
7355+
7356+
// Caspian: EPSG:5706 (Caspian depth) -> EPSG:5705 (Baltic 1977 height)
7357+
// axisswap + inverse of EPSG:5438 (dh=-28)
7358+
checkPipeline("5706", "5705",
7359+
"+proj=pipeline +step +proj=axisswap +order=1,2,-3 "
7360+
"+step +proj=geogoffset +dh=-28");
7361+
7362+
// AIOC95: EPSG:5734 (AIOC95 depth) -> EPSG:5705 (Baltic 1977 height)
7363+
// axisswap + inverse of EPSG:5443 (dh=-26.3)
7364+
checkPipeline("5734", "5705",
7365+
"+proj=pipeline +step +proj=axisswap +order=1,2,-3 "
7366+
"+step +proj=geogoffset +dh=-26.3");
7367+
7368+
// KOC: EPSG:5789 (KOC WD depth) -> EPSG:5790 (KOC CD height)
7369+
// axisswap + inverse of EPSG:7986 (dh=4.74)
7370+
checkPipeline("5789", "5790",
7371+
"+proj=pipeline +step +proj=axisswap +order=1,2,-3 "
7372+
"+step +proj=geogoffset +dh=4.74");
7373+
7374+
// Kuwait PWD: EPSG:5789 (KOC WD depth) -> EPSG:5788 (Kuwait PWD height)
7375+
// axisswap + inverse of EPSG:7981 (dh=4.25)
7376+
checkPipeline("5789", "5788",
7377+
"+proj=pipeline +step +proj=axisswap +order=1,2,-3 "
7378+
"+step +proj=geogoffset +dh=4.25");
7379+
7380+
// KOC ft: EPSG:5614 (KOC WD depth ft) -> EPSG:5790 (KOC CD height)
7381+
// unit conversion ft->m + axisswap + inverse of EPSG:7987 (dh=4.74)
7382+
checkPipeline("5614", "5790",
7383+
"+proj=pipeline "
7384+
"+step +proj=unitconvert +z_in=ft +z_out=m "
7385+
"+step +proj=axisswap +order=1,2,-3 "
7386+
"+step +proj=geogoffset +dh=4.74");
7387+
}
7388+
7389+
// ---------------------------------------------------------------------------
7390+
7391+
TEST(operation, vertCRS_to_vertCRS_height_depth_pivot_blacksea_context) {
7392+
// Test that the intermediate-vert pivot finds a registered cross-datum
7393+
// height-to-height transformation and composes it with an axis reversal.
72827394
//
7283-
// EPSG:5705 (Baltic 1977 height) to EPSG:5706 (Caspian depth)
7284-
// should compose: EPSG:5438 (5705 to 5611, geogoffset +dh=28)
7285-
// + height-to-depth (axisswap order=1,2,-3)
7395+
// EPSG:5705 (Baltic 1977 height) to EPSG:5336 (Black Sea depth)
7396+
// Strategy 1 composes: EPSG:5447 (5705 to 5735, geogoffset +dh=0.4)
7397+
// + height-to-depth (axisswap order=1,2,-3)
7398+
// PARTIAL_INTERSECTION is needed because EPSG:5447's area
7399+
// is slightly smaller than EPSG:5336's extent
72867400
auto authFactory =
72877401
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
72887402
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
7403+
ctxt->setSpatialCriterion(
7404+
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
72897405
auto list = CoordinateOperationFactory::create()->createOperations(
72907406
// Baltic 1977 height
72917407
authFactory->createCoordinateReferenceSystem("5705"),
7292-
// Caspian depth (same datum as Caspian height 5611, but axis: down)
7293-
authFactory->createCoordinateReferenceSystem("5706"), ctxt);
7408+
// Black Sea depth
7409+
authFactory->createCoordinateReferenceSystem("5336"), ctxt);
72947410
ASSERT_GE(list.size(), 1U);
7295-
// Must NOT be a ballpark transformation
72967411
EXPECT_FALSE(list[0]->hasBallparkTransformation());
7297-
// The pipeline should match the registered operation followed by the
7298-
// height-to-depth axis conversion.
72997412
auto projStr =
73007413
list[0]->exportToPROJString(PROJStringFormatter::create().get());
7301-
EXPECT_EQ(projStr, "+proj=pipeline +step +proj=geogoffset +dh=28 +step "
7414+
EXPECT_EQ(projStr, "+proj=pipeline +step +proj=geogoffset +dh=0.4 +step "
73027415
"+proj=axisswap +order=1,2,-3");
73037416
}
73047417

73057418
// ---------------------------------------------------------------------------
73067419

7307-
TEST(operation, vertCRS_to_vertCRS_depth_height_pivot_strategy2_context) {
7308-
// Test Strategy 2: pivot through a candidate sharing the SOURCE's datum.
7420+
TEST(operation, vertCRS_to_vertCRS_depth_height_pivot_blacksea_context) {
7421+
// Test the reverse direction for the Black Sea case.
73097422
//
7310-
// EPSG:5706 (Caspian depth) to EPSG:5705 (Baltic 1977 height)
7311-
// The registered operation EPSG:5438 goes 5705 to 5611 (Caspian height).
7312-
// Strategy 2 finds 5611 (shares Caspian datum with source 5706), discovers
7313-
// the reverse of EPSG:5438 from 5611 to 5705, and composes:
7314-
// 5706 -(depth-to-height conversion)-> 5611 -(inverse of EPSG:5438)->
7315-
// 5705
7423+
// EPSG:5336 (Black Sea depth) to EPSG:5705 (Baltic 1977 height)
7424+
// Strategy 2 composes: depth-to-height (5336 to 5735, axisswap)
7425+
// + inverse of EPSG:5447 (5735 to 5705, dh=-0.4)
7426+
73167427
auto authFactory =
73177428
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
73187429
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
7430+
ctxt->setSpatialCriterion(
7431+
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
73197432
auto list = CoordinateOperationFactory::create()->createOperations(
7320-
// Caspian depth
7321-
authFactory->createCoordinateReferenceSystem("5706"),
7433+
// Black Sea depth
7434+
authFactory->createCoordinateReferenceSystem("5336"),
73227435
// Baltic 1977 height
73237436
authFactory->createCoordinateReferenceSystem("5705"), ctxt);
73247437
ASSERT_GE(list.size(), 1U);
7325-
// Must NOT be a ballpark transformation
73267438
EXPECT_FALSE(list[0]->hasBallparkTransformation());
7327-
// The exported pipeline is normalized to the same PROJ steps as the
7328-
// forward case.
73297439
auto projStr =
73307440
list[0]->exportToPROJString(PROJStringFormatter::create().get());
7331-
EXPECT_EQ(projStr, "+proj=pipeline +step +proj=geogoffset +dh=28 +step "
7332-
"+proj=axisswap +order=1,2,-3");
7441+
EXPECT_EQ(projStr, "+proj=pipeline +step +proj=axisswap +order=1,2,-3 "
7442+
"+step +proj=geogoffset +dh=-0.4");
73337443
}
73347444

73357445
// ---------------------------------------------------------------------------

0 commit comments

Comments
 (0)