Merge remote-tracking branch 'osgeo/master'

6 files changed, 1069 insertions, 295 deletions

diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 7582f37f..6bc1f166 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -79,10 +79,10 @@ static void PROJ_NO_INLINE proj_log_error(PJ_CONTEXT *ctx, const char *function,
     msg += ": ";
     msg += text;
     ctx->logger(ctx->logger_app_data, PJ_LOG_ERROR, msg.c_str());
-    auto previous_errno = pj_ctx_get_errno(ctx);
+    auto previous_errno = proj_context_errno(ctx);
     if (previous_errno == 0) {
         // only set errno if it wasn't set deeper down the call stack
-        pj_ctx_set_errno(ctx, PJD_ERR_GENERIC_ERROR);
+        proj_context_errno_set(ctx, PJD_ERR_GENERIC_ERROR);
     }
 }
 
@@ -675,6 +675,9 @@ PJ *proj_create_from_database(PJ_CONTEXT *ctx, const char *auth_name,
                           codeStr, usePROJAlternativeGridNames != 0)
                       .as_nullable();
             break;
+        case PJ_CATEGORY_DATUM_ENSEMBLE:
+            obj = factory->createDatumEnsemble(codeStr).as_nullable();
+            break;
         }
         return pj_obj_create(ctx, NN_NO_CHECK(obj));
     } catch (const std::exception &e) {
@@ -916,7 +919,7 @@ convertPJObjectTypeToObjectType(PJ_TYPE type, bool &valid) {
         break;
 
     case PJ_TYPE_DATUM_ENSEMBLE:
-        cppType = AuthorityFactory::ObjectType::DATUM;
+        cppType = AuthorityFactory::ObjectType::DATUM_ENSEMBLE;
         break;
 
     case PJ_TYPE_TEMPORAL_DATUM:
@@ -1418,6 +1421,13 @@ const char *proj_get_id_code(const PJ *obj, int index) {
  * variants, for WKT1_GDAL for ProjectedCRS with easting/northing ordering
  * (otherwise stripped), but not for WKT1_ESRI. Setting to YES will output
  * them unconditionally, and to NO will omit them unconditionally.</li>
+ * <li>STRICT=YES/NO. Default is YES. If NO, a Geographic 3D CRS can be for
+ * example exported as WKT1_GDAL with 3 axes, whereas this is normally not
+ * allowed.</li>
+ * <li>ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS=YES/NO. Default is NO. If set
+ * to YES and type == PJ_WKT1_GDAL, a Geographic 3D CRS or a Projected 3D CRS
+ * will be exported as a compound CRS whose vertical part represents an
+ * ellipsoidal height (for example for use with LAS 1.4 WKT1).</li>
  * </ul>
  * @return a string, or NULL in case of error.
  */
@@ -1468,6 +1478,11 @@ const char *proj_as_wkt(PJ_CONTEXT *ctx, const PJ *obj, PJ_WKT_TYPE type,
                 }
             } else if ((value = getOptionValue(*iter, "STRICT="))) {
                 formatter->setStrict(ci_equal(value, "YES"));
+            } else if ((value = getOptionValue(
+                            *iter,
+                            "ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS="))) {
+                formatter->setAllowEllipsoidalHeightAsVerticalCRS(
+                    ci_equal(value, "YES"));
             } else {
                 std::string msg("Unknown option :");
                 msg += *iter;
@@ -2050,6 +2065,8 @@ PJ *proj_get_ellipsoid(PJ_CONTEXT *ctx, const PJ *obj) {
 
 /** \brief Get the horizontal datum from a CRS
  *
+ * This function may return a Datum or DatumEnsemble object.
+ *
  * The returned object must be unreferenced with proj_destroy() after
  * use.
  * It should be used by at most one thread at a time.
@@ -2685,11 +2702,11 @@ proj_get_crs_info_list_from_database(PJ_CONTEXT *ctx, const char *auth_name,
 void proj_crs_info_list_destroy(PROJ_CRS_INFO **list) {
     if (list) {
         for (int i = 0; list[i] != nullptr; i++) {
-            pj_dalloc(list[i]->auth_name);
-            pj_dalloc(list[i]->code);
-            pj_dalloc(list[i]->name);
-            pj_dalloc(list[i]->area_name);
-            pj_dalloc(list[i]->projection_method_name);
+            free(list[i]->auth_name);
+            free(list[i]->code);
+            free(list[i]->name);
+            free(list[i]->area_name);
+            free(list[i]->projection_method_name);
             delete list[i];
         }
         delete[] list;
@@ -2779,11 +2796,11 @@ PROJ_UNIT_INFO **proj_get_units_from_database(PJ_CONTEXT *ctx,
 void proj_unit_list_destroy(PROJ_UNIT_INFO **list) {
     if (list) {
         for (int i = 0; list[i] != nullptr; i++) {
-            pj_dalloc(list[i]->auth_name);
-            pj_dalloc(list[i]->code);
-            pj_dalloc(list[i]->name);
-            pj_dalloc(list[i]->category);
-            pj_dalloc(list[i]->proj_short_name);
+            free(list[i]->auth_name);
+            free(list[i]->code);
+            free(list[i]->name);
+            free(list[i]->category);
+            free(list[i]->proj_short_name);
             delete list[i];
         }
         delete[] list;
@@ -4438,8 +4455,9 @@ PJ *proj_create_cartesian_2D_cs(PJ_CONTEXT *ctx, PJ_CARTESIAN_CS_2D_TYPE type,
  *
  * @param ctx PROJ context, or NULL for default context
  * @param type Coordinate system type.
- * @param unit_name Unit name.
- * @param unit_conv_factor Unit conversion factor to SI.
+ * @param unit_name Name of the angular units. Or NULL for Degree
+ * @param unit_conv_factor Conversion factor from the angular unit to radian.
+ * Or 0 for Degree if unit_name == NULL. Otherwise should be not NULL
  *
  * @return Object that must be unreferenced with
  * proj_destroy(), or NULL in case of error.
@@ -4478,13 +4496,17 @@ PJ *proj_create_ellipsoidal_2D_cs(PJ_CONTEXT *ctx,
  *
  * @param ctx PROJ context, or NULL for default context
  * @param type Coordinate system type.
- * @param horizontal_angular_unit_name Horizontal angular unit name.
- * @param horizontal_angular_unit_conv_factor Horizontal angular unit conversion
- * factor to SI.
- * @param vertical_linear_unit_name Vertical linear unit name.
+ * @param horizontal_angular_unit_name Name of the angular units. Or NULL for
+ * Degree.
+ * @param horizontal_angular_unit_conv_factor Conversion factor from the angular
+ * unit to radian. Or 0 for Degree if horizontal_angular_unit_name == NULL.
+ * Otherwise should be not NULL
+ * @param vertical_linear_unit_name Vertical linear unit name. Or NULL for
+ * Metre.
  * @param vertical_linear_unit_conv_factor Vertical linear unit conversion
- * factor to SI.
- *
+ * factor to metre. Or 0 for Metre if vertical_linear_unit_name == NULL.
+ * Otherwise should be not NULL
+
  * @return Object that must be unreferenced with
  * proj_destroy(), or NULL in case of error.
  * @since 6.3
@@ -7987,6 +8009,9 @@ double proj_coordoperation_get_accuracy(PJ_CONTEXT *ctx,
 
 /** \brief Returns the datum of a SingleCRS.
  *
+ * If that function returns NULL, @see proj_crs_get_datum_ensemble() to
+ * potentially get a DatumEnsemble instead.
+ *
  * The returned object must be unreferenced with proj_destroy() after
  * use.
  * It should be used by at most one thread at a time.
@@ -8018,6 +8043,9 @@ PJ *proj_crs_get_datum(PJ_CONTEXT *ctx, const PJ *crs) {
 
 /** \brief Returns the datum ensemble of a SingleCRS.
  *
+ * If that function returns NULL, @see proj_crs_get_datum() to
+ * potentially get a Datum instead.
+ *
  * The returned object must be unreferenced with proj_destroy() after
  * use.
  * It should be used by at most one thread at a time.
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index be15b3e0..83b626b3 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -64,6 +64,30 @@
 // #define DEBUG_CONCATENATED_OPERATION
 #if defined(DEBUG_SORT) || defined(DEBUG_CONCATENATED_OPERATION)
 #include <iostream>
+
+void dumpWKT(const NS_PROJ::crs::CRS *crs);
+void dumpWKT(const NS_PROJ::crs::CRS *crs) {
+    auto f(NS_PROJ::io::WKTFormatter::create(
+        NS_PROJ::io::WKTFormatter::Convention::WKT2_2019));
+    std::cerr << crs->exportToWKT(f.get()) << std::endl;
+}
+
+void dumpWKT(const NS_PROJ::crs::CRSPtr &crs);
+void dumpWKT(const NS_PROJ::crs::CRSPtr &crs) { dumpWKT(crs.get()); }
+
+void dumpWKT(const NS_PROJ::crs::CRSNNPtr &crs);
+void dumpWKT(const NS_PROJ::crs::CRSNNPtr &crs) {
+    dumpWKT(crs.as_nullable().get());
+}
+
+void dumpWKT(const NS_PROJ::crs::GeographicCRSPtr &crs);
+void dumpWKT(const NS_PROJ::crs::GeographicCRSPtr &crs) { dumpWKT(crs.get()); }
+
+void dumpWKT(const NS_PROJ::crs::GeographicCRSNNPtr &crs);
+void dumpWKT(const NS_PROJ::crs::GeographicCRSNNPtr &crs) {
+    dumpWKT(crs.as_nullable().get());
+}
+
 #endif
 
 using namespace NS_PROJ::internal;
@@ -6010,28 +6034,39 @@ void Conversion::_exportToPROJString(
         !isHeightDepthReversal;
     bool applyTargetCRSModifiers = applySourceCRSModifiers;
 
+    if (formatter->getCRSExport()) {
+        if (methodEPSGCode == EPSG_CODE_METHOD_GEOCENTRIC_TOPOCENTRIC ||
+            methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC) {
+            throw io::FormattingException("Transformation cannot be exported "
+                                          "as a PROJ.4 string (but can be part "
+                                          "of a PROJ pipeline)");
+        }
+    }
+
     auto l_sourceCRS = sourceCRS();
+    crs::GeographicCRSPtr srcGeogCRS;
     if (!formatter->getCRSExport() && l_sourceCRS && applySourceCRSModifiers) {
 
-        crs::CRS *horiz = l_sourceCRS.get();
-        const auto compound = dynamic_cast<const crs::CompoundCRS *>(horiz);
+        crs::CRSPtr horiz = l_sourceCRS;
+        const auto compound =
+            dynamic_cast<const crs::CompoundCRS *>(l_sourceCRS.get());
         if (compound) {
             const auto &components = compound->componentReferenceSystems();
             if (!components.empty()) {
-                horiz = components.front().get();
+                horiz = components.front().as_nullable();
             }
         }
 
-        auto geogCRS = dynamic_cast<const crs::GeographicCRS *>(horiz);
-        if (geogCRS) {
+        srcGeogCRS = std::dynamic_pointer_cast<crs::GeographicCRS>(horiz);
+        if (srcGeogCRS) {
             formatter->setOmitProjLongLatIfPossible(true);
             formatter->startInversion();
-            geogCRS->_exportToPROJString(formatter);
+            srcGeogCRS->_exportToPROJString(formatter);
             formatter->stopInversion();
             formatter->setOmitProjLongLatIfPossible(false);
         }
 
-        auto projCRS = dynamic_cast<const crs::ProjectedCRS *>(horiz);
+        auto projCRS = dynamic_cast<const crs::ProjectedCRS *>(horiz.get());
         if (projCRS) {
             formatter->startInversion();
             formatter->pushOmitZUnitConversion();
@@ -6301,6 +6336,30 @@ void Conversion::_exportToPROJString(
         }
         bConversionDone = true;
         bEllipsoidParametersDone = true;
+    } else if (methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC) {
+        if (!srcGeogCRS) {
+            throw io::FormattingException(
+                "Export of Geographic/Topocentric conversion to a PROJ string "
+                "requires an input geographic CRS");
+        }
+
+        formatter->addStep("cart");
+        srcGeogCRS->ellipsoid()->_exportToPROJString(formatter);
+
+        formatter->addStep("topocentric");
+        const auto latOrigin = parameterValueNumeric(
+            EPSG_CODE_PARAMETER_LATITUDE_TOPOGRAPHIC_ORIGIN,
+            common::UnitOfMeasure::DEGREE);
+        const auto lonOrigin = parameterValueNumeric(
+            EPSG_CODE_PARAMETER_LONGITUDE_TOPOGRAPHIC_ORIGIN,
+            common::UnitOfMeasure::DEGREE);
+        const auto heightOrigin = parameterValueNumeric(
+            EPSG_CODE_PARAMETER_ELLIPSOIDAL_HEIGHT_TOPOCENTRIC_ORIGIN,
+            common::UnitOfMeasure::METRE);
+        formatter->addParam("lat_0", latOrigin);
+        formatter->addParam("lon_0", lonOrigin);
+        formatter->addParam("h_0", heightOrigin);
+        bConversionDone = true;
     }
 
     auto l_targetCRS = targetCRS();
@@ -6425,7 +6484,9 @@ void Conversion::_exportToPROJString(
         }
 
         if (!bEllipsoidParametersDone) {
-            auto targetGeogCRS = horiz->extractGeographicCRS();
+            auto targetGeodCRS = horiz->extractGeodeticCRS();
+            auto targetGeogCRS =
+                std::dynamic_pointer_cast<crs::GeographicCRS>(targetGeodCRS);
             if (targetGeogCRS) {
                 if (formatter->getCRSExport()) {
                     targetGeogCRS->addDatumInfoToPROJString(formatter);
@@ -6434,6 +6495,8 @@ void Conversion::_exportToPROJString(
                     targetGeogCRS->primeMeridian()->_exportToPROJString(
                         formatter);
                 }
+            } else if (targetGeodCRS) {
+                targetGeodCRS->ellipsoid()->_exportToPROJString(formatter);
             }
         }
 
@@ -8836,11 +8899,16 @@ createSimilarPropertiesTransformation(TransformationNNPtr obj) {
     // The domain(s) are unchanged
     addDomains(map, obj.get());
 
-    std::string forwardName = obj->nameStr();
+    const std::string &forwardName = obj->nameStr();
     if (!forwardName.empty()) {
         map.set(common::IdentifiedObject::NAME_KEY, forwardName);
     }
 
+    const std::string &remarks = obj->remarks();
+    if (!remarks.empty()) {
+        map.set(common::IdentifiedObject::REMARKS_KEY, remarks);
+    }
+
     addModifiedIdentifier(map, obj.get(), false, true);
 
     return map;
@@ -9179,6 +9247,14 @@ static void setupPROJGeodeticSourceCRS(io::PROJStringFormatter *formatter,
         formatter->startInversion();
         sourceCRSGeog->_exportToPROJString(formatter);
         formatter->stopInversion();
+        if (util::isOfExactType<crs::DerivedGeographicCRS>(
+                *(sourceCRSGeog.get()))) {
+            // The export of a DerivedGeographicCRS in non-CRS mode adds
+            // unit conversion and axis swapping. We must compensate for that
+            formatter->startInversion();
+            sourceCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
+            formatter->stopInversion();
+        }
 
         if (addPushV3) {
             formatter->addStep("push");
@@ -9212,7 +9288,12 @@ static void setupPROJGeodeticTargetCRS(io::PROJStringFormatter *formatter,
             formatter->addStep("pop");
             formatter->addParam("v_3");
         }
-
+        if (util::isOfExactType<crs::DerivedGeographicCRS>(
+                *(targetCRSGeog.get()))) {
+            // The export of a DerivedGeographicCRS in non-CRS mode adds
+            // unit conversion and axis swapping. We must compensate for that
+            targetCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
+        }
         targetCRSGeog->_exportToPROJString(formatter);
     } else {
         auto targetCRSGeod = dynamic_cast<const crs::GeodeticCRS *>(crs.get());
@@ -11191,6 +11272,12 @@ struct CoordinateOperationFactory::Private {
         const crs::GeographicCRS *geogDst,
         std::vector<CoordinateOperationNNPtr> &res);
 
+    static void createOperationsVertToGeogBallpark(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::VerticalCRS *vertSrc,
+        const crs::GeographicCRS *geogDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
     static void createOperationsBoundToBound(
         const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
         Private::Context &context, const crs::BoundCRS *boundSrc,
@@ -11292,6 +11379,7 @@ struct PrecomputedOpCharacteristics {
     bool gridsKnown_ = false;
     size_t stepCount_ = 0;
     bool isApprox_ = false;
+    bool hasBallparkVertical_ = false;
     bool isNullTransformation_ = false;
 
     PrecomputedOpCharacteristics() = default;
@@ -11299,10 +11387,12 @@ struct PrecomputedOpCharacteristics {
                                  bool isPROJExportable, bool hasGrids,
                                  bool gridsAvailable, bool gridsKnown,
                                  size_t stepCount, bool isApprox,
+                                 bool hasBallparkVertical,
                                  bool isNullTransformation)
         : area_(area), accuracy_(accuracy), isPROJExportable_(isPROJExportable),
           hasGrids_(hasGrids), gridsAvailable_(gridsAvailable),
           gridsKnown_(gridsKnown), stepCount_(stepCount), isApprox_(isApprox),
+          hasBallparkVertical_(hasBallparkVertical),
           isNullTransformation_(isNullTransformation) {}
 };
 
@@ -11346,6 +11436,15 @@ struct SortFunction {
             return false;
         }
 
+        if (!iterA->second.hasBallparkVertical_ &&
+            iterB->second.hasBallparkVertical_) {
+            return true;
+        }
+        if (iterA->second.hasBallparkVertical_ &&
+            !iterB->second.hasBallparkVertical_) {
+            return false;
+        }
+
         if (!iterA->second.isNullTransformation_ &&
             iterB->second.isNullTransformation_) {
             return true;
@@ -11612,7 +11711,9 @@ struct FilterResults {
                 ? CoordinateOperationContext::SpatialCriterion::
                       STRICT_CONTAINMENT
                 : context->getSpatialCriterion();
-        bool hasFoundOpWithExtent = false;
+        bool hasOnlyBallpark = true;
+        bool hasNonBallparkWithoutExtent = false;
+        bool hasNonBallparkOpWithExtent = false;
         const bool allowBallpark = context->getAllowBallparkTransformations();
         for (const auto &op : sourceList) {
             if (desiredAccuracy != 0) {
@@ -11627,9 +11728,15 @@ struct FilterResults {
             if (areaOfInterest) {
                 bool emptyIntersection = false;
                 auto extent = getExtent(op, true, emptyIntersection);
-                if (!extent)
+                if (!extent) {
+                    if (!op->hasBallparkTransformation()) {
+                        hasNonBallparkWithoutExtent = true;
+                    }
                     continue;
-                hasFoundOpWithExtent = true;
+                }
+                if (!op->hasBallparkTransformation()) {
+                    hasNonBallparkOpWithExtent = true;
+                }
                 bool extentContains =
                     extent->contains(NN_NO_CHECK(areaOfInterest));
                 if (!hasOpThatContainsAreaOfInterestAndNoGrid &&
@@ -11656,9 +11763,15 @@ struct FilterResults {
                            BOTH) {
                 bool emptyIntersection = false;
                 auto extent = getExtent(op, true, emptyIntersection);
-                if (!extent)
+                if (!extent) {
+                    if (!op->hasBallparkTransformation()) {
+                        hasNonBallparkWithoutExtent = true;
+                    }
                     continue;
-                hasFoundOpWithExtent = true;
+                }
+                if (!op->hasBallparkTransformation()) {
+                    hasNonBallparkOpWithExtent = true;
+                }
                 bool extentContainsExtent1 =
                     !extent1 || extent->contains(NN_NO_CHECK(extent1));
                 bool extentContainsExtent2 =
@@ -11688,12 +11801,16 @@ struct FilterResults {
                     }
                 }
             }
+            if (!op->hasBallparkTransformation()) {
+                hasOnlyBallpark = false;
+            }
             res.emplace_back(op);
         }
 
         // In case no operation has an extent and no result is found,
         // retain all initial operations that match accuracy criterion.
-        if (res.empty() && !hasFoundOpWithExtent) {
+        if ((res.empty() && !hasNonBallparkOpWithExtent) ||
+            (hasOnlyBallpark && hasNonBallparkWithoutExtent)) {
             for (const auto &op : sourceList) {
                 if (desiredAccuracy != 0) {
                     const double accuracy = getAccuracy(op);
@@ -11797,6 +11914,8 @@ struct FilterResults {
                 area, getAccuracy(op), isPROJExportable, hasGrids,
                 gridsAvailable, gridsKnown, stepCount,
                 op->hasBallparkTransformation(),
+                op->nameStr().find("ballpark vertical transformation") !=
+                    std::string::npos,
                 isNullTransformation(op->nameStr()));
         }
 
@@ -12615,6 +12734,55 @@ createGeodToGeodPROJBased(const crs::CRSNNPtr &geodSrc,
 
 // ---------------------------------------------------------------------------
 
+static std::string
+getRemarks(const std::vector<operation::CoordinateOperationNNPtr> &ops) {
+    std::string remarks;
+    for (const auto &op : ops) {
+        const auto &opRemarks = op->remarks();
+        if (!opRemarks.empty()) {
+            if (!remarks.empty()) {
+                remarks += '\n';
+            }
+
+            std::string opName(op->nameStr());
+            if (starts_with(opName, INVERSE_OF)) {
+                opName = opName.substr(INVERSE_OF.size());
+            }
+
+            remarks += "For ";
+            remarks += opName;
+
+            const auto &ids = op->identifiers();
+            if (!ids.empty()) {
+                std::string authority(*ids.front()->codeSpace());
+                if (starts_with(authority, "INVERSE(") &&
+                    authority.back() == ')') {
+                    authority = authority.substr(strlen("INVERSE("),
+                                                 authority.size() - 1 -
+                                                     strlen("INVERSE("));
+                }
+                if (starts_with(authority, "DERIVED_FROM(") &&
+                    authority.back() == ')') {
+                    authority = authority.substr(strlen("DERIVED_FROM("),
+                                                 authority.size() - 1 -
+                                                     strlen("DERIVED_FROM("));
+                }
+
+                remarks += " (";
+                remarks += authority;
+                remarks += ':';
+                remarks += ids.front()->code();
+                remarks += ')';
+            }
+            remarks += ": ";
+            remarks += opRemarks;
+        }
+    }
+    return remarks;
+}
+
+// ---------------------------------------------------------------------------
+
 static CoordinateOperationNNPtr createHorizVerticalPROJBased(
     const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
     const operation::CoordinateOperationNNPtr &horizTransform,
@@ -12640,6 +12808,10 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased(
             properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
                            NN_NO_CHECK(extent));
         }
+        const auto &remarks = verticalTransform->remarks();
+        if (!remarks.empty()) {
+            properties.set(common::IdentifiedObject::REMARKS_KEY, remarks);
+        }
         return createPROJBased(
             properties, exportable, sourceCRS, targetCRS, nullptr,
             verticalTransform->coordinateOperationAccuracies(),
@@ -12664,6 +12836,11 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased(
                            NN_NO_CHECK(extent));
         }
 
+        const auto remarks = getRemarks(ops);
+        if (!remarks.empty()) {
+            properties.set(common::IdentifiedObject::REMARKS_KEY, remarks);
+        }
+
         std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
         const double accuracy = getAccuracy(ops);
         if (accuracy >= 0.0) {
@@ -12724,6 +12901,11 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased(
                        NN_NO_CHECK(extent));
     }
 
+    const auto remarks = getRemarks(ops);
+    if (!remarks.empty()) {
+        properties.set(common::IdentifiedObject::REMARKS_KEY, remarks);
+    }
+
     std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
     const double accuracy = getAccuracy(ops);
     if (accuracy >= 0.0) {
@@ -13542,11 +13724,17 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
     ENTER_FUNCTION();
 
     if (geogSrc && vertDst) {
-        res = createOperationsGeogToVertFromGeoid(sourceCRS, targetCRS, vertDst,
-                                                  context);
+        createOperationsFromDatabase(targetCRS, sourceCRS, context, geodDst,
+                                     geodSrc, geogDst, geogSrc, vertDst,
+                                     vertSrc, res);
+        res = applyInverse(res);
     } else if (geogDst && vertSrc) {
         res = applyInverse(createOperationsGeogToVertFromGeoid(
             targetCRS, sourceCRS, vertSrc, context));
+        if (!res.empty()) {
+            createOperationsVertToGeogBallpark(sourceCRS, targetCRS, context,
+                                               vertSrc, geogDst, res);
+        }
     }
 
     if (!res.empty()) {
@@ -14222,6 +14410,20 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog(
     const auto &hubSrc = boundSrc->hubCRS();
     auto hubSrcGeog = dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
     auto geogCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractGeographicCRS();
+    {
+        // If geogCRSOfBaseOfBoundSrc is a DerivedGeographicCRS, use its base
+        // instead (if it is a GeographicCRS)
+        auto derivedGeogCRS =
+            std::dynamic_pointer_cast<crs::DerivedGeographicCRS>(
+                geogCRSOfBaseOfBoundSrc);
+        if (derivedGeogCRS) {
+            auto baseCRS = std::dynamic_pointer_cast<crs::GeographicCRS>(
+                derivedGeogCRS->baseCRS().as_nullable());
+            if (baseCRS) {
+                geogCRSOfBaseOfBoundSrc = baseCRS;
+            }
+        }
+    }
 
     const auto &authFactory = context.context->getAuthorityFactory();
     const auto dbContext =
@@ -14643,23 +14845,38 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeog(
                         match.get(),
                         util::IComparable::Criterion::EQUIVALENT) &&
                     !match->identifiers().empty()) {
-                    res = createOperations(
+                    auto resTmp = createOperations(
                         NN_NO_CHECK(
                             util::nn_dynamic_pointer_cast<crs::VerticalCRS>(
                                 match)),
                         targetCRS, context);
+                    res.insert(res.end(), resTmp.begin(), resTmp.end());
                     return;
                 }
             }
         }
     }
 
+    createOperationsVertToGeogBallpark(sourceCRS, targetCRS, context, vertSrc,
+                                       geogDst, res);
+}
+
+// ---------------------------------------------------------------------------
+
+void CoordinateOperationFactory::Private::createOperationsVertToGeogBallpark(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    Private::Context &, const crs::VerticalCRS *vertSrc,
+    const crs::GeographicCRS *geogDst,
+    std::vector<CoordinateOperationNNPtr> &res) {
+
+    ENTER_FUNCTION();
+
     const auto &srcAxis = vertSrc->coordinateSystem()->axisList()[0];
     const double convSrc = srcAxis->unit().conversionToSI();
     double convDst = 1.0;
     const auto &geogAxis = geogDst->coordinateSystem()->axisList();
     bool dstIsUp = true;
-    bool dstIsDown = true;
+    bool dstIsDown = false;
     if (geogAxis.size() == 3) {
         const auto &dstAxis = geogAxis[2];
         convDst = dstAxis->unit().conversionToSI();
@@ -14672,12 +14889,24 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeog(
         ((srcIsUp && dstIsDown) || (srcIsDown && dstIsUp));
 
     const double factor = convSrc / convDst;
-    auto conv = Transformation::createChangeVerticalUnit(
-        util::PropertyMap().set(
-            common::IdentifiedObject::NAME_KEY,
+
+    const auto &sourceCRSExtent = getExtent(sourceCRS);
+    const auto &targetCRSExtent = getExtent(targetCRS);
+    const bool sameExtent =
+        sourceCRSExtent && targetCRSExtent &&
+        sourceCRSExtent->_isEquivalentTo(
+            targetCRSExtent.get(), util::IComparable::Criterion::EQUIVALENT);
+
+    util::PropertyMap map;
+    map.set(common::IdentifiedObject::NAME_KEY,
             buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) +
-                BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT),
-        sourceCRS, targetCRS,
+                BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT)
+        .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
+             sameExtent ? NN_NO_CHECK(sourceCRSExtent)
+                        : metadata::Extent::WORLD);
+
+    auto conv = Transformation::createChangeVerticalUnit(
+        map, sourceCRS, targetCRS,
         common::Scale(heightDepthReversal ? -factor : factor), {});
     conv->setHasBallparkTransformation(true);
     res.push_back(conv);
@@ -15477,149 +15706,6 @@ void CoordinateOperationFactory::Private::createOperationsBoundToCompound(
 
 // ---------------------------------------------------------------------------
 
-static crs::CRSNNPtr
-getResolvedCRS(const crs::CRSNNPtr &crs,
-               const CoordinateOperationContextNNPtr &context,
-               metadata::ExtentPtr &extentOut) {
-    const auto &authFactory = context->getAuthorityFactory();
-    const auto &ids = crs->identifiers();
-    const auto &name = crs->nameStr();
-
-    bool approxExtent;
-    extentOut = getExtentPossiblySynthetized(crs, approxExtent);
-
-    // We try to "identify" the provided CRS with the ones of the database,
-    // but in a more restricted way that what identify() does.
-    // If we get a match from id in priority, and from name as a fallback, and
-    // that they are equivalent to the input CRS, then use the identified CRS.
-    // Even if they aren't equivalent, we update extentOut with the one of the
-    // identified CRS if our input one is absent/not reliable.
-
-    const auto tryToIdentifyByName = [&crs, &name, &authFactory, approxExtent,
-                                      &extentOut](
-        io::AuthorityFactory::ObjectType objectType) {
-        if (name != "unknown" && name != "unnamed") {
-            auto matches = authFactory->createObjectsFromName(
-                name, {objectType}, false, 2);
-            if (matches.size() == 1) {
-                const auto match =
-                    util::nn_static_pointer_cast<crs::CRS>(matches.front());
-                if (approxExtent || !extentOut) {
-                    extentOut = getExtent(match);
-                }
-                if (match->isEquivalentTo(
-                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
-                    return match;
-                }
-            }
-        }
-        return crs;
-    };
-
-    auto geogCRS = dynamic_cast<crs::GeographicCRS *>(crs.get());
-    if (geogCRS && authFactory) {
-        if (!ids.empty()) {
-            const auto tmpAuthFactory = io::AuthorityFactory::create(
-                authFactory->databaseContext(), *ids.front()->codeSpace());
-            try {
-                auto resolvedCrs(
-                    tmpAuthFactory->createGeographicCRS(ids.front()->code()));
-                if (approxExtent || !extentOut) {
-                    extentOut = getExtent(resolvedCrs);
-                }
-                if (resolvedCrs->isEquivalentTo(
-                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
-                    return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
-                }
-            } catch (const std::exception &) {
-            }
-        } else {
-            return tryToIdentifyByName(
-                geogCRS->coordinateSystem()->axisList().size() == 2
-                    ? io::AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS
-                    : io::AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS);
-        }
-    }
-
-    auto projectedCrs = dynamic_cast<crs::ProjectedCRS *>(crs.get());
-    if (projectedCrs && authFactory) {
-        if (!ids.empty()) {
-            const auto tmpAuthFactory = io::AuthorityFactory::create(
-                authFactory->databaseContext(), *ids.front()->codeSpace());
-            try {
-                auto resolvedCrs(
-                    tmpAuthFactory->createProjectedCRS(ids.front()->code()));
-                if (approxExtent || !extentOut) {
-                    extentOut = getExtent(resolvedCrs);
-                }
-                if (resolvedCrs->isEquivalentTo(
-                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
-                    return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
-                }
-            } catch (const std::exception &) {
-            }
-        } else {
-            return tryToIdentifyByName(
-                io::AuthorityFactory::ObjectType::PROJECTED_CRS);
-        }
-    }
-
-    auto compoundCrs = dynamic_cast<crs::CompoundCRS *>(crs.get());
-    if (compoundCrs && authFactory) {
-        if (!ids.empty()) {
-            const auto tmpAuthFactory = io::AuthorityFactory::create(
-                authFactory->databaseContext(), *ids.front()->codeSpace());
-            try {
-                auto resolvedCrs(
-                    tmpAuthFactory->createCompoundCRS(ids.front()->code()));
-                if (approxExtent || !extentOut) {
-                    extentOut = getExtent(resolvedCrs);
-                }
-                if (resolvedCrs->isEquivalentTo(
-                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
-                    return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
-                }
-            } catch (const std::exception &) {
-            }
-        } else {
-            auto outCrs = tryToIdentifyByName(
-                io::AuthorityFactory::ObjectType::COMPOUND_CRS);
-            const auto &components = compoundCrs->componentReferenceSystems();
-            if (outCrs.get() != crs.get()) {
-                bool hasGeoid = false;
-                if (components.size() == 2) {
-                    auto vertCRS =
-                        dynamic_cast<crs::VerticalCRS *>(components[1].get());
-                    if (vertCRS && !vertCRS->geoidModel().empty()) {
-                        hasGeoid = true;
-                    }
-                }
-                if (!hasGeoid) {
-                    return outCrs;
-                }
-            }
-            if (approxExtent || !extentOut) {
-                // If we still did not get a reliable extent, then try to
-                // resolve the components of the compoundCRS, and take the
-                // intersection of their extent.
-                extentOut = metadata::ExtentPtr();
-                for (const auto &component : components) {
-                    metadata::ExtentPtr componentExtent;
-                    getResolvedCRS(component, context, componentExtent);
-                    if (extentOut && componentExtent)
-                        extentOut = extentOut->intersection(
-                            NN_NO_CHECK(componentExtent));
-                    else if (componentExtent)
-                        extentOut = componentExtent;
-                }
-            }
-        }
-    }
-    return crs;
-}
-
-// ---------------------------------------------------------------------------
-
 /** \brief Find a list of CoordinateOperation from sourceCRS to targetCRS.
  *
  * The operations are sorted with the most relevant ones first: by
@@ -15655,13 +15741,14 @@ CoordinateOperationFactory::createOperations(
     const auto &targetBoundCRS = targetCRS->canonicalBoundCRS();
     auto l_sourceCRS = srcBoundCRS ? NN_NO_CHECK(srcBoundCRS) : sourceCRS;
     auto l_targetCRS = targetBoundCRS ? NN_NO_CHECK(targetBoundCRS) : targetCRS;
+    const auto &authFactory = context->getAuthorityFactory();
 
     metadata::ExtentPtr sourceCRSExtent;
     auto l_resolvedSourceCRS =
-        getResolvedCRS(l_sourceCRS, context, sourceCRSExtent);
+        crs::CRS::getResolvedCRS(l_sourceCRS, authFactory, sourceCRSExtent);
     metadata::ExtentPtr targetCRSExtent;
     auto l_resolvedTargetCRS =
-        getResolvedCRS(l_targetCRS, context, targetCRSExtent);
+        crs::CRS::getResolvedCRS(l_targetCRS, authFactory, targetCRSExtent);
     Private::Context contextPrivate(sourceCRSExtent, targetCRSExtent, context);
 
     if (context->getSourceAndTargetCRSExtentUse() ==
@@ -15986,4 +16073,152 @@ PROJBasedOperation::gridsNeeded(const io::DatabaseContextPtr &databaseContext,
 // ---------------------------------------------------------------------------
 
 } // namespace operation
+
+namespace crs {
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
+
+crs::CRSNNPtr CRS::getResolvedCRS(const crs::CRSNNPtr &crs,
+                                  const io::AuthorityFactoryPtr &authFactory,
+                                  metadata::ExtentPtr &extentOut) {
+    const auto &ids = crs->identifiers();
+    const auto &name = crs->nameStr();
+
+    bool approxExtent;
+    extentOut = getExtentPossiblySynthetized(crs, approxExtent);
+
+    // We try to "identify" the provided CRS with the ones of the database,
+    // but in a more restricted way that what identify() does.
+    // If we get a match from id in priority, and from name as a fallback, and
+    // that they are equivalent to the input CRS, then use the identified CRS.
+    // Even if they aren't equivalent, we update extentOut with the one of the
+    // identified CRS if our input one is absent/not reliable.
+
+    const auto tryToIdentifyByName = [&crs, &name, &authFactory, approxExtent,
+                                      &extentOut](
+        io::AuthorityFactory::ObjectType objectType) {
+        if (name != "unknown" && name != "unnamed") {
+            auto matches = authFactory->createObjectsFromName(
+                name, {objectType}, false, 2);
+            if (matches.size() == 1) {
+                const auto match =
+                    util::nn_static_pointer_cast<crs::CRS>(matches.front());
+                if (approxExtent || !extentOut) {
+                    extentOut = getExtent(match);
+                }
+                if (match->isEquivalentTo(
+                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
+                    return match;
+                }
+            }
+        }
+        return crs;
+    };
+
+    auto geogCRS = dynamic_cast<crs::GeographicCRS *>(crs.get());
+    if (geogCRS && authFactory) {
+        if (!ids.empty()) {
+            const auto tmpAuthFactory = io::AuthorityFactory::create(
+                authFactory->databaseContext(), *ids.front()->codeSpace());
+            try {
+                auto resolvedCrs(
+                    tmpAuthFactory->createGeographicCRS(ids.front()->code()));
+                if (approxExtent || !extentOut) {
+                    extentOut = getExtent(resolvedCrs);
+                }
+                if (resolvedCrs->isEquivalentTo(
+                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
+                    return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
+                }
+            } catch (const std::exception &) {
+            }
+        } else {
+            return tryToIdentifyByName(
+                geogCRS->coordinateSystem()->axisList().size() == 2
+                    ? io::AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS
+                    : io::AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS);
+        }
+    }
+
+    auto projectedCrs = dynamic_cast<crs::ProjectedCRS *>(crs.get());
+    if (projectedCrs && authFactory) {
+        if (!ids.empty()) {
+            const auto tmpAuthFactory = io::AuthorityFactory::create(
+                authFactory->databaseContext(), *ids.front()->codeSpace());
+            try {
+                auto resolvedCrs(
+                    tmpAuthFactory->createProjectedCRS(ids.front()->code()));
+                if (approxExtent || !extentOut) {
+                    extentOut = getExtent(resolvedCrs);
+                }
+                if (resolvedCrs->isEquivalentTo(
+                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
+                    return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
+                }
+            } catch (const std::exception &) {
+            }
+        } else {
+            return tryToIdentifyByName(
+                io::AuthorityFactory::ObjectType::PROJECTED_CRS);
+        }
+    }
+
+    auto compoundCrs = dynamic_cast<crs::CompoundCRS *>(crs.get());
+    if (compoundCrs && authFactory) {
+        if (!ids.empty()) {
+            const auto tmpAuthFactory = io::AuthorityFactory::create(
+                authFactory->databaseContext(), *ids.front()->codeSpace());
+            try {
+                auto resolvedCrs(
+                    tmpAuthFactory->createCompoundCRS(ids.front()->code()));
+                if (approxExtent || !extentOut) {
+                    extentOut = getExtent(resolvedCrs);
+                }
+                if (resolvedCrs->isEquivalentTo(
+                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
+                    return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
+                }
+            } catch (const std::exception &) {
+            }
+        } else {
+            auto outCrs = tryToIdentifyByName(
+                io::AuthorityFactory::ObjectType::COMPOUND_CRS);
+            const auto &components = compoundCrs->componentReferenceSystems();
+            if (outCrs.get() != crs.get()) {
+                bool hasGeoid = false;
+                if (components.size() == 2) {
+                    auto vertCRS =
+                        dynamic_cast<crs::VerticalCRS *>(components[1].get());
+                    if (vertCRS && !vertCRS->geoidModel().empty()) {
+                        hasGeoid = true;
+                    }
+                }
+                if (!hasGeoid) {
+                    return outCrs;
+                }
+            }
+            if (approxExtent || !extentOut) {
+                // If we still did not get a reliable extent, then try to
+                // resolve the components of the compoundCRS, and take the
+                // intersection of their extent.
+                extentOut = metadata::ExtentPtr();
+                for (const auto &component : components) {
+                    metadata::ExtentPtr componentExtent;
+                    getResolvedCRS(component, authFactory, componentExtent);
+                    if (extentOut && componentExtent)
+                        extentOut = extentOut->intersection(
+                            NN_NO_CHECK(componentExtent));
+                    else if (componentExtent)
+                        extentOut = componentExtent;
+                }
+            }
+        }
+    }
+    return crs;
+}
+
+//! @endcond
+
+} // namespace crs
 NS_PROJ_END
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index edc8a71f..573dd6db 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -533,8 +533,12 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
             auto authFactory = io::AuthorityFactory::create(
                 NN_NO_CHECK(dbContext),
                 authority == "any" ? std::string() : authority);
+            metadata::ExtentPtr extentResolved(extent);
+            if (!extent) {
+                getResolvedCRS(thisAsCRS, authFactory, extentResolved);
+            }
             auto ctxt = operation::CoordinateOperationContext::create(
-                authFactory, extent, 0.0);
+                authFactory, extentResolved, 0.0);
             ctxt->setAllowUseIntermediateCRS(allowIntermediateCRSUse);
             // ctxt->setSpatialCriterion(
             //    operation::CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
@@ -1616,6 +1620,34 @@ static bool exportAsESRIWktCompoundCRSWithEllipsoidalHeight(
     vertCRSList.front()->_exportToWKT(formatter);
     return true;
 }
+
+// ---------------------------------------------------------------------------
+
+// Try to format a Geographic/ProjectedCRS 3D CRS as a
+// GEOGCS[]/PROJCS[],VERTCS["Ellipsoid (metre)",DATUM["Ellipsoid",2002],...]
+static bool exportAsWKT1CompoundCRSWithEllipsoidalHeight(
+    const CRSNNPtr &base2DCRS,
+    const cs::CoordinateSystemAxisNNPtr &verticalAxis,
+    io::WKTFormatter *formatter) {
+    std::string verticalCRSName = "Ellipsoid (";
+    verticalCRSName += verticalAxis->unit().name();
+    verticalCRSName += ')';
+    auto vertDatum = datum::VerticalReferenceFrame::create(
+        util::PropertyMap()
+            .set(common::IdentifiedObject::NAME_KEY, "Ellipsoid")
+            .set("VERT_DATUM_TYPE", "2002"));
+    auto vertCRS = VerticalCRS::create(
+        util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                verticalCRSName),
+        vertDatum.as_nullable(), nullptr,
+        cs::VerticalCS::create(util::PropertyMap(), verticalAxis));
+    formatter->startNode(io::WKTConstants::COMPD_CS, false);
+    formatter->addQuotedString(base2DCRS->nameStr() + " + " + verticalCRSName);
+    base2DCRS->_exportToWKT(formatter);
+    vertCRS->_exportToWKT(formatter);
+    formatter->endNode();
+    return true;
+}
 //! @endcond
 
 // ---------------------------------------------------------------------------
@@ -1683,6 +1715,13 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const {
             return;
         }
 
+        if (formatter->isAllowedEllipsoidalHeightAsVerticalCRS()) {
+            if (exportAsWKT1CompoundCRSWithEllipsoidalHeight(
+                    geogCRS2D, axisList[2], formatter)) {
+                return;
+            }
+        }
+
         io::FormattingException::Throw(
             "WKT1 does not support Geographic 3D CRS.");
     }
@@ -1922,11 +1961,19 @@ getStandardCriterion(util::IComparable::Criterion criterion) {
 bool GeodeticCRS::_isEquivalentTo(
     const util::IComparable *other, util::IComparable::Criterion criterion,
     const io::DatabaseContextPtr &dbContext) const {
+    if (other == nullptr || !util::isOfExactType<GeodeticCRS>(*other)) {
+        return false;
+    }
+    return _isEquivalentToNoTypeCheck(other, criterion, dbContext);
+}
+
+bool GeodeticCRS::_isEquivalentToNoTypeCheck(
+    const util::IComparable *other, util::IComparable::Criterion criterion,
+    const io::DatabaseContextPtr &dbContext) const {
     const auto standardCriterion = getStandardCriterion(criterion);
-    auto otherGeodCRS = dynamic_cast<const GeodeticCRS *>(other);
+
     // TODO test velocityModel
-    return otherGeodCRS != nullptr &&
-           SingleCRS::baseIsEquivalentTo(other, standardCriterion, dbContext);
+    return SingleCRS::baseIsEquivalentTo(other, standardCriterion, dbContext);
 }
 //! @endcond
 
@@ -2482,12 +2529,13 @@ bool GeographicCRS::is2DPartOf3D(util::nn<const GeographicCRS *> other,
 bool GeographicCRS::_isEquivalentTo(
     const util::IComparable *other, util::IComparable::Criterion criterion,
     const io::DatabaseContextPtr &dbContext) const {
-    auto otherGeogCRS = dynamic_cast<const GeographicCRS *>(other);
-    if (otherGeogCRS == nullptr) {
+    if (other == nullptr || !util::isOfExactType<GeographicCRS>(*other)) {
         return false;
     }
+
     const auto standardCriterion = getStandardCriterion(criterion);
-    if (GeodeticCRS::_isEquivalentTo(other, standardCriterion, dbContext)) {
+    if (GeodeticCRS::_isEquivalentToNoTypeCheck(other, standardCriterion,
+                                                dbContext)) {
         return true;
     }
     if (criterion !=
@@ -2506,7 +2554,29 @@ bool GeographicCRS::_isEquivalentTo(
                            cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH
                        ? cs::EllipsoidalCS::createLatitudeLongitude(unit)
                        : cs::EllipsoidalCS::createLongitudeLatitude(unit))
-            ->GeodeticCRS::_isEquivalentTo(other, standardCriterion, dbContext);
+            ->GeodeticCRS::_isEquivalentToNoTypeCheck(other, standardCriterion,
+                                                      dbContext);
+    }
+    if (axisOrder ==
+            cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH_HEIGHT_UP ||
+        axisOrder ==
+            cs::EllipsoidalCS::AxisOrder::LAT_NORTH_LONG_EAST_HEIGHT_UP) {
+        const auto &angularUnit = coordinateSystem()->axisList()[0]->unit();
+        const auto &linearUnit = coordinateSystem()->axisList()[2]->unit();
+        return GeographicCRS::create(
+                   util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                           nameStr()),
+                   datum(), datumEnsemble(),
+                   axisOrder == cs::EllipsoidalCS::AxisOrder::
+                                    LONG_EAST_LAT_NORTH_HEIGHT_UP
+                       ? cs::EllipsoidalCS::
+                             createLatitudeLongitudeEllipsoidalHeight(
+                                 angularUnit, linearUnit)
+                       : cs::EllipsoidalCS::
+                             createLongitudeLatitudeEllipsoidalHeight(
+                                 angularUnit, linearUnit))
+            ->GeodeticCRS::_isEquivalentToNoTypeCheck(other, standardCriterion,
+                                                      dbContext);
     }
     return false;
 }
@@ -3605,6 +3675,14 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const {
             return;
         }
 
+        if (!formatter->useESRIDialect() &&
+            formatter->isAllowedEllipsoidalHeightAsVerticalCRS()) {
+            if (exportAsWKT1CompoundCRSWithEllipsoidalHeight(
+                    projCRS2D, axisList[2], formatter)) {
+                return;
+            }
+        }
+
         io::FormattingException::Throw(
             "Projected 3D CRS can only be exported since WKT2:2019");
     }
@@ -3885,8 +3963,7 @@ ProjectedCRS::create(const util::PropertyMap &properties,
 bool ProjectedCRS::_isEquivalentTo(
     const util::IComparable *other, util::IComparable::Criterion criterion,
     const io::DatabaseContextPtr &dbContext) const {
-    auto otherProjCRS = dynamic_cast<const ProjectedCRS *>(other);
-    return otherProjCRS != nullptr &&
+    return other != nullptr && util::isOfExactType<ProjectedCRS>(*other) &&
            DerivedCRS::_isEquivalentTo(other, criterion, dbContext);
 }
 
diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp
index 5bc8074c..8e1b8257 100644
--- a/src/iso19111/datum.cpp
+++ b/src/iso19111/datum.cpp
@@ -43,7 +43,6 @@
 // clang-format off
 #include "proj.h"
 #include "proj_internal.h"
-#include "proj_api.h"
 // clang-format on
 
 #include "proj_json_streaming_writer.hpp"
@@ -352,6 +351,23 @@ void PrimeMeridian::_exportToWKT(
     if (!(isWKT2 && formatter->primeMeridianOmittedIfGreenwich() &&
           l_name == "Greenwich")) {
         formatter->startNode(io::WKTConstants::PRIMEM, !identifiers().empty());
+
+        if (formatter->useESRIDialect()) {
+            bool aliasFound = false;
+            const auto &dbContext = formatter->databaseContext();
+            if (dbContext) {
+                auto l_alias = dbContext->getAliasFromOfficialName(
+                    l_name, "prime_meridian", "ESRI");
+                if (!l_alias.empty()) {
+                    l_name = l_alias;
+                    aliasFound = true;
+                }
+            }
+            if (!aliasFound) {
+                l_name = io::WKTFormatter::morphNameToESRI(l_name);
+            }
+        }
+
         formatter->addQuotedString(l_name);
         const auto &l_long = longitude();
         if (formatter->primeMeridianInDegree()) {
@@ -418,7 +434,7 @@ std::string
 PrimeMeridian::getPROJStringWellKnownName(const common::Angle &angle) {
     const double valRad = angle.getSIValue();
     std::string projPMName;
-    projCtx ctxt = pj_ctx_alloc();
+    PJ_CONTEXT *ctxt = proj_context_create();
     auto proj_pm = proj_list_prime_meridians();
     for (int i = 0; proj_pm[i].id != nullptr; ++i) {
         double valRefRad = dmstor_ctx(ctxt, proj_pm[i].defn, nullptr);
@@ -427,7 +443,7 @@ PrimeMeridian::getPROJStringWellKnownName(const common::Angle &angle) {
             break;
         }
     }
-    pj_ctx_free(ctxt);
+    proj_context_destroy(ctxt);
     return projPMName;
 }
 //! @endcond
@@ -1731,6 +1747,14 @@ void DatumEnsemble::_exportToWKT(
     formatter->startNode(io::WKTConstants::ENSEMBLEACCURACY, false);
     formatter->add(positionalAccuracy()->value());
     formatter->endNode();
+
+    // In theory, we should do the following, but currently the WKT grammar
+    // doesn't allow this
+    // ObjectUsage::baseExportToWKT(formatter);
+    if (formatter->outputId()) {
+        formatID(formatter);
+    }
+
     formatter->endNode();
 }
 //! @endcond
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index 5d02aeea..2a03fd4e 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -65,7 +65,6 @@
 // clang-format off
 #include "proj.h"
 #include "proj_internal.h"
-#include "proj_api.h"
 // clang-format on
 
 #include <sqlite3.h>
@@ -94,6 +93,12 @@ namespace io {
 #define GEOG_3D_SINGLE_QUOTED "'geographic 3D'"
 #define GEOCENTRIC_SINGLE_QUOTED "'geocentric'"
 
+// See data/sql/metadata.sql for the semantics of those constants
+constexpr int DATABASE_LAYOUT_VERSION_MAJOR = 1;
+// If the code depends on the new additions, then DATABASE_LAYOUT_VERSION_MINOR
+// must be incremented.
+constexpr int DATABASE_LAYOUT_VERSION_MINOR = 0;
+
 // ---------------------------------------------------------------------------
 
 struct SQLValues {
@@ -192,6 +197,13 @@ struct DatabaseContext::Private {
     void cache(const std::string &code,
                const datum::GeodeticReferenceFrameNNPtr &datum);
 
+    datum::DatumEnsemblePtr
+        // cppcheck-suppress functionStatic
+        getDatumEnsembleFromCache(const std::string &code);
+    // cppcheck-suppress functionStatic
+    void cache(const std::string &code,
+               const datum::DatumEnsembleNNPtr &datumEnsemble);
+
     datum::EllipsoidPtr
         // cppcheck-suppress functionStatic
         getEllipsoidFromCache(const std::string &code);
@@ -258,6 +270,7 @@ struct DatabaseContext::Private {
     LRUCacheOfObjects cacheCRS_{CACHE_SIZE};
     LRUCacheOfObjects cacheEllipsoid_{CACHE_SIZE};
     LRUCacheOfObjects cacheGeodeticDatum_{CACHE_SIZE};
+    LRUCacheOfObjects cacheDatumEnsemble_{CACHE_SIZE};
     LRUCacheOfObjects cachePrimeMeridian_{CACHE_SIZE};
     LRUCacheOfObjects cacheCS_{CACHE_SIZE};
     LRUCacheOfObjects cacheExtent_{CACHE_SIZE};
@@ -270,6 +283,8 @@ struct DatabaseContext::Private {
     lru11::Cache<std::string, std::list<std::string>> cacheAliasNames_{
         CACHE_SIZE};
 
+    void checkDatabaseLayout();
+
     static void insertIntoCache(LRUCacheOfObjects &cache,
                                 const std::string &code,
                                 const util::BaseObjectPtr &obj);
@@ -417,6 +432,22 @@ void DatabaseContext::Private::cache(
 
 // ---------------------------------------------------------------------------
 
+datum::DatumEnsemblePtr
+DatabaseContext::Private::getDatumEnsembleFromCache(const std::string &code) {
+    util::BaseObjectPtr obj;
+    getFromCache(cacheDatumEnsemble_, code, obj);
+    return std::static_pointer_cast<datum::DatumEnsemble>(obj);
+}
+
+// ---------------------------------------------------------------------------
+
+void DatabaseContext::Private::cache(
+    const std::string &code, const datum::DatumEnsembleNNPtr &datumEnsemble) {
+    insertIntoCache(cacheDatumEnsemble_, code, datumEnsemble.as_nullable());
+}
+
+// ---------------------------------------------------------------------------
+
 datum::EllipsoidPtr
 DatabaseContext::Private::getEllipsoidFromCache(const std::string &code) {
     util::BaseObjectPtr obj;
@@ -548,6 +579,61 @@ void DatabaseContext::Private::open(const std::string &databasePath,
 
 // ---------------------------------------------------------------------------
 
+void DatabaseContext::Private::checkDatabaseLayout() {
+    auto res = run("SELECT key, value FROM metadata WHERE key IN "
+                   "('DATABASE.LAYOUT.VERSION.MAJOR', "
+                   "'DATABASE.LAYOUT.VERSION.MINOR')");
+    if (res.size() != 2) {
+        // The database layout of PROJ 7.2 that shipped with EPSG v10.003 is
+        // at the time of writing still compatible of the one we support.
+        static_assert(
+            // cppcheck-suppress knownConditionTrueFalse
+            DATABASE_LAYOUT_VERSION_MAJOR == 1 &&
+                // cppcheck-suppress knownConditionTrueFalse
+                DATABASE_LAYOUT_VERSION_MINOR == 0,
+            "remove that assertion and below lines next time we upgrade "
+            "database structure");
+        res = run("SELECT 1 FROM metadata WHERE key = 'EPSG.VERSION' AND "
+                  "value = 'v10.003'");
+        if (!res.empty()) {
+            return;
+        }
+
+        throw FactoryException(
+            databasePath_ +
+            " lacks DATABASE.LAYOUT.VERSION.MAJOR / "
+            "DATABASE.LAYOUT.VERSION.MINOR "
+            "metadata. It comes from another PROJ installation.");
+    }
+    int nMajor = 0;
+    int nMinor = 0;
+    for (const auto &row : res) {
+        if (row[0] == "DATABASE.LAYOUT.VERSION.MAJOR") {
+            nMajor = atoi(row[1].c_str());
+        } else if (row[0] == "DATABASE.LAYOUT.VERSION.MINOR") {
+            nMinor = atoi(row[1].c_str());
+        }
+    }
+    if (nMajor != DATABASE_LAYOUT_VERSION_MAJOR) {
+        throw FactoryException(databasePath_ +
+                               " contains DATABASE.LAYOUT.VERSION.MAJOR = " +
+                               toString(nMajor) + " whereas " +
+                               toString(DATABASE_LAYOUT_VERSION_MAJOR) +
+                               " is expected. "
+                               "It comes from another PROJ installation.");
+    }
+    if (nMinor < DATABASE_LAYOUT_VERSION_MINOR) {
+        throw FactoryException(databasePath_ +
+                               " contains DATABASE.LAYOUT.VERSION.MINOR = " +
+                               toString(nMinor) + " whereas a number >= " +
+                               toString(DATABASE_LAYOUT_VERSION_MINOR) +
+                               " is expected. "
+                               "It comes from another PROJ installation.");
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 void DatabaseContext::Private::setHandle(sqlite3 *sqlite_handle) {
 
     assert(sqlite_handle);
@@ -865,6 +951,7 @@ DatabaseContext::create(const std::string &databasePath,
     if (!auxiliaryDatabasePaths.empty()) {
         dbCtx->getPrivate()->attachExtraDatabases(auxiliaryDatabasePaths);
     }
+    dbCtx->getPrivate()->checkDatabaseLayout();
     return dbCtx;
 }
 
@@ -1064,18 +1151,6 @@ std::string DatabaseContext::getOldProjGridName(const std::string &gridName) {
 
 // ---------------------------------------------------------------------------
 
-// FIXME: as we don't support datum ensemble yet, add it from name
-static std::string removeEnsembleSuffix(const std::string &name) {
-    if (name == "World Geodetic System 1984 ensemble") {
-        return "World Geodetic System 1984";
-    } else if (name == "European Terrestrial Reference System 1989 ensemble") {
-        return "European Terrestrial Reference System 1989";
-    }
-    return name;
-}
-
-// ---------------------------------------------------------------------------
-
 /** \brief Gets the alias name from an official name.
  *
  * @param officialName Official name. Mandatory
@@ -2012,18 +2087,38 @@ AuthorityFactory::createEllipsoid(const std::string &code) const {
 
 datum::GeodeticReferenceFrameNNPtr
 AuthorityFactory::createGeodeticDatum(const std::string &code) const {
+
+    datum::GeodeticReferenceFramePtr datum;
+    datum::DatumEnsemblePtr datumEnsemble;
+    constexpr bool turnEnsembleAsDatum = true;
+    createGeodeticDatumOrEnsemble(code, datum, datumEnsemble,
+                                  turnEnsembleAsDatum);
+    return NN_NO_CHECK(datum);
+}
+
+// ---------------------------------------------------------------------------
+
+void AuthorityFactory::createGeodeticDatumOrEnsemble(
+    const std::string &code, datum::GeodeticReferenceFramePtr &outDatum,
+    datum::DatumEnsemblePtr &outDatumEnsemble, bool turnEnsembleAsDatum) const {
     const auto cacheKey(d->authority() + code);
     {
-        auto datum = d->context()->d->getGeodeticDatumFromCache(cacheKey);
-        if (datum) {
-            return NN_NO_CHECK(datum);
+        outDatumEnsemble = d->context()->d->getDatumEnsembleFromCache(cacheKey);
+        if (outDatumEnsemble) {
+            if (!turnEnsembleAsDatum)
+                return;
+            outDatumEnsemble = nullptr;
+        }
+        outDatum = d->context()->d->getGeodeticDatumFromCache(cacheKey);
+        if (outDatum) {
+            return;
         }
     }
     auto res =
         d->runWithCodeParam("SELECT name, ellipsoid_auth_name, ellipsoid_code, "
                             "prime_meridian_auth_name, prime_meridian_code, "
                             "publication_date, frame_reference_epoch, "
-                            "deprecated FROM geodetic_datum "
+                            "ensemble_accuracy, deprecated FROM geodetic_datum "
                             "WHERE "
                             "auth_name = ? AND code = ?",
                             code);
@@ -2040,29 +2135,63 @@ AuthorityFactory::createGeodeticDatum(const std::string &code) const {
         const auto &prime_meridian_code = row[4];
         const auto &publication_date = row[5];
         const auto &frame_reference_epoch = row[6];
-        const bool deprecated = row[7] == "1";
-        auto ellipsoid = d->createFactory(ellipsoid_auth_name)
-                             ->createEllipsoid(ellipsoid_code);
-        auto pm = d->createFactory(prime_meridian_auth_name)
-                      ->createPrimeMeridian(prime_meridian_code);
-        auto props = d->createPropertiesSearchUsages(
-            "geodetic_datum", code, removeEnsembleSuffix(name), deprecated);
-        auto anchor = util::optional<std::string>();
-        if (!publication_date.empty()) {
-            props.set("PUBLICATION_DATE", publication_date);
-        }
-        auto datum =
-            frame_reference_epoch.empty()
-                ? datum::GeodeticReferenceFrame::create(props, ellipsoid,
-                                                        anchor, pm)
-                : util::nn_static_pointer_cast<datum::GeodeticReferenceFrame>(
-                      datum::DynamicGeodeticReferenceFrame::create(
-                          props, ellipsoid, anchor, pm,
-                          common::Measure(c_locale_stod(frame_reference_epoch),
-                                          common::UnitOfMeasure::YEAR),
-                          util::optional<std::string>()));
-        d->context()->d->cache(cacheKey, datum);
-        return datum;
+        const auto &ensemble_accuracy = row[7];
+        const bool deprecated = row[8] == "1";
+
+        std::string massagedName = name;
+        if (turnEnsembleAsDatum) {
+            if (name == "World Geodetic System 1984 ensemble") {
+                massagedName = "World Geodetic System 1984";
+            } else if (name ==
+                       "European Terrestrial Reference System 1989 ensemble") {
+                massagedName = "European Terrestrial Reference System 1989";
+            }
+        }
+        auto props = d->createPropertiesSearchUsages("geodetic_datum", code,
+                                                     massagedName, deprecated);
+
+        if (!turnEnsembleAsDatum && !ensemble_accuracy.empty()) {
+            auto resMembers =
+                d->run("SELECT member_auth_name, member_code FROM "
+                       "geodetic_datum_ensemble_member WHERE "
+                       "ensemble_auth_name = ? AND ensemble_code = ? "
+                       "ORDER BY sequence",
+                       {d->authority(), code});
+
+            std::vector<datum::DatumNNPtr> members;
+            for (const auto &memberRow : resMembers) {
+                members.push_back(
+                    d->createFactory(memberRow[0])->createDatum(memberRow[1]));
+            }
+            auto datumEnsemble = datum::DatumEnsemble::create(
+                props, std::move(members),
+                metadata::PositionalAccuracy::create(ensemble_accuracy));
+            d->context()->d->cache(cacheKey, datumEnsemble);
+            outDatumEnsemble = datumEnsemble.as_nullable();
+        } else {
+            auto ellipsoid = d->createFactory(ellipsoid_auth_name)
+                                 ->createEllipsoid(ellipsoid_code);
+            auto pm = d->createFactory(prime_meridian_auth_name)
+                          ->createPrimeMeridian(prime_meridian_code);
+
+            auto anchor = util::optional<std::string>();
+            if (!publication_date.empty()) {
+                props.set("PUBLICATION_DATE", publication_date);
+            }
+            auto datum = frame_reference_epoch.empty()
+                             ? datum::GeodeticReferenceFrame::create(
+                                   props, ellipsoid, anchor, pm)
+                             : util::nn_static_pointer_cast<
+                                   datum::GeodeticReferenceFrame>(
+                                   datum::DynamicGeodeticReferenceFrame::create(
+                                       props, ellipsoid, anchor, pm,
+                                       common::Measure(
+                                           c_locale_stod(frame_reference_epoch),
+                                           common::UnitOfMeasure::YEAR),
+                                       util::optional<std::string>()));
+            d->context()->d->cache(cacheKey, datum);
+            outDatum = datum.as_nullable();
+        }
     } catch (const std::exception &ex) {
         throw buildFactoryException("geodetic reference frame", code, ex);
     }
@@ -2080,9 +2209,23 @@ AuthorityFactory::createGeodeticDatum(const std::string &code) const {
 
 datum::VerticalReferenceFrameNNPtr
 AuthorityFactory::createVerticalDatum(const std::string &code) const {
+    datum::VerticalReferenceFramePtr datum;
+    datum::DatumEnsemblePtr datumEnsemble;
+    constexpr bool turnEnsembleAsDatum = true;
+    createVerticalDatumOrEnsemble(code, datum, datumEnsemble,
+                                  turnEnsembleAsDatum);
+    return NN_NO_CHECK(datum);
+}
+
+// ---------------------------------------------------------------------------
+
+void AuthorityFactory::createVerticalDatumOrEnsemble(
+    const std::string &code, datum::VerticalReferenceFramePtr &outDatum,
+    datum::DatumEnsemblePtr &outDatumEnsemble, bool turnEnsembleAsDatum) const {
     auto res =
         d->runWithCodeParam("SELECT name, publication_date, "
-                            "frame_reference_epoch, deprecated FROM "
+                            "frame_reference_epoch, ensemble_accuracy, "
+                            "deprecated FROM "
                             "vertical_datum WHERE auth_name = ? AND code = ?",
                             code);
     if (res.empty()) {
@@ -2094,24 +2237,49 @@ AuthorityFactory::createVerticalDatum(const std::string &code) const {
         const auto &name = row[0];
         const auto &publication_date = row[1];
         const auto &frame_reference_epoch = row[2];
-        const bool deprecated = row[3] == "1";
+        const auto &ensemble_accuracy = row[3];
+        const bool deprecated = row[4] == "1";
         auto props = d->createPropertiesSearchUsages("vertical_datum", code,
                                                      name, deprecated);
-        if (!publication_date.empty()) {
-            props.set("PUBLICATION_DATE", publication_date);
-        }
-        if (d->authority() == "ESRI" && starts_with(code, "from_geogdatum_")) {
-            props.set("VERT_DATUM_TYPE", "2002");
-        }
-        auto anchor = util::optional<std::string>();
-        if (frame_reference_epoch.empty()) {
-            return datum::VerticalReferenceFrame::create(props, anchor);
+        if (!turnEnsembleAsDatum && !ensemble_accuracy.empty()) {
+            auto resMembers =
+                d->run("SELECT member_auth_name, member_code FROM "
+                       "vertical_datum_ensemble_member WHERE "
+                       "ensemble_auth_name = ? AND ensemble_code = ? "
+                       "ORDER BY sequence",
+                       {d->authority(), code});
+
+            std::vector<datum::DatumNNPtr> members;
+            for (const auto &memberRow : resMembers) {
+                members.push_back(
+                    d->createFactory(memberRow[0])->createDatum(memberRow[1]));
+            }
+            auto datumEnsemble = datum::DatumEnsemble::create(
+                props, std::move(members),
+                metadata::PositionalAccuracy::create(ensemble_accuracy));
+            outDatumEnsemble = datumEnsemble.as_nullable();
         } else {
-            return datum::DynamicVerticalReferenceFrame::create(
-                props, anchor, util::optional<datum::RealizationMethod>(),
-                common::Measure(c_locale_stod(frame_reference_epoch),
-                                common::UnitOfMeasure::YEAR),
-                util::optional<std::string>());
+            if (!publication_date.empty()) {
+                props.set("PUBLICATION_DATE", publication_date);
+            }
+            if (d->authority() == "ESRI" &&
+                starts_with(code, "from_geogdatum_")) {
+                props.set("VERT_DATUM_TYPE", "2002");
+            }
+            auto anchor = util::optional<std::string>();
+            if (frame_reference_epoch.empty()) {
+                outDatum = datum::VerticalReferenceFrame::create(props, anchor)
+                               .as_nullable();
+            } else {
+                outDatum =
+                    datum::DynamicVerticalReferenceFrame::create(
+                        props, anchor,
+                        util::optional<datum::RealizationMethod>(),
+                        common::Measure(c_locale_stod(frame_reference_epoch),
+                                        common::UnitOfMeasure::YEAR),
+                        util::optional<std::string>())
+                        .as_nullable();
+            }
         }
     } catch (const std::exception &ex) {
         throw buildFactoryException("vertical reference frame", code, ex);
@@ -2472,20 +2640,24 @@ AuthorityFactory::createGeodeticCRS(const std::string &code,
 
         auto cs =
             d->createFactory(cs_auth_name)->createCoordinateSystem(cs_code);
-        auto datum =
-            d->createFactory(datum_auth_name)->createGeodeticDatum(datum_code);
+        datum::GeodeticReferenceFramePtr datum;
+        datum::DatumEnsemblePtr datumEnsemble;
+        constexpr bool turnEnsembleAsDatum = false;
+        d->createFactory(datum_auth_name)
+            ->createGeodeticDatumOrEnsemble(datum_code, datum, datumEnsemble,
+                                            turnEnsembleAsDatum);
 
         auto ellipsoidalCS =
             util::nn_dynamic_pointer_cast<cs::EllipsoidalCS>(cs);
         if ((type == GEOG_2D || type == GEOG_3D) && ellipsoidalCS) {
             auto crsRet = crs::GeographicCRS::create(
-                props, datum, NN_NO_CHECK(ellipsoidalCS));
+                props, datum, datumEnsemble, NN_NO_CHECK(ellipsoidalCS));
             d->context()->d->cache(cacheKey, crsRet);
             return crsRet;
         }
         auto geocentricCS = util::nn_dynamic_pointer_cast<cs::CartesianCS>(cs);
         if (type == GEOCENTRIC && geocentricCS) {
-            auto crsRet = crs::GeodeticCRS::create(props, datum,
+            auto crsRet = crs::GeodeticCRS::create(props, datum, datumEnsemble,
                                                    NN_NO_CHECK(geocentricCS));
             d->context()->d->cache(cacheKey, crsRet);
             return crsRet;
@@ -2539,16 +2711,19 @@ AuthorityFactory::createVerticalCRS(const std::string &code) const {
         const bool deprecated = row[5] == "1";
         auto cs =
             d->createFactory(cs_auth_name)->createCoordinateSystem(cs_code);
-        auto datum =
-            d->createFactory(datum_auth_name)->createVerticalDatum(datum_code);
-
+        datum::VerticalReferenceFramePtr datum;
+        datum::DatumEnsemblePtr datumEnsemble;
+        constexpr bool turnEnsembleAsDatum = false;
+        d->createFactory(datum_auth_name)
+            ->createVerticalDatumOrEnsemble(datum_code, datum, datumEnsemble,
+                                            turnEnsembleAsDatum);
         auto props = d->createPropertiesSearchUsages("vertical_crs", code, name,
                                                      deprecated);
 
         auto verticalCS = util::nn_dynamic_pointer_cast<cs::VerticalCS>(cs);
         if (verticalCS) {
-            auto crsRet =
-                crs::VerticalCRS::create(props, datum, NN_NO_CHECK(verticalCS));
+            auto crsRet = crs::VerticalCRS::create(props, datum, datumEnsemble,
+                                                   NN_NO_CHECK(verticalCS));
             d->context()->d->cache(cacheKey, crsRet);
             return crsRet;
         }
@@ -5340,6 +5515,11 @@ AuthorityFactory::getAuthorityCodes(const ObjectType &type,
     case ObjectType::CONCATENATED_OPERATION:
         sql = "SELECT code FROM concatenated_operation WHERE ";
         break;
+    case ObjectType::DATUM_ENSEMBLE:
+        sql = "SELECT code FROM object_view WHERE table_name IN "
+              "('geodetic_datum', 'vertical_datum') AND "
+              "type = 'ensemble' AND ";
+        break;
     }
 
     sql += "auth_name = ?";
@@ -5617,7 +5797,7 @@ std::string AuthorityFactory::getOfficialNameFromAlias(
                 if (res.empty()) { // shouldn't happen normally
                     return std::string();
                 }
-                return removeEnsembleSuffix(res.front()[0]);
+                return res.front()[0];
             }
         }
         return std::string();
@@ -5667,7 +5847,7 @@ std::string AuthorityFactory::getOfficialNameFromAlias(
         outTableName = row[1];
         outAuthName = row[2];
         outCode = row[3];
-        return removeEnsembleSuffix(row[0]);
+        return row[0];
     }
 }
 
@@ -5849,12 +6029,28 @@ AuthorityFactory::createObjectsFromNameEx(
                     res.emplace_back(
                         TableType("concatenated_operation", std::string()));
                     break;
+                case ObjectType::DATUM_ENSEMBLE:
+                    res.emplace_back(TableType("geodetic_datum", "ensemble"));
+                    res.emplace_back(TableType("vertical_datum", "ensemble"));
+                    break;
                 }
             }
         }
         return res;
     };
 
+    bool datumEnsembleAllowed = false;
+    if (allowedObjectTypes.empty()) {
+        datumEnsembleAllowed = true;
+    } else {
+        for (const auto type : allowedObjectTypes) {
+            if (type == ObjectType::DATUM_ENSEMBLE) {
+                datumEnsembleAllowed = true;
+                break;
+            }
+        }
+    }
+
     const auto listTableNameType = getTableAndTypeConstraints();
     bool first = true;
     ListOfParams params;
@@ -5872,6 +6068,8 @@ AuthorityFactory::createObjectsFromNameEx(
         if (!tableNameTypePair.second.empty()) {
             if (tableNameTypePair.second == "frame_reference_epoch") {
                 sql += "AND frame_reference_epoch IS NOT NULL ";
+            } else if (tableNameTypePair.second == "ensemble") {
+                sql += "AND ensemble_accuracy IS NOT NULL ";
             } else {
                 sql += "AND type = '";
                 sql += tableNameTypePair.second;
@@ -5906,6 +6104,8 @@ AuthorityFactory::createObjectsFromNameEx(
         if (!tableNameTypePair.second.empty()) {
             if (tableNameTypePair.second == "frame_reference_epoch") {
                 sql += "AND ov.frame_reference_epoch IS NOT NULL ";
+            } else if (tableNameTypePair.second == "ensemble") {
+                sql += "AND ov.ensemble_accuracy IS NOT NULL ";
             } else {
                 sql += "AND ov.type = '";
                 sql += tableNameTypePair.second;
@@ -6019,6 +6219,8 @@ AuthorityFactory::createObjectsFromNameEx(
         auto sqlRes = d->run(sql, params);
         bool isFirst = true;
         bool firstIsDeprecated = false;
+        bool foundExactMatch = false;
+        std::size_t hashCodeFirstMatch = 0;
         for (const auto &row : sqlRes) {
             const auto &name = row[3];
             if (approximateMatch) {
@@ -6053,7 +6255,7 @@ AuthorityFactory::createObjectsFromNameEx(
                 break;
             }
             auto factory = d->createFactory(auth_name);
-            auto getObject = [&factory](
+            auto getObject = [&factory, datumEnsembleAllowed](
                 const std::string &l_table_name,
                 const std::string &l_code) -> common::IdentifiedObjectNNPtr {
                 if (l_table_name == "prime_meridian") {
@@ -6061,8 +6263,32 @@ AuthorityFactory::createObjectsFromNameEx(
                 } else if (l_table_name == "ellipsoid") {
                     return factory->createEllipsoid(l_code);
                 } else if (l_table_name == "geodetic_datum") {
+                    if (datumEnsembleAllowed) {
+                        datum::GeodeticReferenceFramePtr datum;
+                        datum::DatumEnsemblePtr datumEnsemble;
+                        constexpr bool turnEnsembleAsDatum = false;
+                        factory->createGeodeticDatumOrEnsemble(
+                            l_code, datum, datumEnsemble, turnEnsembleAsDatum);
+                        if (datum) {
+                            return NN_NO_CHECK(datum);
+                        }
+                        assert(datumEnsemble);
+                        return NN_NO_CHECK(datumEnsemble);
+                    }
                     return factory->createGeodeticDatum(l_code);
                 } else if (l_table_name == "vertical_datum") {
+                    if (datumEnsembleAllowed) {
+                        datum::VerticalReferenceFramePtr datum;
+                        datum::DatumEnsemblePtr datumEnsemble;
+                        constexpr bool turnEnsembleAsDatum = false;
+                        factory->createVerticalDatumOrEnsemble(
+                            l_code, datum, datumEnsemble, turnEnsembleAsDatum);
+                        if (datum) {
+                            return NN_NO_CHECK(datum);
+                        }
+                        assert(datumEnsemble);
+                        return NN_NO_CHECK(datumEnsemble);
+                    }
                     return factory->createVerticalDatum(l_code);
                 } else if (l_table_name == "geodetic_crs") {
                     return factory->createGeodeticCRS(l_code);
@@ -6082,11 +6308,38 @@ AuthorityFactory::createObjectsFromNameEx(
                 }
                 throw std::runtime_error("Unsupported table_name");
             };
-            res.emplace_back(PairObjectName(getObject(table_name, code), name));
+            const auto obj = getObject(table_name, code);
+            if (metadata::Identifier::canonicalizeName(obj->nameStr()) ==
+                canonicalizedSearchedName) {
+                foundExactMatch = true;
+            }
+
+            const auto objPtr = obj.get();
+            if (res.empty()) {
+                hashCodeFirstMatch = typeid(*objPtr).hash_code();
+            } else if (hashCodeFirstMatch != typeid(*objPtr).hash_code()) {
+                hashCodeFirstMatch = 0;
+            }
+
+            res.emplace_back(PairObjectName(obj, name));
             if (limitResultCount > 0 && res.size() == limitResultCount) {
                 break;
             }
         }
+
+        // If we found a name that is an exact match, and all objects have the
+        // same type, and we are not in approximate mode, only keep the objet(s)
+        // with the exact name match.
+        if (foundExactMatch && hashCodeFirstMatch != 0 && !approximateMatch) {
+            std::list<PairObjectName> resTmp;
+            for (const auto &pair : res) {
+                if (metadata::Identifier::canonicalizeName(
+                        pair.first->nameStr()) == canonicalizedSearchedName) {
+                    resTmp.emplace_back(pair);
+                }
+            }
+            res = std::move(resTmp);
+        }
     }
 
     auto sortLambda = [](const PairObjectName &a, const PairObjectName &b) {
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index f4ec7740..8a42e7ee 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -70,7 +70,6 @@
 // clang-format off
 #include "proj.h"
 #include "proj_internal.h"
-#include "proj_api.h"
 // clang-format on
 
 using namespace NS_PROJ::common;
@@ -140,6 +139,7 @@ struct WKTFormatter::Private {
         bool primeMeridianInDegree_ = false;
         bool use2019Keywords_ = false;
         bool useESRIDialect_ = false;
+        bool allowEllipsoidalHeightAsVerticalCRS_ = false;
         OutputAxisRule outputAxis_ = WKTFormatter::OutputAxisRule::YES;
     };
     Params params_{};
@@ -251,6 +251,8 @@ WKTFormatter::setOutputAxis(OutputAxisRule outputAxisIn) noexcept {
  *
  * The default is strict mode, in which case a FormattingException can be
  * thrown.
+ * In non-strict mode, a Geographic 3D CRS can be for example exported as
+ * WKT1_GDAL with 3 axes, whereas this is normally not allowed.
  */
 WKTFormatter &WKTFormatter::setStrict(bool strictIn) noexcept {
     d->params_.strict_ = strictIn;
@@ -264,6 +266,28 @@ bool WKTFormatter::isStrict() const noexcept { return d->params_.strict_; }
 
 // ---------------------------------------------------------------------------
 
+/** \brief Set whether the formatter should export, in WKT1, a Geographic or
+ * Projected 3D CRS as a compound CRS whose vertical part represents an
+ * ellipsoidal height.
+ */
+WKTFormatter &
+WKTFormatter::setAllowEllipsoidalHeightAsVerticalCRS(bool allow) noexcept {
+    d->params_.allowEllipsoidalHeightAsVerticalCRS_ = allow;
+    return *this;
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Return whether the formatter should export, in WKT1, a Geographic or
+ * Projected 3D CRS as a compound CRS whose vertical part represents an
+ * ellipsoidal height.
+ */
+bool WKTFormatter::isAllowedEllipsoidalHeightAsVerticalCRS() const noexcept {
+    return d->params_.allowEllipsoidalHeightAsVerticalCRS_;
+}
+
+// ---------------------------------------------------------------------------
+
 /** Returns the WKT string from the formatter. */
 const std::string &WKTFormatter::toString() const {
     if (d->indentLevel_ > 0 || d->level_ > 0) {
@@ -1979,14 +2003,80 @@ PrimeMeridianNNPtr WKTParser::Private::buildPrimeMeridian(
     try {
         double angleValue = asDouble(children[1]);
 
-        // Correct for GDAL WKT1 departure
+        // Correct for GDAL WKT1 and WKT1-ESRI departure
         if (name == "Paris" && std::fabs(angleValue - 2.33722917) < 1e-8 &&
-            unit == UnitOfMeasure::GRAD) {
+            unit._isEquivalentTo(UnitOfMeasure::GRAD,
+                                 util::IComparable::Criterion::EQUIVALENT)) {
             angleValue = 2.5969213;
+        } else {
+            static const struct {
+                const char *name;
+                int deg;
+                int min;
+                double sec;
+            } primeMeridiansDMS[] = {
+                {"Lisbon", -9, 7, 54.862},  {"Bogota", -74, 4, 51.3},
+                {"Madrid", -3, 41, 14.55},  {"Rome", 12, 27, 8.4},
+                {"Bern", 7, 26, 22.5},      {"Jakarta", 106, 48, 27.79},
+                {"Ferro", -17, 40, 0},      {"Brussels", 4, 22, 4.71},
+                {"Stockholm", 18, 3, 29.8}, {"Athens", 23, 42, 58.815},
+                {"Oslo", 10, 43, 22.5},     {"Paris RGS", 2, 20, 13.95},
+                {"Paris_RGS", 2, 20, 13.95}};
+
+            // Current epsg.org output may use the EPSG:9110 "sexagesimal DMS"
+            // unit and a DD.MMSSsss value, but this will likely be changed to
+            // use decimal degree.
+            // Or WKT1 may for example use the Paris RGS decimal degree value
+            // but with a GEOGCS with UNIT["Grad"]
+            for (const auto &pmDef : primeMeridiansDMS) {
+                if (name == pmDef.name) {
+                    double dmsAsDecimalValue =
+                        (pmDef.deg >= 0 ? 1 : -1) *
+                        (std::abs(pmDef.deg) + pmDef.min / 100. +
+                         pmDef.sec / 10000.);
+                    double dmsAsDecimalDegreeValue =
+                        (pmDef.deg >= 0 ? 1 : -1) *
+                        (std::abs(pmDef.deg) + pmDef.min / 60. +
+                         pmDef.sec / 3600.);
+                    if (std::fabs(angleValue - dmsAsDecimalValue) < 1e-8 ||
+                        std::fabs(angleValue - dmsAsDecimalDegreeValue) <
+                            1e-8) {
+                        angleValue = dmsAsDecimalDegreeValue;
+                        unit = UnitOfMeasure::DEGREE;
+                    }
+                    break;
+                }
+            }
+        }
+
+        auto &properties = buildProperties(node);
+        if (dbContext_ && esriStyle_) {
+            std::string outTableName;
+            std::string codeFromAlias;
+            std::string authNameFromAlias;
+            auto authFactory = AuthorityFactory::create(NN_NO_CHECK(dbContext_),
+                                                        std::string());
+            auto officialName = authFactory->getOfficialNameFromAlias(
+                name, "prime_meridian", "ESRI", false, outTableName,
+                authNameFromAlias, codeFromAlias);
+            if (!officialName.empty()) {
+                properties.set(IdentifiedObject::NAME_KEY, officialName);
+                if (!authNameFromAlias.empty()) {
+                    auto identifiers = ArrayOfBaseObject::create();
+                    identifiers->add(Identifier::create(
+                        codeFromAlias,
+                        PropertyMap()
+                            .set(Identifier::CODESPACE_KEY, authNameFromAlias)
+                            .set(Identifier::AUTHORITY_KEY,
+                                 authNameFromAlias)));
+                    properties.set(IdentifiedObject::IDENTIFIERS_KEY,
+                                   identifiers);
+                }
+            }
         }
 
         Angle angle(angleValue, unit);
-        return PrimeMeridian::create(buildProperties(node), angle);
+        return PrimeMeridian::create(properties, angle);
     } catch (const std::exception &e) {
         throw buildRethrow(__FUNCTION__, e);
     }
@@ -2098,9 +2188,15 @@ GeodeticReferenceFrameNNPtr WKTParser::Private::buildGeodeticReferenceFrame(
         return false;
     };
 
-    if (name == "WGS_1984") {
+    // Remap GDAL WGS_1984 to EPSG v9 "World Geodetic System 1984" official
+    // name.
+    // Also remap EPSG v10 datum ensemble names to non-ensemble EPSG v9
+    if (name == "WGS_1984" || name == "World Geodetic System 1984 ensemble") {
         properties.set(IdentifiedObject::NAME_KEY,
                        GeodeticReferenceFrame::EPSG_6326->nameStr());
+    } else if (name == "European Terrestrial Reference System 1989 ensemble") {
+        properties.set(IdentifiedObject::NAME_KEY,
+                       "European Terrestrial Reference System 1989");
     } else if (starts_with(name, "D_")) {
         esriStyle_ = true;
         const char *tableNameForAlias = nullptr;
@@ -2110,6 +2206,10 @@ GeodeticReferenceFrameNNPtr WKTParser::Private::buildGeodeticReferenceFrame(
             name = "World Geodetic System 1984";
             authNameFromAlias = Identifier::EPSG;
             codeFromAlias = "6326";
+        } else if (name == "D_ETRS_1989") {
+            name = "European Terrestrial Reference System 1989";
+            authNameFromAlias = Identifier::EPSG;
+            codeFromAlias = "6258";
         } else {
             tableNameForAlias = "geodetic_datum";
         }
@@ -2734,6 +2834,9 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) {
         throw ParsingException("Missing DATUM or ENSEMBLE node");
     }
 
+    // Do that now so that esriStyle_ can be set before buildPrimeMeridian()
+    auto props = buildProperties(node);
+
     auto &dynamicNode = nodeP->lookForChild(WKTConstants::DYNAMIC);
 
     auto &csNode = nodeP->lookForChild(WKTConstants::CS_);
@@ -2771,7 +2874,6 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) {
         angularUnit = primeMeridian->longitude().unit();
     }
 
-    auto props = buildProperties(node);
     addExtensionProj4ToProp(nodeP, props);
 
     // No explicit AXIS node ? (WKT1)
@@ -2790,6 +2892,32 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) {
                   .as_nullable()
             : nullptr;
     auto cs = buildCS(csNode, node, angularUnit);
+
+    // If there's no CS[] node, typically for a BASEGEODCRS of a projected CRS,
+    // in a few rare cases, this might be a Geocentric CRS, and thus a
+    // Cartesian CS, and not the ellipsoidalCS we assumed above. The only way
+    // to figure that is to resolve the CRS from its code...
+    if (isNull(csNode) && dbContext_ &&
+        ci_equal(nodeName, WKTConstants::BASEGEODCRS)) {
+        const auto &nodeChildren = nodeP->children();
+        for (const auto &subNode : nodeChildren) {
+            const auto &subNodeName(subNode->GP()->value());
+            if (ci_equal(subNodeName, WKTConstants::ID) ||
+                ci_equal(subNodeName, WKTConstants::AUTHORITY)) {
+                auto id = buildId(subNode, true, false);
+                if (id) {
+                    try {
+                        auto authFactory = AuthorityFactory::create(
+                            NN_NO_CHECK(dbContext_), *id->codeSpace());
+                        auto dbCRS = authFactory->createGeodeticCRS(id->code());
+                        cs = dbCRS->coordinateSystem();
+                    } catch (const util::Exception &) {
+                    }
+                }
+            }
+        }
+    }
+
     auto ellipsoidalCS = nn_dynamic_pointer_cast<EllipsoidalCS>(cs);
     if (ellipsoidalCS) {
         if (ci_equal(nodeName, WKTConstants::GEOCCS)) {
@@ -2940,7 +3068,8 @@ UnitOfMeasure WKTParser::Private::guessUnitForParameter(
     const UnitOfMeasure &defaultAngularUnit) {
     UnitOfMeasure unit;
     // scale must be first because of 'Scale factor on pseudo standard parallel'
-    if (ci_find(paramName, "scale") != std::string::npos) {
+    if (ci_find(paramName, "scale") != std::string::npos ||
+        ci_find(paramName, "scaling factor") != std::string::npos) {
         unit = UnitOfMeasure::SCALE_UNITY;
     } else if (ci_find(paramName, "latitude") != std::string::npos ||
                ci_find(paramName, "longitude") != std::string::npos ||
@@ -3872,8 +4001,16 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
         return createPseudoMercator(props, NN_NO_CHECK(cartesianCS));
     }
 
-    auto linearUnit = buildUnitInSubNode(node, UnitOfMeasure::Type::LINEAR);
-    auto angularUnit = baseGeodCRS->coordinateSystem()->axisList()[0]->unit();
+    // For WKT2, if there is no explicit parameter unit, use metre for linear
+    // units and degree for angular units
+    auto linearUnit =
+        !isNull(conversionNode)
+            ? UnitOfMeasure::METRE
+            : buildUnitInSubNode(node, UnitOfMeasure::Type::LINEAR);
+    auto angularUnit =
+        !isNull(conversionNode)
+            ? UnitOfMeasure::DEGREE
+            : baseGeodCRS->coordinateSystem()->axisList()[0]->unit();
 
     auto conversion =
         !isNull(conversionNode)
@@ -4941,6 +5078,10 @@ class JSONParser {
     TransformationNNPtr buildTransformation(const json &j);
     ConcatenatedOperationNNPtr buildConcatenatedOperation(const json &j);
 
+    void buildGeodeticDatumOrDatumEnsemble(const json &j,
+                                           GeodeticReferenceFramePtr &datum,
+                                           DatumEnsemblePtr &datumEnsemble);
+
     static util::optional<std::string> getAnchor(const json &j) {
         util::optional<std::string> anchor;
         if (j.contains("anchor")) {
@@ -5400,9 +5541,9 @@ BaseObjectNNPtr JSONParser::create(const json &j)
 
 // ---------------------------------------------------------------------------
 
-GeographicCRSNNPtr JSONParser::buildGeographicCRS(const json &j) {
-    GeodeticReferenceFramePtr datum;
-    DatumEnsemblePtr datumEnsemble;
+void JSONParser::buildGeodeticDatumOrDatumEnsemble(
+    const json &j, GeodeticReferenceFramePtr &datum,
+    DatumEnsemblePtr &datumEnsemble) {
     if (j.contains("datum")) {
         auto datumJ = getObject(j, "datum");
         datum = util::nn_dynamic_pointer_cast<GeodeticReferenceFrame>(
@@ -5415,6 +5556,14 @@ GeographicCRSNNPtr JSONParser::buildGeographicCRS(const json &j) {
         datumEnsemble =
             buildDatumEnsemble(getObject(j, "datum_ensemble")).as_nullable();
     }
+}
+
+// ---------------------------------------------------------------------------
+
+GeographicCRSNNPtr JSONParser::buildGeographicCRS(const json &j) {
+    GeodeticReferenceFramePtr datum;
+    DatumEnsemblePtr datumEnsemble;
+    buildGeodeticDatumOrDatumEnsemble(j, datum, datumEnsemble);
     auto csJ = getObject(j, "coordinate_system");
     auto ellipsoidalCS =
         util::nn_dynamic_pointer_cast<EllipsoidalCS>(buildCS(csJ));
@@ -5428,12 +5577,9 @@ GeographicCRSNNPtr JSONParser::buildGeographicCRS(const json &j) {
 // ---------------------------------------------------------------------------
 
 GeodeticCRSNNPtr JSONParser::buildGeodeticCRS(const json &j) {
-    auto datumJ = getObject(j, "datum");
-    if (getType(datumJ) != "GeodeticReferenceFrame") {
-        throw ParsingException("Unsupported type for datum.");
-    }
-    auto datum = buildGeodeticReferenceFrame(datumJ);
+    GeodeticReferenceFramePtr datum;
     DatumEnsemblePtr datumEnsemble;
+    buildGeodeticDatumOrDatumEnsemble(j, datum, datumEnsemble);
     auto csJ = getObject(j, "coordinate_system");
     auto cs = buildCS(csJ);
     auto props = buildProperties(j);
@@ -5468,7 +5614,13 @@ GeodeticCRSNNPtr JSONParser::buildGeodeticCRS(const json &j) {
 // ---------------------------------------------------------------------------
 
 ProjectedCRSNNPtr JSONParser::buildProjectedCRS(const json &j) {
-    auto baseCRS = buildGeographicCRS(getObject(j, "base_crs"));
+    auto jBaseCRS = getObject(j, "base_crs");
+    auto jBaseCS = getObject(jBaseCRS, "coordinate_system");
+    auto baseCS = buildCS(jBaseCS);
+    auto baseCRS = dynamic_cast<EllipsoidalCS *>(baseCS.get()) != nullptr
+                       ? util::nn_static_pointer_cast<GeodeticCRS>(
+                             buildGeographicCRS(jBaseCRS))
+                       : buildGeodeticCRS(jBaseCRS);
     auto csJ = getObject(j, "coordinate_system");
     auto cartesianCS = util::nn_dynamic_pointer_cast<CartesianCS>(buildCS(csJ));
     if (!cartesianCS) {
@@ -6263,6 +6415,9 @@ static BaseObjectNNPtr createFromUserInput(const std::string &text,
         if (type == "datum") {
             return factory->createDatum(code);
         }
+        if (type == "ensemble") {
+            return factory->createDatumEnsemble(code);
+        }
         if (type == "ellipsoid") {
             return factory->createEllipsoid(code);
         }
@@ -6381,6 +6536,8 @@ static BaseObjectNNPtr createFromUserInput(const std::string &text,
                                               AuthorityFactory::ObjectType::
                                                   DATUM,
                                               AuthorityFactory::ObjectType::
+                                                  DATUM_ENSEMBLE,
+                                              AuthorityFactory::ObjectType::
                                                   COORDINATE_OPERATION},
                     goOn);
             } catch (const std::exception &) {
@@ -9460,8 +9617,8 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
         std::string methodName = "PROJ " + step.name;
         for (const auto &param : step.paramValues) {
             if (is_in_stringlist(param.key,
-                                 "wktext,no_defs,datum,ellps,a,b,R,towgs84,"
-                                 "nadgrids,geoidgrids,"
+                                 "wktext,no_defs,datum,ellps,a,b,R,f,rf,"
+                                 "towgs84,nadgrids,geoidgrids,"
                                  "units,to_meter,vunits,vto_meter,type")) {
                 continue;
             }
@@ -9761,7 +9918,7 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
         paralist *list = pj_expand_init(ctx, init);
         ctx->projStringParserCreateFromPROJStringRecursionCounter--;
         if (!list) {
-            pj_dealloc(init);
+            free(init);
             throw ParsingException("cannot expand " + projString);
         }
         std::string expanded;
@@ -9784,7 +9941,7 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
             }
 
             auto n = t->next;
-            pj_dealloc(t);
+            free(t);
             t = n;
         }
         for (const auto &pair : d->steps_[0].paramValues) {