API: add CRS::promoteTo3D(), proj_crs_promote_to_3D() and proj_crs_create_projected_3D_crs_from_2D() (fixes #1587)

Also add a --3d switch to projinfo

11 files changed, 485 insertions, 19 deletions

diff --git a/docs/source/apps/projinfo.rst b/docs/source/apps/projinfo.rst
index 6681881c..8fbf0546 100644
--- a/docs/source/apps/projinfo.rst
+++ b/docs/source/apps/projinfo.rst
@@ -24,7 +24,7 @@ Synopsis
     |    [--pivot-crs always|if_no_direct_transformation|never|{auth:code[,auth:code]*}]
     |    [--boundcrs-to-wgs84]
     |    [--main-db-path path] [--aux-db-path path]*
-    |    [--identify]
+    |    [--identify] [--3d]
     |    [--c-ify] [--single-line]
     |    {object_definition} | {object_reference} | (-s {srs_def} -t {srs_def})
     |
@@ -220,6 +220,14 @@ The following control parameters can appear in any order:
     For example, `+proj=utm +zone=31 +datum=WGS84 +type=crs` will be identified
     with a likelihood of 70% to EPSG:32631
 
+.. option:: --3d
+
+    .. versionadded:: 7.0
+
+    "Promote" the CRS(s) to their 3D version. Useful for example when wanting
+    to transform between a 2D projected CRS with elevations as ellipsoidal
+    height to a 3D geographic CRS or a compoundCRS.
+
 .. option:: --c-ify
 
     For developers only. Modify the string output of the utility so that it
diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp
index 7062885d..65459058 100644
--- a/include/proj/crs.hpp
+++ b/include/proj/crs.hpp
@@ -107,6 +107,10 @@ class PROJ_GCC_DLL CRS : public common::ObjectUsage,
     PROJ_DLL std::list<CRSNNPtr>
     getNonDeprecated(const io::DatabaseContextNNPtr &dbContext) const;
 
+    PROJ_DLL CRSNNPtr
+    promoteTo3D(const std::string &newName,
+                const io::DatabaseContextPtr &dbContext) const;
+
     PROJ_PRIVATE :
         //! @cond Doxygen_Suppress
         PROJ_INTERNAL const GeodeticCRS *
diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index c6937d18..3206c3c8 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -114,6 +114,7 @@ osgeo::proj::crs::CRS::extractGeographicCRS() const
 osgeo::proj::crs::CRS::extractVerticalCRS() const
 osgeo::proj::crs::CRS::getNonDeprecated(dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::io::DatabaseContext> > const&) const
 osgeo::proj::crs::CRS::identify(std::shared_ptr<osgeo::proj::io::AuthorityFactory> const&) const
+osgeo::proj::crs::CRS::promoteTo3D(std::string const&, std::shared_ptr<osgeo::proj::io::DatabaseContext> const&) const
 osgeo::proj::crs::CRS::shallowClone() const
 osgeo::proj::crs::CRS::stripVerticalComponent() const
 osgeo::proj::crs::DerivedCRS::baseCRS() const
@@ -899,6 +900,7 @@ proj_crs_alter_geodetic_crs
 proj_crs_alter_parameters_linear_unit
 proj_crs_create_bound_crs
 proj_crs_create_bound_crs_to_WGS84
+proj_crs_create_projected_3D_crs_from_2D
 proj_crs_get_coordinate_system
 proj_crs_get_coordoperation
 proj_crs_get_datum
@@ -906,6 +908,7 @@ proj_crs_get_geodetic_crs
 proj_crs_get_horizontal_datum
 proj_crs_get_sub_crs
 proj_crs_info_list_destroy
+proj_crs_promote_to_3D
 proj_cs_get_axis_count
 proj_cs_get_axis_info
 proj_cs_get_type
diff --git a/src/apps/projinfo.cpp b/src/apps/projinfo.cpp
index ba9e8350..1760ae29 100644
--- a/src/apps/projinfo.cpp
+++ b/src/apps/projinfo.cpp
@@ -95,7 +95,7 @@ static void usage() {
         << "                [--boundcrs-to-wgs84]" << std::endl
         << "                [--main-db-path path] [--aux-db-path path]*"
         << std::endl
-        << "                [--identify]" << std::endl
+        << "                [--identify] [--3d]" << std::endl
         << "                [--c-ify] [--single-line]" << std::endl
         << "                {object_definition} | (-s {srs_def} -t {srs_def})"
         << std::endl;
@@ -143,7 +143,7 @@ static BaseObjectNNPtr buildObject(
     const std::string &kind, const std::string &context,
     bool buildBoundCRSToWGS84,
     CoordinateOperationContext::IntermediateCRSUse allowUseIntermediateCRS,
-    bool quiet) {
+    bool promoteTo3D, bool quiet) {
     BaseObjectPtr obj;
 
     std::string l_user_string(user_string);
@@ -225,6 +225,13 @@ static BaseObjectNNPtr buildObject(
         }
     }
 
+    if (promoteTo3D) {
+        auto crs = std::dynamic_pointer_cast<CRS>(obj);
+        if (crs) {
+            obj = crs->promoteTo3D(std::string(), dbContext).as_nullable();
+        }
+    }
+
     return NN_NO_CHECK(obj);
 }
 
@@ -594,6 +601,37 @@ static void outputOperationSummary(const CoordinateOperationNNPtr &op,
 
 // ---------------------------------------------------------------------------
 
+static size_t getAxisCount(const CRSNNPtr &crs) {
+    const auto singleCRS = dynamic_cast<const SingleCRS *>(crs.get());
+    if (singleCRS) {
+        return singleCRS->coordinateSystem()->axisList().size();
+    }
+    const auto compoundCRS = dynamic_cast<const CompoundCRS *>(crs.get());
+    if (compoundCRS) {
+        size_t axisCount = 0;
+        const auto &components = compoundCRS->componentReferenceSystems();
+        for (const auto &subCRS : components) {
+            axisCount += getAxisCount(subCRS);
+        }
+        return axisCount;
+    }
+    const auto boundCRS = dynamic_cast<const BoundCRS *>(crs.get());
+    if (boundCRS) {
+        return getAxisCount(boundCRS->baseCRS());
+    }
+    return 0;
+}
+
+// ---------------------------------------------------------------------------
+
+static bool is2D(const CRSNNPtr &crs) { return getAxisCount(crs) == 2; }
+
+// ---------------------------------------------------------------------------
+
+static bool is3D(const CRSNNPtr &crs) { return getAxisCount(crs) == 3; }
+
+// ---------------------------------------------------------------------------
+
 static void outputOperations(
     DatabaseContextPtr dbContext, const std::string &sourceCRSStr,
     const std::string &targetCRSStr, const ExtentPtr &bboxFilter,
@@ -604,24 +642,38 @@ static void outputOperations(
     CoordinateOperationContext::IntermediateCRSUse allowUseIntermediateCRS,
     const std::vector<std::pair<std::string, std::string>> &pivots,
     const std::string &authority, bool usePROJGridAlternatives,
-    bool showSuperseded, const OutputOptions &outputOpt, bool summary) {
-    auto sourceObj = buildObject(
-        dbContext, sourceCRSStr, "crs", "source CRS", false,
-        CoordinateOperationContext::IntermediateCRSUse::NEVER, outputOpt.quiet);
+    bool showSuperseded, bool promoteTo3D, const OutputOptions &outputOpt,
+    bool summary) {
+    auto sourceObj =
+        buildObject(dbContext, sourceCRSStr, "crs", "source CRS", false,
+                    CoordinateOperationContext::IntermediateCRSUse::NEVER,
+                    promoteTo3D, outputOpt.quiet);
     auto sourceCRS = nn_dynamic_pointer_cast<CRS>(sourceObj);
     if (!sourceCRS) {
         std::cerr << "source CRS string is not a CRS" << std::endl;
         std::exit(1);
     }
+    auto nnSourceCRS = NN_NO_CHECK(sourceCRS);
 
-    auto targetObj = buildObject(
-        dbContext, targetCRSStr, "crs", "target CRS", false,
-        CoordinateOperationContext::IntermediateCRSUse::NEVER, outputOpt.quiet);
+    auto targetObj =
+        buildObject(dbContext, targetCRSStr, "crs", "target CRS", false,
+                    CoordinateOperationContext::IntermediateCRSUse::NEVER,
+                    promoteTo3D, outputOpt.quiet);
     auto targetCRS = nn_dynamic_pointer_cast<CRS>(targetObj);
     if (!targetCRS) {
         std::cerr << "target CRS string is not a CRS" << std::endl;
         std::exit(1);
     }
+    auto nnTargetCRS = NN_NO_CHECK(targetCRS);
+
+    if (!promoteTo3D && !outputOpt.quiet &&
+        ((is2D(nnSourceCRS) && is3D(nnTargetCRS)) ||
+         (is3D(nnSourceCRS) && is2D(nnTargetCRS)))) {
+        std::cerr << "Warning: mix of 2D and 3D CRS. Vertical transformations, "
+                     "if available, will not be applied. Consider using 3D "
+                     "version of the CRS, or the --3d switch"
+                  << std::endl;
+    }
 
     std::vector<CoordinateOperationNNPtr> list;
     size_t spatialCriterionPartialIntersectionResultCount = 0;
@@ -641,7 +693,7 @@ static void outputOperations(
         ctxt->setUsePROJAlternativeGridNames(usePROJGridAlternatives);
         ctxt->setDiscardSuperseded(!showSuperseded);
         list = CoordinateOperationFactory::create()->createOperations(
-            NN_NO_CHECK(sourceCRS), NN_NO_CHECK(targetCRS), ctxt);
+            nnSourceCRS, nnTargetCRS, ctxt);
         if (!spatialCriterionExplicitlySpecified &&
             spatialCriterion == CoordinateOperationContext::SpatialCriterion::
                                     STRICT_CONTAINMENT) {
@@ -651,8 +703,7 @@ static void outputOperations(
                         PARTIAL_INTERSECTION);
                 spatialCriterionPartialIntersectionResultCount =
                     CoordinateOperationFactory::create()
-                        ->createOperations(NN_NO_CHECK(sourceCRS),
-                                           NN_NO_CHECK(targetCRS), ctxt)
+                        ->createOperations(nnSourceCRS, nnTargetCRS, ctxt)
                         .size();
             } catch (const std::exception &) {
             }
@@ -735,6 +786,7 @@ int main(int argc, char **argv) {
     std::string authority;
     bool identify = false;
     bool showSuperseded = false;
+    bool promoteTo3D = false;
 
     for (int i = 1; i < argc; i++) {
         std::string arg(argv[i]);
@@ -984,6 +1036,8 @@ int main(int argc, char **argv) {
             showSuperseded = true;
         } else if (arg == "--lax") {
             outputOpt.strict = false;
+        } else if (ci_equal(arg, "--3d")) {
+            promoteTo3D = true;
         } else if (arg == "-?" || arg == "--help") {
             usage();
         } else if (arg[0] == '-') {
@@ -1054,7 +1108,8 @@ int main(int argc, char **argv) {
         try {
             auto obj(buildObject(dbContext, user_string, objectKind,
                                  "input string", buildBoundCRSToWGS84,
-                                 allowUseIntermediateCRS, outputOpt.quiet));
+                                 allowUseIntermediateCRS, promoteTo3D,
+                                 outputOpt.quiet));
             if (guessDialect) {
                 auto dialect = WKTParser().guessDialect(user_string);
                 std::cout << "Guessed WKT dialect: ";
@@ -1193,7 +1248,7 @@ int main(int argc, char **argv) {
                              spatialCriterionExplicitlySpecified, crsExtentUse,
                              gridAvailabilityUse, allowUseIntermediateCRS,
                              pivots, authority, usePROJGridAlternatives,
-                             showSuperseded, outputOpt, summary);
+                             showSuperseded, promoteTo3D, outputOpt, summary);
         } catch (const std::exception &e) {
             std::cerr << "outputOperations() failed with: " << e.what()
                       << std::endl;
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 45eb16d1..3a1b66af 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -3148,7 +3148,7 @@ PJ *proj_crs_alter_cs_linear_unit(PJ_CONTEXT *ctx, const PJ *obj,
 
 // ---------------------------------------------------------------------------
 
-/** \brief Return a copy of the CRS with the lineaer units of the parameters
+/** \brief Return a copy of the CRS with the linear units of the parameters
  * of its conversion modified.
  *
  * The CRS must be or contain a ProjectedCRS, VerticalCRS or a GeocentricCRS.
@@ -3197,6 +3197,151 @@ PJ *proj_crs_alter_parameters_linear_unit(PJ_CONTEXT *ctx, const PJ *obj,
 
 // ---------------------------------------------------------------------------
 
+/** \brief Create a 3D CRS from an existing 2D CRS.
+ *
+ * The new axis will be ellipsoidal height, oriented upwards, and with metre
+ * units.
+ *
+ * See osgeo::proj::crs::CRS::promoteTo3D().
+ *
+ * The returned object must be unreferenced with proj_destroy() after
+ * use.
+ * It should be used by at most one thread at a time.
+ *
+ * @param ctx PROJ context, or NULL for default context
+ * @param crs_3D_name CRS name. Or NULL (in which case the name of crs_2D
+ * will be used)
+ * @param crs_2D 2D CRS to be "promoted" to 3D. Must not be NULL.
+ *
+ * @return Object that must be unreferenced with
+ * proj_destroy(), or NULL in case of error.
+ * @since 7.0
+ */
+PJ *proj_crs_promote_to_3D(PJ_CONTEXT *ctx, const char *crs_3D_name,
+                           const PJ *crs_2D) {
+    SANITIZE_CTX(ctx);
+    auto cpp_2D_crs = dynamic_cast<const CRS *>(crs_2D->iso_obj.get());
+    if (!cpp_2D_crs) {
+        proj_log_error(ctx, __FUNCTION__, "crs_2D is not a CRS");
+        return nullptr;
+    }
+    try {
+        auto dbContext = getDBcontextNoException(ctx, __FUNCTION__);
+        return pj_obj_create(
+            ctx, cpp_2D_crs->promoteTo3D(crs_3D_name ? std::string(crs_3D_name)
+                                                     : cpp_2D_crs->nameStr(),
+                                         dbContext));
+    } catch (const std::exception &e) {
+        proj_log_error(ctx, __FUNCTION__, e.what());
+        if (ctx->cpp_context) {
+            ctx->cpp_context->autoCloseDbIfNeeded();
+        }
+        return nullptr;
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Create a projected 3D CRS from an existing projected 2D CRS.
+ *
+ * The passed projected_2D_crs is used so that its name is replaced by
+ * crs_name and its base geographic CRS is replaced by geog_3D_crs. The vertical
+ * axis of geog_3D_crs (ellipsoidal height) will be added as the 3rd axis of
+ * the resulting projected 3D CRS.
+ * Normally, the passed geog_3D_crs should be the 3D counterpart of the original
+ * 2D base geographic CRS of projected_2D_crs, but such no check is done.
+ *
+ * It is also possible to invoke this function with a NULL geog_3D_crs. In which
+ * case, the existing base geographic 2D CRS of projected_2D_crs will be
+ * automatically promoted to 3D by assuming a 3rd axis being an ellipsoidal
+ * height, oriented upwards, and with metre units. This is equivalent to using
+ * proj_crs_promote_to_3D().
+ *
+ * The returned object must be unreferenced with proj_destroy() after
+ * use.
+ * It should be used by at most one thread at a time.
+ *
+ * @param ctx PROJ context, or NULL for default context
+ * @param crs_name CRS name. Or NULL (in which case the name of projected_2D_crs
+ * will be used)
+ * @param projected_2D_crs Projected 2D CRS to be "promoted" to 3D. Must not be
+ * NULL.
+ * @param geog_3D_crs Base geographic 3D CRS for the new CRS. May be NULL.
+ *
+ * @return Object that must be unreferenced with
+ * proj_destroy(), or NULL in case of error.
+ * @since 7.0
+ */
+PJ *proj_crs_create_projected_3D_crs_from_2D(PJ_CONTEXT *ctx,
+                                             const char *crs_name,
+                                             const PJ *projected_2D_crs,
+                                             const PJ *geog_3D_crs) {
+    SANITIZE_CTX(ctx);
+    auto cpp_projected_2D_crs =
+        dynamic_cast<const ProjectedCRS *>(projected_2D_crs->iso_obj.get());
+    if (!cpp_projected_2D_crs) {
+        proj_log_error(ctx, __FUNCTION__,
+                       "projected_2D_crs is not a Projected CRS");
+        return nullptr;
+    }
+    const auto &oldCS = cpp_projected_2D_crs->coordinateSystem();
+    const auto &oldCSAxisList = oldCS->axisList();
+
+    if (geog_3D_crs && geog_3D_crs->iso_obj) {
+        auto cpp_geog_3D_CRS =
+            std::dynamic_pointer_cast<GeographicCRS>(geog_3D_crs->iso_obj);
+        if (!cpp_geog_3D_CRS) {
+            proj_log_error(ctx, __FUNCTION__,
+                           "geog_3D_crs is not a Geographic CRS");
+            return nullptr;
+        }
+
+        const auto &geogCS = cpp_geog_3D_CRS->coordinateSystem();
+        const auto &geogCSAxisList = geogCS->axisList();
+        if (geogCSAxisList.size() != 3) {
+            proj_log_error(ctx, __FUNCTION__,
+                           "geog_3D_crs is not a Geographic 3D CRS");
+            return nullptr;
+        }
+        try {
+            auto newCS =
+                cs::CartesianCS::create(PropertyMap(), oldCSAxisList[0],
+                                        oldCSAxisList[1], geogCSAxisList[2]);
+            return pj_obj_create(
+                ctx,
+                ProjectedCRS::create(
+                    createPropertyMapName(
+                        crs_name ? crs_name
+                                 : cpp_projected_2D_crs->nameStr().c_str()),
+                    NN_NO_CHECK(cpp_geog_3D_CRS),
+                    cpp_projected_2D_crs->derivingConversionRef(), newCS));
+        } catch (const std::exception &e) {
+            proj_log_error(ctx, __FUNCTION__, e.what());
+            if (ctx->cpp_context) {
+                ctx->cpp_context->autoCloseDbIfNeeded();
+            }
+            return nullptr;
+        }
+    } else {
+        try {
+            auto dbContext = getDBcontextNoException(ctx, __FUNCTION__);
+            return pj_obj_create(ctx,
+                                 cpp_projected_2D_crs->promoteTo3D(
+                                     crs_name ? std::string(crs_name)
+                                              : cpp_projected_2D_crs->nameStr(),
+                                     dbContext));
+        } catch (const std::exception &e) {
+            proj_log_error(ctx, __FUNCTION__, e.what());
+            if (ctx->cpp_context) {
+                ctx->cpp_context->autoCloseDbIfNeeded();
+            }
+            return nullptr;
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Instantiate a EngineeringCRS with just a name
  *
  * The returned object must be unreferenced with proj_destroy() after
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 07f4eb69..1b7078c3 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -756,6 +756,96 @@ CRS::getNonDeprecated(const io::DatabaseContextNNPtr &dbContext) const {
 
 // ---------------------------------------------------------------------------
 
+/** \brief Return a variant of this CRS "promoted" to a 3D one, if not already
+ * the case.
+ *
+ * The new axis will be ellipsoidal height, oriented upwards, and with metre
+ * units.
+ *
+ * @param newName Name of the new CRS. If empty, nameStr() will be used.
+ * @param dbContext Database context to look for potentially already registered
+ *                  3D CRS. May be nullptr.
+ * @return a new CRS promoted to 3D, or the current one if already 3D or not
+ * applicable.
+ * @since 7.0
+ */
+CRSNNPtr CRS::promoteTo3D(const std::string &newName,
+                          const io::DatabaseContextPtr &dbContext) const {
+
+    const auto geogCRS = dynamic_cast<const GeographicCRS *>(this);
+    if (geogCRS) {
+        const auto &axisList = geogCRS->coordinateSystem()->axisList();
+        if (axisList.size() == 2) {
+            const auto &l_identifiers = identifiers();
+            // First check if there is a Geographic 3D CRS in the database
+            // of the same name.
+            // This is the common practice in the EPSG dataset.
+            if (dbContext && l_identifiers.size() == 1) {
+                auto authFactory = io::AuthorityFactory::create(
+                    NN_NO_CHECK(dbContext), *(l_identifiers[0]->codeSpace()));
+                auto res = authFactory->createObjectsFromName(
+                    nameStr(),
+                    {io::AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS},
+                    false);
+                if (!res.empty()) {
+                    const auto &firstRes = res.front();
+                    if (geogCRS->is2DPartOf3D(NN_NO_CHECK(
+                            dynamic_cast<GeographicCRS *>(firstRes.get())))) {
+                        return NN_NO_CHECK(
+                            util::nn_dynamic_pointer_cast<CRS>(firstRes));
+                    }
+                }
+            }
+
+            auto upAxis = cs::CoordinateSystemAxis::create(
+                util::PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                        cs::AxisName::Ellipsoidal_height),
+                cs::AxisAbbreviation::h, cs::AxisDirection::UP,
+                common::UnitOfMeasure::METRE);
+            auto cs = cs::EllipsoidalCS::create(
+                util::PropertyMap(), axisList[0], axisList[1], upAxis);
+            return util::nn_static_pointer_cast<CRS>(GeographicCRS::create(
+                util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                        !newName.empty() ? newName : nameStr()),
+                geogCRS->datum(), geogCRS->datumEnsemble(), cs));
+        }
+    }
+
+    const auto projCRS = dynamic_cast<const ProjectedCRS *>(this);
+    if (projCRS) {
+        const auto &axisList = projCRS->coordinateSystem()->axisList();
+        if (axisList.size() == 2) {
+            auto base3DCRS =
+                projCRS->baseCRS()->promoteTo3D(std::string(), dbContext);
+            auto upAxis = cs::CoordinateSystemAxis::create(
+                util::PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                        cs::AxisName::Ellipsoidal_height),
+                cs::AxisAbbreviation::h, cs::AxisDirection::UP,
+                common::UnitOfMeasure::METRE);
+            auto cs = cs::CartesianCS::create(util::PropertyMap(), axisList[0],
+                                              axisList[1], upAxis);
+            return util::nn_static_pointer_cast<CRS>(ProjectedCRS::create(
+                util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                        !newName.empty() ? newName : nameStr()),
+                NN_NO_CHECK(
+                    util::nn_dynamic_pointer_cast<GeodeticCRS>(base3DCRS)),
+                projCRS->derivingConversionRef(), cs));
+        }
+    }
+
+    const auto boundCRS = dynamic_cast<const BoundCRS *>(this);
+    if (boundCRS) {
+        return BoundCRS::create(
+            boundCRS->baseCRS()->promoteTo3D(newName, dbContext),
+            boundCRS->hubCRS(), boundCRS->transformation());
+    }
+
+    return NN_NO_CHECK(
+        std::static_pointer_cast<CRS>(shared_from_this().as_nullable()));
+}
+
+// ---------------------------------------------------------------------------
+
 //! @cond Doxygen_Suppress
 
 std::list<std::pair<CRSNNPtr, int>>
@@ -1898,14 +1988,18 @@ bool GeographicCRS::is2DPartOf3D(util::nn<const GeographicCRS *> other)
     const auto &secondAxis = axis[1];
     const auto &otherFirstAxis = otherAxis[0];
     const auto &otherSecondAxis = otherAxis[1];
-    if (!(firstAxis->_isEquivalentTo(otherFirstAxis.get()) &&
-          secondAxis->_isEquivalentTo(otherSecondAxis.get()))) {
+    if (!(firstAxis->_isEquivalentTo(
+              otherFirstAxis.get(), util::IComparable::Criterion::EQUIVALENT) &&
+          secondAxis->_isEquivalentTo(
+              otherSecondAxis.get(),
+              util::IComparable::Criterion::EQUIVALENT))) {
         return false;
     }
     const auto &thisDatum = GeodeticCRS::getPrivate()->datum_;
     const auto &otherDatum = other->GeodeticCRS::getPrivate()->datum_;
     if (thisDatum && otherDatum) {
-        return thisDatum->_isEquivalentTo(otherDatum.get());
+        return thisDatum->_isEquivalentTo(
+            otherDatum.get(), util::IComparable::Criterion::EQUIVALENT);
     }
     return false;
 }
diff --git a/src/proj_experimental.h b/src/proj_experimental.h
index 5a96203c..5e6bbb13 100644
--- a/src/proj_experimental.h
+++ b/src/proj_experimental.h
@@ -226,6 +226,15 @@ PJ PROJ_DLL *proj_crs_alter_parameters_linear_unit(
                                                     const char *unit_code,
                                                     int convert_to_new_unit);
 
+PJ PROJ_DLL *proj_crs_promote_to_3D(PJ_CONTEXT *ctx,
+                                    const char* crs_3D_name,
+                                    const PJ* crs_2D);
+
+PJ PROJ_DLL *proj_crs_create_projected_3D_crs_from_2D(PJ_CONTEXT *ctx,
+                                                      const char* crs_name,
+                                                      const PJ* projected_2D_crs,
+                                                      const PJ* geog_3D_crs);
+
 PJ PROJ_DLL *proj_create_engineering_crs(PJ_CONTEXT *ctx,
                                                  const char *crsName);
 
diff --git a/test/cli/testprojinfo b/test/cli/testprojinfo
index f070c999..2b3432a9 100755
--- a/test/cli/testprojinfo
+++ b/test/cli/testprojinfo
@@ -115,6 +115,14 @@ echo "Testing NGF IGN69 height to RGF93: projinfo -s EPSG:5720 -t EPSG:4965 -o P
 $EXE -s EPSG:5720 -t EPSG:4965 -o PROJ >>${OUT} 2>&1
 echo "" >>${OUT}
 
+echo "Testing -s EPSG:32631 -t EPSG:4326+3855 --summary" >> ${OUT}
+$EXE -s EPSG:32631 -t EPSG:4326+3855  --summary 2>>${OUT}
+echo "" >>${OUT}
+
+echo "Testing -s EPSG:32631 -t EPSG:4326+3855 --3d --summary" >> ${OUT}
+$EXE -s EPSG:32631 -t EPSG:4326+3855 --3d --summary >>${OUT} 2>&1
+echo "" >>${OUT}
+
 # do 'diff' with distribution results
 echo "diff ${OUT} with testprojinfo_out.dist"
 diff -u ${OUT} ${TEST_CLI_DIR}/testprojinfo_out.dist
diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist
index 8cb37cf7..9ec74fde 100644
--- a/test/cli/testprojinfo_out.dist
+++ b/test/cli/testprojinfo_out.dist
@@ -821,3 +821,12 @@ INVERSE(EPSG):10000, Inverse of RGF93 to NGF IGN69 height (1), 0.5 m, France - m
 PROJ string:
 +proj=pipeline +step +inv +proj=vgridshift +grids=ggf97a.txt +multiplier=1
 
+Testing -s EPSG:32631 -t EPSG:4326+3855 --summary
+Warning: mix of 2D and 3D CRS. Vertical transformations, if available, will not be applied. Consider using 3D version of the CRS, or the --3d switch
+
+Testing -s EPSG:32631 -t EPSG:4326+3855 --3d --summary
+Candidate operations found: 3
+unknown id, Inverse of UTM zone 31N + WGS 84 to EGM2008 height (1), 1 m, World
+unknown id, Inverse of UTM zone 31N + WGS 84 to EGM2008 height (2), 0.5 m, World
+unknown id, Inverse of UTM zone 31N + Inverse of Transformation from EGM2008 height to WGS 84 (ballpark vertical transformation, without ellipsoid height to vertical height correction), unknown accuracy, World, has ballpark transformation
+
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index fd129c80..b8dde430 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -3892,4 +3892,78 @@ TEST_F(CApi, proj_create_ellipsoidal_3D_cs) {
     }
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_crs_promote_to_3D) {
+
+    auto crs2D =
+        proj_create(m_ctxt, GeographicCRS::EPSG_4326
+                                ->exportToWKT(WKTFormatter::create().get())
+                                .c_str());
+    ObjectKeeper keeper_crs2D(crs2D);
+    EXPECT_NE(crs2D, nullptr);
+
+    auto crs3D = proj_crs_promote_to_3D(m_ctxt, nullptr, crs2D);
+    ObjectKeeper keeper_crs3D(crs3D);
+    EXPECT_NE(crs3D, nullptr);
+
+    auto cs = proj_crs_get_coordinate_system(m_ctxt, crs3D);
+    ASSERT_NE(cs, nullptr);
+    ObjectKeeper keeperCs(cs);
+    EXPECT_EQ(proj_cs_get_axis_count(m_ctxt, cs), 3);
+
+    auto code = proj_get_id_code(crs3D, 0);
+    ASSERT_TRUE(code != nullptr);
+    EXPECT_EQ(code, std::string("4979"));
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_crs_create_projected_3D_crs_from_2D) {
+
+    auto projCRS = proj_create_from_database(m_ctxt, "EPSG", "32631",
+                                             PJ_CATEGORY_CRS, false, nullptr);
+    ASSERT_NE(projCRS, nullptr);
+    ObjectKeeper keeper_projCRS(projCRS);
+
+    {
+        auto geog3DCRS = proj_create_from_database(
+            m_ctxt, "EPSG", "4979", PJ_CATEGORY_CRS, false, nullptr);
+        ASSERT_NE(geog3DCRS, nullptr);
+        ObjectKeeper keeper_geog3DCRS(geog3DCRS);
+
+        auto crs3D = proj_crs_create_projected_3D_crs_from_2D(
+            m_ctxt, nullptr, projCRS, geog3DCRS);
+        ObjectKeeper keeper_crs3D(crs3D);
+        EXPECT_NE(crs3D, nullptr);
+
+        EXPECT_EQ(proj_get_type(crs3D), PJ_TYPE_PROJECTED_CRS);
+
+        EXPECT_EQ(std::string(proj_get_name(crs3D)),
+                  std::string(proj_get_name(projCRS)));
+
+        auto cs = proj_crs_get_coordinate_system(m_ctxt, crs3D);
+        ASSERT_NE(cs, nullptr);
+        ObjectKeeper keeperCs(cs);
+        EXPECT_EQ(proj_cs_get_axis_count(m_ctxt, cs), 3);
+    }
+
+    {
+        auto crs3D = proj_crs_create_projected_3D_crs_from_2D(m_ctxt, nullptr,
+                                                              projCRS, nullptr);
+        ObjectKeeper keeper_crs3D(crs3D);
+        EXPECT_NE(crs3D, nullptr);
+
+        EXPECT_EQ(proj_get_type(crs3D), PJ_TYPE_PROJECTED_CRS);
+
+        EXPECT_EQ(std::string(proj_get_name(crs3D)),
+                  std::string(proj_get_name(projCRS)));
+
+        auto cs = proj_crs_get_coordinate_system(m_ctxt, crs3D);
+        ASSERT_NE(cs, nullptr);
+        ObjectKeeper keeperCs(cs);
+        EXPECT_EQ(proj_cs_get_axis_count(m_ctxt, cs), 3);
+    }
+}
+
 } // namespace
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index 5ee62ce4..605df714 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -5224,3 +5224,60 @@ TEST(crs, getNonDeprecated) {
         ASSERT_EQ(list.size(), 1U);
     }
 }
+
+// ---------------------------------------------------------------------------
+
+TEST(crs, promoteTo3D) {
+    auto dbContext = DatabaseContext::create();
+    {
+        auto crs = GeographicCRS::EPSG_4326;
+        auto crs3D = crs->promoteTo3D(std::string(), nullptr);
+        auto crs3DAsGeog = nn_dynamic_pointer_cast<GeographicCRS>(crs3D);
+        ASSERT_TRUE(crs3DAsGeog != nullptr);
+        EXPECT_EQ(crs3DAsGeog->coordinateSystem()->axisList().size(), 3U);
+        EXPECT_TRUE(crs3D->promoteTo3D(std::string(), nullptr)
+                        ->isEquivalentTo(crs3D.get()));
+    }
+    {
+        auto crs = GeographicCRS::EPSG_4326;
+        auto crs3D = crs->promoteTo3D(std::string(), dbContext);
+        auto crs3DAsGeog = nn_dynamic_pointer_cast<GeographicCRS>(crs3D);
+        ASSERT_TRUE(crs3DAsGeog != nullptr);
+        EXPECT_EQ(crs3DAsGeog->coordinateSystem()->axisList().size(), 3U);
+        EXPECT_TRUE(!crs3DAsGeog->identifiers().empty());
+    }
+    {
+        auto crs = createProjected();
+        auto crs3D = crs->promoteTo3D(std::string(), nullptr);
+        auto crs3DAsProjected = nn_dynamic_pointer_cast<ProjectedCRS>(crs3D);
+        ASSERT_TRUE(crs3DAsProjected != nullptr);
+        EXPECT_EQ(crs3DAsProjected->coordinateSystem()->axisList().size(), 3U);
+        EXPECT_EQ(
+            crs3DAsProjected->baseCRS()->coordinateSystem()->axisList().size(),
+            3U);
+        EXPECT_TRUE(crs3D->promoteTo3D(std::string(), nullptr)
+                        ->isEquivalentTo(crs3D.get()));
+    }
+    {
+        auto crs = createProjected();
+        auto crs3D = crs->promoteTo3D(std::string(), dbContext);
+        auto crs3DAsProjected = nn_dynamic_pointer_cast<ProjectedCRS>(crs3D);
+        ASSERT_TRUE(crs3DAsProjected != nullptr);
+        EXPECT_EQ(crs3DAsProjected->coordinateSystem()->axisList().size(), 3U);
+        EXPECT_EQ(
+            crs3DAsProjected->baseCRS()->coordinateSystem()->axisList().size(),
+            3U);
+        EXPECT_TRUE(!crs3DAsProjected->baseCRS()->identifiers().empty());
+    }
+    {
+        auto crs = BoundCRS::createFromTOWGS84(
+            createProjected(), std::vector<double>{1, 2, 3, 4, 5, 6, 7});
+        auto crs3D = crs->promoteTo3D(std::string(), dbContext);
+        auto crs3DAsBound = nn_dynamic_pointer_cast<BoundCRS>(crs3D);
+        ASSERT_TRUE(crs3DAsBound != nullptr);
+        auto baseCRS =
+            nn_dynamic_pointer_cast<ProjectedCRS>(crs3DAsBound->baseCRS());
+        ASSERT_TRUE(baseCRS != nullptr);
+        EXPECT_EQ(baseCRS->coordinateSystem()->axisList().size(), 3U);
+    }
+}