Skip to content

Commit

Permalink
Merge pull request #3365 from rouault/fix_findsOpsInRegistryWithInter…
Browse files Browse the repository at this point in the history 
…mediate

findsOpsInRegistryWithIntermediate(): make it work when source/target geodetic CRS has no known id
  • Loading branch information
@rouault @github-actions
rouault authored and github-actions[bot] committed Oct 7, 2022
1 parent 08bcb4f commit 428e815
Show file tree
Hide file tree
Showing 2 changed files with 182 additions and 51 deletions.
152 changes: 101 additions & 51 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1814,6 +1814,56 @@ CoordinateOperationFactory::Private::findOpsInRegistryDirectTo(

// ---------------------------------------------------------------------------

static std::vector<crs::GeodeticCRSNNPtr>
findCandidateGeodCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
const crs::GeodeticCRS *crs,
const datum::GeodeticReferenceFrame *datum) {
std::vector<crs::GeodeticCRSNNPtr> candidates;
assert(datum);
const auto &ids = datum->identifiers();
const auto &datumName = datum->nameStr();
if (!ids.empty()) {
for (const auto &id : ids) {
const auto &authName = *(id->codeSpace());
const auto &code = id->code();
if (!authName.empty()) {
const auto crsIds = crs->identifiers();
const auto tmpFactory =
(crsIds.size() == 1 &&
*(crsIds.front()->codeSpace()) == authName)
? io::AuthorityFactory::create(
authFactory->databaseContext(), authName)
.as_nullable()
: authFactory;
auto l_candidates = tmpFactory->createGeodeticCRSFromDatum(
authName, code, std::string());
for (const auto &candidate : l_candidates) {
candidates.emplace_back(candidate);
}
}
}
} else if (datumName != "unknown" && datumName != "unnamed") {
auto matches = authFactory->createObjectsFromName(
datumName,
{io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME}, false,
2);
if (matches.size() == 1) {
const auto &match = matches.front();
if (datum->_isEquivalentTo(
match.get(), util::IComparable::Criterion::EQUIVALENT) &&
!match->identifiers().empty()) {
return findCandidateGeodCRSForDatum(
authFactory, crs,
dynamic_cast<const datum::GeodeticReferenceFrame *>(
match.get()));
}
}
}
return candidates;
}

// ---------------------------------------------------------------------------

//! @cond Doxygen_Suppress

// Look in the authority registry for operations from sourceCRS to targetCRS
Expand All @@ -1834,9 +1884,59 @@ CoordinateOperationFactory::Private::findsOpsInRegistryWithIntermediate(
assert(authFactory);

std::list<std::pair<std::string, std::string>> sourceIds;
std::list<std::pair<std::string, std::string>> targetIds;
buildCRSIds(sourceCRS, context, sourceIds);
if (sourceIds.empty()) {
auto geodSrc = dynamic_cast<crs::GeodeticCRS *>(sourceCRS.get());
if (geodSrc) {
const auto candidatesSrcGeod(findCandidateGeodCRSForDatum(
authFactory, geodSrc,
geodSrc
->datumNonNull(authFactory->databaseContext().as_nullable())
.get()));
std::vector<CoordinateOperationNNPtr> res;
for (const auto &candidateSrcGeod : candidatesSrcGeod) {
if (candidateSrcGeod->coordinateSystem()->axisList().size() ==
geodSrc->coordinateSystem()->axisList().size() &&
((dynamic_cast<crs::GeographicCRS *>(sourceCRS.get()) !=
nullptr) ==
(dynamic_cast<crs::GeographicCRS *>(
candidateSrcGeod.get()) != nullptr))) {
const auto opsWithIntermediate =
findsOpsInRegistryWithIntermediate(
candidateSrcGeod, targetCRS, context,
useCreateBetweenGeodeticCRSWithDatumBasedIntermediates);
if (!opsWithIntermediate.empty()) {
const auto opsFirst = createOperations(
sourceCRS, candidateSrcGeod, context);
for (const auto &opFirst : opsFirst) {
for (const auto &opSecond : opsWithIntermediate) {
try {
res.emplace_back(
ConcatenatedOperation::
createComputeMetadata(
{opFirst, opSecond},
disallowEmptyIntersection));
} catch (
const InvalidOperationEmptyIntersection &) {
}
}
}
if (!res.empty())
return res;
}
}
}
}
return std::vector<CoordinateOperationNNPtr>();
}

std::list<std::pair<std::string, std::string>> targetIds;
buildCRSIds(targetCRS, context, targetIds);
if (targetIds.empty()) {
return applyInverse(findsOpsInRegistryWithIntermediate(
targetCRS, sourceCRS, context,
useCreateBetweenGeodeticCRSWithDatumBasedIntermediates));
}

const auto gridAvailabilityUse = context.context->getGridAvailabilityUse();
for (const auto &idSrc : sourceIds) {
Expand Down Expand Up @@ -2659,56 +2759,6 @@ static bool hasIdentifiers(const CoordinateOperationNNPtr &op) {

// ---------------------------------------------------------------------------

static std::vector<crs::CRSNNPtr>
findCandidateGeodCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
const crs::GeodeticCRS *crs,
const datum::GeodeticReferenceFrame *datum) {
std::vector<crs::CRSNNPtr> candidates;
assert(datum);
const auto &ids = datum->identifiers();
const auto &datumName = datum->nameStr();
if (!ids.empty()) {
for (const auto &id : ids) {
const auto &authName = *(id->codeSpace());
const auto &code = id->code();
if (!authName.empty()) {
const auto crsIds = crs->identifiers();
const auto tmpFactory =
(crsIds.size() == 1 &&
*(crsIds.front()->codeSpace()) == authName)
? io::AuthorityFactory::create(
authFactory->databaseContext(), authName)
.as_nullable()
: authFactory;
auto l_candidates = tmpFactory->createGeodeticCRSFromDatum(
authName, code, std::string());
for (const auto &candidate : l_candidates) {
candidates.emplace_back(candidate);
}
}
}
} else if (datumName != "unknown" && datumName != "unnamed") {
auto matches = authFactory->createObjectsFromName(
datumName,
{io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME}, false,
2);
if (matches.size() == 1) {
const auto &match = matches.front();
if (datum->_isEquivalentTo(
match.get(), util::IComparable::Criterion::EQUIVALENT) &&
!match->identifiers().empty()) {
return findCandidateGeodCRSForDatum(
authFactory, crs,
dynamic_cast<const datum::GeodeticReferenceFrame *>(
match.get()));
}
}
}
return candidates;
}

// ---------------------------------------------------------------------------

void CoordinateOperationFactory::Private::setCRSs(
CoordinateOperation *co, const crs::CRSNNPtr &sourceCRS,
const crs::CRSNNPtr &targetCRS) {
Expand Down
81 changes: 81 additions & 0 deletions test/unit/test_operationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1767,6 +1767,87 @@ TEST(operation, geogCRS_to_geogCRS_WGS84_to_GDA2020) {

// ---------------------------------------------------------------------------

TEST(operation, geogCRS_to_geogCRS_with_intermediate_no_ids) {

auto dbContext = DatabaseContext::create();
auto authFactory = AuthorityFactory::create(dbContext, std::string());
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
ctxt->setGridAvailabilityUse(
CoordinateOperationContext::GridAvailabilityUse::
IGNORE_GRID_AVAILABILITY);

auto objSrc = WKTParser().createFromWKT(
"GEOGCRS[\"input\",\n"
" DATUM[\"International Terrestrial Reference Frame 2014\",\n"
" ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n"
" LENGTHUNIT[\"metre\",1]],\n"
" ID[\"EPSG\",1165]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8901]],\n"
" CS[ellipsoidal,3],\n"
" AXIS[\"longitude\",east,\n"
" ORDER[1],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433,\n"
" ID[\"EPSG\",9122]]],\n"
" AXIS[\"latitude\",north,\n"
" ORDER[2],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433,\n"
" ID[\"EPSG\",9122]]],\n"
" AXIS[\"ellipsoidal height (h)\",up,\n"
" ORDER[3],\n"
" LENGTHUNIT[\"metre\",1]]]");
auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
ASSERT_TRUE(src != nullptr);

auto objDest = WKTParser().createFromWKT(
"GEOGCRS[\"output\",\n"
" DATUM[\"Estonia 1997\",\n"
" ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n"
" LENGTHUNIT[\"metre\",1]],\n"
" ID[\"EPSG\",6180]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8901]],\n"
" CS[ellipsoidal,3],\n"
" AXIS[\"longitude\",east,\n"
" ORDER[1],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433,\n"
" ID[\"EPSG\",9122]]],\n"
" AXIS[\"latitude\",north,\n"
" ORDER[2],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433,\n"
" ID[\"EPSG\",9122]]],\n"
" AXIS[\"ellipsoidal height (h)\",up,\n"
" ORDER[3],\n"
" LENGTHUNIT[\"metre\",1]]]");
auto dest = nn_dynamic_pointer_cast<CRS>(objDest);
ASSERT_TRUE(dest != nullptr);

{
auto list = CoordinateOperationFactory::create()->createOperations(
NN_NO_CHECK(src), NN_NO_CHECK(dest), ctxt);
ASSERT_GE(list.size(), 1U);
// Test that a non-noop operation is returned
EXPECT_EQ(
list[0]->nameStr(),
"axis order change (geographic3D horizontal) + "
"Conversion from ITRF2014 (geog3D) to ITRF2014 (geocentric) + "
"ITRF2014 to ETRF2014 (1) + "
"Inverse of NKG_ETRF14 to ETRF2014 + "
"NKG_ETRF14 to ETRF96@2000.0 + "
"ETRF96@2000.0 to ETRF96@1997.56 + "
"Conversion from ETRS89 (geocentric) to ETRS89 (geog2D) + "
"Inverse of EST97 to ETRS89 (1) + "
"Null geographic offset from EST97 (geog2D) to EST97 (geog3D) + "
"axis order change (geographic3D horizontal)");
}
}

// ---------------------------------------------------------------------------

static ProjectedCRSNNPtr createUTM31_WGS84() {
return ProjectedCRS::create(
PropertyMap(), GeographicCRS::EPSG_4326,
Expand Down

0 comments on commit 428e815

Please sign in to comment.