From 4ee2ed65f728985e53ba0c76590f07b253b03587 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 19 May 2020 01:31:25 +0200 Subject: createOperations(): fix 'caching' bugs causing potential exception about Inconsistent chainging of CRS (fixes #2232) --- src/iso19111/coordinateoperation.cpp | 4 ++ test/unit/test_operation.cpp | 79 ++++++++++++++++++++++++++++++++++++ 2 files changed, 83 insertions(+) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 9815a22a..d1bc8eb5 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -6587,6 +6587,10 @@ TransformationNNPtr Transformation::shallowClone() const { auto transf = Transformation::nn_make_shared(*this); transf->assignSelf(transf); transf->setCRSs(this, false); + if (transf->d->forwardOperation_) { + transf->d->forwardOperation_ = + transf->d->forwardOperation_->shallowClone().as_nullable(); + } return transf; } diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 92b5f923..077ee1c0 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -7566,6 +7566,85 @@ TEST( // --------------------------------------------------------------------------- +TEST(operation, compoundCRS_to_compoundCRS_issue_2232) { + auto objSrc = WKTParser().createFromWKT( + "COMPD_CS[\"NAD83 / Alabama West + NAVD88 height - Geoid12B " + "(Meters)\",\n" + " PROJCS[\"NAD83 / Alabama West\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " PROJECTION[\"Transverse_Mercator\"],\n" + " PARAMETER[\"latitude_of_origin\",30],\n" + " PARAMETER[\"central_meridian\",-87.5],\n" + " PARAMETER[\"scale_factor\",0.999933333],\n" + " PARAMETER[\"false_easting\",600000],\n" + " PARAMETER[\"false_northing\",0],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"X\",EAST],\n" + " AXIS[\"Y\",NORTH],\n" + " AUTHORITY[\"EPSG\",\"26930\"]],\n" + " VERT_CS[\"NAVD88 height - Geoid12B (Meters)\",\n" + " VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n" + " EXTENSION[\"PROJ4_GRIDS\",\"foo.gtx\"],\n" + " AUTHORITY[\"EPSG\",\"5103\"]],\n" + " UNIT[\"metre\",1.0,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Gravity-related height\",UP],\n" + " AUTHORITY[\"EPSG\",\"5703\"]]]"); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + + auto objDst = WKTParser().createFromWKT( + "COMPD_CS[\"NAD83 + some CRS (US Feet)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"some CRS (US Feet)\",\n" + " VERT_DATUM[\"some datum\",2005],\n" + " UNIT[\"US survey foot\",0.3048006096012192,\n" + " AUTHORITY[\"EPSG\",\"9003\"]],\n" + " AXIS[\"Up\",UP]]]"); + auto dst = nn_dynamic_pointer_cast(objDst); + ASSERT_TRUE(dst != nullptr); + + auto dbContext = DatabaseContext::create(); + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + + auto list = CoordinateOperationFactory::create()->createOperations( + NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst), ctxt); + EXPECT_GE(list.size(), 1U); + + auto list2 = CoordinateOperationFactory::create()->createOperations( + NN_CHECK_ASSERT(dst), NN_CHECK_ASSERT(src), ctxt); + EXPECT_EQ(list2.size(), list.size()); +} + +// --------------------------------------------------------------------------- + TEST(operation, compoundCRS_to_compoundCRS_context) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); -- cgit v1.2.3 From 166a1c681bba659095b1e1296cbb66e9f1637bbd Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 19 May 2020 01:33:28 +0200 Subject: createOperationsWithDatumPivot(): remove useless code --- src/iso19111/coordinateoperation.cpp | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index d1bc8eb5..4bf5c41e 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -12946,7 +12946,6 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot( (isNullFirst || isNullThird || sourceAndTargetAre3D) ? opSecond->shallowClone() : opSecond); - CoordinateOperation *invCOForward = nullptr; if (isNullFirst || isNullThird) { if (opSecondCloned->identifiers().size() == 1 && (*opSecondCloned->identifiers()[0]->codeSpace()) @@ -12960,7 +12959,7 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot( auto invCO = dynamic_cast( opSecondCloned.get()); if (invCO) { - invCOForward = invCO->forwardOperation().get(); + auto invCOForward = invCO->forwardOperation().get(); if (invCOForward->identifiers().size() == 1 && (*invCOForward->identifiers()[0]->codeSpace()) .find("DERIVED_FROM") == @@ -12978,25 +12977,19 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot( auto invCO = dynamic_cast( opSecondCloned.get()); if (invCO) { - invCOForward = invCO->forwardOperation().get(); + auto invCOForward = invCO->forwardOperation().get(); invCOForward->getPrivate()->use3DHelmert_ = true; } } if (isNullFirst) { auto oldTarget(NN_CHECK_ASSERT(opSecondCloned->targetCRS())); setCRSs(opSecondCloned.get(), sourceCRS, oldTarget); - if (invCOForward) { - setCRSs(invCOForward, oldTarget, sourceCRS); - } } else { subOps.emplace_back(opFirst); } if (isNullThird) { auto oldSource(NN_CHECK_ASSERT(opSecondCloned->sourceCRS())); setCRSs(opSecondCloned.get(), oldSource, targetCRS); - if (invCOForward) { - setCRSs(invCOForward, targetCRS, oldSource); - } subOps.emplace_back(opSecondCloned); } else { subOps.emplace_back(opSecondCloned); -- cgit v1.2.3 From 1a715234754146ebe224fb849a87ca6575fdc88f Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 19 May 2020 01:46:28 +0200 Subject: createOperations(): speed optimizations for transforming between a BoundCRS of a datum and the same datum (relates to #2232) --- src/iso19111/coordinateoperation.cpp | 22 ++++++++++++++++++---- test/unit/test_operation.cpp | 22 +++++++++++++++++++++- 2 files changed, 39 insertions(+), 5 deletions(-) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 4bf5c41e..231d31a0 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -13997,6 +13997,20 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( const auto &hubSrc = boundSrc->hubCRS(); auto hubSrcGeog = dynamic_cast(hubSrc.get()); auto geogCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractGeographicCRS(); + auto geogDstDatum = geogDst->datum(); + + // If the underlying datum of the source is the same as the target, do + // not consider the boundCRS at all, but just its base + if (geogCRSOfBaseOfBoundSrc && geogDstDatum) { + auto geogCRSOfBaseOfBoundSrcDatum = geogCRSOfBaseOfBoundSrc->datum(); + if (geogCRSOfBaseOfBoundSrcDatum && + geogCRSOfBaseOfBoundSrcDatum->_isEquivalentTo( + geogDstDatum.get(), util::IComparable::Criterion::EQUIVALENT)) { + res = createOperations(boundSrc->baseCRS(), targetCRS, context); + return; + } + } + bool triedBoundCrsToGeogCRSSameAsHubCRS = false; // Is it: boundCRS to a geogCRS that is the same as the hubCRS ? if (hubSrcGeog && geogCRSOfBaseOfBoundSrc && @@ -14043,9 +14057,9 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( } // If the datum are equivalent, this is also fine } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog && hubSrcGeog->datum() && - geogDst->datum() && + geogDstDatum && hubSrcGeog->datum()->_isEquivalentTo( - geogDst->datum().get(), + geogDstDatum.get(), util::IComparable::Criterion::EQUIVALENT)) { auto opsFirst = createOperations( boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context); @@ -14088,14 +14102,14 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( // +nadgrids=ntv1_can.dat,conus" // to "+proj=latlong +datum=NAD83" } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog && hubSrcGeog->datum() && - geogDst->datum() && + geogDstDatum && geogCRSOfBaseOfBoundSrc->ellipsoid()->_isEquivalentTo( datum::Ellipsoid::CLARKE_1866.get(), util::IComparable::Criterion::EQUIVALENT) && hubSrcGeog->datum()->_isEquivalentTo( datum::GeodeticReferenceFrame::EPSG_6326.get(), util::IComparable::Criterion::EQUIVALENT) && - geogDst->datum()->_isEquivalentTo( + geogDstDatum->_isEquivalentTo( datum::GeodeticReferenceFrame::EPSG_6269.get(), util::IComparable::Criterion::EQUIVALENT)) { auto nnGeogCRSOfBaseOfBoundSrc = NN_NO_CHECK(geogCRSOfBaseOfBoundSrc); diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 077ee1c0..bdf648e4 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -6379,6 +6379,26 @@ TEST(operation, geogCRS_to_boundCRS_of_geogCRS) { // --------------------------------------------------------------------------- +TEST(operation, boundCRS_to_geogCRS_same_datum_context) { + auto boundCRS = BoundCRS::createFromTOWGS84( + GeographicCRS::EPSG_4269, std::vector{1, 2, 3, 4, 5, 6, 7}); + auto dbContext = DatabaseContext::create(); + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + boundCRS, GeographicCRS::EPSG_4269, ctxt); + ASSERT_EQ(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=noop"); +} + +// --------------------------------------------------------------------------- + TEST(operation, boundCRS_to_boundCRS) { auto utm31 = ProjectedCRS::create( PropertyMap(), GeographicCRS::EPSG_4807, @@ -7545,7 +7565,7 @@ TEST( CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); auto list = CoordinateOperationFactory::create()->createOperations( NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst), ctxt); - ASSERT_GE(list.size(), 1U); + ASSERT_EQ(list.size(), 1U); EXPECT_EQ(list[0]->nameStr(), "Inverse of unnamed + " "Transformation from NAD83 to WGS84 + " -- cgit v1.2.3 From ec3fdd00f133736560f807765dd73367c85f4bdc Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 19 May 2020 12:42:50 +0200 Subject: WKT1 ingestion: fix ingestion of COMPD_CS with ellipsoidal vertical datum and non metre units (contributes to fixes #2232) --- include/proj/crs.hpp | 19 ++++++---- src/iso19111/crs.cpp | 100 ++++++++++++++++++++++++++++---------------------- test/unit/test_io.cpp | 38 +++++++++++++++++++ 3 files changed, 107 insertions(+), 50 deletions(-) diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp index 7aa74c41..bbdc9565 100644 --- a/include/proj/crs.hpp +++ b/include/proj/crs.hpp @@ -68,6 +68,12 @@ using BoundCRSPtr = std::shared_ptr; /** Non-null shared pointer of BoundCRS */ using BoundCRSNNPtr = util::nn; +class CompoundCRS; +/** Shared pointer of CompoundCRS */ +using CompoundCRSPtr = std::shared_ptr; +/** Non-null shared pointer of CompoundCRS */ +using CompoundCRSNNPtr = util::nn; + // --------------------------------------------------------------------------- class CRS; @@ -141,7 +147,12 @@ class PROJ_GCC_DLL CRS : public common::ObjectUsage, PROJ_INTERNAL CRSNNPtr allowNonConformantWKT1Export() const; PROJ_INTERNAL CRSNNPtr - attachOriginalVertCRS(const VerticalCRSNNPtr &vertCRS) const; + attachOriginalCompoundCRS(const CompoundCRSNNPtr &compoundCRS) const; + + PROJ_INTERNAL CRSNNPtr promoteTo3D( + const std::string &newName, const io::DatabaseContextPtr &dbContext, + const cs::CoordinateSystemAxisNNPtr &verticalAxisIfNotAlreadyPresent) + const; //! @endcond @@ -855,12 +866,6 @@ class PROJ_GCC_DLL InvalidCompoundCRSException : public util::Exception { // --------------------------------------------------------------------------- -class CompoundCRS; -/** Shared pointer of CompoundCRS */ -using CompoundCRSPtr = std::shared_ptr; -/** Non-null shared pointer of CompoundCRS */ -using CompoundCRSNNPtr = util::nn; - /** \brief A coordinate reference system describing the position of points * through two or more independent single coordinate reference systems. * diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index bad7deea..35e11d84 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -98,7 +98,7 @@ struct CRS::Private { bool allowNonConformantWKT1Export_ = false; // for what was initially a COMPD_CS with a VERT_CS with a datum type == // ellipsoidal height / 2002 - VerticalCRSPtr originalVertCRS_{}; + CompoundCRSPtr originalCompoundCRS_{}; void setImplicitCS(const util::PropertyMap &properties) { const auto pVal = properties.get("IMPLICIT_CS"); @@ -583,17 +583,18 @@ CRSNNPtr CRS::allowNonConformantWKT1Export() const { //! @cond Doxygen_Suppress -CRSNNPtr CRS::attachOriginalVertCRS(const VerticalCRSNNPtr &vertCRS) const { +CRSNNPtr +CRS::attachOriginalCompoundCRS(const CompoundCRSNNPtr &compoundCRS) const { const auto boundCRS = dynamic_cast(this); if (boundCRS) { return BoundCRS::create( - boundCRS->baseCRS()->attachOriginalVertCRS(vertCRS), + boundCRS->baseCRS()->attachOriginalCompoundCRS(compoundCRS), boundCRS->hubCRS(), boundCRS->transformation()); } auto crs(shallowClone()); - crs->d->originalVertCRS_ = vertCRS.as_nullable(); + crs->d->originalCompoundCRS_ = compoundCRS.as_nullable(); return crs; } @@ -868,7 +869,22 @@ CRS::getNonDeprecated(const io::DatabaseContextNNPtr &dbContext) const { */ CRSNNPtr CRS::promoteTo3D(const std::string &newName, const io::DatabaseContextPtr &dbContext) const { + auto upAxis = cs::CoordinateSystemAxis::create( + util::PropertyMap().set(IdentifiedObject::NAME_KEY, + cs::AxisName::Ellipsoidal_height), + cs::AxisAbbreviation::h, cs::AxisDirection::UP, + common::UnitOfMeasure::METRE); + return promoteTo3D(newName, dbContext, upAxis); +} + +// --------------------------------------------------------------------------- + +//! @cond Doxygen_Suppress +CRSNNPtr CRS::promoteTo3D(const std::string &newName, + const io::DatabaseContextPtr &dbContext, + const cs::CoordinateSystemAxisNNPtr + &verticalAxisIfNotAlreadyPresent) const { const auto geogCRS = dynamic_cast(this); if (geogCRS) { const auto &axisList = geogCRS->coordinateSystem()->axisList(); @@ -886,21 +902,23 @@ CRSNNPtr CRS::promoteTo3D(const std::string &newName, false); if (!res.empty()) { const auto &firstRes = res.front(); - if (geogCRS->is2DPartOf3D(NN_NO_CHECK( - dynamic_cast(firstRes.get())))) { + const auto firstResGeog = + dynamic_cast(firstRes.get()); + const auto &firstResAxisList = + firstResGeog->coordinateSystem()->axisList(); + if (firstResAxisList[2]->_isEquivalentTo( + verticalAxisIfNotAlreadyPresent.get(), + util::IComparable::Criterion::EQUIVALENT) && + geogCRS->is2DPartOf3D(NN_NO_CHECK(firstResGeog))) { return NN_NO_CHECK( util::nn_dynamic_pointer_cast(firstRes)); } } } - auto upAxis = cs::CoordinateSystemAxis::create( - util::PropertyMap().set(IdentifiedObject::NAME_KEY, - cs::AxisName::Ellipsoidal_height), - cs::AxisAbbreviation::h, cs::AxisDirection::UP, - common::UnitOfMeasure::METRE); auto cs = cs::EllipsoidalCS::create( - util::PropertyMap(), axisList[0], axisList[1], upAxis); + util::PropertyMap(), axisList[0], axisList[1], + verticalAxisIfNotAlreadyPresent); return util::nn_static_pointer_cast(GeographicCRS::create( util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, !newName.empty() ? newName : nameStr()), @@ -914,13 +932,9 @@ CRSNNPtr CRS::promoteTo3D(const std::string &newName, if (axisList.size() == 2) { auto base3DCRS = projCRS->baseCRS()->promoteTo3D(std::string(), dbContext); - auto upAxis = cs::CoordinateSystemAxis::create( - util::PropertyMap().set(IdentifiedObject::NAME_KEY, - cs::AxisName::Ellipsoidal_height), - cs::AxisAbbreviation::h, cs::AxisDirection::UP, - common::UnitOfMeasure::METRE); auto cs = cs::CartesianCS::create(util::PropertyMap(), axisList[0], - axisList[1], upAxis); + axisList[1], + verticalAxisIfNotAlreadyPresent); return util::nn_static_pointer_cast(ProjectedCRS::create( util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, !newName.empty() ? newName : nameStr()), @@ -932,7 +946,8 @@ CRSNNPtr CRS::promoteTo3D(const std::string &newName, const auto boundCRS = dynamic_cast(this); if (boundCRS) { - auto base3DCRS = boundCRS->baseCRS()->promoteTo3D(newName, dbContext); + auto base3DCRS = boundCRS->baseCRS()->promoteTo3D( + newName, dbContext, verticalAxisIfNotAlreadyPresent); auto transf = boundCRS->transformation(); try { transf->getTOWGS84Parameters(); @@ -948,6 +963,9 @@ CRSNNPtr CRS::promoteTo3D(const std::string &newName, return NN_NO_CHECK( std::static_pointer_cast(shared_from_this().as_nullable())); } + +//! @endcond + // --------------------------------------------------------------------------- /** \brief Return a variant of this CRS "demoted" to a 2D one, if not already @@ -1436,14 +1454,9 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const { return; } - auto &originalVertCRS = CRS::getPrivate()->originalVertCRS_; - if (originalVertCRS) { - formatter->startNode(io::WKTConstants::COMPD_CS, false); - formatter->addQuotedString(l_name + " + " + - originalVertCRS->nameStr()); - geogCRS2D->_exportToWKT(formatter); - originalVertCRS->_exportToWKT(formatter); - formatter->endNode(); + auto &originalCompoundCRS = CRS::getPrivate()->originalCompoundCRS_; + if (originalCompoundCRS) { + originalCompoundCRS->_exportToWKT(formatter); return; } @@ -3292,14 +3305,9 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const { return; } - auto &originalVertCRS = CRS::getPrivate()->originalVertCRS_; - if (!formatter->useESRIDialect() && originalVertCRS) { - formatter->startNode(io::WKTConstants::COMPD_CS, false); - formatter->addQuotedString(l_name + " + " + - originalVertCRS->nameStr()); - projCRS2D->_exportToWKT(formatter); - originalVertCRS->_exportToWKT(formatter); - formatter->endNode(); + auto &originalCompoundCRS = CRS::getPrivate()->originalCompoundCRS_; + if (!formatter->useESRIDialect() && originalCompoundCRS) { + originalCompoundCRS->_exportToWKT(formatter); return; } @@ -4254,8 +4262,8 @@ CRSNNPtr CompoundCRS::createLax(const util::PropertyMap &properties, auto comp1 = components[1].get(); auto comp0Geog = dynamic_cast(comp0); auto comp0Proj = dynamic_cast(comp0); + auto comp0Bound = dynamic_cast(comp0); if (comp0Geog == nullptr && comp0Proj == nullptr) { - auto comp0Bound = dynamic_cast(comp0); if (comp0Bound) { const auto *baseCRS = comp0Bound->baseCRS().get(); comp0Geog = dynamic_cast(baseCRS); @@ -4286,14 +4294,20 @@ CRSNNPtr CompoundCRS::createLax(const util::PropertyMap &properties, if (comp1Vert != nullptr && comp1Vert->datum() && comp1Vert->datum()->getWKT1DatumType() == "2002") { const auto &axis = comp1Vert->coordinateSystem()->axisList()[0]; - if (axis->unit()._isEquivalentTo( - common::UnitOfMeasure::METRE, - util::IComparable::Criterion::EQUIVALENT) && - &(axis->direction()) == &(cs::AxisDirection::UP)) { - return components[0] - ->promoteTo3D(std::string(), dbContext) - ->attachOriginalVertCRS(NN_NO_CHECK(comp1Vert)); + std::string name(components[0]->nameStr()); + if (!(axis->unit()._isEquivalentTo( + common::UnitOfMeasure::METRE, + util::IComparable::Criterion::EQUIVALENT) && + &(axis->direction()) == &(cs::AxisDirection::UP))) { + name += " (" + comp1Vert->nameStr() + ')'; } + return components[0] + ->promoteTo3D(name, dbContext, axis) + ->attachOriginalCompoundCRS(create( + properties, + comp0Bound ? std::vector{comp0Bound->baseCRS(), + components[1]} + : components)); } } diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index d62e5c59..17fddaaf 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -2614,6 +2614,44 @@ TEST(wkt_parse, // --------------------------------------------------------------------------- +TEST(wkt_parse, + COMPD_CS_horizontal_geog_plus_vertical_ellipsoidal_height_non_metre) { + // See https://github.com/OSGeo/PROJ/issues/2232 + const char *wkt = + "COMPD_CS[\"NAD83 + Ellipsoid (US Feet)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"Ellipsoid (US Feet)\",\n" + " VERT_DATUM[\"Ellipsoid\",2002],\n" + " UNIT[\"US survey foot\",0.304800609601219,\n" + " AUTHORITY[\"EPSG\",\"9003\"]],\n" + " AXIS[\"Up\",UP]]]"; + auto dbContext = DatabaseContext::create(); + auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ(crs->nameStr(), "NAD83 (Ellipsoid (US Feet))"); + EXPECT_EQ(crs->coordinateSystem()->axisList().size(), 3U); + EXPECT_NEAR(crs->coordinateSystem()->axisList()[2]->unit().conversionToSI(), + 0.304800609601219, 1e-15); + + EXPECT_EQ( + crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL, dbContext) + .get()), + wkt); +} + +// --------------------------------------------------------------------------- + TEST(wkt_parse, COORDINATEOPERATION) { std::string src_wkt; -- cgit v1.2.3 From 8dc8e1de06dcb71da11facaaff6af9e46294deb0 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 19 May 2020 12:57:22 +0200 Subject: createOperations(): fix wrong optimization in some instances of transforming a BoundCRS to a GeographicCRS (contributes to fixes #2232) --- src/iso19111/coordinateoperation.cpp | 39 +++++++++++++++++++++++------------ test/unit/test_operation.cpp | 40 ++++++++++++++++++++++++++++++++++++ 2 files changed, 66 insertions(+), 13 deletions(-) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 231d31a0..32a4b9c7 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -14018,25 +14018,38 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( geogDst, util::IComparable::Criterion::EQUIVALENT) || hubSrcGeog->is2DPartOf3D(NN_NO_CHECK(geogDst)))) { triedBoundCrsToGeogCRSSameAsHubCRS = true; + + CoordinateOperationPtr opIntermediate; + if (!geogCRSOfBaseOfBoundSrc->_isEquivalentTo( + boundSrc->transformation()->sourceCRS().get(), + util::IComparable::Criterion::EQUIVALENT)) { + auto opsIntermediate = createOperations( + NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), + boundSrc->transformation()->sourceCRS(), context); + assert(!opsIntermediate.empty()); + opIntermediate = opsIntermediate.front(); + } + if (boundSrc->baseCRS() == geogCRSOfBaseOfBoundSrc) { - // Optimization to avoid creating a useless concatenated - // operation - res.emplace_back(boundSrc->transformation()); + if (opIntermediate) { + try { + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + {NN_NO_CHECK(opIntermediate), + boundSrc->transformation()}, + !allowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { + } + } else { + // Optimization to avoid creating a useless concatenated + // operation + res.emplace_back(boundSrc->transformation()); + } return; } auto opsFirst = createOperations( boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context); if (!opsFirst.empty()) { - CoordinateOperationPtr opIntermediate; - if (!geogCRSOfBaseOfBoundSrc->_isEquivalentTo( - boundSrc->transformation()->sourceCRS().get(), - util::IComparable::Criterion::EQUIVALENT)) { - auto opsIntermediate = createOperations( - NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), - boundSrc->transformation()->sourceCRS(), context); - assert(!opsIntermediate.empty()); - opIntermediate = opsIntermediate.front(); - } for (const auto &opFirst : opsFirst) { try { std::vector subops; diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index bdf648e4..0cd3b3f5 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -6399,6 +6399,46 @@ TEST(operation, boundCRS_to_geogCRS_same_datum_context) { // --------------------------------------------------------------------------- +TEST(operation, boundCRS_to_geogCRS_hubCRS_and_targetCRS_same_but_baseCRS_not) { + const char *wkt = + "COMPD_CS[\"NAD83 + Ellipsoid (US Feet)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"Ellipsoid (US Feet)\",\n" + " VERT_DATUM[\"Ellipsoid\",2002],\n" + " UNIT[\"US survey foot\",0.304800609601219,\n" + " AUTHORITY[\"EPSG\",\"9003\"]],\n" + " AXIS[\"Up\",UP]]]"; + + auto dbContext = DatabaseContext::create(); + auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt); + auto boundCRS = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(boundCRS != nullptr); + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(boundCRS), GeographicCRS::EPSG_4979, ctxt); + ASSERT_EQ(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=unitconvert +z_in=us-ft +z_out=m"); +} + +// --------------------------------------------------------------------------- + TEST(operation, boundCRS_to_boundCRS) { auto utm31 = ProjectedCRS::create( PropertyMap(), GeographicCRS::EPSG_4807, -- cgit v1.2.3 From 5aa99d497433377fc08cb569b6826473dc0c46d1 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 19 May 2020 13:15:40 +0200 Subject: createOperations(): take correctly into account vertical unit change in a case of transformation between Compound of BoundVerticalCRS to GeographicCRS (contributes to fixes #2232) --- src/iso19111/coordinateoperation.cpp | 53 +++++++++++++++++++++++++- test/unit/test_operation.cpp | 74 ++++++++++++++++++++++++++++++++++++ 2 files changed, 125 insertions(+), 2 deletions(-) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 32a4b9c7..6dff57c1 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -14221,8 +14221,57 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( if (vertCRSOfBaseOfBoundSrc && hubSrcGeog) { auto opsFirst = createOperations(sourceCRS, hubSrc, context); if (context.skipHorizontalTransformation) { - if (!opsFirst.empty()) - res = opsFirst; + if (!opsFirst.empty()) { + const auto &hubAxisList = + hubSrcGeog->coordinateSystem()->axisList(); + const auto &targetAxisList = + geogDst->coordinateSystem()->axisList(); + if (hubAxisList.size() == 3 && hubAxisList.size() == 3 && + !hubAxisList[2]->_isEquivalentTo( + targetAxisList[2].get(), + util::IComparable::Criterion::EQUIVALENT)) { + + const auto &srcAxis = hubAxisList[2]; + const double convSrc = srcAxis->unit().conversionToSI(); + const auto &dstAxis = targetAxisList[2]; + const double convDst = dstAxis->unit().conversionToSI(); + const bool srcIsUp = + srcAxis->direction() == cs::AxisDirection::UP; + const bool srcIsDown = + srcAxis->direction() == cs::AxisDirection::DOWN; + const bool dstIsUp = + dstAxis->direction() == cs::AxisDirection::UP; + const bool dstIsDown = + dstAxis->direction() == cs::AxisDirection::DOWN; + const bool heightDepthReversal = + ((srcIsUp && dstIsDown) || (srcIsDown && dstIsUp)); + + const double factor = convSrc / convDst; + auto conv = Conversion::createChangeVerticalUnit( + util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + "Change of vertical unit"), + common::Scale(heightDepthReversal ? -factor : factor)); + auto dbContext = context.context->getAuthorityFactory() + ->databaseContext(); + conv->setCRSs( + hubSrc, + hubSrc->demoteTo2D(std::string(), dbContext) + ->promoteTo3D(std::string(), dbContext, dstAxis), + nullptr); + + for (const auto &op : opsFirst) { + try { + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + {op, conv}, !allowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { + } + } + } else { + res = opsFirst; + } + } return; } else { auto opsSecond = createOperations(hubSrc, targetCRS, context); diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 0cd3b3f5..69bc2a93 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -7255,6 +7255,80 @@ TEST(operation, // --------------------------------------------------------------------------- +TEST(operation, + compoundCRS_with_boundVerticalCRS_from_grids_to_geogCRS_with_ftus_ctxt) { + + auto dbContext = DatabaseContext::create(); + + const char *wktSrc = + "COMPD_CS[\"NAD83 + NAVD88 height - Geoid12B (Meters)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"NAVD88 height - Geoid12B (Meters)\",\n" + " VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n" + " EXTENSION[\"PROJ4_GRIDS\",\"@foo.gtx\"],\n" + " AUTHORITY[\"EPSG\",\"5103\"]],\n" + " UNIT[\"metre\",1.0,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Gravity-related height\",UP],\n" + " AUTHORITY[\"EPSG\",\"5703\"]]]"; + auto objSrc = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktSrc); + auto srcCRS = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(srcCRS != nullptr); + + const char *wktDst = + "COMPD_CS[\"NAD83 + Ellipsoid (US Feet)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"Ellipsoid (US Feet)\",\n" + " VERT_DATUM[\"Ellipsoid\",2002],\n" + " UNIT[\"US survey foot\",0.304800609601219,\n" + " AUTHORITY[\"EPSG\",\"9003\"]],\n" + " AXIS[\"Up\",UP]]]"; + auto objDst = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktDst); + auto dstCRS = nn_dynamic_pointer_cast(objDst); + ASSERT_TRUE(dstCRS != nullptr); + + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(srcCRS), NN_NO_CHECK(dstCRS), ctxt); + ASSERT_GE(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 " + "+step +proj=unitconvert +xy_in=rad +z_in=m " + "+xy_out=deg +z_out=us-ft " + "+step +proj=axisswap +order=2,1"); +} + +// --------------------------------------------------------------------------- + TEST(operation, geocent_to_compoundCRS) { auto objSrc = PROJStringParser().createFromPROJString( "+proj=geocent +datum=WGS84 +units=m +type=crs"); -- cgit v1.2.3 From b6fc4159252b100fa04b808f97269177836f0acf Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 19 May 2020 13:34:26 +0200 Subject: createOperations(): optimization in generated pipeline in a case involving vertical unit change --- src/iso19111/io.cpp | 32 ++++++++++++++++++ test/unit/test_operation.cpp | 79 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 111 insertions(+) diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index ebec053a..e7f8769c 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -6963,6 +6963,38 @@ const std::string &PROJStringFormatter::toString() const { break; } + // +step +proj=unitconvert +xy_in=X1 +z_in=Z1 +xy_out=X2 +z_out=Z2 + // +step +proj=unitconvert +z_in=Z2 +z_out=Z3 + // ==> step +proj=unitconvert +xy_in=X1 +z_in=Z1 +xy_out=X2 + // +z_out=Z3 + if (prevStep.name == "unitconvert" && + curStep.name == "unitconvert" && !prevStep.inverted && + !curStep.inverted && prevStep.paramValues.size() == 4 && + curStep.paramValues.size() == 2 && + prevStep.paramValues[0].keyEquals("xy_in") && + prevStep.paramValues[1].keyEquals("z_in") && + prevStep.paramValues[2].keyEquals("xy_out") && + prevStep.paramValues[3].keyEquals("z_out") && + curStep.paramValues[0].keyEquals("z_in") && + curStep.paramValues[1].keyEquals("z_out") && + prevStep.paramValues[3].value == curStep.paramValues[0].value) { + auto xy_in = prevStep.paramValues[0].value; + auto z_in = prevStep.paramValues[1].value; + auto xy_out = prevStep.paramValues[2].value; + auto z_out = curStep.paramValues[1].value; + d->steps_.erase(iterPrev, iterCur); + iterCur->paramValues.clear(); + iterCur->paramValues.emplace_back( + Step::KeyValue("xy_in", xy_in)); + iterCur->paramValues.emplace_back(Step::KeyValue("z_in", z_in)); + iterCur->paramValues.emplace_back( + Step::KeyValue("xy_out", xy_out)); + iterCur->paramValues.emplace_back( + Step::KeyValue("z_out", z_out)); + changeDone = true; + break; + } + // unitconvert (1), axisswap order=2,1, unitconvert(2) ==> // axisswap order=2,1, unitconvert (1), unitconvert(2) which // will get further optimized by previous case diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 69bc2a93..e1d1e5e2 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -7329,6 +7329,85 @@ TEST(operation, // --------------------------------------------------------------------------- +TEST( + operation, + compoundCRS_with_boundGeogCRS_boundVerticalCRS_from_grids_to_boundGeogCRS_with_ftus_ctxt) { + + // Variant of above but with TOWGS84 in source & target CRS + + auto dbContext = DatabaseContext::create(); + + const char *wktSrc = + "COMPD_CS[\"NAD83 + NAVD88 height - Geoid12B (Meters)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"NAVD88 height - Geoid12B (Meters)\",\n" + " VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n" + " EXTENSION[\"PROJ4_GRIDS\",\"@foo.gtx\"],\n" + " AUTHORITY[\"EPSG\",\"5103\"]],\n" + " UNIT[\"metre\",1.0,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Gravity-related height\",UP],\n" + " AUTHORITY[\"EPSG\",\"5703\"]]]"; + auto objSrc = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktSrc); + auto srcCRS = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(srcCRS != nullptr); + + const char *wktDst = + "COMPD_CS[\"NAD83 + Ellipsoid (US Feet)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"Ellipsoid (US Feet)\",\n" + " VERT_DATUM[\"Ellipsoid\",2002],\n" + " UNIT[\"US survey foot\",0.304800609601219,\n" + " AUTHORITY[\"EPSG\",\"9003\"]],\n" + " AXIS[\"Up\",UP]]]"; + auto objDst = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktDst); + auto dstCRS = nn_dynamic_pointer_cast(objDst); + ASSERT_TRUE(dstCRS != nullptr); + + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(srcCRS), NN_NO_CHECK(dstCRS), ctxt); + ASSERT_EQ(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=rad +z_in=m " + "+xy_out=deg +z_out=us-ft"); +} + +// --------------------------------------------------------------------------- + TEST(operation, geocent_to_compoundCRS) { auto objSrc = PROJStringParser().createFromPROJString( "+proj=geocent +datum=WGS84 +units=m +type=crs"); -- cgit v1.2.3 From b0b33b8447972ac6e60d68213d6c24b0a4989421 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 19 May 2020 13:38:23 +0200 Subject: createOperations(): fix bug when transforming between CompoundCRS and BoundCRS in the general case let to 0 result (contributes to fixes #2232) --- src/iso19111/coordinateoperation.cpp | 5 +++ test/unit/test_operation.cpp | 78 ++++++++++++++++++++++++++++++++++++ 2 files changed, 83 insertions(+) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 6dff57c1..bdb2ad2e 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -15218,10 +15218,15 @@ void CoordinateOperationFactory::Private::createOperationsBoundToCompound( } } } + return; } } } } + + // There might be better things to do, but for now just ignore the + // transformation of the bound CRS + res = createOperations(boundSrc->baseCRS(), targetCRS, context); } //! @endcond diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index e1d1e5e2..51bbf50b 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -7408,6 +7408,84 @@ TEST( // --------------------------------------------------------------------------- +TEST( + operation, + compoundCRS_with_boundVerticalCRS_from_grids_to_boundGeogCRS_with_ftus_ctxt) { + + // Variant of above but with TOWGS84 in target CRS only + + auto dbContext = DatabaseContext::create(); + + const char *wktSrc = + "COMPD_CS[\"NAD83 + NAVD88 height - Geoid12B (Meters)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"NAVD88 height - Geoid12B (Meters)\",\n" + " VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n" + " EXTENSION[\"PROJ4_GRIDS\",\"@foo.gtx\"],\n" + " AUTHORITY[\"EPSG\",\"5103\"]],\n" + " UNIT[\"metre\",1.0,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Gravity-related height\",UP],\n" + " AUTHORITY[\"EPSG\",\"5703\"]]]"; + auto objSrc = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktSrc); + auto srcCRS = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(srcCRS != nullptr); + + const char *wktDst = + "COMPD_CS[\"NAD83 + Ellipsoid (US Feet)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"Ellipsoid (US Feet)\",\n" + " VERT_DATUM[\"Ellipsoid\",2002],\n" + " UNIT[\"US survey foot\",0.304800609601219,\n" + " AUTHORITY[\"EPSG\",\"9003\"]],\n" + " AXIS[\"Up\",UP]]]"; + auto objDst = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktDst); + auto dstCRS = nn_dynamic_pointer_cast(objDst); + ASSERT_TRUE(dstCRS != nullptr); + + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(srcCRS), NN_NO_CHECK(dstCRS), ctxt); + ASSERT_GE(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 " + "+step +proj=unitconvert +xy_in=rad +z_in=m " + "+xy_out=deg +z_out=us-ft " + "+step +proj=axisswap +order=2,1"); +} + +// --------------------------------------------------------------------------- + TEST(operation, geocent_to_compoundCRS) { auto objSrc = PROJStringParser().createFromPROJString( "+proj=geocent +datum=WGS84 +units=m +type=crs"); -- cgit v1.2.3