Skip to content

Commit

Permalink
createOperations(): do not stop at the first operation in the PROJ na…
Browse files Browse the repository at this point in the history
…mespace for vertical transformations

In particular helps with transformation between "NAD83 + NAVD88 height" and WGS 84 that
have regressed in 8.2.0

Fixes #2936
  • Loading branch information
rouault committed Nov 11, 2021
1 parent f7e9db5 commit b0d8717
Show file tree
Hide file tree
Showing 2 changed files with 71 additions and 12 deletions.
36 changes: 24 additions & 12 deletions src/iso19111/operation/coordinateoperationfactory.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1624,6 +1624,8 @@ CoordinateOperationFactory::Private::findOpsInRegistryDirect(
context.extent1, context.extent2);
res.insert(res.end(), resTmp.begin(), resTmp.end());
if (authName == "PROJ") {
// Do not stop at the first transformations available in
// the PROJ namespace, but allow the next authority too
continue;
}
if (!res.empty()) {
Expand Down Expand Up @@ -1669,24 +1671,34 @@ CoordinateOperationFactory::Private::findOpsInRegistryDirectTo(

const auto authorities(getCandidateAuthorities(
authFactory, targetAuthName, targetAuthName));
std::vector<CoordinateOperationNNPtr> res;
for (const auto &authority : authorities) {
const auto authName =
authority == "any" ? std::string() : authority;
const auto tmpAuthFactory = io::AuthorityFactory::create(
authFactory->databaseContext(),
authority == "any" ? std::string() : authority);
auto res = tmpAuthFactory->createFromCoordinateReferenceSystemCodes(
std::string(), std::string(), targetAuthName, targetCode,
context.context->getUsePROJAlternativeGridNames(),
authFactory->databaseContext(), authName);
auto resTmp =
tmpAuthFactory->createFromCoordinateReferenceSystemCodes(
std::string(), std::string(), targetAuthName, targetCode,
context.context->getUsePROJAlternativeGridNames(),

gridAvailabilityUse ==
CoordinateOperationContext::GridAvailabilityUse::
DISCARD_OPERATION_IF_MISSING_GRID ||
gridAvailabilityUse ==
CoordinateOperationContext::GridAvailabilityUse::
DISCARD_OPERATION_IF_MISSING_GRID ||
gridAvailabilityUse ==
CoordinateOperationContext::GridAvailabilityUse::
KNOWN_AVAILABLE,
gridAvailabilityUse ==
CoordinateOperationContext::GridAvailabilityUse::
KNOWN_AVAILABLE,
gridAvailabilityUse == CoordinateOperationContext::
GridAvailabilityUse::KNOWN_AVAILABLE,
context.context->getDiscardSuperseded(), true, true,
context.extent1, context.extent2);
context.context->getDiscardSuperseded(), true, true,
context.extent1, context.extent2);
res.insert(res.end(), resTmp.begin(), resTmp.end());
if (authName == "PROJ") {
// Do not stop at the first transformations available in
// the PROJ namespace, but allow the next authority too
continue;
}
if (!res.empty()) {
auto resFiltered =
FilterResults(res, context.context, context.extent1,
Expand Down
47 changes: 47 additions & 0 deletions test/unit/test_operationfactory.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4922,6 +4922,53 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}

// Check that we can handle vertical transformations where there is a
// mix of available ones in the PROJ namespace (mx_inegi_ggm10) and in
// in the EPSG namespace (us_noaa_g2018u0)
// This test might no longer test this scenario if mx_inegi_ggm10 is
// referenced one day by EPSG, but at least this tests a common use case.
{
auto authFactoryAll =
AuthorityFactory::create(DatabaseContext::create(), std::string());
auto ctxt =
CoordinateOperationContext::create(authFactoryAll, nullptr, 0.0);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
ctxt->setGridAvailabilityUse(
CoordinateOperationContext::GridAvailabilityUse::
IGNORE_GRID_AVAILABILITY);
// NAD83(2011) + NAVD88 height
auto srcObj = createFromUserInput(
"EPSG:6318+5703", authFactory->databaseContext(), false);
auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
ASSERT_TRUE(src != nullptr);
auto nnSrc = NN_NO_CHECK(src);

auto list = CoordinateOperationFactory::create()->createOperations(
nnSrc,
authFactory->createCoordinateReferenceSystem("4979"), // WGS 84 3D
ctxt);
bool foundGeoid2018 = false;
bool foundGGM10 = false;
for (const auto &op : list) {
try {
const auto projString = op->exportToPROJString(
PROJStringFormatter::create(
PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get());
if (projString.find("us_noaa_g2018u0.tif") != std::string::npos)
foundGeoid2018 = true;
else if (projString.find("mx_inegi_ggm10.tif") !=
std::string::npos)
foundGGM10 = true;
} catch (const std::exception &) {
}
}
EXPECT_TRUE(foundGeoid2018);
EXPECT_TRUE(foundGGM10);
}
}

// ---------------------------------------------------------------------------
Expand Down

0 comments on commit b0d8717

Please sign in to comment.