diff options
| -rw-r--r-- | include/proj/coordinateoperation.hpp | 2 | ||||
| -rw-r--r-- | include/proj/crs.hpp | 4 | ||||
| -rw-r--r-- | src/iso19111/c_api.cpp | 32 | ||||
| -rw-r--r-- | src/iso19111/coordinateoperation.cpp | 40 | ||||
| -rw-r--r-- | src/iso19111/crs.cpp | 98 | ||||
| -rw-r--r-- | src/proj.h | 1 | ||||
| -rw-r--r-- | test/unit/test_c_api.cpp | 27 | ||||
| -rw-r--r-- | test/unit/test_operation.cpp | 120 |
8 files changed, 324 insertions, 0 deletions
diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp index 4397ed60..9e7812c7 100644 --- a/include/proj/coordinateoperation.hpp +++ b/include/proj/coordinateoperation.hpp @@ -152,6 +152,8 @@ class PROJ_GCC_DLL CoordinateOperation : public common::ObjectUsage, PROJ_DLL static const std::string OPERATION_VERSION_KEY; + PROJ_DLL CoordinateOperationNNPtr normalizeForVisualization() const; + protected: PROJ_INTERNAL CoordinateOperation(); PROJ_INTERNAL CoordinateOperation(const CoordinateOperation &other); diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp index ed3463cd..33bf3b61 100644 --- a/include/proj/crs.hpp +++ b/include/proj/crs.hpp @@ -126,6 +126,10 @@ class PROJ_GCC_DLL CRS : public common::ObjectUsage { PROJ_FOR_TEST CRSNNPtr alterCSLinearUnit(const common::UnitOfMeasure &unit) const; + PROJ_INTERNAL bool mustAxisOrderBeSwitchedForVisualization() const; + + PROJ_INTERNAL CRSNNPtr normalizeForVisualization() const; + //! @endcond protected: diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index e62d3f56..1f4cecde 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -6767,3 +6767,35 @@ int proj_cs_get_axis_info(PJ_CONTEXT *ctx, const PJ *cs, int index, } return true; } + +// --------------------------------------------------------------------------- + +/** \brief Returns a PJ* object whose axis order is the one expected for + * visualization purposes. + * + * The input object must be a coordinate operation, that has been created with + * proj_create_crs_to_crs(). + * If the axis order of its source or target CRS is northing,easting, then an + * axis swap operation will be inserted. + * + * @param ctx PROJ context, or NULL for default context + * @param obj Object of type CoordinateOperation, + * created with proj_create_crs_to_crs() (must not be NULL) + * @return a new PJ* object to free with proj_destroy() in case of success, or + * nullptr in case of error + */ +PJ *proj_normalize_for_visualization(PJ_CONTEXT *ctx, const PJ *obj) { + auto co = dynamic_cast<const CoordinateOperation *>(obj->iso_obj.get()); + if (!co) { + proj_log_error(ctx, __FUNCTION__, "Object is not a CoordinateOperation " + "created with " + "proj_create_crs_to_crs"); + return nullptr; + } + try { + return pj_obj_create(ctx, co->normalizeForVisualization()); + } catch (const std::exception &e) { + proj_log_debug(ctx, __FUNCTION__, e.what()); + return nullptr; + } +} diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 1abb74b3..58e97272 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -765,6 +765,46 @@ void CoordinateOperation::setProperties( // --------------------------------------------------------------------------- +/** \brief Return a variation of the current coordinate operation whose axis + * order is the one expected for visualization purposes. + */ +CoordinateOperationNNPtr +CoordinateOperation::normalizeForVisualization() const { + auto l_sourceCRS = sourceCRS(); + auto l_targetCRS = targetCRS(); + if (!l_sourceCRS || !l_targetCRS) { + throw util::UnsupportedOperationException( + "Cannot retrieve source or target CRS"); + } + const bool swapSource = + l_sourceCRS->mustAxisOrderBeSwitchedForVisualization(); + const bool swapTarget = + l_targetCRS->mustAxisOrderBeSwitchedForVisualization(); + auto l_this = NN_NO_CHECK(std::dynamic_pointer_cast<CoordinateOperation>( + shared_from_this().as_nullable())); + if (!swapSource && !swapTarget) { + return l_this; + } + std::vector<CoordinateOperationNNPtr> subOps; + if (swapSource) { + auto op = Conversion::createAxisOrderReversal(false); + op->setCRSs(l_sourceCRS->normalizeForVisualization(), + NN_NO_CHECK(l_sourceCRS), nullptr); + subOps.emplace_back(op); + } + subOps.emplace_back(l_this); + if (swapTarget) { + auto op = Conversion::createAxisOrderReversal(false); + op->setCRSs(NN_NO_CHECK(l_targetCRS), + l_targetCRS->normalizeForVisualization(), nullptr); + subOps.emplace_back(op); + } + return util::nn_static_pointer_cast<CoordinateOperation>( + ConcatenatedOperation::createComputeMetadata(subOps, true)); +} + +// --------------------------------------------------------------------------- + //! @cond Doxygen_Suppress struct OperationMethod::Private { util::optional<std::string> formula_{}; diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index b51d03c9..1ef7dcda 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -591,6 +591,104 @@ CRSNNPtr CRS::alterId(const std::string &authName, // --------------------------------------------------------------------------- +//! @cond Doxygen_Suppress + +static bool isAxisListNorthEast( + const std::vector<cs::CoordinateSystemAxisNNPtr> &axisList) { + const auto &dir0 = axisList[0]->direction(); + const auto &dir1 = axisList[1]->direction(); + return (&dir0 == &cs::AxisDirection::NORTH && + &dir1 == &cs::AxisDirection::EAST); +} +// --------------------------------------------------------------------------- + +bool CRS::mustAxisOrderBeSwitchedForVisualization() const { + + const CompoundCRS *compoundCRS = dynamic_cast<const CompoundCRS *>(this); + if (compoundCRS) { + const auto &comps = compoundCRS->componentReferenceSystems(); + if (!comps.empty()) { + return comps[0]->mustAxisOrderBeSwitchedForVisualization(); + } + } + + const GeographicCRS *geogCRS = dynamic_cast<const GeographicCRS *>(this); + if (geogCRS) { + return isAxisListNorthEast(geogCRS->coordinateSystem()->axisList()); + } + + const ProjectedCRS *projCRS = dynamic_cast<const ProjectedCRS *>(this); + if (projCRS) { + return isAxisListNorthEast(projCRS->coordinateSystem()->axisList()); + } + + return false; +} + +//! @endcond + +// --------------------------------------------------------------------------- + +//! @cond Doxygen_Suppress + +CRSNNPtr CRS::normalizeForVisualization() const { + auto props = util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + nameStr() + " (with axis order normalized for visualization)"); + + const CompoundCRS *compoundCRS = dynamic_cast<const CompoundCRS *>(this); + if (compoundCRS) { + const auto &comps = compoundCRS->componentReferenceSystems(); + if (!comps.empty()) { + std::vector<CRSNNPtr> newComps; + newComps.emplace_back(comps[0]->normalizeForVisualization()); + for (size_t i = 1; i < comps.size(); i++) { + newComps.emplace_back(comps[i]); + } + return util::nn_static_pointer_cast<CRS>( + CompoundCRS::create(props, newComps)); + } + } + + const GeographicCRS *geogCRS = dynamic_cast<const GeographicCRS *>(this); + if (geogCRS) { + const auto &axisList = geogCRS->coordinateSystem()->axisList(); + if (isAxisListNorthEast(axisList)) { + auto cs = axisList.size() == 2 + ? cs::EllipsoidalCS::create(util::PropertyMap(), + axisList[1], axisList[0]) + : cs::EllipsoidalCS::create(util::PropertyMap(), + axisList[1], axisList[0], + axisList[2]); + return util::nn_static_pointer_cast<CRS>(GeographicCRS::create( + props, geogCRS->datum(), geogCRS->datumEnsemble(), cs)); + } + } + + const ProjectedCRS *projCRS = dynamic_cast<const ProjectedCRS *>(this); + if (projCRS) { + const auto &axisList = projCRS->coordinateSystem()->axisList(); + if (isAxisListNorthEast(axisList)) { + auto cs = + axisList.size() == 2 + ? cs::CartesianCS::create(util::PropertyMap(), axisList[1], + axisList[0]) + : cs::CartesianCS::create(util::PropertyMap(), axisList[1], + axisList[0], axisList[2]); + return util::nn_static_pointer_cast<CRS>( + ProjectedCRS::create(props, projCRS->baseCRS(), + projCRS->derivingConversionRef(), cs)); + } + } + + return NN_NO_CHECK( + std::static_pointer_cast<CRS>(shared_from_this().as_nullable())); +} + +//! @endcond + +// --------------------------------------------------------------------------- + /** \brief Identify the CRS with reference CRSs. * * The candidate CRSs are either hard-coded, or looked in the database when @@ -356,6 +356,7 @@ int PROJ_DLL proj_context_get_use_proj4_init_rules(PJ_CONTEXT *ctx, int from_leg PJ PROJ_DLL *proj_create (PJ_CONTEXT *ctx, const char *definition); PJ PROJ_DLL *proj_create_argv (PJ_CONTEXT *ctx, int argc, char **argv); PJ PROJ_DLL *proj_create_crs_to_crs(PJ_CONTEXT *ctx, const char *source_crs, const char *target_crs, PJ_AREA *area); +PJ PROJ_DLL *proj_normalize_for_visualization(PJ_CONTEXT *ctx, const PJ* obj); PJ PROJ_DLL *proj_destroy (PJ *P); diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index 486ab0c7..12e98d7e 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -3220,4 +3220,31 @@ TEST_F(CApi, proj_get_crs_info_list_from_database) { proj_crs_info_list_destroy(list); } } + +// --------------------------------------------------------------------------- + +TEST_F(CApi, proj_normalize_for_visualization) { + + { + auto P = proj_create(m_ctxt, "+proj=utm +zone=31 +ellps=WGS84"); + ObjectKeeper keeper_P(P); + ASSERT_NE(P, nullptr); + auto Pnormalized = proj_normalize_for_visualization(m_ctxt, P); + ObjectKeeper keeper_Pnormalized(Pnormalized); + EXPECT_EQ(Pnormalized, nullptr); + } + + auto P = proj_create_crs_to_crs(m_ctxt, "EPSG:4326", "EPSG:32631", nullptr); + ObjectKeeper keeper_P(P); + ASSERT_NE(P, nullptr); + auto Pnormalized = proj_normalize_for_visualization(m_ctxt, P); + ObjectKeeper keeper_Pnormalized(Pnormalized); + ASSERT_NE(Pnormalized, nullptr); + auto projstr = proj_as_proj_string(m_ctxt, Pnormalized, PJ_PROJ_5, nullptr); + ASSERT_NE(projstr, nullptr); + EXPECT_EQ(std::string(projstr), + "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=utm +zone=31 +ellps=WGS84"); +} + } // namespace diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 0f719a6d..a1774416 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -7648,3 +7648,123 @@ TEST(operation, validateParameters) { EXPECT_EQ(validation, expected); } } + +// --------------------------------------------------------------------------- + +TEST(operation, normalizeForVisualization) { + + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + + // Source(geographic) must be inverted + { + auto src = authFactory->createCoordinateReferenceSystem("4326"); + auto dst = authFactory->createCoordinateReferenceSystem("32631"); + auto op = + CoordinateOperationFactory::create()->createOperation(src, dst); + ASSERT_TRUE(op != nullptr); + auto opNormalized = op->normalizeForVisualization(); + EXPECT_FALSE(opNormalized->_isEquivalentTo(op.get())); + EXPECT_EQ(opNormalized->exportToPROJString( + PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=utm +zone=31 +ellps=WGS84"); + } + + // Target(geographic) must be inverted + { + auto src = authFactory->createCoordinateReferenceSystem("32631"); + auto dst = authFactory->createCoordinateReferenceSystem("4326"); + auto op = + CoordinateOperationFactory::create()->createOperation(src, dst); + ASSERT_TRUE(op != nullptr); + auto opNormalized = op->normalizeForVisualization(); + EXPECT_FALSE(opNormalized->_isEquivalentTo(op.get())); + EXPECT_EQ(opNormalized->exportToPROJString( + PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +inv +proj=utm +zone=31 +ellps=WGS84 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); + } + + // Source(geographic) and target(projected) must be inverted + { + auto src = authFactory->createCoordinateReferenceSystem("4326"); + auto dst = authFactory->createCoordinateReferenceSystem("3040"); + auto op = + CoordinateOperationFactory::create()->createOperation(src, dst); + ASSERT_TRUE(op != nullptr); + auto opNormalized = op->normalizeForVisualization(); + EXPECT_FALSE(opNormalized->_isEquivalentTo(op.get())); + EXPECT_EQ(opNormalized->exportToPROJString( + PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=utm +zone=28 +ellps=GRS80"); + } + + // No inversion + { + auto src = authFactory->createCoordinateReferenceSystem("32631"); + auto dst = authFactory->createCoordinateReferenceSystem("32632"); + auto op = + CoordinateOperationFactory::create()->createOperation(src, dst); + ASSERT_TRUE(op != nullptr); + auto opNormalized = op->normalizeForVisualization(); + EXPECT_TRUE(opNormalized->_isEquivalentTo(op.get())); + } + + // Source(compoundCRS) and target(geographic 3D) must be inverted + { + auto ctxt = + CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setUsePROJAlternativeGridNames(false); + auto src = CompoundCRS::create( + PropertyMap(), + std::vector<CRSNNPtr>{ + authFactory->createCoordinateReferenceSystem("4326"), + // EGM2008 height + authFactory->createCoordinateReferenceSystem("3855")}); + auto list = CoordinateOperationFactory::create()->createOperations( + src, + authFactory->createCoordinateReferenceSystem("4979"), // WGS 84 3D + ctxt); + ASSERT_EQ(list.size(), 2U); + auto op = list[1]; + auto opNormalized = op->normalizeForVisualization(); + EXPECT_FALSE(opNormalized->_isEquivalentTo(op.get())); + EXPECT_EQ(opNormalized->exportToPROJString( + PROJStringFormatter::create( + PROJStringFormatter::Convention::PROJ_5, + authFactory->databaseContext()) + .get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=vgridshift +grids=egm08_25.gtx +multiplier=1 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); + } + + // Source(boundCRS) and target(geographic) must be inverted + { + auto src = BoundCRS::createFromTOWGS84( + GeographicCRS::EPSG_4269, std::vector<double>{1, 2, 3, 4, 5, 6, 7}); + auto dst = authFactory->createCoordinateReferenceSystem("4326"); + auto op = + CoordinateOperationFactory::create()->createOperation(src, dst); + ASSERT_TRUE(op != nullptr); + auto opNormalized = op->normalizeForVisualization(); + EXPECT_FALSE(opNormalized->_isEquivalentTo(op.get())); + EXPECT_EQ(opNormalized->exportToPROJString( + PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=push +v_3 " + "+step +proj=cart +ellps=GRS80 " + "+step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 " + "+convention=position_vector " + "+step +inv +proj=cart +ellps=WGS84 " + "+step +proj=pop +v_3 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); + } +} |