Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

createOperations(): improvement for "NAD83(CSRS) + CGVD28 height" to "NAD83(CSRS) + CGVD2013(CGG2013) height" #2976

Merged
merged 1 commit into from
Dec 9, 2021
Merged

createOperations(): improvement for "NAD83(CSRS) + CGVD28 height" to "NAD83(CSRS) + CGVD2013(CGG2013) height" #2976

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
152 changes: 120 additions & 32 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -729,10 +729,12 @@ struct PrecomputedOpCharacteristics {
struct SortFunction {

const std::map<CoordinateOperation *, PrecomputedOpCharacteristics> &map;
const std::string BALLPARK_GEOGRAPHIC_OFFSET_FROM;

explicit SortFunction(const std::map<CoordinateOperation *,
PrecomputedOpCharacteristics> &mapIn)
: map(mapIn) {}
: map(mapIn), BALLPARK_GEOGRAPHIC_OFFSET_FROM(
std::string(BALLPARK_GEOGRAPHIC_OFFSET) + " from ") {}

// Sorting function
// Return true if a < b
Expand Down Expand Up @@ -860,6 +862,48 @@ struct SortFunction {

const auto &a_name = a->nameStr();
const auto &b_name = b->nameStr();

// Make sure that
// "Ballpark geographic offset from NAD83(CSRS)v6 to NAD83(CSRS)"
// has more priority than
// "Ballpark geographic offset from ITRF2008 to NAD83(CSRS)"
const auto posA = a_name.find(BALLPARK_GEOGRAPHIC_OFFSET_FROM);
const auto posB = b_name.find(BALLPARK_GEOGRAPHIC_OFFSET_FROM);
if (posA != std::string::npos && posB != std::string::npos) {
const auto pos2A = a_name.find(" to ", posA);
const auto pos2B = b_name.find(" to ", posB);
if (pos2A != std::string::npos && pos2B != std::string::npos) {
const auto pos3A = a_name.find(" + ", pos2A);
const auto pos3B = b_name.find(" + ", pos2B);
const std::string fromA = a_name.substr(
posA + BALLPARK_GEOGRAPHIC_OFFSET_FROM.size(),
pos2A - (posA + BALLPARK_GEOGRAPHIC_OFFSET_FROM.size()));
const std::string toA =
a_name.substr(pos2A + strlen(" to "),
pos3A == std::string::npos
? pos3A
: pos3A - (pos2A + strlen(" to ")));
const std::string fromB = b_name.substr(
posB + BALLPARK_GEOGRAPHIC_OFFSET_FROM.size(),
pos2B - (posB + BALLPARK_GEOGRAPHIC_OFFSET_FROM.size()));
const std::string toB =
b_name.substr(pos2B + strlen(" to "),
pos3B == std::string::npos
? pos3B
: pos3B - (pos2B + strlen(" to ")));
const bool similarCRSInA =
(fromA.find(toA) == 0 || toA.find(fromA) == 0);
const bool similarCRSInB =
(fromB.find(toB) == 0 || toB.find(fromB) == 0);
if (similarCRSInA && !similarCRSInB) {
return true;
}
if (!similarCRSInA && similarCRSInB) {
return false;
}
}
}

// The shorter name, the better ?
if (a_name.size() < b_name.size()) {
return true;
Expand Down Expand Up @@ -1240,7 +1284,7 @@ struct FilterResults {
area, getAccuracy(op), isPROJExportable, hasGrids,
gridsAvailable, gridsKnown, stepCount,
op->hasBallparkTransformation(),
op->nameStr().find("ballpark vertical transformation") !=
op->nameStr().find(BALLPARK_VERTICAL_TRANSFORMATION) !=
std::string::npos,
isNullTransformation(op->nameStr()));
}
Expand Down Expand Up @@ -4482,7 +4526,9 @@ void CoordinateOperationFactory::Private::createOperationsVertToVert(
const double factor = convSrc / convDst;
if (!equivalentVDatum) {
auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
name += " (";
name += BALLPARK_VERTICAL_TRANSFORMATION;
name += ')';
auto conv = Transformation::createChangeVerticalUnit(
util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name),
sourceCRS, targetCRS,
Expand Down Expand Up @@ -4588,9 +4634,12 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeogBallpark(
targetCRSExtent.get(), util::IComparable::Criterion::EQUIVALENT);

util::PropertyMap map;
map.set(common::IdentifiedObject::NAME_KEY,
buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) +
BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT)
std::string transfName(
buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()));
transfName += " (";
transfName += BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT;
transfName += ')';
map.set(common::IdentifiedObject::NAME_KEY, transfName)
.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
sameExtent ? NN_NO_CHECK(sourceCRSExtent)
: metadata::Extent::WORLD);
Expand Down Expand Up @@ -5280,45 +5329,84 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
createOperations(sourceCRS, intermGeogSrc, context);
const auto opsGeogToTarget =
createOperations(intermGeogDst, targetCRS, context);

// We preferably accept using an intermediate if the operations
// to it do not include ballpark operations, but if they do include
// grids, using such operations result will be better than nothing.
// This will help for example for "NAD83(CSRS) + CGVD28 height" to
// "NAD83(CSRS) + CGVD2013(CGG2013) height" where the transformation
// between "CGVD2013(CGG2013) height" and "NAD83(CSRS)" actually
// involves going through NAD83(CSRS)v6
const auto hasKnownGrid =
[&dbContext](const CoordinateOperationNNPtr &op) {
const auto grids = op->gridsNeeded(dbContext, true);
if (grids.empty()) {
return false;
}
for (const auto &grid : grids) {
if (grid.available) {
return true;
}
}
return false;
};

const bool hasNonTrivialSrcTransf =
!opsSrcToGeog.empty() &&
!opsSrcToGeog.front()->hasBallparkTransformation();
(!opsSrcToGeog.front()->hasBallparkTransformation() ||
hasKnownGrid(opsSrcToGeog.front()));
const bool hasNonTrivialTargetTransf =
!opsGeogToTarget.empty() &&
!opsGeogToTarget.front()->hasBallparkTransformation();
(!opsGeogToTarget.front()->hasBallparkTransformation() ||
hasKnownGrid(opsGeogToTarget.front()));
if (hasNonTrivialSrcTransf && hasNonTrivialTargetTransf) {
const auto opsGeogSrcToGeogDst =
createOperations(intermGeogSrc, intermGeogDst, context);
for (const auto &op1 : opsSrcToGeog) {
if (op1->hasBallparkTransformation()) {
// std::cerr << "excluded " << op1->nameStr() << std::endl;
continue;
}
for (const auto &op2 : opsGeogSrcToGeogDst) {
for (const auto &op3 : opsGeogToTarget) {
if (op3->hasBallparkTransformation()) {
// std::cerr << "excluded " << op3->nameStr() <<
// std::endl;
continue;
}
try {
res.emplace_back(
ConcatenatedOperation::createComputeMetadata(
intermGeogSrcIsSameAsIntermGeogDst
? std::vector<
CoordinateOperationNNPtr>{op1,
op3}
: std::vector<
CoordinateOperationNNPtr>{op1,
op2,
op3},
disallowEmptyIntersection));
} catch (const std::exception &) {
// In first pass, exclude (horizontal) ballpark operations, but
// accept them in second pass.
for (int pass = 0; pass <= 1 && res.empty(); pass++) {
for (const auto &op1 : opsSrcToGeog) {
if (pass == 0 && op1->hasBallparkTransformation()) {
// std::cerr << "excluded " << op1->nameStr() <<
// std::endl;
continue;
}
if (op1->nameStr().find(BALLPARK_VERTICAL_TRANSFORMATION) !=
std::string::npos) {
continue;
}
for (const auto &op2 : opsGeogSrcToGeogDst) {
for (const auto &op3 : opsGeogToTarget) {
if (pass == 0 && op3->hasBallparkTransformation()) {
// std::cerr << "excluded " << op3->nameStr() <<
// std::endl;
continue;
}
if (op3->nameStr().find(
BALLPARK_VERTICAL_TRANSFORMATION) !=
std::string::npos) {
continue;
}
try {
res.emplace_back(
ConcatenatedOperation::createComputeMetadata(
intermGeogSrcIsSameAsIntermGeogDst
? std::vector<
CoordinateOperationNNPtr>{op1,
op3}
: std::vector<
CoordinateOperationNNPtr>{op1,
op2,
op3},
disallowEmptyIntersection));
} catch (const std::exception &) {
}
}
}
}
}
}

if (!res.empty()) {
return;
}
Expand Down
6 changes: 3 additions & 3 deletions src/iso19111/operation/oputils.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,10 +61,10 @@ const char *NULL_GEOGRAPHIC_OFFSET = "Null geographic offset";
const char *NULL_GEOCENTRIC_TRANSLATION = "Null geocentric translation";
const char *BALLPARK_GEOGRAPHIC_OFFSET = "Ballpark geographic offset";
const char *BALLPARK_VERTICAL_TRANSFORMATION =
" (ballpark vertical transformation)";
"ballpark vertical transformation";
const char *BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT =
" (ballpark vertical transformation, without ellipsoid height to vertical "
"height correction)";
"ballpark vertical transformation, without ellipsoid height to vertical "
"height correction";

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

Expand Down
57 changes: 57 additions & 0 deletions test/unit/test_operationfactory.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4518,6 +4518,63 @@ TEST(operation, compoundCRS_to_compoundCRS_issue_2720) {

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

TEST(
operation,
compoundCRS_to_compoundCRS_concatenated_operation_with_two_vert_transformation_and_ballpark_geog) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
// "NAD83(CSRS) + CGVD28 height"
auto srcObj = createFromUserInput("EPSG:4617+5713",
authFactory->databaseContext(), false);
auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
ASSERT_TRUE(src != nullptr);

// "NAD83(CSRS) + CGVD2013(CGG2013) height"
auto dstObj = createFromUserInput("EPSG:4617+6647",
authFactory->databaseContext(), false);
auto dst = nn_dynamic_pointer_cast<CRS>(dstObj);
ASSERT_TRUE(dst != nullptr);

// That transformation involves doing CGVD28 height to CGVD2013(CGG2013)
// height by doing:
// - CGVD28 height to NAD83(CSRS): EPSG registered operation
// - NAD83(CSRS) to CGVD2013(CGG2013) height by doing:
// * NAD83(CSRS) to NAD83(CSRS)v6: ballpark
// * NAD83(CSRS)v6 to CGVD2013(CGG2013): EPSG registered operation
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
{
auto list = CoordinateOperationFactory::create()->createOperations(
NN_NO_CHECK(src), NN_NO_CHECK(dst), ctxt);
ASSERT_GE(list.size(), 1U);
// Check that we have the transformation using NAD83(CSRS)v6 first
// (as well as the one between NAD83(CSRS) to CGVD28 height)
EXPECT_EQ(list[0]->nameStr(),
"Inverse of NAD83(CSRS) to CGVD28 height (1) + "
"Inverse of Ballpark geographic offset from NAD83(CSRS)v6 to "
"NAD83(CSRS) + "
"NAD83(CSRS)v6 to CGVD2013(CGG2013) height (1) + "
"Inverse of Ballpark geographic offset from NAD83(CSRS) to "
"NAD83(CSRS)v6");
}
{
auto list = CoordinateOperationFactory::create()->createOperations(
NN_NO_CHECK(dst), NN_NO_CHECK(src), ctxt);
ASSERT_GE(list.size(), 1U);
// Check that we have the transformation using NAD83(CSRS)v6 first
// (as well as the one between NAD83(CSRS) to CGVD28 height)
EXPECT_EQ(
list[0]->nameStr(),
"Ballpark geographic offset from NAD83(CSRS) to NAD83(CSRS)v6 + "
"Inverse of NAD83(CSRS)v6 to CGVD2013(CGG2013) height (1) + "
"Ballpark geographic offset from NAD83(CSRS)v6 to NAD83(CSRS) + "
"NAD83(CSRS) to CGVD28 height (1)");
}
}

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

TEST(operation, vertCRS_to_vertCRS) {

auto vertcrs_m_obj = PROJStringParser().createFromPROJString("+vunits=m");
Expand Down