From 601eb6a8058328efccfc792539b2e245614f4676 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 4 Aug 2022 17:06:01 +0200 Subject: [PATCH] createOperations(): take into account axis unit and inversion of target DerivedProjectedCRS Fixes https://lists.osgeo.org/pipermail/proj/2022-August/010697.html --- include/proj/crs.hpp | 9 +++- src/iso19111/crs.cpp | 33 ++++++++++---- src/iso19111/operation/conversion.cpp | 12 +++++- test/unit/test_operationfactory.cpp | 62 +++++++++++++++++++++++++++ 4 files changed, 106 insertions(+), 10 deletions(-) diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp index 1eef5ee2db..6d10a21794 100644 --- a/include/proj/crs.hpp +++ b/include/proj/crs.hpp @@ -661,6 +661,10 @@ class PROJ_GCC_DLL ProjectedCRS final : public DerivedCRS, addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter, bool axisSpecFound) const; + PROJ_INTERNAL static void addUnitConvertAndAxisSwap( + const std::vector &axisListIn, + io::PROJStringFormatter *formatter, bool axisSpecFound); + PROJ_INTERNAL void _exportToWKT(io::WKTFormatter *formatter) const override; // throw(io::FormattingException) @@ -1272,7 +1276,10 @@ class PROJ_GCC_DLL DerivedProjectedCRS final : public DerivedCRS { //! @cond Doxygen_Suppress PROJ_INTERNAL void _exportToWKT(io::WKTFormatter *formatter) const override; // throw(io::FormattingException) - //! @endcond + + PROJ_INTERNAL void + addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter) const; + //! @endcond protected: PROJ_INTERNAL diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index d19cd6cc8b..9fcc51acec 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -4474,9 +4474,16 @@ ProjectedCRS::alterParametersLinearUnit(const common::UnitOfMeasure &unit, //! @cond Doxygen_Suppress void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter, bool axisSpecFound) const { - const auto &axisList = d->coordinateSystem()->axisList(); - const auto &unit = axisList[0]->unit(); - const auto *zUnit = axisList.size() == 3 ? &(axisList[2]->unit()) : nullptr; + ProjectedCRS::addUnitConvertAndAxisSwap(d->coordinateSystem()->axisList(), + formatter, axisSpecFound); +} + +void ProjectedCRS::addUnitConvertAndAxisSwap( + const std::vector &axisListIn, + io::PROJStringFormatter *formatter, bool axisSpecFound) { + const auto &unit = axisListIn[0]->unit(); + const auto *zUnit = + axisListIn.size() == 3 ? &(axisListIn[2]->unit()) : nullptr; if (!unit._isEquivalentTo(common::UnitOfMeasure::METRE, util::IComparable::Criterion::EQUIVALENT) || (zUnit && @@ -4518,8 +4525,8 @@ void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter, if (!axisSpecFound && (!formatter->getCRSExport() || formatter->getLegacyCRSToCRSContext())) { - const auto &dir0 = axisList[0]->direction(); - const auto &dir1 = axisList[1]->direction(); + const auto &dir0 = axisListIn[0]->direction(); + const auto &dir1 = axisListIn[1]->direction(); if (!(&dir0 == &cs::AxisDirection::EAST && &dir1 == &cs::AxisDirection::NORTH) && // For polar projections, that have south+south direction, @@ -4528,7 +4535,7 @@ void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter, const char *order[2] = {nullptr, nullptr}; for (int i = 0; i < 2; i++) { - const auto &dir = axisList[i]->direction(); + const auto &dir = axisListIn[i]->direction(); if (&dir == &cs::AxisDirection::WEST) order[i] = "-1"; else if (&dir == &cs::AxisDirection::EAST) @@ -4546,8 +4553,8 @@ void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter, formatter->addParam("order", orderStr); } } else { - const auto &name0 = axisList[0]->nameStr(); - const auto &name1 = axisList[1]->nameStr(); + const auto &name0 = axisListIn[0]->nameStr(); + const auto &name1 = axisListIn[1]->nameStr(); const bool northingEasting = ci_starts_with(name0, "northing") && ci_starts_with(name1, "easting"); // case of EPSG:32661 ["WGS 84 / UPS North (N,E)]" @@ -6548,6 +6555,16 @@ bool DerivedProjectedCRS::_isEquivalentTo( // --------------------------------------------------------------------------- +//! @cond Doxygen_Suppress +void DerivedProjectedCRS::addUnitConvertAndAxisSwap( + io::PROJStringFormatter *formatter) const { + ProjectedCRS::addUnitConvertAndAxisSwap(coordinateSystem()->axisList(), + formatter, false); +} +//! @endcond + +// --------------------------------------------------------------------------- + //! @cond Doxygen_Suppress struct TemporalCRS::Private {}; //! @endcond diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp index 37ca29a8e4..1e9d8df2b0 100644 --- a/src/iso19111/operation/conversion.cpp +++ b/src/iso19111/operation/conversion.cpp @@ -4086,8 +4086,12 @@ void Conversion::_exportToPROJString( } } + auto derivedProjCRS = + dynamic_cast(horiz); + // horiz != nullptr: only to make clang static analyzer happy - if (!bEllipsoidParametersDone && horiz != nullptr) { + if (!bEllipsoidParametersDone && horiz != nullptr && + derivedProjCRS == nullptr) { auto targetGeodCRS = horiz->extractGeodeticCRS(); auto targetGeogCRS = std::dynamic_pointer_cast(targetGeodCRS); @@ -4119,6 +4123,12 @@ void Conversion::_exportToPROJString( if (projCRS->hasOver()) { formatter->addParam("over"); } + } else { + if (derivedProjCRS) { + formatter->pushOmitZUnitConversion(); + derivedProjCRS->addUnitConvertAndAxisSwap(formatter); + formatter->popOmitZUnitConversion(); + } } auto derivedGeographicCRS = diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 40db0726bc..21fd3068ea 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -8411,3 +8411,65 @@ TEST(operation, createOperation_ossfuzz_47873_simplified_if_i_might_say) { } catch (const std::exception &) { } } + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_derived_projected_crs) { + + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto src = authFactory->createCoordinateReferenceSystem("6507"); + + auto wkt = + "DERIVEDPROJCRS[\"Custom Site Calibrated CRS\",\n" + " BASEPROJCRS[\"NAD83(2011) / Mississippi East (ftUS)\",\n" + " BASEGEOGCRS[\"NAD83(2011)\",\n" + " DATUM[\"NAD83 (National Spatial Reference System " + "2011)\",\n" + " ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n" + " LENGTHUNIT[\"metre\",1]]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]]],\n" + " CONVERSION[\"SPCS83 Mississippi East zone (US Survey " + "feet)\",\n" + " METHOD[\"Transverse Mercator\",\n" + " ID[\"EPSG\",9807]],\n" + " PARAMETER[\"Latitude of natural origin\",29.5,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433],\n" + " ID[\"EPSG\",8801]],\n" + " PARAMETER[\"Longitude of natural " + "origin\",-88.8333333333333,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433],\n" + " ID[\"EPSG\",8802]],\n" + " PARAMETER[\"Scale factor at natural origin\",0.99995,\n" + " SCALEUNIT[\"unity\",1],\n" + " ID[\"EPSG\",8805]],\n" + " PARAMETER[\"False easting\",984250,\n" + " LENGTHUNIT[\"US survey foot\",0.304800609601219],\n" + " ID[\"EPSG\",8806]],\n" + " PARAMETER[\"False northing\",0,\n" + " LENGTHUNIT[\"US survey foot\",0.304800609601219],\n" + " ID[\"EPSG\",8807]]]],\n" + " DERIVINGCONVERSION[\"Affine transformation as PROJ-based\",\n" + " METHOD[\"PROJ-based operation method: " + "+proj=pipeline +step +proj=unitconvert +xy_in=m +xy_out=us-ft " + "+step +proj=affine +xoff=20 " + "+step +proj=unitconvert +xy_in=us-ft +xy_out=m\"]],\n" + " CS[Cartesian,2],\n" + " AXIS[\"northing (Y)\",north,\n" + " LENGTHUNIT[\"US survey foot\",0.304800609601219]],\n" + " AXIS[\"easting (X)\",east,\n" + " LENGTHUNIT[\"US survey foot\",0.304800609601219]],\n" + " REMARK[\"EPSG:6507 with 20 feet offset and axis inversion\"]]"; + auto objDst = WKTParser().createFromWKT(wkt); + auto dst = nn_dynamic_pointer_cast(objDst); + ASSERT_TRUE(dst != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + src, NN_NO_CHECK(dst)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=affine +xoff=20 " + "+step +proj=axisswap +order=2,1"); +}