Add proj_normalize_for_visualization()

Fixes #1301

This function takes the output PJ from proj_create_crs_to_crs(),
and add (or undo) the needed axis swap operations so that the
object returned by proj_normalize_for_visualization() has the usual
GIS axis order.
In this implementation, this does something only if the coordinate
system of the source or target CRS, geographic or projected, has
NORTH, EAST ordering.
CompoundCRS wrapping those objects are also handled.

4 files changed, 171 insertions, 0 deletions

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
diff --git a/src/proj.h b/src/proj.h
index 6f53bb07..a412f266 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -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);