createOperations(): use an explicit conversion operation for geodetic <--> geocentric latitude

4 files changed, 138 insertions, 9 deletions

diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp
index c1e5cb44..f7bb8ec8 100644
--- a/src/iso19111/operation/conversion.cpp
+++ b/src/iso19111/operation/conversion.cpp
@@ -79,6 +79,7 @@ constexpr double UTM_SCALE_FACTOR = 0.9996;
 constexpr double UTM_FALSE_EASTING = 500000.0;
 constexpr double UTM_NORTH_FALSE_NORTHING = 0.0;
 constexpr double UTM_SOUTH_FALSE_NORTHING = 10000000.0;
+
 //! @endcond
 
 // ---------------------------------------------------------------------------
@@ -2403,6 +2404,28 @@ Conversion::createGeographicGeocentric(const crs::CRSNNPtr &sourceCRS,
 
 // ---------------------------------------------------------------------------
 
+/** \brief Instantiate a conversion between a GeographicCRS and a spherical
+ * planetocentric GeodeticCRS
+ *
+ * This method peforms conversion between geodetic latitude and geocentric
+ * latitude
+ *
+ * @return a new Conversion.
+ */
+ConversionNNPtr
+Conversion::createGeographicGeocentricLatitude(const crs::CRSNNPtr &sourceCRS,
+                                               const crs::CRSNNPtr &targetCRS) {
+    auto properties = util::PropertyMap().set(
+        common::IdentifiedObject::NAME_KEY,
+        buildOpName("Conversion", sourceCRS, targetCRS));
+    auto conv = create(
+        properties, PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE, {});
+    conv->setCRSs(sourceCRS, targetCRS, nullptr);
+    return conv;
+}
+
+// ---------------------------------------------------------------------------
+
 InverseConversion::InverseConversion(const ConversionNNPtr &forward)
     : Conversion(
           OperationMethod::create(createPropertiesForInverse(forward->method()),
@@ -2509,6 +2532,15 @@ CoordinateOperationNNPtr Conversion::inverse() const {
         return conv;
     }
 
+    if (method()->nameStr() ==
+        PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE) {
+        auto conv =
+            create(createPropertiesForInverse(this, false, false),
+                   PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE, {});
+        conv->setCRSs(this, true);
+        return conv;
+    }
+
     return InverseConversion::create(NN_NO_CHECK(
         util::nn_dynamic_pointer_cast<Conversion>(shared_from_this())));
 }
@@ -3447,6 +3479,77 @@ void Conversion::_exportToPROJString(
     auto l_sourceCRS = sourceCRS();
     auto l_targetCRS = targetCRS();
 
+    if (methodName == PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE) {
+
+        const auto extractGeodeticCRSIfGeodeticCRSOrEquivalent =
+            [](const crs::CRSPtr &crs) {
+                auto geodCRS = std::dynamic_pointer_cast<crs::GeodeticCRS>(crs);
+                if (!geodCRS) {
+                    auto compoundCRS =
+                        std::dynamic_pointer_cast<crs::CompoundCRS>(crs);
+                    if (compoundCRS) {
+                        const auto &components =
+                            compoundCRS->componentReferenceSystems();
+                        if (!components.empty()) {
+                            geodCRS =
+                                util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
+                                    components[0]);
+                            if (!geodCRS) {
+                                auto boundCRS = util::nn_dynamic_pointer_cast<
+                                    crs::BoundCRS>(components[0]);
+                                if (boundCRS) {
+                                    geodCRS = util::nn_dynamic_pointer_cast<
+                                        crs::GeodeticCRS>(boundCRS->baseCRS());
+                                }
+                            }
+                        }
+                    } else {
+                        auto boundCRS =
+                            std::dynamic_pointer_cast<crs::BoundCRS>(crs);
+                        if (boundCRS) {
+                            geodCRS =
+                                util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
+                                    boundCRS->baseCRS());
+                        }
+                    }
+                }
+                return geodCRS;
+            };
+
+        auto sourceCRSGeod = dynamic_cast<const crs::GeodeticCRS *>(
+            extractGeodeticCRSIfGeodeticCRSOrEquivalent(l_sourceCRS).get());
+        auto targetCRSGeod = dynamic_cast<const crs::GeodeticCRS *>(
+            extractGeodeticCRSIfGeodeticCRSOrEquivalent(l_targetCRS).get());
+        if (sourceCRSGeod && targetCRSGeod) {
+            auto sourceCRSGeog =
+                dynamic_cast<const crs::GeographicCRS *>(sourceCRSGeod);
+            auto targetCRSGeog =
+                dynamic_cast<const crs::GeographicCRS *>(targetCRSGeod);
+            bool isSrcGeocentricLat =
+                sourceCRSGeod->isSphericalPlanetocentric();
+            bool isSrcGeographic = sourceCRSGeog != nullptr;
+            bool isTargetGeocentricLat =
+                targetCRSGeod->isSphericalPlanetocentric();
+            bool isTargetGeographic = targetCRSGeog != nullptr;
+            if ((isSrcGeocentricLat && isTargetGeographic) ||
+                (isSrcGeographic && isTargetGeocentricLat)) {
+
+                formatter->startInversion();
+                sourceCRSGeod->_exportToPROJString(formatter);
+                formatter->stopInversion();
+
+                targetCRSGeod->_exportToPROJString(formatter);
+
+                return;
+            }
+        }
+
+        throw io::FormattingException("Invalid nature of source and/or "
+                                      "targetCRS for Geographic latitude / "
+                                      "Geocentric latitude"
+                                      "conversion");
+    }
+
     crs::GeographicCRSPtr srcGeogCRS;
     if (!formatter->getCRSExport() && l_sourceCRS && applySourceCRSModifiers) {
 
diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
index 60365713..1b1cae9b 100644
--- a/src/iso19111/operation/coordinateoperationfactory.cpp
+++ b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -2993,13 +2993,6 @@ CoordinateOperationFactory::Private::createOperations(
         }
     }
 
-    // Special case if both CRS are geodetic
-    if (geodSrc && geodDst && !derivedSrc && !derivedDst) {
-        createOperationsGeodToGeod(sourceCRS, targetCRS, context, geodSrc,
-                                   geodDst, res);
-        return res;
-    }
-
     if (geodSrc && geodSrc->isSphericalPlanetocentric()) {
         createOperationsFromSphericalPlanetocentric(sourceCRS, targetCRS,
                                                     context, geodSrc, res);
@@ -3008,6 +3001,13 @@ CoordinateOperationFactory::Private::createOperations(
         return applyInverse(createOperations(targetCRS, sourceCRS, context));
     }
 
+    // Special case if both CRS are geodetic
+    if (geodSrc && geodDst && !derivedSrc && !derivedDst) {
+        createOperationsGeodToGeod(sourceCRS, targetCRS, context, geodSrc,
+                                   geodDst, res);
+        return res;
+    }
+
     if (boundSrc) {
         auto geodSrcBase = util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
             boundSrc->baseCRS());
@@ -3940,6 +3940,24 @@ void CoordinateOperationFactory::Private::
 
     ENTER_FUNCTION();
 
+    const auto IsSameDatum = [&context,
+                              &geodSrc](const crs::GeodeticCRS *geodDst) {
+        const auto &authFactory = context.context->getAuthorityFactory();
+        const auto dbContext =
+            authFactory ? authFactory->databaseContext().as_nullable()
+                        : nullptr;
+
+        return geodSrc->datumNonNull(dbContext)->_isEquivalentTo(
+            geodDst->datumNonNull(dbContext).get(),
+            util::IComparable::Criterion::EQUIVALENT);
+    };
+    auto geogDst = dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
+    if (geogDst && IsSameDatum(geogDst)) {
+        res.emplace_back(Conversion::createGeographicGeocentricLatitude(
+            sourceCRS, targetCRS));
+        return;
+    }
+
     // Create an intermediate geographic CRS with the same datum as the
     // source spherical planetocentric one
     std::string interm_crs_name(geodSrc->nameStr());
@@ -3953,7 +3971,8 @@ void CoordinateOperationFactory::Private::
             cs::EllipsoidalCS::createLatitudeLongitude(
                 common::UnitOfMeasure::DEGREE)));
 
-    auto opFirst = createGeodToGeodPROJBased(sourceCRS, interm_crs);
+    auto opFirst =
+        Conversion::createGeographicGeocentricLatitude(sourceCRS, interm_crs);
     auto opsSecond = createOperations(interm_crs, targetCRS, context);
     for (const auto &opSecond : opsSecond) {
         try {
@@ -3999,7 +4018,8 @@ void CoordinateOperationFactory::Private::
     auto intermBoundCRS =
         crs::BoundCRS::create(intermGeog, boundSrc->hubCRS(), transf);
 
-    auto opFirst = createGeodToGeodPROJBased(geodSrcBase, intermGeog);
+    auto opFirst =
+        Conversion::createGeographicGeocentricLatitude(geodSrcBase, intermGeog);
     setCRSs(opFirst.get(), sourceCRS, intermBoundCRS);
     auto opsSecond = createOperations(intermBoundCRS, targetCRS, context);
     for (const auto &opSecond : opsSecond) {
diff --git a/src/iso19111/operation/parammappings.cpp b/src/iso19111/operation/parammappings.cpp
index 9acc9e93..df5ae7c5 100644
--- a/src/iso19111/operation/parammappings.cpp
+++ b/src/iso19111/operation/parammappings.cpp
@@ -1340,6 +1340,8 @@ static const MethodMapping otherMethodMappings[] = {
     {EPSG_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC,
      EPSG_CODE_METHOD_GEOGRAPHIC_GEOCENTRIC, nullptr, nullptr, nullptr,
      nullptr},
+    {PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE, 0, nullptr, nullptr,
+     nullptr, nullptr},
     {EPSG_NAME_METHOD_LONGITUDE_ROTATION, EPSG_CODE_METHOD_LONGITUDE_ROTATION,
      nullptr, nullptr, nullptr, paramsLongitudeRotation},
     {EPSG_NAME_METHOD_AFFINE_PARAMETRIC_TRANSFORMATION,
diff --git a/src/proj_constants.h b/src/proj_constants.h
index f67a0583..5985f1b1 100644
--- a/src/proj_constants.h
+++ b/src/proj_constants.h
@@ -693,5 +693,9 @@
 #define EPSG_NAME_METHOD_GEOGRAPHIC_TOPOCENTRIC "Geographic/topocentric conversions"
 #define EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC 9837
 
+/* ------------------------------------------------------------------------ */
+
+#define PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE                   \
+    "Geographic latitude / Geocentric latitude"
 
 #endif /* PROJ_CONSTANTS_INCLUDED */