diff options
| -rw-r--r-- | src/iso19111/coordinateoperation.cpp | 51 | ||||
| -rw-r--r-- | src/iso19111/io.cpp | 5 | ||||
| -rw-r--r-- | test/unit/test_operation.cpp | 47 |
3 files changed, 93 insertions, 10 deletions
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index c51b4d3d..287611f8 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -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()) @@ -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; } @@ -13590,7 +13614,7 @@ 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 @@ -13598,7 +13622,7 @@ void CoordinateOperationFactory::Private::createOperationsVertToVert( 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); @@ -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); } @@ -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 &) { } diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 645bec0b..a0a87f65 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -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); diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 07b0daea..ce01f967 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -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"); |
