cs2cs: upgrade to use proj_create_crs_to_crs()

9 files changed, 935 insertions, 283 deletions

diff --git a/src/c_api.cpp b/src/c_api.cpp
index ba1b9534..3b51d905 100644
--- a/src/c_api.cpp
+++ b/src/c_api.cpp
@@ -39,6 +39,7 @@
 
 #include "proj/common.hpp"
 #include "proj/coordinateoperation.hpp"
+#include "proj/coordinatesystem.hpp"
 #include "proj/crs.hpp"
 #include "proj/datum.hpp"
 #include "proj/io.hpp"
@@ -56,6 +57,7 @@
 
 using namespace NS_PROJ::common;
 using namespace NS_PROJ::crs;
+using namespace NS_PROJ::cs;
 using namespace NS_PROJ::datum;
 using namespace NS_PROJ::io;
 using namespace NS_PROJ::internal;
@@ -1354,8 +1356,8 @@ int proj_obj_prime_meridian_get_parameters(PJ_OBJ *prime_meridian,
 
 // ---------------------------------------------------------------------------
 
-/** \brief Return the base CRS of a BoundCRS or the source CRS of a
- * CoordinateOperation.
+/** \brief Return the base CRS of a BoundCRS or a DerivedCRS/ProjectedCRS, or
+ * the source CRS of a CoordinateOperation.
  *
  * The returned object must be unreferenced with proj_obj_unref() after
  * use.
@@ -1372,6 +1374,10 @@ PJ_OBJ *proj_obj_get_source_crs(PJ_OBJ *obj) {
     if (boundCRS) {
         return PJ_OBJ::create(obj->ctx, boundCRS->baseCRS());
     }
+    auto derivedCRS = dynamic_cast<const DerivedCRS *>(ptr);
+    if (derivedCRS) {
+        return PJ_OBJ::create(obj->ctx, derivedCRS->baseCRS());
+    }
     auto co = dynamic_cast<const CoordinateOperation *>(ptr);
     if (co) {
         auto sourceCRS = co->sourceCRS();
@@ -4071,3 +4077,140 @@ double proj_coordoperation_get_accuracy(PJ_OBJ *coordoperation) {
     }
     return -1;
 }
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns the coordinate system of a SingleCRS.
+ *
+ * The returned object must be unreferenced with proj_obj_unref() after
+ * use.
+ * It should be used by at most one thread at a time.
+ *
+ * @param crs Objet of type SingleCRS (must not be NULL)
+ * @return Object that must be unreferenced with proj_obj_unref(), or NULL
+ * in case of error.
+ */
+PJ_OBJ *proj_obj_crs_get_coordinate_system(PJ_OBJ *crs) {
+    assert(crs);
+    auto l_crs = dynamic_cast<const SingleCRS *>(crs->obj.get());
+    if (!l_crs) {
+        proj_log_error(crs->ctx, __FUNCTION__, "Object is not a SingleCRS");
+        return nullptr;
+    }
+    return PJ_OBJ::create(crs->ctx, l_crs->coordinateSystem());
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns the type of the coordinate system.
+ *
+ * @param cs Objet of type CoordinateSystem (must not be NULL)
+ * @return type, or NULL in case of error.
+ */
+const char *proj_obj_cs_get_type(PJ_OBJ *cs) {
+    assert(cs);
+    auto l_cs = dynamic_cast<const CoordinateSystem *>(cs->obj.get());
+    if (!l_cs) {
+        proj_log_error(cs->ctx, __FUNCTION__,
+                       "Object is not a CoordinateSystem");
+        return nullptr;
+    }
+    if (dynamic_cast<const CartesianCS *>(l_cs)) {
+        return "Cartesian";
+    }
+    if (dynamic_cast<const EllipsoidalCS *>(l_cs)) {
+        return "Ellipsoidal";
+    }
+    if (dynamic_cast<const VerticalCS *>(l_cs)) {
+        return "Vertical";
+    }
+    if (dynamic_cast<const SphericalCS *>(l_cs)) {
+        return "Spherical";
+    }
+    if (dynamic_cast<const OrdinalCS *>(l_cs)) {
+        return "Ordinal";
+    }
+    if (dynamic_cast<const ParametricCS *>(l_cs)) {
+        return "Parametric";
+    }
+    if (dynamic_cast<const DateTimeTemporalCS *>(l_cs)) {
+        return "DateTimeTemporal";
+    }
+    if (dynamic_cast<const TemporalCountCS *>(l_cs)) {
+        return "TemporalCount";
+    }
+    if (dynamic_cast<const TemporalMeasureCS *>(l_cs)) {
+        return "TemporalMeasure";
+    }
+    return "unknown";
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns the number of axis of the coordinate system.
+ *
+ * @param cs Objet of type CoordinateSystem (must not be NULL)
+ * @return number of axis, or -1 in case of error.
+ */
+int proj_obj_cs_get_axis_count(PJ_OBJ *cs) {
+    assert(cs);
+    auto l_cs = dynamic_cast<const CoordinateSystem *>(cs->obj.get());
+    if (!l_cs) {
+        proj_log_error(cs->ctx, __FUNCTION__,
+                       "Object is not a CoordinateSystem");
+        return -1;
+    }
+    return static_cast<int>(l_cs->axisList().size());
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns information on an axis
+ *
+ * @param cs Objet of type CoordinateSystem (must not be NULL)
+ * @param index Index of the coordinate system (between 0 and
+ * proj_obj_cs_get_axis_count() - 1)
+ * @param pName Pointer to a string value to store the axis name. or NULL
+ * @param pAbbrev Pointer to a string value to store the axis abbreviation. or
+ * NULL
+ * @param pDirection Pointer to a string value to store the axis direction. or
+ * NULL
+ * @param pUnitConvFactor Pointer to a double value to store the axis
+ * unit conversion factor. or NULL
+ * @param pUnitName Pointer to a string value to store the axis
+ * unit name. or NULL
+ * @return TRUE in case of success
+ */
+int proj_obj_cs_get_axis_info(PJ_OBJ *cs, int index, const char **pName,
+                              const char **pAbbrev, const char **pDirection,
+                              double *pUnitConvFactor, const char **pUnitName) {
+    assert(cs);
+    auto l_cs = dynamic_cast<const CoordinateSystem *>(cs->obj.get());
+    if (!l_cs) {
+        proj_log_error(cs->ctx, __FUNCTION__,
+                       "Object is not a CoordinateSystem");
+        return false;
+    }
+    const auto &axisList = l_cs->axisList();
+    if (index < 0 || static_cast<size_t>(index) >= axisList.size()) {
+        proj_log_error(cs->ctx, __FUNCTION__, "Invalid index");
+        return false;
+    }
+    const auto &axis = axisList[index];
+    if (pName) {
+        *pName = axis->nameStr().c_str();
+    }
+    if (pAbbrev) {
+        *pAbbrev = axis->abbreviation().c_str();
+    }
+    if (pDirection) {
+        *pDirection = axis->direction().toString().c_str();
+    }
+    if (pUnitConvFactor) {
+        *pUnitConvFactor = axis->unit().conversionToSI();
+    }
+    if (pUnitName) {
+        *pUnitName = axis->unit().name().c_str();
+    }
+    return true;
+}
diff --git a/src/coordinateoperation.cpp b/src/coordinateoperation.cpp
index e0e02931..893b52d3 100644
--- a/src/coordinateoperation.cpp
+++ b/src/coordinateoperation.cpp
@@ -5471,6 +5471,12 @@ Transformation::~Transformation() = default;
 
 // ---------------------------------------------------------------------------
 
+Transformation::Transformation(const Transformation &other)
+    : CoordinateOperation(other), SingleOperation(other),
+      d(internal::make_unique<Private>(*other.d)) {}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Return the source crs::CRS of the transformation.
  *
  * @return the source CRS.
@@ -5492,6 +5498,17 @@ const crs::CRSNNPtr &Transformation::targetCRS() PROJ_CONST_DEFN {
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+TransformationNNPtr Transformation::shallowClone() const {
+    auto conv = Transformation::nn_make_shared<Transformation>(*this);
+    conv->assignSelf(conv);
+    conv->setCRSs(this, false);
+    return conv;
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 /** \brief Return the TOWGS84 parameters of the transformation.
  *
  * If this transformation uses Coordinate Frame Rotation, Position Vector
@@ -7595,6 +7612,8 @@ void Transformation::_exportToPROJString(
             "Transformation cannot be exported as a PROJ.4 string");
     }
 
+    formatter->setCoordinateOperationOptimizations(true);
+
     bool positionVectorConvention = true;
     bool sevenParamsTransform = false;
     bool threeParamsTransform = false;
@@ -10156,6 +10175,72 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
     const auto candidatesDstGeod(
         findCandidateGeodCRSForDatum(authFactory, geodDst->datum()));
 
+    auto createTransformations = [&](const crs::CRSNNPtr &candidateSrcGeod,
+                                     const crs::CRSNNPtr &candidateDstGeod,
+                                     const CoordinateOperationNNPtr &opFirst,
+                                     bool isNullFirst) {
+        const auto opsSecond =
+            createOperations(candidateSrcGeod, candidateDstGeod, context);
+        const auto opsThird =
+            createOperations(candidateDstGeod, targetCRS, context);
+        assert(!opsThird.empty());
+
+        for (auto &opSecond : opsSecond) {
+            // Check that it is not a transformation synthetized by
+            // ourselves
+            if (!hasIdentifiers(opSecond)) {
+                continue;
+            }
+            // And even if it is a referenced transformation, check that
+            // it is not a trivial one
+            auto so = dynamic_cast<const SingleOperation *>(opSecond.get());
+            if (so && isAxisOrderReversal(so->method()->getEPSGCode())) {
+                continue;
+            }
+
+            std::vector<CoordinateOperationNNPtr> subOps;
+            if (isNullFirst) {
+                opSecond->setCRSs(
+                    sourceCRS, NN_CHECK_ASSERT(opSecond->targetCRS()), nullptr);
+            } else {
+                subOps.emplace_back(opFirst);
+            }
+            if (isNullTransformation(opsThird[0]->nameStr())) {
+                opSecond->setCRSs(NN_CHECK_ASSERT(opSecond->sourceCRS()),
+                                  targetCRS, nullptr);
+                subOps.emplace_back(opSecond);
+            } else {
+                subOps.emplace_back(opSecond);
+                subOps.emplace_back(opsThird[0]);
+            }
+            res.emplace_back(ConcatenatedOperation::createComputeMetadata(
+                subOps, !allowEmptyIntersection));
+        }
+    };
+
+    // Start in priority with candidates that have exactly the same name as
+    // the sourcCRS and targetCRS. Typically for the case of init=IGNF:XXXX
+    for (const auto &candidateSrcGeod : candidatesSrcGeod) {
+        if (candidateSrcGeod->nameStr() == sourceCRS->nameStr()) {
+            for (const auto &candidateDstGeod : candidatesDstGeod) {
+                if (candidateDstGeod->nameStr() == targetCRS->nameStr()) {
+                    const auto opsFirst =
+                        createOperations(sourceCRS, candidateSrcGeod, context);
+                    assert(!opsFirst.empty());
+                    const bool isNullFirst =
+                        isNullTransformation(opsFirst[0]->nameStr());
+                    createTransformations(candidateSrcGeod, candidateDstGeod,
+                                          opsFirst[0], isNullFirst);
+                    if (!res.empty()) {
+                        return;
+                    }
+                    break;
+                }
+            }
+            break;
+        }
+    }
+
     for (const auto &candidateSrcGeod : candidatesSrcGeod) {
         const auto opsFirst =
             createOperations(sourceCRS, candidateSrcGeod, context);
@@ -10163,44 +10248,8 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
         const bool isNullFirst = isNullTransformation(opsFirst[0]->nameStr());
 
         for (const auto &candidateDstGeod : candidatesDstGeod) {
-            const auto opsSecond =
-                createOperations(candidateSrcGeod, candidateDstGeod, context);
-            const auto opsThird =
-                createOperations(candidateDstGeod, targetCRS, context);
-            assert(!opsThird.empty());
-
-            for (auto &opSecond : opsSecond) {
-                // Check that it is not a transformation synthetized by
-                // ourselves
-                if (!hasIdentifiers(opSecond)) {
-                    continue;
-                }
-                // And even if it is a referenced transformation, check that
-                // it is not a trivial one
-                auto so = dynamic_cast<const SingleOperation *>(opSecond.get());
-                if (so && isAxisOrderReversal(so->method()->getEPSGCode())) {
-                    continue;
-                }
-
-                std::vector<CoordinateOperationNNPtr> subOps;
-                if (isNullFirst) {
-                    opSecond->setCRSs(sourceCRS,
-                                      NN_CHECK_ASSERT(opSecond->targetCRS()),
-                                      nullptr);
-                } else {
-                    subOps.emplace_back(opsFirst[0]);
-                }
-                if (isNullTransformation(opsThird[0]->nameStr())) {
-                    opSecond->setCRSs(NN_CHECK_ASSERT(opSecond->sourceCRS()),
-                                      targetCRS, nullptr);
-                    subOps.emplace_back(opSecond);
-                } else {
-                    subOps.emplace_back(opSecond);
-                    subOps.emplace_back(opsThird[0]);
-                }
-                res.emplace_back(ConcatenatedOperation::createComputeMetadata(
-                    subOps, !allowEmptyIntersection));
-            }
+            createTransformations(candidateSrcGeod, candidateDstGeod,
+                                  opsFirst[0], isNullFirst);
         }
         if (!res.empty()) {
             return;
@@ -10261,6 +10310,56 @@ CoordinateOperationFactory::Private::createOperations(
     std::vector<CoordinateOperationNNPtr> res;
     const bool allowEmptyIntersection = true;
 
+    const auto &sourceProj4Ext = sourceCRS->getExtensionProj4();
+    const auto &targetProj4Ext = targetCRS->getExtensionProj4();
+    if (!sourceProj4Ext.empty() || !targetProj4Ext.empty()) {
+
+        auto sourceProjExportable =
+            dynamic_cast<io::IPROJStringExportable *>(sourceCRS.get());
+        auto targetProjExportable =
+            dynamic_cast<io::IPROJStringExportable *>(targetCRS.get());
+        if (!sourceProjExportable) {
+            throw InvalidOperation("Source CRS is not PROJ exportable");
+        }
+        if (!targetProjExportable) {
+            throw InvalidOperation("Target CRS is not PROJ exportable");
+        }
+        auto projFormatter = io::PROJStringFormatter::create(
+            io::PROJStringFormatter::Convention::PROJ_4);
+        projFormatter->startInversion();
+        sourceProjExportable->_exportToPROJString(projFormatter.get());
+
+        auto geogSrc =
+            dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
+        if (geogSrc) {
+            auto proj5Formatter = io::PROJStringFormatter::create(
+                io::PROJStringFormatter::Convention::PROJ_5);
+            geogSrc->addAngularUnitConvertAndAxisSwap(proj5Formatter.get());
+            projFormatter->ingestPROJString(proj5Formatter->toString());
+        }
+
+        projFormatter->stopInversion();
+
+        targetProjExportable->_exportToPROJString(projFormatter.get());
+
+        auto geogDst =
+            dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
+        if (geogDst) {
+            auto proj5Formatter = io::PROJStringFormatter::create(
+                io::PROJStringFormatter::Convention::PROJ_5);
+            geogDst->addAngularUnitConvertAndAxisSwap(proj5Formatter.get());
+            projFormatter->ingestPROJString(proj5Formatter->toString());
+        }
+
+        const auto PROJString = projFormatter->toString();
+        auto properties = util::PropertyMap().set(
+            common::IdentifiedObject::NAME_KEY,
+            buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()));
+        res.emplace_back(SingleOperation::createPROJBased(
+            properties, PROJString, sourceCRS, targetCRS, {}));
+        return res;
+    }
+
     auto geodSrc = dynamic_cast<const crs::GeodeticCRS *>(sourceCRS.get());
     auto geodDst = dynamic_cast<const crs::GeodeticCRS *>(targetCRS.get());
 
@@ -10317,7 +10416,9 @@ CoordinateOperationFactory::Private::createOperations(
             const auto &dstDatum = geodDst->datum();
             const bool dstHasDatumWithId =
                 dstDatum && !dstDatum->identifiers().empty();
-            if (srcHasDatumWithId && dstHasDatumWithId) {
+            if (srcHasDatumWithId && dstHasDatumWithId &&
+                !srcDatum->_isEquivalentTo(
+                    dstDatum.get(), util::IComparable::Criterion::EQUIVALENT)) {
                 createOperationsWithDatumPivot(res, sourceCRS, targetCRS,
                                                geodSrc, geodDst, context);
                 doFilterAndCheckPerfectOp = !res.empty();
@@ -10443,11 +10544,10 @@ CoordinateOperationFactory::Private::createOperations(
             dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
         auto geogCRSOfBaseOfBoundSrc =
             boundSrc->baseCRS()->extractGeographicCRS();
-        if (hubSrcGeog &&
+        if (hubSrcGeog && geogCRSOfBaseOfBoundSrc &&
             (hubSrcGeog->_isEquivalentTo(
                  geogDst, util::IComparable::Criterion::EQUIVALENT) ||
-             hubSrcGeog->is2DPartOf3D(NN_NO_CHECK(geogDst))) &&
-            geogCRSOfBaseOfBoundSrc) {
+             hubSrcGeog->is2DPartOf3D(NN_NO_CHECK(geogDst)))) {
             if (boundSrc->baseCRS() == geogCRSOfBaseOfBoundSrc) {
                 // Optimization to avoid creating a useless concatenated
                 // operation
@@ -10471,6 +10571,86 @@ CoordinateOperationFactory::Private::createOperations(
                     return res;
                 }
             }
+            // If the datum are equivalent, this is also fine
+        } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() &&
+                   geogDst->datum() &&
+                   hubSrcGeog->datum()->_isEquivalentTo(
+                       geogDst->datum().get(),
+                       util::IComparable::Criterion::EQUIVALENT)) {
+            auto opsFirst =
+                createOperations(boundSrc->baseCRS(),
+                                 NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
+            auto opsLast = createOperations(hubSrc, targetCRS, context);
+            if (!opsFirst.empty() && !opsLast.empty()) {
+                for (const auto &opFirst : opsFirst) {
+                    for (const auto &opLast : opsLast) {
+                        try {
+                            res.emplace_back(
+                                ConcatenatedOperation::createComputeMetadata(
+                                    {opFirst, boundSrc->transformation(),
+                                     opLast},
+                                    !allowEmptyIntersection));
+                        } catch (const InvalidOperationEmptyIntersection &) {
+                        }
+                    }
+                }
+                if (!res.empty()) {
+                    return res;
+                }
+            }
+            // Consider WGS 84 and NAD83 as equivalent in that context if the
+            // geogCRSOfBaseOfBoundSrc ellipsoid is Clarke66 (for NAD27)
+            // Case of "+proj=latlong +ellps=clrk66
+            // +nadgrids=ntv1_can.dat,conus"
+            // to "+proj=latlong +datum=NAD83"
+        } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() &&
+                   geogDst->datum() &&
+                   geogCRSOfBaseOfBoundSrc->ellipsoid()->_isEquivalentTo(
+                       datum::Ellipsoid::CLARKE_1866.get(),
+                       util::IComparable::Criterion::EQUIVALENT) &&
+                   hubSrcGeog->datum()->_isEquivalentTo(
+                       datum::GeodeticReferenceFrame::EPSG_6326.get(),
+                       util::IComparable::Criterion::EQUIVALENT) &&
+                   geogDst->datum()->_isEquivalentTo(
+                       datum::GeodeticReferenceFrame::EPSG_6269.get(),
+                       util::IComparable::Criterion::EQUIVALENT)) {
+            auto nnGeogCRSOfBaseOfBoundSrc =
+                NN_NO_CHECK(geogCRSOfBaseOfBoundSrc);
+            if (boundSrc->baseCRS()->_isEquivalentTo(
+                    nnGeogCRSOfBaseOfBoundSrc.get(),
+                    util::IComparable::Criterion::EQUIVALENT)) {
+                auto transf = boundSrc->transformation()->shallowClone();
+                transf->setProperties(util::PropertyMap().set(
+                    common::IdentifiedObject::NAME_KEY,
+                    buildTransfName(boundSrc->baseCRS()->nameStr(),
+                                    targetCRS->nameStr())));
+                transf->setCRSs(boundSrc->baseCRS(), targetCRS, nullptr);
+                res.emplace_back(transf);
+                return res;
+            } else {
+                auto opsFirst = createOperations(
+                    boundSrc->baseCRS(), nnGeogCRSOfBaseOfBoundSrc, context);
+                auto transf = boundSrc->transformation()->shallowClone();
+                transf->setProperties(util::PropertyMap().set(
+                    common::IdentifiedObject::NAME_KEY,
+                    buildTransfName(nnGeogCRSOfBaseOfBoundSrc->nameStr(),
+                                    targetCRS->nameStr())));
+                transf->setCRSs(nnGeogCRSOfBaseOfBoundSrc, targetCRS, nullptr);
+                if (!opsFirst.empty()) {
+                    for (const auto &opFirst : opsFirst) {
+                        try {
+                            res.emplace_back(
+                                ConcatenatedOperation::createComputeMetadata(
+                                    {opFirst, transf},
+                                    !allowEmptyIntersection));
+                        } catch (const InvalidOperationEmptyIntersection &) {
+                        }
+                    }
+                    if (!res.empty()) {
+                        return res;
+                    }
+                }
+            }
         }
 
         if (hubSrcGeog &&
diff --git a/src/crs.cpp b/src/crs.cpp
index a204a037..eec7a926 100644
--- a/src/crs.cpp
+++ b/src/crs.cpp
@@ -170,6 +170,14 @@ const GeodeticCRS *CRS::extractGeodeticCRSRaw() const {
 
 // ---------------------------------------------------------------------------
 
+//! @cond Doxygen_Suppress
+const std::string &CRS::getExtensionProj4() const noexcept {
+    return d->extensionProj4_;
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
 /** \brief Return the GeographicCRS of the CRS.
  *
  * Returns the GeographicCRS contained in a CRS. This works currently with
diff --git a/src/cs2cs.cpp b/src/cs2cs.cpp
index 6f302ae3..6f4c4a55 100644
--- a/src/cs2cs.cpp
+++ b/src/cs2cs.cpp
@@ -26,6 +26,8 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
+#define FROM_PROJ_CPP
+
 #include <ctype.h>
 #include <locale.h>
 #include <math.h>
@@ -33,6 +35,11 @@
 #include <stdlib.h>
 #include <string.h>
 
+#include <cassert>
+#include <string>
+
+#include <proj/internal/internal.hpp>
+
 // PROJ include order is sensitive
 // clang-format off
 #include "proj.h"
@@ -41,9 +48,15 @@
 // clang-format on
 
 #define MAX_LINE 1000
-#define MAX_PARGS 100
 
-static projPJ fromProj, toProj;
+static PJ *transformation = nullptr;
+
+static bool srcIsGeog = false;
+static double srcToRadians = 0.0;
+
+static bool destIsGeog = false;
+static double destToRadians = 0.0;
+static bool destIsLatLong = false;
 
 static int reversein = 0, /* != 0 reverse input arguments */
     reverseout = 0,       /* != 0 reverse output arguments */
@@ -116,19 +129,45 @@ static void process(FILE *fid)
         }
 
         if (data.u != HUGE_VAL) {
-            if (pj_transform(fromProj, toProj, 1, 0, &(data.u), &(data.v),
-                             &z) != 0) {
-                data.u = HUGE_VAL;
-                data.v = HUGE_VAL;
-                emess(-3, "pj_transform(): %s", pj_strerrno(pj_errno));
+
+            if (srcIsGeog) {
+                /* dmstor gives values to radians. Convert now to the SRS unit
+                 */
+                data.u /= srcToRadians;
+                data.v /= srcToRadians;
             }
+
+            PJ_COORD coord;
+            coord.xyzt.x = data.u;
+            coord.xyzt.y = data.v;
+            coord.xyzt.z = z;
+            coord.xyzt.t = HUGE_VAL;
+            coord = proj_trans(transformation, PJ_FWD, coord);
+            data.u = coord.xyz.x;
+            data.v = coord.xyz.y;
+            z = coord.xyz.z;
         }
 
         if (data.u == HUGE_VAL) /* error output */
             fputs(oterr, stdout);
 
-        else if (pj_is_latlong(toProj) && !oform) { /*ascii DMS output */
-            if (reverseout) {
+        else if (destIsGeog && !oform) { /*ascii DMS output */
+
+            // rtodms() expect radians: convert from the output SRS unit
+            data.u *= destToRadians;
+            data.v *= destToRadians;
+
+            if (destIsLatLong) {
+                if (reverseout) {
+                    fputs(rtodms(pline, data.v, 'E', 'W'), stdout);
+                    putchar('\t');
+                    fputs(rtodms(pline, data.u, 'N', 'S'), stdout);
+                } else {
+                    fputs(rtodms(pline, data.u, 'N', 'S'), stdout);
+                    putchar('\t');
+                    fputs(rtodms(pline, data.v, 'E', 'W'), stdout);
+                }
+            } else if (reverseout) {
                 fputs(rtodms(pline, data.v, 'N', 'S'), stdout);
                 putchar('\t');
                 fputs(rtodms(pline, data.u, 'E', 'W'), stdout);
@@ -139,9 +178,9 @@ static void process(FILE *fid)
             }
 
         } else { /* x-y or decimal degree ascii output */
-            if (proj_angular_output(toProj, PJ_FWD)) {
-                data.v *= RAD_TO_DEG;
-                data.u *= RAD_TO_DEG;
+            if (destIsGeog) {
+                data.v *= destToRadians * RAD_TO_DEG;
+                data.u *= destToRadians * RAD_TO_DEG;
             }
             if (reverseout) {
                 printf(oform, data.v);
@@ -167,17 +206,116 @@ static void process(FILE *fid)
 }
 
 /************************************************************************/
+/*                          instanciate_crs()                           */
+/************************************************************************/
+
+static PJ_OBJ *instanciate_crs(const std::string &definition,
+                               const char *const *optionsImportCRS,
+                               bool &isGeog, double &toRadians,
+                               bool &isLatFirst) {
+    PJ_OBJ *crs = proj_obj_create_from_user_input(nullptr, definition.c_str(),
+                                                  optionsImportCRS);
+    if (!crs) {
+        return nullptr;
+    }
+
+    isGeog = false;
+    toRadians = 0.0;
+    isLatFirst = false;
+
+    auto type = proj_obj_get_type(crs);
+    if (type == PJ_OBJ_TYPE_BOUND_CRS) {
+        auto base = proj_obj_get_source_crs(crs);
+        proj_obj_unref(crs);
+        crs = base;
+        type = proj_obj_get_type(crs);
+    }
+    if (type == PJ_OBJ_TYPE_GEOGRAPHIC_2D_CRS ||
+        type == PJ_OBJ_TYPE_GEOGRAPHIC_3D_CRS) {
+        auto cs = proj_obj_crs_get_coordinate_system(crs);
+        assert(cs);
+
+        isGeog = true;
+        const char *axisName = "";
+        proj_obj_cs_get_axis_info(cs, 0,
+                                  &axisName, // name,
+                                  nullptr,   // abbrev
+                                  nullptr,   // direction
+                                  &toRadians,
+                                  nullptr // unit name
+                                  );
+        isLatFirst =
+            NS_PROJ::internal::ci_find(std::string(axisName), "latitude") !=
+            std::string::npos;
+
+        proj_obj_unref(cs);
+    }
+
+    return crs;
+}
+
+/************************************************************************/
+/*               get_geog_crs_proj_string_from_proj_crs()               */
+/************************************************************************/
+
+static std::string get_geog_crs_proj_string_from_proj_crs(PJ_OBJ *src,
+                                                          double &toRadians,
+                                                          bool &isLatFirst) {
+    auto srcType = proj_obj_get_type(src);
+    if (srcType == PJ_OBJ_TYPE_BOUND_CRS) {
+        auto base = proj_obj_get_source_crs(src);
+        assert(base);
+        proj_obj_unref(src);
+        src = base;
+        srcType = proj_obj_get_type(src);
+    }
+    if (srcType != PJ_OBJ_TYPE_PROJECTED_CRS) {
+        return std::string();
+    }
+
+    auto base = proj_obj_get_source_crs(src);
+    assert(base);
+    auto baseType = proj_obj_get_type(base);
+    if (baseType != PJ_OBJ_TYPE_GEOGRAPHIC_2D_CRS &&
+        baseType != PJ_OBJ_TYPE_GEOGRAPHIC_3D_CRS) {
+        proj_obj_unref(base);
+        return std::string();
+    }
+
+    auto cs = proj_obj_crs_get_coordinate_system(base);
+    assert(cs);
+
+    const char *axisName = "";
+    proj_obj_cs_get_axis_info(cs, 0,
+                              &axisName, // name,
+                              nullptr,   // abbrev
+                              nullptr,   // direction
+                              &toRadians,
+                              nullptr // unit name
+                              );
+    isLatFirst = NS_PROJ::internal::ci_find(std::string(axisName),
+                                            "latitude") != std::string::npos;
+
+    proj_obj_unref(cs);
+
+    auto retCStr = proj_obj_as_proj_string(base, PJ_PROJ_5, nullptr);
+    std::string ret(retCStr ? retCStr : "");
+    proj_obj_unref(base);
+    return ret;
+}
+
+/************************************************************************/
 /*                                main()                                */
 /************************************************************************/
 
 int main(int argc, char **argv) {
     char *arg;
     char **eargv = argv;
-    char *from_argv[MAX_PARGS];
-    char *to_argv[MAX_PARGS];
+    std::string fromStr;
+    std::string toStr;
     FILE *fid;
-    int from_argc = 0, to_argc = 0, eargc = 0, mon = 0;
-    int have_to_flag = 0, inverse = 0, i;
+    int eargc = 0, mon = 0;
+    int have_to_flag = 0, inverse = 0;
     int use_env_locale = 0;
 
     /* This is just to check that pj_init() is locale-safe */
@@ -185,6 +323,11 @@ int main(int argc, char **argv) {
     if (getenv("PROJ_USE_ENV_LOCALE") != nullptr)
         use_env_locale = 1;
 
+    /* Enable compatibility mode for init=epsg:XXXX by default */
+    if (getenv("PROJ_USE_PROJ4_INIT_RULES") == nullptr) {
+        proj_context_use_proj4_init_rules(nullptr, true);
+    }
+
     if ((emess_dat.Prog_name = strrchr(*argv, DIR_CHAR)) != nullptr)
         ++emess_dat.Prog_name;
     else
@@ -194,9 +337,27 @@ int main(int argc, char **argv) {
         (void)fprintf(stderr, usage, pj_get_release(), emess_dat.Prog_name);
         exit(0);
     }
+
+    // First pass to check if we have "cs2cs [-bla]* <SRC> <DEST>" syntax
+    int countNonOptionArg = 0;
+    for (int i = 1; i < argc; i++) {
+        if (argv[i][0] == '-') {
+            if (argv[i][1] == '\0') {
+                countNonOptionArg++;
+            }
+        } else {
+            if (strcmp(argv[i], "+to") == 0) {
+                countNonOptionArg = -1;
+                break;
+            }
+            countNonOptionArg++;
+        }
+    }
+    const bool isSrcDestSyntax = (countNonOptionArg == 2);
+
     /* process run line arguments */
     while (--argc > 0) { /* collect run line arguments */
-        if (**++argv == '-')
+        if (**++argv == '-') {
             for (arg = *argv;;) {
                 switch (*++arg) {
                 case '\0': /* position of "stdin" */
@@ -325,21 +486,28 @@ int main(int argc, char **argv) {
                 }
                 break;
             }
-        else if (strcmp(*argv, "+to") == 0) {
+        } else if (isSrcDestSyntax) {
+            if (fromStr.empty())
+                fromStr = *argv;
+            else
+                toStr = *argv;
+        } else if (strcmp(*argv, "+to") == 0) {
             have_to_flag = 1;
 
         } else if (**argv == '+') { /* + argument */
             if (have_to_flag) {
-                if (to_argc < MAX_PARGS)
-                    to_argv[to_argc++] = *argv + 1;
-                else
-                    emess(1, "overflowed + argument table");
+                if (!toStr.empty())
+                    toStr += ' ';
+                toStr += *argv;
             } else {
-                if (from_argc < MAX_PARGS)
-                    from_argv[from_argc++] = *argv + 1;
-                else
-                    emess(1, "overflowed + argument table");
+                if (!fromStr.empty())
+                    fromStr += ' ';
+                fromStr += *argv;
             }
+        } else if (!have_to_flag) {
+            fromStr = *argv;
+        } else if (toStr.empty()) {
+            toStr = *argv;
         } else /* assumed to be input file name(s) */
             eargv[eargc++] = *argv;
     }
@@ -351,17 +519,7 @@ int main(int argc, char **argv) {
      * coordinate systems.
      */
     if (inverse) {
-        int argcount;
-
-        for (i = 0; i < MAX_PARGS; i++) {
-            arg = from_argv[i];
-            from_argv[i] = to_argv[i];
-            to_argv[i] = arg;
-        }
-
-        argcount = from_argc;
-        from_argc = to_argc;
-        to_argc = argcount;
+        std::swap(fromStr, toStr);
     }
 
     if (use_env_locale) {
@@ -369,48 +527,64 @@ int main(int argc, char **argv) {
         setlocale(LC_ALL, "");
     }
 
-    if (from_argc == 0 && to_argc != 0) {
-        /* we will generate the from proj as the latlong of the +to in a bit */
-    } else if (!(fromProj = pj_init(from_argc, from_argv))) {
-        printf("Using from definition: ");
-        for (i = 0; i < from_argc; i++)
-            printf("%s ", from_argv[i]);
-        printf("\n");
-
-        emess(3, "projection initialization failure\ncause: %s",
-              pj_strerrno(pj_errno));
+    if (fromStr.empty() && toStr.empty()) {
+        emess(3, "missing source and target coordinate systems");
     }
 
-    if (to_argc == 0) {
-        if (!(toProj = pj_latlong_from_proj(fromProj))) {
-            printf("Using to definition: ");
-            for (i = 0; i < to_argc; i++)
-                printf("%s ", to_argv[i]);
-            printf("\n");
+    const char *const optionsProj4Mode[] = {"USE_PROJ4_INIT_RULES=YES",
+                                            nullptr};
+    const char *const *optionsImportCRS =
+        proj_context_get_use_proj4_init_rules(nullptr) ? optionsProj4Mode
+                                                       : nullptr;
+
+    PJ_OBJ *src = nullptr;
+    if (!fromStr.empty()) {
+        bool ignored;
+        src = instanciate_crs(fromStr, optionsImportCRS, srcIsGeog,
+                              srcToRadians, ignored);
+        if (!src) {
+            emess(3, "cannot instanciate source coordinate system");
+        }
+    }
 
-            emess(3, "projection initialization failure\ncause: %s",
-                  pj_strerrno(pj_errno));
+    PJ_OBJ *dst = nullptr;
+    if (!toStr.empty()) {
+        dst = instanciate_crs(toStr, optionsImportCRS, destIsGeog,
+                              destToRadians, destIsLatLong);
+        if (!dst) {
+            emess(3, "cannot instanciate target coordinate system");
         }
-    } else if (!(toProj = pj_init(to_argc, to_argv))) {
-        printf("Using to definition: ");
-        for (i = 0; i < to_argc; i++)
-            printf("%s ", to_argv[i]);
-        printf("\n");
+    }
 
-        emess(3, "projection initialization failure\ncause: %s",
-              pj_strerrno(pj_errno));
+    if (toStr.empty()) {
+        assert(src);
+        toStr = get_geog_crs_proj_string_from_proj_crs(src, destToRadians,
+                                                       destIsLatLong);
+        if (toStr.empty()) {
+            emess(3,
+                  "missing target CRS and source CRS is not a projected CRS");
+        }
+        destIsGeog = true;
+    } else if (fromStr.empty()) {
+        assert(dst);
+        bool ignored;
+        fromStr =
+            get_geog_crs_proj_string_from_proj_crs(dst, srcToRadians, ignored);
+        if (fromStr.empty()) {
+            emess(3,
+                  "missing source CRS and target CRS is not a projected CRS");
+        }
+        srcIsGeog = true;
     }
 
-    if (from_argc == 0 && toProj != nullptr) {
-        if (!(fromProj = pj_latlong_from_proj(toProj))) {
-            printf("Using to definition: ");
-            for (i = 0; i < to_argc; i++)
-                printf("%s ", to_argv[i]);
-            printf("\n");
+    proj_obj_unref(src);
+    proj_obj_unref(dst);
 
-            emess(3, "projection initialization failure\ncause: %s",
-                  pj_strerrno(pj_errno));
-        }
+    transformation = proj_create_crs_to_crs(nullptr, fromStr.c_str(),
+                                            toStr.c_str(), nullptr);
+    if (!transformation) {
+        emess(3, "cannot initialize transformation\ncause: %s",
+              pj_strerrno(pj_errno));
     }
 
     if (use_env_locale) {
@@ -420,19 +594,19 @@ int main(int argc, char **argv) {
 
     if (mon) {
         printf("%c ---- From Coordinate System ----\n", tag);
-        pj_pr_list(fromProj);
+        printf("%s\n", fromStr.c_str());
         printf("%c ---- To Coordinate System ----\n", tag);
-        pj_pr_list(toProj);
+        printf("%s\n", toStr.c_str());
     }
 
     /* set input formatting control */
-    if (!fromProj->is_latlong)
+    if (!srcIsGeog)
         informat = strtod;
     else {
         informat = dmstor;
     }
 
-    if (!toProj->is_latlong && !oform)
+    if (!destIsGeog && !oform)
         oform = "%.2f";
 
     /* process input file list */
@@ -454,8 +628,7 @@ int main(int argc, char **argv) {
         emess_dat.File_name = nullptr;
     }
 
-    pj_free(fromProj);
-    pj_free(toProj);
+    proj_destroy(transformation);
 
     pj_deallocate_grids();
 
diff --git a/src/io.cpp b/src/io.cpp
index 0d220e13..15df312b 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -4422,6 +4422,7 @@ struct PROJStringFormatter::Private {
     bool useETMercForTMerc_ = false;
     bool useETMercForTMercSet_ = false;
     bool addNoDefs_ = true;
+    bool coordOperationOptimizations_ = false;
 
     std::string result_{};
 
@@ -4504,6 +4505,12 @@ const std::string &PROJStringFormatter::toString() const {
               step.paramValues[6].equals("s", "0") &&
               step.paramValues[7].keyEquals("convention")))) {
             iter = d->steps_.erase(iter);
+        } else if (d->coordOperationOptimizations_ &&
+                   step.name == "unitconvert" && paramCount == 2 &&
+                   step.paramValues[0].keyEquals("xy_in") &&
+                   step.paramValues[1].keyEquals("xy_out") &&
+                   step.paramValues[0].value == step.paramValues[1].value) {
+            iter = d->steps_.erase(iter);
         } else {
             ++iter;
         }
@@ -4925,6 +4932,12 @@ bool PROJStringFormatter::getUseETMercForTMerc(bool &settingSetOut) const {
 
 // ---------------------------------------------------------------------------
 
+void PROJStringFormatter::setCoordinateOperationOptimizations(bool enable) {
+    d->coordOperationOptimizations_ = enable;
+}
+
+// ---------------------------------------------------------------------------
+
 void PROJStringFormatter::Private::appendToResult(const char *str) {
     if (!result_.empty()) {
         result_ += " ";
@@ -4984,11 +4997,13 @@ PROJStringSyntaxParser(const std::string &projString, std::vector<Step> &steps,
             }
             prevWasStep = false;
         } else if (starts_with(word, "proj=")) {
+            auto stepName = word.substr(strlen("proj="));
             if (prevWasInit) {
-                throw ParsingException("+init= found at unexpected place");
+                steps.back() = Step();
+                prevWasInit = false;
+            } else {
+                steps.push_back(Step());
             }
-            auto stepName = word.substr(strlen("proj="));
-            steps.push_back(Step());
             steps.back().name = stepName;
             steps.back().inverted = inverted;
             prevWasStep = false;
@@ -5041,7 +5056,7 @@ void PROJStringFormatter::ingestPROJString(
     std::string vto_meter;
     PROJStringSyntaxParser(str, steps, d->globalParamValues_, title, vunits,
                            vto_meter);
-    d->steps_.insert(d->steps_.begin(), steps.begin(), steps.end());
+    d->steps_.insert(d->steps_.end(), steps.begin(), steps.end());
 }
 
 // ---------------------------------------------------------------------------
@@ -5633,11 +5648,12 @@ PROJStringParser::Private::buildPrimeMeridian(const Step &step) {
     PrimeMeridianNNPtr pm = PrimeMeridian::GREENWICH;
     const auto &pmStr = getParamValue(step, "pm");
     if (!pmStr.empty()) {
-        try {
-            double pmValue = c_locale_stod(pmStr);
+        char *end;
+        double pmValue = dmstor(pmStr.c_str(), &end) * RAD_TO_DEG;
+        if (pmValue != HUGE_VAL && *end == '\0') {
             pm = PrimeMeridian::create(createMapWithUnknownName(),
                                        Angle(pmValue));
-        } catch (const std::invalid_argument &) {
+        } else {
             bool found = false;
             if (pmStr == "paris") {
                 found = true;
@@ -5650,9 +5666,8 @@ PROJStringParser::Private::buildPrimeMeridian(const Step &step) {
                     found = true;
                     std::string name = static_cast<char>(::toupper(pmStr[0])) +
                                        pmStr.substr(1);
-                    double pmValue =
-                        dmstor(proj_prime_meridians[i].defn, nullptr) *
-                        RAD_TO_DEG;
+                    pmValue = dmstor(proj_prime_meridians[i].defn, nullptr) *
+                              RAD_TO_DEG;
                     pm = PrimeMeridian::create(
                         PropertyMap().set(IdentifiedObject::NAME_KEY, name),
                         Angle(pmValue));
@@ -5681,12 +5696,20 @@ PROJStringParser::Private::buildDatum(const Step &step,
 
     const auto &ellpsStr = getParamValue(step, "ellps");
     const auto &datumStr = getParamValue(step, "datum");
+    const auto &RStr = getParamValue(step, "R");
     const auto &aStr = getParamValue(step, "a");
     const auto &bStr = getParamValue(step, "b");
     const auto &rfStr = getParamValue(step, "rf");
     const auto &fStr = getParamValue(step, "f");
-    const auto &RStr = getParamValue(step, "R");
+    const auto &esStr = getParamValue(step, "es");
+    const auto &eStr = getParamValue(step, "e");
+    double a = -1.0;
+    double b = -1.0;
+    double rf = -1.0;
     const util::optional<std::string> optionalEmptyString{};
+    const bool numericParamPresent =
+        !RStr.empty() || !aStr.empty() || !bStr.empty() || !rfStr.empty() ||
+        !fStr.empty() || !esStr.empty() || !eStr.empty();
 
     PrimeMeridianNNPtr pm(buildPrimeMeridian(step));
     PropertyMap grfMap;
@@ -5709,104 +5732,151 @@ PROJStringParser::Private::buildDatum(const Step &step,
             }
         };
 
+    // R take precedence
+    if (!RStr.empty()) {
+        double R;
+        try {
+            R = c_locale_stod(RStr);
+        } catch (const std::invalid_argument &) {
+            throw ParsingException("Invalid R value");
+        }
+        auto ellipsoid = Ellipsoid::createSphere(createMapWithUnknownName(),
+                                                 Length(R), guessBodyName(R));
+        return GeodeticReferenceFrame::create(
+            grfMap.set(IdentifiedObject::NAME_KEY,
+                       title.empty() ? "unknown" : title.c_str()),
+            ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
+    }
+
     if (!datumStr.empty()) {
-        if (datumStr == "WGS84") {
-            return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6326);
-        } else if (datumStr == "NAD83") {
-            return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6269);
-        } else if (datumStr == "NAD27") {
-            return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6267);
-        } else {
+        auto l_datum = [&datumStr, &overridePmIfNeeded, &grfMap,
+                        &optionalEmptyString, &pm]() {
+            if (datumStr == "WGS84") {
+                return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6326);
+            } else if (datumStr == "NAD83") {
+                return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6269);
+            } else if (datumStr == "NAD27") {
+                return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6267);
+            } else {
 
-            for (const auto &datumDesc : datumDescs) {
-                if (datumStr == datumDesc.projName) {
-                    (void)datumDesc.gcsName; // to please cppcheck
-                    (void)datumDesc.gcsCode; // to please cppcheck
-                    auto ellipsoid = Ellipsoid::createFlattenedSphere(
-                        grfMap
-                            .set(IdentifiedObject::NAME_KEY,
-                                 datumDesc.ellipsoidName)
-                            .set(Identifier::CODESPACE_KEY, Identifier::EPSG)
-                            .set(Identifier::CODE_KEY, datumDesc.ellipsoidCode),
-                        Length(datumDesc.a), Scale(datumDesc.rf));
-                    return GeodeticReferenceFrame::create(
-                        grfMap
-                            .set(IdentifiedObject::NAME_KEY,
-                                 datumDesc.datumName)
-                            .set(Identifier::CODESPACE_KEY, Identifier::EPSG)
-                            .set(Identifier::CODE_KEY, datumDesc.datumCode),
-                        ellipsoid, optionalEmptyString, pm);
+                for (const auto &datumDesc : datumDescs) {
+                    if (datumStr == datumDesc.projName) {
+                        (void)datumDesc.gcsName; // to please cppcheck
+                        (void)datumDesc.gcsCode; // to please cppcheck
+                        auto ellipsoid = Ellipsoid::createFlattenedSphere(
+                            grfMap
+                                .set(IdentifiedObject::NAME_KEY,
+                                     datumDesc.ellipsoidName)
+                                .set(Identifier::CODESPACE_KEY,
+                                     Identifier::EPSG)
+                                .set(Identifier::CODE_KEY,
+                                     datumDesc.ellipsoidCode),
+                            Length(datumDesc.a), Scale(datumDesc.rf));
+                        return GeodeticReferenceFrame::create(
+                            grfMap
+                                .set(IdentifiedObject::NAME_KEY,
+                                     datumDesc.datumName)
+                                .set(Identifier::CODESPACE_KEY,
+                                     Identifier::EPSG)
+                                .set(Identifier::CODE_KEY, datumDesc.datumCode),
+                            ellipsoid, optionalEmptyString, pm);
+                    }
                 }
             }
+            throw ParsingException("unknown datum " + datumStr);
+        }();
+        if (!numericParamPresent) {
+            return l_datum;
         }
-        throw ParsingException("unknown datum " + datumStr);
+        a = l_datum->ellipsoid()->semiMajorAxis().getSIValue();
+        rf = l_datum->ellipsoid()->computedInverseFlattening();
     }
 
     else if (!ellpsStr.empty()) {
-        if (ellpsStr == "WGS84") {
-            return GeodeticReferenceFrame::create(
-                grfMap.set(IdentifiedObject::NAME_KEY,
-                           title.empty() ? "Unknown based on WGS84 ellipsoid"
-                                         : title.c_str()),
-                Ellipsoid::WGS84, optionalEmptyString, pm);
-        } else if (ellpsStr == "GRS80") {
-            return GeodeticReferenceFrame::create(
-                grfMap.set(IdentifiedObject::NAME_KEY,
-                           title.empty() ? "Unknown based on GRS80 ellipsoid"
-                                         : title.c_str()),
-                Ellipsoid::GRS1980, optionalEmptyString, pm);
-        } else {
-            auto proj_ellps = proj_list_ellps();
-            for (int i = 0; proj_ellps[i].id != nullptr; i++) {
-                if (ellpsStr == proj_ellps[i].id) {
-                    assert(strncmp(proj_ellps[i].major, "a=", 2) == 0);
-                    const double a_iter =
-                        c_locale_stod(proj_ellps[i].major + 2);
-                    EllipsoidPtr ellipsoid;
-                    PropertyMap ellpsMap;
-                    if (strncmp(proj_ellps[i].ell, "b=", 2) == 0) {
-                        const double b_iter =
-                            c_locale_stod(proj_ellps[i].ell + 2);
-                        ellipsoid = Ellipsoid::createTwoAxis(
-                                        ellpsMap.set(IdentifiedObject::NAME_KEY,
-                                                     proj_ellps[i].name),
-                                        Length(a_iter), Length(b_iter))
-                                        .as_nullable();
-                    } else {
-                        assert(strncmp(proj_ellps[i].ell, "rf=", 3) == 0);
-                        const double rf_iter =
-                            c_locale_stod(proj_ellps[i].ell + 3);
-                        ellipsoid = Ellipsoid::createFlattenedSphere(
-                                        ellpsMap.set(IdentifiedObject::NAME_KEY,
-                                                     proj_ellps[i].name),
-                                        Length(a_iter), Scale(rf_iter))
-                                        .as_nullable();
+        auto l_datum = [&ellpsStr, &title, &grfMap, &optionalEmptyString,
+                        &pm]() {
+            if (ellpsStr == "WGS84") {
+                return GeodeticReferenceFrame::create(
+                    grfMap.set(IdentifiedObject::NAME_KEY,
+                               title.empty()
+                                   ? "Unknown based on WGS84 ellipsoid"
+                                   : title.c_str()),
+                    Ellipsoid::WGS84, optionalEmptyString, pm);
+            } else if (ellpsStr == "GRS80") {
+                return GeodeticReferenceFrame::create(
+                    grfMap.set(IdentifiedObject::NAME_KEY,
+                               title.empty()
+                                   ? "Unknown based on GRS80 ellipsoid"
+                                   : title.c_str()),
+                    Ellipsoid::GRS1980, optionalEmptyString, pm);
+            } else {
+                auto proj_ellps = proj_list_ellps();
+                for (int i = 0; proj_ellps[i].id != nullptr; i++) {
+                    if (ellpsStr == proj_ellps[i].id) {
+                        assert(strncmp(proj_ellps[i].major, "a=", 2) == 0);
+                        const double a_iter =
+                            c_locale_stod(proj_ellps[i].major + 2);
+                        EllipsoidPtr ellipsoid;
+                        PropertyMap ellpsMap;
+                        if (strncmp(proj_ellps[i].ell, "b=", 2) == 0) {
+                            const double b_iter =
+                                c_locale_stod(proj_ellps[i].ell + 2);
+                            ellipsoid =
+                                Ellipsoid::createTwoAxis(
+                                    ellpsMap.set(IdentifiedObject::NAME_KEY,
+                                                 proj_ellps[i].name),
+                                    Length(a_iter), Length(b_iter))
+                                    .as_nullable();
+                        } else {
+                            assert(strncmp(proj_ellps[i].ell, "rf=", 3) == 0);
+                            const double rf_iter =
+                                c_locale_stod(proj_ellps[i].ell + 3);
+                            ellipsoid =
+                                Ellipsoid::createFlattenedSphere(
+                                    ellpsMap.set(IdentifiedObject::NAME_KEY,
+                                                 proj_ellps[i].name),
+                                    Length(a_iter), Scale(rf_iter))
+                                    .as_nullable();
+                        }
+                        return GeodeticReferenceFrame::create(
+                            grfMap.set(IdentifiedObject::NAME_KEY,
+                                       title.empty()
+                                           ? std::string("Unknown based on ") +
+                                                 proj_ellps[i].name +
+                                                 " ellipsoid"
+                                           : title),
+                            NN_NO_CHECK(ellipsoid), optionalEmptyString, pm);
                     }
-                    return GeodeticReferenceFrame::create(
-                        grfMap.set(IdentifiedObject::NAME_KEY,
-                                   title.empty()
-                                       ? std::string("Unknown based on ") +
-                                             proj_ellps[i].name + " ellipsoid"
-                                       : title),
-                        NN_NO_CHECK(ellipsoid), optionalEmptyString, pm);
                 }
+                throw ParsingException("unknown ellipsoid " + ellpsStr);
             }
-            throw ParsingException("unknown ellipsoid " + ellpsStr);
+        }();
+        if (!numericParamPresent) {
+            return l_datum;
+        }
+        a = l_datum->ellipsoid()->semiMajorAxis().getSIValue();
+        if (l_datum->ellipsoid()->semiMinorAxis().has_value()) {
+            b = l_datum->ellipsoid()->semiMinorAxis()->getSIValue();
+        } else {
+            rf = l_datum->ellipsoid()->computedInverseFlattening();
         }
     }
 
-    else if (!aStr.empty() && !bStr.empty()) {
-        double a;
+    if (!aStr.empty()) {
         try {
             a = c_locale_stod(aStr);
         } catch (const std::invalid_argument &) {
             throw ParsingException("Invalid a value");
         }
-        double b;
-        try {
-            b = c_locale_stod(bStr);
-        } catch (const std::invalid_argument &) {
-            throw ParsingException("Invalid b value");
+    }
+
+    if (a > 0 && (b > 0 || !bStr.empty())) {
+        if (!bStr.empty()) {
+            try {
+                b = c_locale_stod(bStr);
+            } catch (const std::invalid_argument &) {
+                throw ParsingException("Invalid b value");
+            }
         }
         auto ellipsoid =
             Ellipsoid::createTwoAxis(createMapWithUnknownName(), Length(a),
@@ -5818,18 +5888,13 @@ PROJStringParser::Private::buildDatum(const Step &step,
             ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
     }
 
-    else if (!aStr.empty() && !rfStr.empty()) {
-        double a;
-        try {
-            a = c_locale_stod(aStr);
-        } catch (const std::invalid_argument &) {
-            throw ParsingException("Invalid a value");
-        }
-        double rf;
-        try {
-            rf = c_locale_stod(rfStr);
-        } catch (const std::invalid_argument &) {
-            throw ParsingException("Invalid rf value");
+    else if (a > 0 && (rf >= 0 || !rfStr.empty())) {
+        if (!rfStr.empty()) {
+            try {
+                rf = c_locale_stod(rfStr);
+            } catch (const std::invalid_argument &) {
+                throw ParsingException("Invalid rf value");
+            }
         }
         auto ellipsoid = Ellipsoid::createFlattenedSphere(
                              createMapWithUnknownName(), Length(a), Scale(rf),
@@ -5841,13 +5906,7 @@ PROJStringParser::Private::buildDatum(const Step &step,
             ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
     }
 
-    else if (!aStr.empty() && !fStr.empty()) {
-        double a;
-        try {
-            a = c_locale_stod(aStr);
-        } catch (const std::invalid_argument &) {
-            throw ParsingException("Invalid a value");
-        }
+    else if (a > 0 && !fStr.empty()) {
         double f;
         try {
             f = c_locale_stod(fStr);
@@ -5864,23 +5923,52 @@ PROJStringParser::Private::buildDatum(const Step &step,
             ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
     }
 
-    else if (!RStr.empty()) {
-        double R;
+    else if (a > 0 && !eStr.empty()) {
+        double e;
         try {
-            R = c_locale_stod(RStr);
+            e = c_locale_stod(eStr);
         } catch (const std::invalid_argument &) {
-            throw ParsingException("Invalid R value");
+            throw ParsingException("Invalid e value");
         }
-        auto ellipsoid = Ellipsoid::createSphere(createMapWithUnknownName(),
-                                                 Length(R), guessBodyName(R));
+        double alpha = asin(e);    /* angular eccentricity */
+        double f = 1 - cos(alpha); /* = 1 - sqrt (1 - es); */
+        auto ellipsoid = Ellipsoid::createFlattenedSphere(
+                             createMapWithUnknownName(), Length(a),
+                             Scale(f != 0.0 ? 1.0 / f : 0.0), guessBodyName(a))
+                             ->identify();
         return GeodeticReferenceFrame::create(
             grfMap.set(IdentifiedObject::NAME_KEY,
                        title.empty() ? "unknown" : title.c_str()),
             ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
     }
 
-    if (!aStr.empty() && bStr.empty() && rfStr.empty()) {
-        throw ParsingException("a found, but b, f or rf missing");
+    else if (a > 0 && !esStr.empty()) {
+        double es;
+        try {
+            es = c_locale_stod(esStr);
+        } catch (const std::invalid_argument &) {
+            throw ParsingException("Invalid es value");
+        }
+        double f = 1 - sqrt(1 - es);
+        auto ellipsoid = Ellipsoid::createFlattenedSphere(
+                             createMapWithUnknownName(), Length(a),
+                             Scale(f != 0.0 ? 1.0 / f : 0.0), guessBodyName(a))
+                             ->identify();
+        return GeodeticReferenceFrame::create(
+            grfMap.set(IdentifiedObject::NAME_KEY,
+                       title.empty() ? "unknown" : title.c_str()),
+            ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
+    }
+
+    // If only a is specified, create a sphere
+    if (a > 0 && bStr.empty() && rfStr.empty() && eStr.empty() &&
+        esStr.empty()) {
+        auto ellipsoid = Ellipsoid::createSphere(createMapWithUnknownName(),
+                                                 Length(a), guessBodyName(a));
+        return GeodeticReferenceFrame::create(
+            grfMap.set(IdentifiedObject::NAME_KEY,
+                       title.empty() ? "unknown" : title.c_str()),
+            ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm));
     }
 
     if (!bStr.empty() && aStr.empty()) {
@@ -5895,6 +5983,14 @@ PROJStringParser::Private::buildDatum(const Step &step,
         throw ParsingException("f found, but a missing");
     }
 
+    if (!eStr.empty() && aStr.empty()) {
+        throw ParsingException("e found, but a missing");
+    }
+
+    if (!esStr.empty() && aStr.empty()) {
+        throw ParsingException("es found, but a missing");
+    }
+
     return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6326);
 }
 
@@ -6064,6 +6160,22 @@ PROJStringParser::Private::buildEllipsoidalCS(int iStep, int iUnitConvert,
 
 // ---------------------------------------------------------------------------
 
+static double getNumericValue(const std::string &paramValue,
+                              bool *pHasError = nullptr) {
+    try {
+        double value = c_locale_stod(paramValue);
+        if (pHasError)
+            *pHasError = false;
+        return value;
+    } catch (const std::invalid_argument &) {
+        if (pHasError)
+            *pHasError = true;
+        return 0.0;
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 GeographicCRSNNPtr
 PROJStringParser::Private::buildGeographicCRS(int iStep, int iUnitConvert,
                                               int iAxisSwap, bool ignoreVUnits,
@@ -6077,7 +6189,10 @@ PROJStringParser::Private::buildGeographicCRS(int iStep, int iUnitConvert,
 
     auto props = PropertyMap().set(IdentifiedObject::NAME_KEY,
                                    title.empty() ? "unknown" : title);
-    if (l_isGeographicStep && hasParamValue(step, "wktext")) {
+    if (l_isGeographicStep &&
+        (hasParamValue(step, "wktext") ||
+         hasParamValue(step, "lon_wrap") | hasParamValue(step, "geoc") ||
+         getNumericValue(getParamValue(step, "lon_0")) != 0.0)) {
         props.set("EXTENSION_PROJ4", projString_);
     }
 
@@ -6214,22 +6329,6 @@ static double getAngularValue(const std::string &paramValue,
 
 // ---------------------------------------------------------------------------
 
-static double getNumericValue(const std::string &paramValue,
-                              bool *pHasError = nullptr) {
-    try {
-        double value = c_locale_stod(paramValue);
-        if (pHasError)
-            *pHasError = false;
-        return value;
-    } catch (const std::invalid_argument &) {
-        if (pHasError)
-            *pHasError = true;
-        return 0.0;
-    }
-}
-
-// ---------------------------------------------------------------------------
-
 CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
     int iStep, GeographicCRSNNPtr geogCRS, int iUnitConvert, int iAxisSwap) {
     auto &step = steps_[iStep];
@@ -6285,7 +6384,10 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
         mapping =
             getMapping(PROJ_WKT2_NAME_METHOD_GEOSTATIONARY_SATELLITE_SWEEP_Y);
     } else if (step.name == "omerc") {
-        if (hasParamValue(step, "no_uoff") || hasParamValue(step, "no_off")) {
+        if (hasParamValue(step, "no_rot")) {
+            mapping = nullptr;
+        } else if (hasParamValue(step, "no_uoff") ||
+                   hasParamValue(step, "no_off")) {
             mapping =
                 getMapping(EPSG_CODE_METHOD_HOTINE_OBLIQUE_MERCATOR_VARIANT_A);
         } else if (hasParamValue(step, "lat_1") &&
@@ -6839,6 +6941,21 @@ PROJStringParser::Private::buildMolodenskyTransformation(
 
 // ---------------------------------------------------------------------------
 
+//! @cond Doxygen_Suppress
+static const metadata::ExtentPtr nullExtent{};
+
+static const metadata::ExtentPtr &getExtent(const crs::CRS *crs) {
+    const auto &domains = crs->domains();
+    if (!domains.empty()) {
+        return domains[0]->domainOfValidity();
+    }
+    return nullExtent;
+}
+
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
 /** \brief Instanciate a sub-class of BaseObject from a PROJ string.
  * @throw ParsingException
  */
@@ -6909,30 +7026,36 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
             }
             auto obj =
                 createFromUserInput(d->steps_[0].name, d->dbContext_, true);
-            auto geogCRS = dynamic_cast<GeographicCRS *>(obj.get());
-            if (geogCRS) {
-                // Override with longitude latitude in radian
-                return GeographicCRS::create(
-                    PropertyMap().set(IdentifiedObject::NAME_KEY,
-                                      geogCRS->nameStr()),
-                    geogCRS->datum(), geogCRS->datumEnsemble(),
-                    EllipsoidalCS::createLongitudeLatitude(
-                        UnitOfMeasure::RADIAN));
-            }
-            auto projCRS = dynamic_cast<ProjectedCRS *>(obj.get());
-            if (projCRS) {
-                // Override with easting northing orer
-                auto conv = projCRS->derivingConversionRef();
-                if (conv->method()->getEPSGCode() !=
-                    EPSG_CODE_METHOD_TRANSVERSE_MERCATOR_SOUTH_ORIENTATED) {
-                    return ProjectedCRS::create(
-                        PropertyMap().set(IdentifiedObject::NAME_KEY,
-                                          projCRS->nameStr()),
-                        projCRS->baseCRS(), conv,
-                        CartesianCS::createEastingNorthing(
-                            projCRS->coordinateSystem()
-                                ->axisList()[0]
-                                ->unit()));
+            auto crs = dynamic_cast<CRS *>(obj.get());
+            if (crs) {
+                PropertyMap properties;
+                properties.set(IdentifiedObject::NAME_KEY, crs->nameStr());
+                const auto &extent = getExtent(crs);
+                if (extent) {
+                    properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
+                                   NN_NO_CHECK(extent));
+                }
+                auto geogCRS = dynamic_cast<GeographicCRS *>(crs);
+                if (geogCRS) {
+                    // Override with longitude latitude in radian
+                    return GeographicCRS::create(
+                        properties, geogCRS->datum(), geogCRS->datumEnsemble(),
+                        EllipsoidalCS::createLongitudeLatitude(
+                            UnitOfMeasure::RADIAN));
+                }
+                auto projCRS = dynamic_cast<ProjectedCRS *>(crs);
+                if (projCRS) {
+                    // Override with easting northing order
+                    const auto &conv = projCRS->derivingConversionRef();
+                    if (conv->method()->getEPSGCode() !=
+                        EPSG_CODE_METHOD_TRANSVERSE_MERCATOR_SOUTH_ORIENTATED) {
+                        return ProjectedCRS::create(
+                            properties, projCRS->baseCRS(), conv,
+                            CartesianCS::createEastingNorthing(
+                                projCRS->coordinateSystem()
+                                    ->axisList()[0]
+                                    ->unit()));
+                    }
                 }
             }
             return obj;
diff --git a/src/pj_apply_gridshift.c b/src/pj_apply_gridshift.c
index 31e0124e..45ce5c8e 100644
--- a/src/pj_apply_gridshift.c
+++ b/src/pj_apply_gridshift.c
@@ -349,7 +349,7 @@ LP proj_hgrid_apply(PJ *P, LP lp, PJ_DIRECTION direction) {
     out = nad_cvt(lp, inverse, ct);
 
     if (out.lam == HUGE_VAL || out.phi == HUGE_VAL)
-        pj_ctx_set_errno(P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);
 
     return out;
 
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index 1a970374..38443f07 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -103,7 +103,6 @@ static PJ_COORD fwd_prepare (PJ *P, PJ_COORD coo) {
 }
 
 
-
 static PJ_COORD fwd_finalize (PJ *P, PJ_COORD coo) {
 
     switch (OUTPUT_UNITS) {
@@ -138,6 +137,14 @@ static PJ_COORD fwd_finalize (PJ *P, PJ_COORD coo) {
 
     case PJ_IO_UNITS_ANGULAR:
         coo.lpz.z = P->vfr_meter * (coo.lpz.z + P->z0);
+
+        if( P->is_long_wrap_set ) {
+            if( coo.lpz.lam != HUGE_VAL ) {
+                coo.lpz.lam  = P->long_wrap_center +
+                               adjlon(coo.lpz.lam - P->long_wrap_center);
+            }
+        }
+
         break;
     }
 
diff --git a/src/proj.h b/src/proj.h
index 4b599eba..31cd730c 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -1381,6 +1381,19 @@ PJ_OBJ PROJ_DLL *proj_obj_crs_get_sub_crs(PJ_OBJ *crs, int index);
 
 PJ_OBJ PROJ_DLL *proj_obj_crs_create_bound_crs_to_WGS84(PJ_OBJ *crs);
 
+PJ_OBJ PROJ_DLL *proj_obj_crs_get_coordinate_system(PJ_OBJ *crs);
+
+const char PROJ_DLL *proj_obj_cs_get_type(PJ_OBJ* cs);
+
+int PROJ_DLL proj_obj_cs_get_axis_count(PJ_OBJ *cs);
+
+int PROJ_DLL proj_obj_cs_get_axis_info(PJ_OBJ *cs, int index,
+                                       const char **pName,
+                                       const char **pAbbrev,
+                                       const char **pDirection,
+                                       double *pUnitConvFactor,
+                                       const char **pUnitName);
+
 PJ_OBJ PROJ_DLL *proj_obj_get_ellipsoid(PJ_OBJ *obj);
 
 int PROJ_DLL proj_obj_ellipsoid_get_parameters(PJ_OBJ *ellipsoid,
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 72e1a2d6..b7b500a7 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -794,7 +794,12 @@ PJ  *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *srid_from, const char
         return NULL;
     }
 
-    P = proj_create(ctx, proj_string);
+    if( proj_string[0] == '\0' ) {
+        /* Null transform ? */
+        P = proj_create(ctx, "proj=affine");
+    } else {
+        P = proj_create(ctx, proj_string);
+    }
 
     proj_obj_unref(op);