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.

6 files changed, 366 insertions, 150 deletions

diff --git a/src/apps/cs2cs.cpp b/src/apps/cs2cs.cpp
index 19f924d4..03fd4dba 100644
--- a/src/apps/cs2cs.cpp
+++ b/src/apps/cs2cs.cpp
@@ -59,10 +59,10 @@
 
 static PJ *transformation = nullptr;
 
-static bool srcIsGeog = false;
+static bool srcIsLongLat = false;
 static double srcToRadians = 0.0;
 
-static bool destIsGeog = false;
+static bool destIsLongLat = false;
 static double destToRadians = 0.0;
 static bool destIsLatLong = false;
 
@@ -158,7 +158,7 @@ static void process(FILE *fid)
 
         if (data.u != HUGE_VAL) {
 
-            if (srcIsGeog) {
+            if (srcIsLongLat) {
                 /* dmstor gives values to radians. Convert now to the SRS unit
                  */
                 data.u /= srcToRadians;
@@ -179,7 +179,7 @@ static void process(FILE *fid)
         if (data.u == HUGE_VAL) /* error output */
             fputs(oterr, stdout);
 
-        else if (destIsGeog && !oform) { /*ascii DMS output */
+        else if (destIsLongLat && !oform) { /*ascii DMS output */
 
             // rtodms() expect radians: convert from the output SRS unit
             data.u *= destToRadians;
@@ -206,7 +206,7 @@ static void process(FILE *fid)
             }
 
         } else { /* x-y or decimal degree ascii output */
-            if (destIsGeog) {
+            if (destIsLongLat) {
                 data.v *= destToRadians * RAD_TO_DEG;
                 data.u *= destToRadians * RAD_TO_DEG;
             }
@@ -239,7 +239,7 @@ static void process(FILE *fid)
 /************************************************************************/
 
 static PJ *instantiate_crs(const std::string &definition,
-                               bool &isGeog, double &toRadians,
+                               bool &isLongLatCS, double &toRadians,
                                bool &isLatFirst) {
     PJ *crs = proj_create(nullptr,
                           pj_add_type_crs_if_needed(definition).c_str());
@@ -247,7 +247,7 @@ static PJ *instantiate_crs(const std::string &definition,
         return nullptr;
     }
 
-    isGeog = false;
+    isLongLatCS = false;
     toRadians = 0.0;
     isLatFirst = false;
 
@@ -259,11 +259,11 @@ static PJ *instantiate_crs(const std::string &definition,
         type = proj_get_type(crs);
     }
     if (type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
-        type == PJ_TYPE_GEOGRAPHIC_3D_CRS) {
+        type == PJ_TYPE_GEOGRAPHIC_3D_CRS ||
+        type == PJ_TYPE_GEODETIC_CRS) {
         auto cs = proj_crs_get_coordinate_system(nullptr, crs);
         assert(cs);
 
-        isGeog = true;
         const char *axisName = "";
         proj_cs_get_axis_info(nullptr, cs, 0,
                                   &axisName, // name,
@@ -277,6 +277,9 @@ static PJ *instantiate_crs(const std::string &definition,
         isLatFirst =
             NS_PROJ::internal::ci_find(std::string(axisName), "latitude") !=
             std::string::npos;
+        isLongLatCS = isLatFirst ||
+            NS_PROJ::internal::ci_find(std::string(axisName), "longitude") !=
+            std::string::npos;
 
         proj_destroy(cs);
     }
@@ -736,7 +739,7 @@ int main(int argc, char **argv) {
     PJ *src = nullptr;
     if (!fromStr.empty()) {
         bool ignored;
-        src = instantiate_crs(fromStr, srcIsGeog,
+        src = instantiate_crs(fromStr, srcIsLongLat,
                               srcToRadians, ignored);
         if (!src) {
             emess(3, "cannot instantiate source coordinate system");
@@ -745,7 +748,7 @@ int main(int argc, char **argv) {
 
     PJ *dst = nullptr;
     if (!toStr.empty()) {
-        dst = instantiate_crs(toStr, destIsGeog,
+        dst = instantiate_crs(toStr, destIsLongLat,
                               destToRadians, destIsLatLong);
         if (!dst) {
             emess(3, "cannot instantiate target coordinate system");
@@ -760,7 +763,7 @@ int main(int argc, char **argv) {
             emess(3,
                   "missing target CRS and source CRS is not a projected CRS");
         }
-        destIsGeog = true;
+        destIsLongLat = true;
     } else if (fromStr.empty()) {
         assert(dst);
         bool ignored;
@@ -770,7 +773,7 @@ int main(int argc, char **argv) {
             emess(3,
                   "missing source CRS and target CRS is not a projected CRS");
         }
-        srcIsGeog = true;
+        srcIsLongLat = true;
     }
     proj_destroy(src);
     proj_destroy(dst);
@@ -835,13 +838,13 @@ int main(int argc, char **argv) {
     }
 
     /* set input formatting control */
-    if (!srcIsGeog)
+    if (!srcIsLongLat)
         informat = strtod;
     else {
         informat = dmstor;
     }
 
-    if (!destIsGeog && !oform)
+    if (!destIsLongLat && !oform)
         oform = "%.2f";
 
     /* process input file list */
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 6f511043..be593878 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -1996,6 +1996,61 @@ void GeodeticCRS::addGeocentricUnitConversionIntoPROJString(
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+void GeodeticCRS::addAngularUnitConvertAndAxisSwap(
+    io::PROJStringFormatter *formatter) const {
+    const auto &axisList = coordinateSystem()->axisList();
+
+    formatter->addStep("unitconvert");
+    formatter->addParam("xy_in", "rad");
+    if (axisList.size() == 3 && !formatter->omitZUnitConversion()) {
+        formatter->addParam("z_in", "m");
+    }
+    {
+        const auto &unitHoriz = axisList[0]->unit();
+        const auto projUnit = unitHoriz.exportToPROJString();
+        if (projUnit.empty()) {
+            formatter->addParam("xy_out", unitHoriz.conversionToSI());
+        } else {
+            formatter->addParam("xy_out", projUnit);
+        }
+    }
+    if (axisList.size() == 3 && !formatter->omitZUnitConversion()) {
+        const auto &unitZ = axisList[2]->unit();
+        auto projVUnit = unitZ.exportToPROJString();
+        if (projVUnit.empty()) {
+            formatter->addParam("z_out", unitZ.conversionToSI());
+        } else {
+            formatter->addParam("z_out", projVUnit);
+        }
+    }
+
+    const char *order[2] = {nullptr, nullptr};
+    const char *one = "1";
+    const char *two = "2";
+    for (int i = 0; i < 2; i++) {
+        const auto &dir = axisList[i]->direction();
+        if (&dir == &cs::AxisDirection::WEST) {
+            order[i] = "-1";
+        } else if (&dir == &cs::AxisDirection::EAST) {
+            order[i] = one;
+        } else if (&dir == &cs::AxisDirection::SOUTH) {
+            order[i] = "-2";
+        } else if (&dir == &cs::AxisDirection::NORTH) {
+            order[i] = two;
+        }
+    }
+    if (order[0] && order[1] && (order[0] != one || order[1] != two)) {
+        formatter->addStep("axisswap");
+        char orderStr[10];
+        sprintf(orderStr, "%.2s,%.2s", order[0], order[1]);
+        formatter->addParam("order", orderStr);
+    }
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 void GeodeticCRS::_exportToPROJString(
     io::PROJStringFormatter *formatter) const // throw(io::FormattingException)
 {
@@ -2007,19 +2062,38 @@ void GeodeticCRS::_exportToPROJString(
         return;
     }
 
-    if (!isGeocentric()) {
-        io::FormattingException::Throw(
-            "GeodeticCRS::exportToPROJString() only "
-            "supports geocentric coordinate systems");
-    }
+    if (isGeocentric()) {
+        if (!formatter->getCRSExport()) {
+            formatter->addStep("cart");
+        } else {
+            formatter->addStep("geocent");
+        }
 
-    if (!formatter->getCRSExport()) {
-        formatter->addStep("cart");
+        addDatumInfoToPROJString(formatter);
+        addGeocentricUnitConversionIntoPROJString(formatter);
+    } else if (isSphericalPlanetocentric()) {
+        if (!formatter->getCRSExport()) {
+            formatter->addStep("geoc");
+            addDatumInfoToPROJString(formatter);
+            addAngularUnitConvertAndAxisSwap(formatter);
+        } else {
+            io::FormattingException::Throw(
+                "GeodeticCRS::exportToPROJString() not supported on spherical "
+                "planetocentric coordinate systems");
+            // The below code now works as input to PROJ, but I'm not sure we
+            // want to propagate this, given that we got cs2cs doing conversion
+            // in the wrong direction in past versions.
+            /*formatter->addStep("longlat");
+            formatter->addParam("geoc");
+
+            addDatumInfoToPROJString(formatter);*/
+        }
     } else {
-        formatter->addStep("geocent");
+        io::FormattingException::Throw(
+            "GeodeticCRS::exportToPROJString() only "
+            "supports geocentric or spherical planetocentric "
+            "coordinate systems");
     }
-    addDatumInfoToPROJString(formatter);
-    addGeocentricUnitConversionIntoPROJString(formatter);
 }
 //! @endcond
 
@@ -2884,61 +2958,6 @@ GeographicCRS::demoteTo2D(const std::string &newName,
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
-void GeographicCRS::addAngularUnitConvertAndAxisSwap(
-    io::PROJStringFormatter *formatter) const {
-    const auto &axisList = coordinateSystem()->axisList();
-
-    formatter->addStep("unitconvert");
-    formatter->addParam("xy_in", "rad");
-    if (axisList.size() == 3 && !formatter->omitZUnitConversion()) {
-        formatter->addParam("z_in", "m");
-    }
-    {
-        const auto &unitHoriz = axisList[0]->unit();
-        const auto projUnit = unitHoriz.exportToPROJString();
-        if (projUnit.empty()) {
-            formatter->addParam("xy_out", unitHoriz.conversionToSI());
-        } else {
-            formatter->addParam("xy_out", projUnit);
-        }
-    }
-    if (axisList.size() == 3 && !formatter->omitZUnitConversion()) {
-        const auto &unitZ = axisList[2]->unit();
-        auto projVUnit = unitZ.exportToPROJString();
-        if (projVUnit.empty()) {
-            formatter->addParam("z_out", unitZ.conversionToSI());
-        } else {
-            formatter->addParam("z_out", projVUnit);
-        }
-    }
-
-    const char *order[2] = {nullptr, nullptr};
-    const char *one = "1";
-    const char *two = "2";
-    for (int i = 0; i < 2; i++) {
-        const auto &dir = axisList[i]->direction();
-        if (&dir == &cs::AxisDirection::WEST) {
-            order[i] = "-1";
-        } else if (&dir == &cs::AxisDirection::EAST) {
-            order[i] = one;
-        } else if (&dir == &cs::AxisDirection::SOUTH) {
-            order[i] = "-2";
-        } else if (&dir == &cs::AxisDirection::NORTH) {
-            order[i] = two;
-        }
-    }
-    if (order[0] && order[1] && (order[0] != one || order[1] != two)) {
-        formatter->addStep("axisswap");
-        char orderStr[10];
-        sprintf(orderStr, "%.2s,%.2s", order[0], order[1]);
-        formatter->addParam("order", orderStr);
-    }
-}
-//! @endcond
-
-// ---------------------------------------------------------------------------
-
-//! @cond Doxygen_Suppress
 void GeographicCRS::_exportToPROJString(
     io::PROJStringFormatter *formatter) const // throw(io::FormattingException)
 {
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 1373c583..5d7e951e 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -8639,10 +8639,10 @@ struct PROJStringParser::Private {
     PrimeMeridianNNPtr buildPrimeMeridian(Step &step);
     GeodeticReferenceFrameNNPtr buildDatum(Step &step,
                                            const std::string &title);
-    GeographicCRSNNPtr buildGeographicCRS(int iStep, int iUnitConvert,
-                                          int iAxisSwap, bool ignorePROJAxis);
+    GeodeticCRSNNPtr buildGeodeticCRS(int iStep, int iUnitConvert,
+                                      int iAxisSwap, bool ignorePROJAxis);
     GeodeticCRSNNPtr buildGeocentricCRS(int iStep, int iUnitConvert);
-    CRSNNPtr buildProjectedCRS(int iStep, GeographicCRSNNPtr geogCRS,
+    CRSNNPtr buildProjectedCRS(int iStep, const GeodeticCRSNNPtr &geogCRS,
                                int iUnitConvert, int iAxisSwap);
     CRSNNPtr buildBoundOrCompoundCRSIfNeeded(int iStep, CRSNNPtr crs);
     UnitOfMeasure buildUnit(Step &step, const std::string &unitsParamName,
@@ -8656,6 +8656,9 @@ struct PROJStringParser::Private {
 
     EllipsoidalCSNNPtr buildEllipsoidalCS(int iStep, int iUnitConvert,
                                           int iAxisSwap, bool ignorePROJAxis);
+
+    SphericalCSNNPtr buildSphericalCS(int iStep, int iUnitConvert,
+                                      int iAxisSwap, bool ignorePROJAxis);
 };
 
 // ---------------------------------------------------------------------------
@@ -9282,10 +9285,13 @@ PROJStringParser::Private::processAxisSwap(Step &step,
     assert(iAxisSwap < 0 || ci_equal(steps_[iAxisSwap].name, "axisswap"));
 
     const bool isGeographic = unit.type() == UnitOfMeasure::Type::ANGULAR;
+    const bool isSpherical = isGeographic && hasParamValue(step, "geoc");
     const auto &eastName =
-        isGeographic ? AxisName::Longitude : AxisName::Easting;
-    const auto &eastAbbev =
-        isGeographic ? AxisAbbreviation::lon : AxisAbbreviation::E;
+        isSpherical ? "Planetocentric longitude"
+                    : isGeographic ? AxisName::Longitude : AxisName::Easting;
+    const auto &eastAbbev = isSpherical ? "V"
+                                        : isGeographic ? AxisAbbreviation::lon
+                                                       : AxisAbbreviation::E;
     const auto &eastDir = isGeographic
                               ? AxisDirection::EAST
                               : (axisType == AxisType::NORTH_POLE)
@@ -9301,9 +9307,11 @@ PROJStringParser::Private::processAxisSwap(Step &step,
             : nullMeridian);
 
     const auto &northName =
-        isGeographic ? AxisName::Latitude : AxisName::Northing;
-    const auto &northAbbev =
-        isGeographic ? AxisAbbreviation::lat : AxisAbbreviation::N;
+        isSpherical ? "Planetocentric latitude"
+                    : isGeographic ? AxisName::Latitude : AxisName::Northing;
+    const auto &northAbbev = isSpherical ? "U"
+                                         : isGeographic ? AxisAbbreviation::lat
+                                                        : AxisAbbreviation::N;
     const auto &northDir = isGeographic
                                ? AxisDirection::NORTH
                                : (axisType == AxisType::NORTH_POLE)
@@ -9320,15 +9328,19 @@ PROJStringParser::Private::processAxisSwap(Step &step,
                         .as_nullable()
                   : nullMeridian);
 
-    CoordinateSystemAxisNNPtr west =
-        createAxis(isGeographic ? AxisName::Longitude : AxisName::Westing,
-                   isGeographic ? AxisAbbreviation::lon : std::string(),
-                   AxisDirection::WEST, unit);
+    CoordinateSystemAxisNNPtr west = createAxis(
+        isSpherical ? "Planetocentric longitude"
+                    : isGeographic ? AxisName::Longitude : AxisName::Westing,
+        isSpherical ? "V"
+                    : isGeographic ? AxisAbbreviation::lon : std::string(),
+        AxisDirection::WEST, unit);
 
-    CoordinateSystemAxisNNPtr south =
-        createAxis(isGeographic ? AxisName::Latitude : AxisName::Southing,
-                   isGeographic ? AxisAbbreviation::lat : std::string(),
-                   AxisDirection::SOUTH, unit);
+    CoordinateSystemAxisNNPtr south = createAxis(
+        isSpherical ? "Planetocentric latitude"
+                    : isGeographic ? AxisName::Latitude : AxisName::Southing,
+        isSpherical ? "U"
+                    : isGeographic ? AxisAbbreviation::lat : std::string(),
+        AxisDirection::SOUTH, unit);
 
     std::vector<CoordinateSystemAxisNNPtr> axis{east, north};
 
@@ -9428,6 +9440,42 @@ EllipsoidalCSNNPtr PROJStringParser::Private::buildEllipsoidalCS(
 
 // ---------------------------------------------------------------------------
 
+SphericalCSNNPtr PROJStringParser::Private::buildSphericalCS(
+    int iStep, int iUnitConvert, int iAxisSwap, bool ignorePROJAxis) {
+    auto &step = steps_[iStep];
+    assert(iUnitConvert < 0 ||
+           ci_equal(steps_[iUnitConvert].name, "unitconvert"));
+
+    UnitOfMeasure angularUnit = UnitOfMeasure::DEGREE;
+    if (iUnitConvert >= 0) {
+        auto &stepUnitConvert = steps_[iUnitConvert];
+        const std::string *xy_in = &getParamValue(stepUnitConvert, "xy_in");
+        const std::string *xy_out = &getParamValue(stepUnitConvert, "xy_out");
+        if (stepUnitConvert.inverted) {
+            std::swap(xy_in, xy_out);
+        }
+        if (iUnitConvert < iStep) {
+            std::swap(xy_in, xy_out);
+        }
+        if (xy_in->empty() || xy_out->empty() || *xy_in != "rad" ||
+            (*xy_out != "rad" && *xy_out != "deg" && *xy_out != "grad")) {
+            throw ParsingException("unhandled values for xy_in and/or xy_out");
+        }
+        if (*xy_out == "rad") {
+            angularUnit = UnitOfMeasure::RADIAN;
+        } else if (*xy_out == "grad") {
+            angularUnit = UnitOfMeasure::GRAD;
+        }
+    }
+
+    std::vector<CoordinateSystemAxisNNPtr> axis = processAxisSwap(
+        step, angularUnit, iAxisSwap, AxisType::REGULAR, ignorePROJAxis);
+
+    return SphericalCS::create(emptyPropertyMap, axis[0], axis[1]);
+}
+
+// ---------------------------------------------------------------------------
+
 static double getNumericValue(const std::string &paramValue,
                               bool *pHasError = nullptr) {
     try {
@@ -9447,7 +9495,7 @@ namespace {
 template <class T> inline void ignoreRetVal(T) {}
 } // namespace
 
-GeographicCRSNNPtr PROJStringParser::Private::buildGeographicCRS(
+GeodeticCRSNNPtr PROJStringParser::Private::buildGeodeticCRS(
     int iStep, int iUnitConvert, int iAxisSwap, bool ignorePROJAxis) {
     auto &step = steps_[iStep];
 
@@ -9462,8 +9510,6 @@ GeographicCRSNNPtr PROJStringParser::Private::buildGeographicCRS(
 
     auto props = PropertyMap().set(IdentifiedObject::NAME_KEY,
                                    title.empty() ? "unknown" : title);
-    auto cs =
-        buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, ignorePROJAxis);
 
     if (l_isGeographicStep &&
         (hasUnusedParameters(step) ||
@@ -9472,7 +9518,17 @@ GeographicCRSNNPtr PROJStringParser::Private::buildGeographicCRS(
     }
     props.set("IMPLICIT_CS", true);
 
-    return GeographicCRS::create(props, datum, cs);
+    if (!hasParamValue(step, "geoc")) {
+        auto cs =
+            buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, ignorePROJAxis);
+
+        return GeographicCRS::create(props, datum, cs);
+    } else {
+        auto cs =
+            buildSphericalCS(iStep, iUnitConvert, iAxisSwap, ignorePROJAxis);
+
+        return GeodeticCRS::create(props, datum, cs);
+    }
 }
 
 // ---------------------------------------------------------------------------
@@ -9623,8 +9679,10 @@ static bool is_in_stringlist(const std::string &str, const char *stringlist) {
 
 // ---------------------------------------------------------------------------
 
-CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
-    int iStep, GeographicCRSNNPtr geogCRS, int iUnitConvert, int iAxisSwap) {
+CRSNNPtr
+PROJStringParser::Private::buildProjectedCRS(int iStep,
+                                             const GeodeticCRSNNPtr &geodCRS,
+                                             int iUnitConvert, int iAxisSwap) {
     auto &step = steps_[iStep];
     const auto mappings = getMappingsFromPROJName(step.name);
     const MethodMapping *mapping = mappings.empty() ? nullptr : mappings[0];
@@ -9650,7 +9708,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
     }
 
     if (mapping) {
-        mapping = selectSphericalOrEllipsoidal(mapping, geogCRS);
+        mapping = selectSphericalOrEllipsoidal(mapping, geodCRS);
     }
 
     assert(isProjectedStep(step.name));
@@ -9660,7 +9718,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
     const auto &title = title_;
 
     if (!buildPrimeMeridian(step)->longitude()._isEquivalentTo(
-            geogCRS->primeMeridian()->longitude(),
+            geodCRS->primeMeridian()->longitude(),
             util::IComparable::Criterion::EQUIVALENT)) {
         throw ParsingException("inconsistent pm values between projectedCRS "
                                "and its base geographicalCRS");
@@ -9926,7 +9984,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                     }
                     if (k >= 0 && k <= 1) {
                         const double es =
-                            geogCRS->ellipsoid()->squaredEccentricity();
+                            geodCRS->ellipsoid()->squaredEccentricity();
                         if (es < 0 || es == 1) {
                             throw ParsingException("Invalid flattening");
                         }
@@ -10037,10 +10095,16 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
              {"PROJ ob_tran o_proj=longlat", "PROJ ob_tran o_proj=lonlat",
               "PROJ ob_tran o_proj=latlon", "PROJ ob_tran o_proj=latlong"}) {
             if (starts_with(methodName, substr)) {
-                return DerivedGeographicCRS::create(
-                    PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"),
-                    geogCRS, NN_NO_CHECK(conv),
-                    buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, false));
+                auto geogCRS =
+                    util::nn_dynamic_pointer_cast<GeographicCRS>(geodCRS);
+                if (geogCRS) {
+                    return DerivedGeographicCRS::create(
+                        PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                          "unnamed"),
+                        NN_NO_CHECK(geogCRS), NN_NO_CHECK(conv),
+                        buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap,
+                                           false));
+                }
             }
         }
     }
@@ -10048,11 +10112,11 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
     std::vector<CoordinateSystemAxisNNPtr> axis =
         processAxisSwap(step, unit, iAxisSwap, axisType, false);
 
-    auto csGeogCRS = geogCRS->coordinateSystem();
-    auto cs = csGeogCRS->axisList().size() == 2
+    auto csGeodCRS = geodCRS->coordinateSystem();
+    auto cs = csGeodCRS->axisList().size() == 2
                   ? CartesianCS::create(emptyPropertyMap, axis[0], axis[1])
                   : CartesianCS::create(emptyPropertyMap, axis[0], axis[1],
-                                        csGeogCRS->axisList()[2]);
+                                        csGeodCRS->axisList()[2]);
 
     auto props = PropertyMap().set(IdentifiedObject::NAME_KEY,
                                    title.empty() ? "unknown" : title);
@@ -10066,7 +10130,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
         bWebMercator
             ? createPseudoMercator(
                   props.set(IdentifiedObject::NAME_KEY, webMercatorName), cs)
-            : ProjectedCRS::create(props, geogCRS, NN_NO_CHECK(conv), cs);
+            : ProjectedCRS::create(props, geodCRS, NN_NO_CHECK(conv), cs);
 
     return crs;
 }
@@ -10537,8 +10601,8 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
             // First run is dry run to mark all recognized/unrecognized tokens
             for (int iter = 0; iter < 2; iter++) {
                 auto obj = d->buildBoundOrCompoundCRSIfNeeded(
-                    0, d->buildGeographicCRS(iFirstGeogStep, iFirstUnitConvert,
-                                             iFirstAxisSwap, false));
+                    0, d->buildGeodeticCRS(iFirstGeogStep, iFirstUnitConvert,
+                                           iFirstAxisSwap, false));
                 if (iter == 1) {
                     return nn_static_pointer_cast<BaseObject>(obj);
                 }
@@ -10555,14 +10619,14 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
                     iProjStep,
                     d->buildProjectedCRS(
                         iProjStep,
-                        d->buildGeographicCRS(iFirstGeogStep,
-                                              iFirstUnitConvert < iFirstGeogStep
-                                                  ? iFirstUnitConvert
-                                                  : -1,
-                                              iFirstAxisSwap < iFirstGeogStep
-                                                  ? iFirstAxisSwap
-                                                  : -1,
-                                              true),
+                        d->buildGeodeticCRS(iFirstGeogStep,
+                                            iFirstUnitConvert < iFirstGeogStep
+                                                ? iFirstUnitConvert
+                                                : -1,
+                                            iFirstAxisSwap < iFirstGeogStep
+                                                ? iFirstAxisSwap
+                                                : -1,
+                                            true),
                         iFirstUnitConvert < iFirstGeogStep ? iSecondUnitConvert
                                                            : iFirstUnitConvert,
                         iFirstAxisSwap < iFirstGeogStep ? iSecondAxisSwap
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], {});
 }
 
 // ---------------------------------------------------------------------------