Skip to content

Commit

Permalink
Merge pull request #1731 from rouault/fix_createoperations_with_geoid…
Browse files Browse the repository at this point in the history 
…grids_and_non_metre_vunits

createOperations(): fix transformation computation from/to a CRS with +geoidgrids and +vunits != m
  • Loading branch information
@rouault
rouault committed Nov 14, 2019
2 parents c42801b + 2cfd4d6 commit 448af96
Show file tree
Hide file tree
Showing 3 changed files with 93 additions and 10 deletions.
51 changes: 42 additions & 9 deletions src/iso19111/coordinateoperation.cpp
View file
Expand Up @@ -13100,10 +13100,16 @@ std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
// which is likely/hopefully the only one.
for (const auto &opFirst : resTmp) {
if (hasIdentifiers(opFirst)) {
res.emplace_back(
ConcatenatedOperation::createComputeMetadata(
{opFirst, opsSecond.front()},
!allowEmptyIntersection));
if (candidateVert->_isEquivalentTo(
targetCRS.get(),
util::IComparable::Criterion::EQUIVALENT)) {
res.emplace_back(opFirst);
} else {
res.emplace_back(
ConcatenatedOperation::createComputeMetadata(
{opFirst, opsSecond.front()},
!allowEmptyIntersection));
}
}
}
if (!res.empty())
Expand Down Expand Up @@ -13456,6 +13462,24 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog(
hubSrcGeog->_isEquivalentTo(geogDst,
util::IComparable::Criterion::EQUIVALENT) &&
dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get())) {
auto transfSrc = boundSrc->transformation()->sourceCRS();
if (dynamic_cast<const crs::VerticalCRS *>(transfSrc.get()) &&
!boundSrc->baseCRS()->_isEquivalentTo(
transfSrc.get(), util::IComparable::Criterion::EQUIVALENT)) {
auto opsFirst =
createOperations(boundSrc->baseCRS(), transfSrc, context);
for (const auto &opFirst : opsFirst) {
try {
res.emplace_back(
ConcatenatedOperation::createComputeMetadata(
{opFirst, boundSrc->transformation()},
!allowEmptyIntersection));
} catch (const InvalidOperationEmptyIntersection &) {
}
}
return;
}

res.emplace_back(boundSrc->transformation());
return;
}
Expand Down Expand Up @@ -13590,15 +13614,15 @@ void CoordinateOperationFactory::Private::createOperationsVertToVert(
common::Scale(heightDepthReversal ? -factor : factor), {});
conv->setHasBallparkTransformation(true);
res.push_back(conv);
} else if (convSrc != convDst) {
} else if (convSrc != convDst || !heightDepthReversal) {
auto conv = Conversion::createChangeVerticalUnit(
util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name),
// In case of a height depth reversal, we should probably have
// 2 steps instead of putting a negative factor...
common::Scale(heightDepthReversal ? -factor : factor));
conv->setCRSs(sourceCRS, targetCRS, nullptr);
res.push_back(conv);
} else if (heightDepthReversal) {
} else {
auto conv = Conversion::createHeightDepthReversal(
util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name));
conv->setCRSs(sourceCRS, targetCRS, nullptr);
Expand Down Expand Up @@ -14131,10 +14155,13 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
if (componentsSrc[0]->extractGeographicCRS() &&
componentsDst[0]->extractGeographicCRS()) {

bool verticalTransfIsNoOp = false;
std::vector<CoordinateOperationNNPtr> verticalTransforms;
if (componentsSrc.size() >= 2 &&
componentsSrc[1]->extractVerticalCRS() &&
componentsDst[1]->extractVerticalCRS()) {
verticalTransfIsNoOp =
componentsSrc[1]->_isEquivalentTo(componentsDst[1].get());
verticalTransforms = createOperations(
componentsSrc[1], componentsDst[1], context);
}
Expand Down Expand Up @@ -14181,9 +14208,15 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
for (const auto &opDst : opGeogCRStoDstCRS) {

try {
auto op = createHorizVerticalHorizPROJBased(
sourceCRS, targetCRS, opSrc, verticalTransform,
opDst, interpolationGeogCRS, true);
auto op = verticalTransfIsNoOp
? ConcatenatedOperation::
createComputeMetadata(
{opSrc, opDst},
!allowEmptyIntersection)
: createHorizVerticalHorizPROJBased(
sourceCRS, targetCRS, opSrc,
verticalTransform, opDst,
interpolationGeogCRS, true);
res.emplace_back(op);
} catch (const InvalidOperationEmptyIntersection &) {
}
Expand Down
5 changes: 4 additions & 1 deletion src/iso19111/io.cpp
View file
Expand Up @@ -8364,7 +8364,10 @@ PROJStringParser::Private::buildBoundOrCompoundCRSIfNeeded(int iStep,
Transformation::createGravityRelatedHeightToGeographic3D(
PropertyMap().set(IdentifiedObject::NAME_KEY,
"unknown to WGS84 ellipsoidal height"),
crs, GeographicCRS::EPSG_4979, nullptr, geoidgrids,
VerticalCRS::create(createMapWithUnknownName(), vdatum,
VerticalCS::createGravityRelatedHeight(
common::UnitOfMeasure::METRE)),
GeographicCRS::EPSG_4979, nullptr, geoidgrids,
std::vector<PositionalAccuracyNNPtr>());
auto boundvcrs =
BoundCRS::create(vcrs, GeographicCRS::EPSG_4979, transformation);
Expand Down
47 changes: 47 additions & 0 deletions test/unit/test_operation.cpp
View file
Expand Up @@ -6638,6 +6638,53 @@ TEST(operation, compoundCRS_with_boundProjCRS_and_boundVerticalCRS_to_geogCRS) {

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

TEST(operation,
compoundCRS_with_boundVerticalCRS_from_geoidgrids_with_m_to_geogCRS) {

auto objSrc = PROJStringParser().createFromPROJString(
"+proj=longlat +datum=WGS84 +geoidgrids=@foo.gtx +type=crs");
auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
ASSERT_TRUE(src != nullptr);
auto op = CoordinateOperationFactory::create()->createOperation(
NN_NO_CHECK(src), GeographicCRS::EPSG_4979);
ASSERT_TRUE(op != nullptr);
EXPECT_EQ(op->nameStr(), "axis order change (2D) + "
"unknown to WGS84 ellipsoidal height");
EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
"+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}

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

TEST(operation,
compoundCRS_with_boundVerticalCRS_from_geoidgrids_with_ftus_to_geogCRS) {

auto objSrc = PROJStringParser().createFromPROJString(
"+proj=longlat +datum=WGS84 +geoidgrids=@foo.gtx +vunits=us-ft "
"+type=crs");
auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
ASSERT_TRUE(src != nullptr);
auto op = CoordinateOperationFactory::create()->createOperation(
NN_NO_CHECK(src), GeographicCRS::EPSG_4979);
ASSERT_TRUE(op != nullptr);
EXPECT_EQ(op->nameStr(), "axis order change (2D) + "
"Transformation from unknown to unknown + "
"unknown to WGS84 ellipsoidal height");
EXPECT_EQ(
op->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=rad +z_out=m "
"+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}

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

TEST(operation, geocent_to_compoundCRS) {
auto objSrc = PROJStringParser().createFromPROJString(
"+proj=geocent +datum=WGS84 +units=m +type=crs");
Expand Down

0 comments on commit 448af96

Please sign in to comment.