createOperations(): deal with spherical planetocentric geodetic CRS

This also fixes conversion between geocentric latlong and geodetic latlong
with cs2cs. This was dealt with in PR 1093, but in the wrong direction
(the geocentric latitude must be <= in absolute value to the geodetic one)
The issue here was linked to the semantics of the +geoc specifier, which
affects the semantics of the input coordinates in the forward direction
(+geoc means that the input coordinate is is a geocentric latitude),
which defeats the logic of doing A to B by using the inverse path of A
and the forward path of B.

3 files changed, 156 insertions, 26 deletions

diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp
index 3d09c9fa..b8d9c0b7 100644
--- a/src/iso19111/operation/conversion.cpp
+++ b/src/iso19111/operation/conversion.cpp
@@ -3459,11 +3459,15 @@ void Conversion::_exportToPROJString(
             }
         }
 
-        srcGeogCRS = std::dynamic_pointer_cast<crs::GeographicCRS>(horiz);
-        if (srcGeogCRS) {
+        auto srcGeodCRS = dynamic_cast<const crs::GeodeticCRS *>(horiz.get());
+        if (srcGeodCRS) {
+            srcGeogCRS = std::dynamic_pointer_cast<crs::GeographicCRS>(horiz);
+        }
+        if (srcGeodCRS &&
+            (srcGeogCRS || srcGeodCRS->isSphericalPlanetocentric())) {
             formatter->setOmitProjLongLatIfPossible(true);
             formatter->startInversion();
-            srcGeogCRS->_exportToPROJString(formatter);
+            srcGeodCRS->_exportToPROJString(formatter);
             formatter->stopInversion();
             formatter->setOmitProjLongLatIfPossible(false);
         }
diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
index fc64bd2e..60365713 100644
--- a/src/iso19111/operation/coordinateoperationfactory.cpp
+++ b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -558,6 +558,17 @@ struct CoordinateOperationFactory::Private {
         const crs::GeodeticCRS *geodDst,
         std::vector<CoordinateOperationNNPtr> &res);
 
+    static void createOperationsFromSphericalPlanetocentric(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::GeodeticCRS *geodSrc,
+        std::vector<CoordinateOperationNNPtr> &res);
+
+    static void createOperationsFromBoundOfSphericalPlanetocentric(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::BoundCRS *boundSrc,
+        const crs::GeodeticCRSNNPtr &geodSrcBase,
+        std::vector<CoordinateOperationNNPtr> &res);
+
     static void createOperationsDerivedTo(
         const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
         Private::Context &context, const crs::DerivedCRS *derivedSrc,
@@ -2989,6 +3000,32 @@ CoordinateOperationFactory::Private::createOperations(
         return res;
     }
 
+    if (geodSrc && geodSrc->isSphericalPlanetocentric()) {
+        createOperationsFromSphericalPlanetocentric(sourceCRS, targetCRS,
+                                                    context, geodSrc, res);
+        return res;
+    } else if (geodDst && geodDst->isSphericalPlanetocentric()) {
+        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+    }
+
+    if (boundSrc) {
+        auto geodSrcBase = util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
+            boundSrc->baseCRS());
+        if (geodSrcBase && geodSrcBase->isSphericalPlanetocentric()) {
+            createOperationsFromBoundOfSphericalPlanetocentric(
+                sourceCRS, targetCRS, context, boundSrc,
+                NN_NO_CHECK(geodSrcBase), res);
+            return res;
+        }
+    } else if (boundDst) {
+        auto geodDstBase = util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
+            boundDst->baseCRS());
+        if (geodDstBase && geodDstBase->isSphericalPlanetocentric()) {
+            return applyInverse(
+                createOperations(targetCRS, sourceCRS, context));
+        }
+    }
+
     // If the source is a derived CRS, then chain the inverse of its
     // deriving conversion, with transforms from its baseCRS to the
     // targetCRS
@@ -3895,6 +3932,92 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod(
 
 // ---------------------------------------------------------------------------
 
+void CoordinateOperationFactory::Private::
+    createOperationsFromSphericalPlanetocentric(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::GeodeticCRS *geodSrc,
+        std::vector<CoordinateOperationNNPtr> &res) {
+
+    ENTER_FUNCTION();
+
+    // Create an intermediate geographic CRS with the same datum as the
+    // source spherical planetocentric one
+    std::string interm_crs_name(geodSrc->nameStr());
+    interm_crs_name += " (geographic)";
+    auto interm_crs =
+        util::nn_static_pointer_cast<crs::CRS>(crs::GeographicCRS::create(
+            addDomains(util::PropertyMap().set(
+                           common::IdentifiedObject::NAME_KEY, interm_crs_name),
+                       geodSrc),
+            geodSrc->datum(), geodSrc->datumEnsemble(),
+            cs::EllipsoidalCS::createLatitudeLongitude(
+                common::UnitOfMeasure::DEGREE)));
+
+    auto opFirst = createGeodToGeodPROJBased(sourceCRS, interm_crs);
+    auto opsSecond = createOperations(interm_crs, targetCRS, context);
+    for (const auto &opSecond : opsSecond) {
+        try {
+            res.emplace_back(ConcatenatedOperation::createComputeMetadata(
+                {opFirst, opSecond}, disallowEmptyIntersection));
+        } catch (const InvalidOperationEmptyIntersection &) {
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+void CoordinateOperationFactory::Private::
+    createOperationsFromBoundOfSphericalPlanetocentric(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::BoundCRS *boundSrc,
+        const crs::GeodeticCRSNNPtr &geodSrcBase,
+        std::vector<CoordinateOperationNNPtr> &res) {
+
+    ENTER_FUNCTION();
+
+    // Create an intermediate geographic CRS with the same datum as the
+    // source spherical planetocentric one
+    std::string interm_crs_name(geodSrcBase->nameStr());
+    interm_crs_name += " (geographic)";
+    auto intermGeog =
+        util::nn_static_pointer_cast<crs::CRS>(crs::GeographicCRS::create(
+            addDomains(util::PropertyMap().set(
+                           common::IdentifiedObject::NAME_KEY, interm_crs_name),
+                       geodSrcBase.get()),
+            geodSrcBase->datum(), geodSrcBase->datumEnsemble(),
+            cs::EllipsoidalCS::createLatitudeLongitude(
+                common::UnitOfMeasure::DEGREE)));
+
+    // Create an intermediate boundCRS wrapping the above intermediate
+    // geographic CRS
+    auto transf = boundSrc->transformation()->shallowClone();
+    // keep a reference to the target before patching it with itself
+    // (this is due to our abuse of passing shared_ptr by reference
+    auto transfTarget = transf->targetCRS();
+    setCRSs(transf.get(), intermGeog, transfTarget);
+
+    auto intermBoundCRS =
+        crs::BoundCRS::create(intermGeog, boundSrc->hubCRS(), transf);
+
+    auto opFirst = createGeodToGeodPROJBased(geodSrcBase, intermGeog);
+    setCRSs(opFirst.get(), sourceCRS, intermBoundCRS);
+    auto opsSecond = createOperations(intermBoundCRS, targetCRS, context);
+    for (const auto &opSecond : opsSecond) {
+        try {
+            auto opSecondClone = opSecond->shallowClone();
+            // In theory, we should not need that setCRSs() forcing, but due
+            // how BoundCRS transformations are implemented currently, we
+            // need it in practice.
+            setCRSs(opSecondClone.get(), intermBoundCRS, targetCRS);
+            res.emplace_back(ConcatenatedOperation::createComputeMetadata(
+                {opFirst, opSecondClone}, disallowEmptyIntersection));
+        } catch (const InvalidOperationEmptyIntersection &) {
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 void CoordinateOperationFactory::Private::createOperationsDerivedTo(
     const crs::CRSNNPtr & /*sourceCRS*/, const crs::CRSNNPtr &targetCRS,
     Private::Context &context, const crs::DerivedCRS *derivedSrc,
diff --git a/src/iso19111/operation/transformation.cpp b/src/iso19111/operation/transformation.cpp
index 273b636f..a30e1248 100644
--- a/src/iso19111/operation/transformation.cpp
+++ b/src/iso19111/operation/transformation.cpp
@@ -474,13 +474,16 @@ static void getTransformationType(const crs::CRSNNPtr &sourceCRSIn,
         dynamic_cast<const crs::GeographicCRS *>(sourceCRSIn.get());
     auto targetCRSGeog =
         dynamic_cast<const crs::GeographicCRS *>(targetCRSIn.get());
-    if (!sourceCRSGeog || !targetCRSGeog) {
+    if (!(sourceCRSGeog ||
+          (sourceCRSGeod && sourceCRSGeod->isSphericalPlanetocentric())) ||
+        !(targetCRSGeog ||
+          (targetCRSGeod && targetCRSGeod->isSphericalPlanetocentric()))) {
         throw InvalidOperation("Inconsistent CRS type");
     }
     const auto nSrcAxisCount =
-        sourceCRSGeog->coordinateSystem()->axisList().size();
+        sourceCRSGeod->coordinateSystem()->axisList().size();
     const auto nTargetAxisCount =
-        targetCRSGeog->coordinateSystem()->axisList().size();
+        targetCRSGeod->coordinateSystem()->axisList().size();
     isGeog2D = nSrcAxisCount == 2 && nTargetAxisCount == 2;
     isGeog3D = !isGeog2D && nSrcAxisCount >= 2 && nTargetAxisCount >= 2;
 }
@@ -1003,36 +1006,36 @@ TransformationNNPtr Transformation::createTOWGS84(
             "Invalid number of elements in TOWGS84Parameters");
     }
 
-    crs::CRSPtr transformSourceCRS = sourceCRSIn->extractGeodeticCRS();
-    if (!transformSourceCRS) {
+    auto transformSourceGeodCRS = sourceCRSIn->extractGeodeticCRS();
+    if (!transformSourceGeodCRS) {
         throw InvalidOperation(
             "Cannot find GeodeticCRS in sourceCRS of TOWGS84 transformation");
     }
 
     util::PropertyMap properties;
     properties.set(common::IdentifiedObject::NAME_KEY,
-                   concat("Transformation from ", transformSourceCRS->nameStr(),
-                          " to WGS84"));
-
-    auto targetCRS =
-        dynamic_cast<const crs::GeographicCRS *>(transformSourceCRS.get())
-            ? util::nn_static_pointer_cast<crs::CRS>(
-                  crs::GeographicCRS::EPSG_4326)
-            : util::nn_static_pointer_cast<crs::CRS>(
-                  crs::GeodeticCRS::EPSG_4978);
-
+                   concat("Transformation from ",
+                          transformSourceGeodCRS->nameStr(), " to WGS84"));
+
+    auto targetCRS = dynamic_cast<const crs::GeographicCRS *>(
+                         transformSourceGeodCRS.get()) ||
+                             transformSourceGeodCRS->isSphericalPlanetocentric()
+                         ? util::nn_static_pointer_cast<crs::CRS>(
+                               crs::GeographicCRS::EPSG_4326)
+                         : util::nn_static_pointer_cast<crs::CRS>(
+                               crs::GeodeticCRS::EPSG_4978);
+
+    crs::CRSNNPtr transformSourceCRS = NN_NO_CHECK(transformSourceGeodCRS);
     if (TOWGS84Parameters.size() == 3) {
         return createGeocentricTranslations(
-            properties, NN_NO_CHECK(transformSourceCRS), targetCRS,
-            TOWGS84Parameters[0], TOWGS84Parameters[1], TOWGS84Parameters[2],
-            {});
+            properties, transformSourceCRS, targetCRS, TOWGS84Parameters[0],
+            TOWGS84Parameters[1], TOWGS84Parameters[2], {});
     }
 
-    return createPositionVector(properties, NN_NO_CHECK(transformSourceCRS),
-                                targetCRS, TOWGS84Parameters[0],
-                                TOWGS84Parameters[1], TOWGS84Parameters[2],
-                                TOWGS84Parameters[3], TOWGS84Parameters[4],
-                                TOWGS84Parameters[5], TOWGS84Parameters[6], {});
+    return createPositionVector(
+        properties, transformSourceCRS, targetCRS, TOWGS84Parameters[0],
+        TOWGS84Parameters[1], TOWGS84Parameters[2], TOWGS84Parameters[3],
+        TOWGS84Parameters[4], TOWGS84Parameters[5], TOWGS84Parameters[6], {});
 }
 
 // ---------------------------------------------------------------------------