From 47cece34a149945e700b05b1c2e22817f00bb879 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 28 Feb 2020 16:39:52 +0100 Subject: [PATCH] createOperations(): fix wrong pipeline generation with CRS that has +nadgrids= and +pm= Fixes issue reported at https://lists.osgeo.org/pipermail/gdal-dev/2020-February/051749.html The generated pipeline assumes that the input coordinates for the grid transformation were related to the non-Greenwich based datum, so we must compensate for that and add logic to go back to Greenwich. --- src/iso19111/coordinateoperation.cpp | 39 +++++++++++-- src/iso19111/crs.cpp | 26 +++++++-- src/iso19111/io.cpp | 86 +++++++++++++++++----------- test/unit/test_operation.cpp | 59 ++++++++++++++++++- 4 files changed, 168 insertions(+), 42 deletions(-) diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 16325423b8..63b11761b3 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -13850,12 +13850,27 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( 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; + subops.emplace_back(opFirst); + if (opIntermediate) { + subops.emplace_back(NN_NO_CHECK(opIntermediate)); + } + subops.emplace_back(boundSrc->transformation()); res.emplace_back( ConcatenatedOperation::createComputeMetadata( - {opFirst, boundSrc->transformation()}, - !allowEmptyIntersection)); + subops, !allowEmptyIntersection)); } catch (const InvalidOperationEmptyIntersection &) { } } @@ -13873,13 +13888,29 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context); auto opsLast = createOperations(hubSrc, targetCRS, context); if (!opsFirst.empty() && !opsLast.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) { for (const auto &opLast : opsLast) { try { + std::vector subops; + subops.emplace_back(opFirst); + if (opIntermediate) { + subops.emplace_back(NN_NO_CHECK(opIntermediate)); + } + subops.emplace_back(boundSrc->transformation()); + subops.emplace_back(opLast); res.emplace_back( ConcatenatedOperation::createComputeMetadata( - {opFirst, boundSrc->transformation(), opLast}, - !allowEmptyIntersection)); + subops, !allowEmptyIntersection)); } catch (const InvalidOperationEmptyIntersection &) { } } diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index dbb68f0baa..45e2309ca0 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -4525,9 +4525,27 @@ BoundCRS::createFromTOWGS84(const CRSNNPtr &baseCRSIn, */ BoundCRSNNPtr BoundCRS::createFromNadgrids(const CRSNNPtr &baseCRSIn, const std::string &filename) { - const CRSPtr sourceGeographicCRS = baseCRSIn->extractGeographicCRS(); + const auto sourceGeographicCRS = baseCRSIn->extractGeographicCRS(); auto transformationSourceCRS = - sourceGeographicCRS ? sourceGeographicCRS : baseCRSIn.as_nullable(); + sourceGeographicCRS + ? NN_NO_CHECK(std::static_pointer_cast(sourceGeographicCRS)) + : baseCRSIn; + if (sourceGeographicCRS != nullptr && + sourceGeographicCRS->datum() != nullptr && + sourceGeographicCRS->primeMeridian()->longitude().getSIValue() != 0.0) { + transformationSourceCRS = GeographicCRS::create( + util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, + sourceGeographicCRS->nameStr() + + " (with Greenwich prime meridian)"), + datum::GeodeticReferenceFrame::create( + util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + sourceGeographicCRS->datum()->nameStr() + + " (with Greenwich prime meridian)"), + sourceGeographicCRS->datum()->ellipsoid(), + util::optional(), datum::PrimeMeridian::GREENWICH), + sourceGeographicCRS->coordinateSystem()); + } std::string transformationName = transformationSourceCRS->nameStr(); transformationName += " to WGS84"; @@ -4536,8 +4554,8 @@ BoundCRSNNPtr BoundCRS::createFromNadgrids(const CRSNNPtr &baseCRSIn, operation::Transformation::createNTv2( util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, transformationName), - NN_NO_CHECK(transformationSourceCRS), GeographicCRS::EPSG_4326, - filename, std::vector())); + transformationSourceCRS, GeographicCRS::EPSG_4326, filename, + std::vector())); } // --------------------------------------------------------------------------- diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 8a8ad52661..5b12270722 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -4060,6 +4060,52 @@ WKTParser::Private::buildCompoundCRS(const WKTNodeNNPtr &node) { // --------------------------------------------------------------------------- +static CRSNNPtr +createBoundCRSSourceTransformationCRS(const crs::CRSPtr &sourceCRS, + const crs::CRSPtr &targetCRS) { + CRSPtr sourceTransformationCRS; + if (dynamic_cast(targetCRS.get())) { + GeographicCRSPtr sourceGeographicCRS = + sourceCRS->extractGeographicCRS(); + sourceTransformationCRS = sourceGeographicCRS; + if (sourceGeographicCRS) { + if (sourceGeographicCRS->datum() != nullptr && + sourceGeographicCRS->primeMeridian() + ->longitude() + .getSIValue() != 0.0) { + sourceTransformationCRS = + GeographicCRS::create( + util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + sourceGeographicCRS->nameStr() + + " (with Greenwich prime meridian)"), + datum::GeodeticReferenceFrame::create( + util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + sourceGeographicCRS->datum()->nameStr() + + " (with Greenwich prime meridian)"), + sourceGeographicCRS->datum()->ellipsoid(), + util::optional(), + datum::PrimeMeridian::GREENWICH), + sourceGeographicCRS->coordinateSystem()) + .as_nullable(); + } + } else { + sourceTransformationCRS = + std::dynamic_pointer_cast(sourceCRS); + if (!sourceTransformationCRS) { + throw ParsingException( + "Cannot find GeographicCRS or VerticalCRS in sourceCRS"); + } + } + } else { + sourceTransformationCRS = sourceCRS; + } + return NN_NO_CHECK(sourceTransformationCRS); +} + +// --------------------------------------------------------------------------- + BoundCRSNNPtr WKTParser::Private::buildBoundCRS(const WKTNodeNNPtr &node) { const auto *nodeP = node->GP(); auto &abridgedNode = @@ -4103,23 +4149,10 @@ BoundCRSNNPtr WKTParser::Private::buildBoundCRS(const WKTNodeNNPtr &node) { consumeParameters(abridgedNode, true, parameters, values, defaultLinearUnit, defaultAngularUnit); - CRSPtr sourceTransformationCRS; - if (dynamic_cast(targetCRS.get())) { - sourceTransformationCRS = sourceCRS->extractGeographicCRS(); - if (!sourceTransformationCRS) { - sourceTransformationCRS = - std::dynamic_pointer_cast(sourceCRS); - if (!sourceTransformationCRS) { - throw ParsingException( - "Cannot find GeographicCRS or VerticalCRS in sourceCRS"); - } - } - } else { - sourceTransformationCRS = sourceCRS; - } - + const auto sourceTransformationCRS( + createBoundCRSSourceTransformationCRS(sourceCRS, targetCRS)); auto transformation = Transformation::create( - buildProperties(abridgedNode), NN_NO_CHECK(sourceTransformationCRS), + buildProperties(abridgedNode), sourceTransformationCRS, NN_NO_CHECK(targetCRS), nullptr, buildProperties(methodNode), parameters, values, std::vector()); @@ -5265,24 +5298,11 @@ BoundCRSNNPtr JSONParser::buildBoundCRS(const json &j) { values.emplace_back(ParameterValue::create(getMeasure(param))); } - CRSPtr sourceTransformationCRS; - if (dynamic_cast(targetCRS.get())) { - sourceTransformationCRS = sourceCRS->extractGeographicCRS(); - if (!sourceTransformationCRS) { - sourceTransformationCRS = - std::dynamic_pointer_cast(sourceCRS.as_nullable()); - if (!sourceTransformationCRS) { - throw ParsingException( - "Cannot find GeographicCRS or VerticalCRS in sourceCRS"); - } - } - } else { - sourceTransformationCRS = sourceCRS; - } - + const auto sourceTransformationCRS( + createBoundCRSSourceTransformationCRS(sourceCRS, targetCRS)); auto transformation = Transformation::create( - buildProperties(transformationJ), NN_NO_CHECK(sourceTransformationCRS), - targetCRS, nullptr, buildProperties(methodJ), parameters, values, + buildProperties(transformationJ), sourceTransformationCRS, targetCRS, + nullptr, buildProperties(methodJ), parameters, values, std::vector()); return BoundCRS::create(sourceCRS, targetCRS, transformation); diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 84415efb05..19b71cee3f 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -6567,6 +6567,64 @@ TEST(operation, ETRS89_3D_to_proj_string_with_geoidgrids_nadgrids) { // --------------------------------------------------------------------------- +TEST(operation, nadgrids_with_pm) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto objSrc = PROJStringParser().createFromPROJString( + "+proj=tmerc +lat_0=39.66666666666666 +lon_0=1 +k=1 +x_0=200000 " + "+y_0=300000 +ellps=intl +nadgrids=foo.gsb +pm=lisbon " + "+units=m +type=crs"); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + auto dst = authFactory->createCoordinateReferenceSystem("4326"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(src), dst, ctxt); + ASSERT_EQ(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +inv +proj=tmerc +lat_0=39.6666666666667 +lon_0=1 " + "+k=1 +x_0=200000 +y_0=300000 +ellps=intl +pm=lisbon " + // Check that there is no extra +step +proj=longlat +pm=lisbon + "+step +proj=hgridshift +grids=foo.gsb " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1"); + + // ETRS89 + dst = authFactory->createCoordinateReferenceSystem("4258"); + list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(src), dst, ctxt); + ASSERT_EQ(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +inv +proj=tmerc +lat_0=39.6666666666667 +lon_0=1 " + "+k=1 +x_0=200000 +y_0=300000 +ellps=intl +pm=lisbon " + // Check that there is no extra +step +proj=longlat +pm=lisbon + "+step +proj=hgridshift +grids=foo.gsb " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1"); + + // From WKT BOUNDCRS + auto formatter = WKTFormatter::create(WKTFormatter::Convention::WKT2_2019); + auto src_wkt = src->exportToWKT(formatter.get()); + auto objFromWkt = WKTParser().createFromWKT(src_wkt); + auto crsFromWkt = nn_dynamic_pointer_cast(objFromWkt); + ASSERT_TRUE(crsFromWkt); + list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(crsFromWkt), dst, ctxt); + ASSERT_EQ(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +inv +proj=tmerc +lat_0=39.6666666666667 +lon_0=1 " + "+k=1 +x_0=200000 +y_0=300000 +ellps=intl +pm=lisbon " + // Check that there is no extra +step +proj=longlat +pm=lisbon + "+step +proj=hgridshift +grids=foo.gsb " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1"); +} + +// --------------------------------------------------------------------------- + TEST(operation, WGS84_G1762_to_compoundCRS_with_bound_vertCRS) { auto authFactoryEPSG = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); @@ -6767,7 +6825,6 @@ TEST(operation, transformation_BEV_AT_to_PROJ_string) { "+multiplier=1"); } - // --------------------------------------------------------------------------- TEST(operation, transformation_longitude_rotation_to_PROJ_string) {