Merge branch 'master' into 6.0

14 files changed, 474 insertions, 289 deletions

diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 1b3374f3..4f13f238 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -1102,8 +1102,8 @@ PJ  *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *source_crs, const char
             double north_lat = 0.0;
 
             const char* name = proj_get_name(op);
-            if( name && (strstr(name, "Null geographic offset") ||
-                         strstr(name, "Null geocentric translation")) )
+            if( name && (strstr(name, "Ballpark geographic offset") ||
+                         strstr(name, "Ballpark geocentric translation")) )
             {
                 // Skip default transformations
             }
diff --git a/src/apply_vgridshift.cpp b/src/apply_vgridshift.cpp
index 61e0c528..ae23e39a 100644
--- a/src/apply_vgridshift.cpp
+++ b/src/apply_vgridshift.cpp
@@ -35,16 +35,16 @@
 #include "proj_internal.h"
 #include "proj_internal.h"
 
-static int is_nodata(float value)
+static int is_nodata(float value, double vmultiplier)
 {
     /* nodata?  */
     /* GTX official nodata value if  -88.88880f, but some grids also */
     /* use other  big values for nodata (e.g naptrans2008.gtx has */
     /* nodata values like -2147479936), so test them too */
-    return value > 1000 || value < -1000 || value == -88.88880f;
+    return value * vmultiplier > 1000 || value * vmultiplier < -1000 || value == -88.88880f;
 }
 
-static double read_vgrid_value( PJ *defn, PJ_LP input, int *gridlist_count_p, PJ_GRIDINFO **tables, struct CTABLE *ct) {
+static double read_vgrid_value( PJ *defn, PJ_LP input, double vmultiplier, int *gridlist_count_p, PJ_GRIDINFO **tables, struct CTABLE *ct) {
     int  itable = 0;
     double value = HUGE_VAL;
     double grid_x, grid_y;
@@ -129,28 +129,28 @@ static double read_vgrid_value( PJ *defn, PJ_LP input, int *gridlist_count_p, PJ
             double total_weight = 0.0;
             int n_weights = 0;
             value = 0.0f;
-            if( !is_nodata(value_a) )
+            if( !is_nodata(value_a, vmultiplier) )
             {
                 double weight = (1.0-grid_x) * (1.0-grid_y);
                 value += value_a * weight;
                 total_weight += weight;
                 n_weights ++;
             }
-            if( !is_nodata(value_b) )
+            if( !is_nodata(value_b, vmultiplier) )
             {
                 double weight = (grid_x) * (1.0-grid_y);
                 value += value_b * weight;
                 total_weight += weight;
                 n_weights ++;
             }
-            if( !is_nodata(value_c) )
+            if( !is_nodata(value_c, vmultiplier) )
             {
                 double weight = (1.0-grid_x) * (grid_y);
                 value += value_c * weight;
                 total_weight += weight;
                 n_weights ++;
             }
-            if( !is_nodata(value_d) )
+            if( !is_nodata(value_d, vmultiplier) )
             {
                 double weight = (grid_x) * (grid_y);
                 value += value_d * weight;
@@ -165,7 +165,7 @@ static double read_vgrid_value( PJ *defn, PJ_LP input, int *gridlist_count_p, PJ
 
     }
 
-    return value;
+    return value * vmultiplier;
 }
 
 /************************************************************************/
@@ -218,7 +218,7 @@ int pj_apply_vgridshift( PJ *defn, const char *listname,
         input.phi = y[io];
         input.lam = x[io];
 
-        value = read_vgrid_value(defn, input, gridlist_count_p, tables, &ct);
+        value = read_vgrid_value(defn, input, 1.0, gridlist_count_p, tables, &ct);
 
         if( inverse )
             z[io] -= value;
@@ -310,7 +310,7 @@ int proj_vgrid_init(PJ* P, const char *grids) {
 }
 
 /***********************************************/
-double proj_vgrid_value(PJ *P, PJ_LP lp){
+double proj_vgrid_value(PJ *P, PJ_LP lp, double vmultiplier){
 /***********************************************
 
   Read grid value at position lp in grids loaded
@@ -324,7 +324,7 @@ double proj_vgrid_value(PJ *P, PJ_LP lp){
     double value;
     memset(&used_grid, 0, sizeof(struct CTABLE));
 
-    value = read_vgrid_value(P, lp, &(P->vgridlist_geoid_count), P->vgridlist_geoid, &used_grid);
+    value = read_vgrid_value(P, lp, vmultiplier, &(P->vgridlist_geoid_count), P->vgridlist_geoid, &used_grid);
     proj_log_trace(P, "proj_vgrid_value: (%f, %f) = %f", lp.lam*RAD_TO_DEG, lp.phi*RAD_TO_DEG, value);
 
     return value;
diff --git a/src/apps/cs2cs.cpp b/src/apps/cs2cs.cpp
index 150548c5..dafd06f8 100644
--- a/src/apps/cs2cs.cpp
+++ b/src/apps/cs2cs.cpp
@@ -206,10 +206,10 @@ static void process(FILE *fid)
 }
 
 /************************************************************************/
-/*                          instanciate_crs()                           */
+/*                          instantiate_crs()                           */
 /************************************************************************/
 
-static PJ *instanciate_crs(const std::string &definition,
+static PJ *instantiate_crs(const std::string &definition,
                                bool &isGeog, double &toRadians,
                                bool &isLatFirst) {
     PJ *crs = proj_create(nullptr,
@@ -541,7 +541,7 @@ int main(int argc, char **argv) {
     PJ *src = nullptr;
     if (!fromStr.empty()) {
         bool ignored;
-        src = instanciate_crs(fromStr, srcIsGeog,
+        src = instantiate_crs(fromStr, srcIsGeog,
                               srcToRadians, ignored);
         if (!src) {
             emess(3, "cannot instantiate source coordinate system");
@@ -550,7 +550,7 @@ int main(int argc, char **argv) {
 
     PJ *dst = nullptr;
     if (!toStr.empty()) {
-        dst = instanciate_crs(toStr, destIsGeog,
+        dst = instantiate_crs(toStr, destIsGeog,
                               destToRadians, destIsLatLong);
         if (!dst) {
             emess(3, "cannot instantiate target coordinate system");
diff --git a/src/apps/projinfo.cpp b/src/apps/projinfo.cpp
index 9f908c8a..094587e2 100644
--- a/src/apps/projinfo.cpp
+++ b/src/apps/projinfo.cpp
@@ -136,11 +136,11 @@ static std::string c_ify_string(const std::string &str) {
 
 // ---------------------------------------------------------------------------
 
-static BaseObjectNNPtr buildObject(DatabaseContextPtr dbContext,
-                                   const std::string &user_string,
-                                   bool kindIsCRS, const std::string &context,
-                                   bool buildBoundCRSToWGS84, CoordinateOperationContext::IntermediateCRSUse allowUseIntermediateCRS,
-                                   bool quiet) {
+static BaseObjectNNPtr buildObject(
+    DatabaseContextPtr dbContext, const std::string &user_string,
+    bool kindIsCRS, const std::string &context, bool buildBoundCRSToWGS84,
+    CoordinateOperationContext::IntermediateCRSUse allowUseIntermediateCRS,
+    bool quiet) {
     BaseObjectPtr obj;
 
     std::string l_user_string(user_string);
@@ -213,7 +213,8 @@ static BaseObjectNNPtr buildObject(DatabaseContextPtr dbContext,
     if (buildBoundCRSToWGS84) {
         auto crs = std::dynamic_pointer_cast<CRS>(obj);
         if (crs) {
-            obj = crs->createBoundCRSToWGS84IfPossible(dbContext, allowUseIntermediateCRS)
+            obj = crs->createBoundCRSToWGS84IfPossible(dbContext,
+                                                       allowUseIntermediateCRS)
                       .as_nullable();
         }
     }
@@ -223,8 +224,10 @@ static BaseObjectNNPtr buildObject(DatabaseContextPtr dbContext,
 
 // ---------------------------------------------------------------------------
 
-static void outputObject(DatabaseContextPtr dbContext, BaseObjectNNPtr obj,
-                         CoordinateOperationContext::IntermediateCRSUse allowUseIntermediateCRS, const OutputOptions &outputOpt) {
+static void outputObject(
+    DatabaseContextPtr dbContext, BaseObjectNNPtr obj,
+    CoordinateOperationContext::IntermediateCRSUse allowUseIntermediateCRS,
+    const OutputOptions &outputOpt) {
 
     auto identified = dynamic_cast<const IdentifiedObject *>(obj.get());
     if (!outputOpt.quiet && identified && identified->isDeprecated()) {
@@ -260,7 +263,7 @@ static void outputObject(DatabaseContextPtr dbContext, BaseObjectNNPtr obj,
                 }
                 auto crs = nn_dynamic_pointer_cast<CRS>(obj);
                 if (!outputOpt.quiet) {
-                    if( crs ) {
+                    if (crs) {
                         std::cout << "PROJ.4 string:" << std::endl;
                     } else {
                         std::cout << "PROJ string:" << std::endl;
@@ -271,8 +274,8 @@ static void outputObject(DatabaseContextPtr dbContext, BaseObjectNNPtr obj,
                 if (crs) {
                     objToExport =
                         nn_dynamic_pointer_cast<IPROJStringExportable>(
-                            crs->createBoundCRSToWGS84IfPossible(dbContext,
-                                                                 allowUseIntermediateCRS));
+                            crs->createBoundCRSToWGS84IfPossible(
+                                dbContext, allowUseIntermediateCRS));
                 }
                 if (!objToExport) {
                     objToExport = projStringExportable;
@@ -411,8 +414,8 @@ static void outputObject(DatabaseContextPtr dbContext, BaseObjectNNPtr obj,
                 std::shared_ptr<IWKTExportable> objToExport;
                 if (crs) {
                     objToExport = nn_dynamic_pointer_cast<IWKTExportable>(
-                        crs->createBoundCRSToWGS84IfPossible(dbContext,
-                                                             allowUseIntermediateCRS));
+                        crs->createBoundCRSToWGS84IfPossible(
+                            dbContext, allowUseIntermediateCRS));
                 }
                 if (!objToExport) {
                     objToExport = wktExportable;
@@ -505,6 +508,10 @@ static void outputOperationSummary(const CoordinateOperationNNPtr &op) {
         std::cout << "unknown domain of validity";
     }
 
+    if (op->hasBallparkTransformation()) {
+        std::cout << ", has ballpark transformation";
+    }
+
     std::cout << std::endl;
 }
 
@@ -514,26 +521,25 @@ static void outputOperations(
     DatabaseContextPtr dbContext, const std::string &sourceCRSStr,
     const std::string &targetCRSStr, const ExtentPtr &bboxFilter,
     CoordinateOperationContext::SpatialCriterion spatialCriterion,
+    bool spatialCriterionExplicitlySpecified,
     CoordinateOperationContext::SourceTargetCRSExtentUse crsExtentUse,
     CoordinateOperationContext::GridAvailabilityUse gridAvailabilityUse,
     CoordinateOperationContext::IntermediateCRSUse allowUseIntermediateCRS,
     const std::vector<std::pair<std::string, std::string>> &pivots,
     const std::string &authority, bool usePROJGridAlternatives,
     bool showSuperseded, const OutputOptions &outputOpt, bool summary) {
-    auto sourceObj = buildObject(dbContext, sourceCRSStr, true, "source CRS",
-                                 false,
-                                 CoordinateOperationContext::IntermediateCRSUse::NEVER,
-                                 outputOpt.quiet);
+    auto sourceObj = buildObject(
+        dbContext, sourceCRSStr, true, "source CRS", false,
+        CoordinateOperationContext::IntermediateCRSUse::NEVER, outputOpt.quiet);
     auto sourceCRS = nn_dynamic_pointer_cast<CRS>(sourceObj);
     if (!sourceCRS) {
         std::cerr << "source CRS string is not a CRS" << std::endl;
         std::exit(1);
     }
 
-    auto targetObj = buildObject(dbContext, targetCRSStr, true, "target CRS",
-                                 false,
-                                 CoordinateOperationContext::IntermediateCRSUse::NEVER,
-                                 outputOpt.quiet);
+    auto targetObj = buildObject(
+        dbContext, targetCRSStr, true, "target CRS", false,
+        CoordinateOperationContext::IntermediateCRSUse::NEVER, outputOpt.quiet);
     auto targetCRS = nn_dynamic_pointer_cast<CRS>(targetObj);
     if (!targetCRS) {
         std::cerr << "target CRS string is not a CRS" << std::endl;
@@ -541,6 +547,7 @@ static void outputOperations(
     }
 
     std::vector<CoordinateOperationNNPtr> list;
+    size_t spatialCriterionPartialIntersectionResultCount = 0;
     try {
         auto authFactory =
             dbContext
@@ -558,6 +565,21 @@ static void outputOperations(
         ctxt->setDiscardSuperseded(!showSuperseded);
         list = CoordinateOperationFactory::create()->createOperations(
             NN_NO_CHECK(sourceCRS), NN_NO_CHECK(targetCRS), ctxt);
+        if (!spatialCriterionExplicitlySpecified &&
+            spatialCriterion == CoordinateOperationContext::SpatialCriterion::
+                                    STRICT_CONTAINMENT) {
+            try {
+                ctxt->setSpatialCriterion(
+                    CoordinateOperationContext::SpatialCriterion::
+                        PARTIAL_INTERSECTION);
+                spatialCriterionPartialIntersectionResultCount =
+                    CoordinateOperationFactory::create()
+                        ->createOperations(NN_NO_CHECK(sourceCRS),
+                                           NN_NO_CHECK(targetCRS), ctxt)
+                        .size();
+            } catch (const std::exception &) {
+            }
+        }
     } catch (const std::exception &e) {
         std::cerr << "createOperations() failed with: " << e.what()
                   << std::endl;
@@ -567,8 +589,14 @@ static void outputOperations(
         outputObject(dbContext, list[0], allowUseIntermediateCRS, outputOpt);
         return;
     }
+    std::cout << "Candidate operations found: " << list.size() << std::endl;
+    if (spatialCriterionPartialIntersectionResultCount > list.size()) {
+        std::cout << "Note: using '--spatial-test intersects' would bring "
+                     "more results ("
+                  << spatialCriterionPartialIntersectionResultCount << ")"
+                  << std::endl;
+    }
     if (summary) {
-        std::cout << "Candidate operations found: " << list.size() << std::endl;
         for (const auto &op : list) {
             outputOperationSummary(op);
         }
@@ -576,18 +604,15 @@ static void outputOperations(
         bool first = true;
         for (size_t i = 0; i < list.size(); ++i) {
             const auto &op = list[i];
-            if (list.size() > 1) {
-                if (!first) {
-                    std::cout << std::endl;
-                }
-                first = false;
-                std::cout << "-------------------------------------"
-                          << std::endl;
-                std::cout << "Operation n"
-                             "\xC2\xB0"
-                          << (i + 1) << ":" << std::endl
-                          << std::endl;
+            if (!first) {
+                std::cout << std::endl;
             }
+            first = false;
+            std::cout << "-------------------------------------" << std::endl;
+            std::cout << "Operation n"
+                         "\xC2\xB0"
+                      << (i + 1) << ":" << std::endl
+                      << std::endl;
             outputOperationSummary(op);
             std::cout << std::endl;
             outputObject(dbContext, op, allowUseIntermediateCRS, outputOpt);
@@ -614,6 +639,7 @@ int main(int argc, char **argv) {
     bool summary = false;
     ExtentPtr bboxFilter = nullptr;
     std::string area;
+    bool spatialCriterionExplicitlySpecified = false;
     CoordinateOperationContext::SpatialCriterion spatialCriterion =
         CoordinateOperationContext::SpatialCriterion::STRICT_CONTAINMENT;
     CoordinateOperationContext::SourceTargetCRSExtentUse crsExtentUse =
@@ -622,7 +648,8 @@ int main(int argc, char **argv) {
     CoordinateOperationContext::GridAvailabilityUse gridAvailabilityUse =
         CoordinateOperationContext::GridAvailabilityUse::USE_FOR_SORTING;
     CoordinateOperationContext::IntermediateCRSUse allowUseIntermediateCRS =
-    CoordinateOperationContext::IntermediateCRSUse::IF_NO_DIRECT_TRANSFORMATION;
+        CoordinateOperationContext::IntermediateCRSUse::
+            IF_NO_DIRECT_TRANSFORMATION;
     std::vector<std::pair<std::string, std::string>> pivots;
     bool usePROJGridAlternatives = true;
     std::string mainDBPath;
@@ -768,6 +795,7 @@ int main(int argc, char **argv) {
         } else if (arg == "--spatial-test" && i + 1 < argc) {
             i++;
             std::string value(argv[i]);
+            spatialCriterionExplicitlySpecified = true;
             if (ci_equal(value, "contains")) {
                 spatialCriterion = CoordinateOperationContext::
                     SpatialCriterion::STRICT_CONTAINMENT;
@@ -822,9 +850,10 @@ int main(int argc, char **argv) {
             if (ci_equal(std::string(value), "always")) {
                 allowUseIntermediateCRS =
                     CoordinateOperationContext::IntermediateCRSUse::ALWAYS;
-            } else if (ci_equal(std::string(value), "if_no_direct_transformation")) {
-                allowUseIntermediateCRS =
-                    CoordinateOperationContext::IntermediateCRSUse::IF_NO_DIRECT_TRANSFORMATION;
+            } else if (ci_equal(std::string(value),
+                                "if_no_direct_transformation")) {
+                allowUseIntermediateCRS = CoordinateOperationContext::
+                    IntermediateCRSUse::IF_NO_DIRECT_TRANSFORMATION;
             } else if (ci_equal(std::string(value), "never")) {
                 allowUseIntermediateCRS =
                     CoordinateOperationContext::IntermediateCRSUse::NEVER;
@@ -918,8 +947,8 @@ int main(int argc, char **argv) {
     }
 
     if (outputOpt.quiet &&
-        (outputOpt.PROJ5 + outputOpt.WKT2_2018 +
-         outputOpt.WKT2_2015 + outputOpt.WKT1_GDAL) != 1) {
+        (outputOpt.PROJ5 + outputOpt.WKT2_2018 + outputOpt.WKT2_2015 +
+         outputOpt.WKT1_GDAL) != 1) {
         std::cerr << "-q can only be used with a single output format"
                   << std::endl;
         usage();
@@ -1055,7 +1084,8 @@ int main(int argc, char **argv) {
 
         outputOperations(
             dbContext, sourceCRSStr, targetCRSStr, bboxFilter, spatialCriterion,
-            crsExtentUse, gridAvailabilityUse, allowUseIntermediateCRS, pivots, authority,
+            spatialCriterionExplicitlySpecified, crsExtentUse,
+            gridAvailabilityUse, allowUseIntermediateCRS, pivots, authority,
             usePROJGridAlternatives, showSuperseded, outputOpt, summary);
     }
 
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index b3f200fe..f7dcd354 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -1950,7 +1950,7 @@ void proj_string_list_destroy(PROJ_STRING_LIST list) {
 
 // ---------------------------------------------------------------------------
 
-/** \brief Instanciate a default set of parameters to be used by
+/** \brief Instantiate a default set of parameters to be used by
  * proj_get_crs_list().
  *
  * @return a new object to free with proj_get_crs_list_parameters_destroy() */
@@ -1987,7 +1987,7 @@ void proj_get_crs_list_parameters_destroy(PROJ_CRS_LIST_PARAMETERS *params) {
  * entry is NULL. This array should be freed with proj_crs_info_list_destroy()
  *
  * When no filter parameters are set, this is functionnaly equivalent to
- * proj_get_crs_info_list_from_database(), instanciating a PJ* object for each
+ * proj_get_crs_info_list_from_database(), instantiating a PJ* object for each
  * of the proj_create_from_database() and retrieving information with the
  * various getters. However this function will be much faster.
  *
@@ -5734,7 +5734,7 @@ PJ *proj_create_conversion_equal_earth(PJ_CONTEXT *ctx, double center_long,
  * @return TRUE or FALSE.
  */
 
-int proj_coordoperation_is_instanciable(PJ_CONTEXT *ctx,
+int proj_coordoperation_is_instantiable(PJ_CONTEXT *ctx,
                                         const PJ *coordoperation) {
     assert(coordoperation);
     auto op = dynamic_cast<const CoordinateOperation *>(
@@ -5746,7 +5746,7 @@ int proj_coordoperation_is_instanciable(PJ_CONTEXT *ctx,
     }
     auto dbContext = getDBcontextNoException(ctx, __FUNCTION__);
     try {
-        return op->isPROJInstanciable(dbContext) ? 1 : 0;
+        return op->isPROJInstantiable(dbContext) ? 1 : 0;
     } catch (const std::exception &) {
         return 0;
     }
@@ -5754,6 +5754,36 @@ int proj_coordoperation_is_instanciable(PJ_CONTEXT *ctx,
 
 // ---------------------------------------------------------------------------
 
+/** \brief Return whether a coordinate operation has a "ballpark"
+ * transformation,
+ * that is a very approximate one, due to lack of more accurate transformations.
+ *
+ * Typically a null geographic offset between two horizontal datum, or a
+ * null vertical offset (or limited to unit changes) between two vertical
+ * datum. Errors of several tens to one hundred meters might be expected,
+ * compared to more accurate transformations.
+ *
+ * @param ctx PROJ context, or NULL for default context
+ * @param coordoperation Objet of type CoordinateOperation or derived classes
+ * (must not be NULL)
+ * @return TRUE or FALSE.
+ */
+
+int proj_coordoperation_has_ballpark_transformation(PJ_CONTEXT *ctx,
+                                                    const PJ *coordoperation) {
+    assert(coordoperation);
+    auto op = dynamic_cast<const CoordinateOperation *>(
+        coordoperation->iso_obj.get());
+    if (!op) {
+        proj_log_error(ctx, __FUNCTION__,
+                       "Object is not a CoordinateOperation");
+        return 0;
+    }
+    return op->hasBallparkTransformation();
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Return the number of parameters of a SingleOperation
  *
  * @param ctx PROJ context, or NULL for default context
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 2128124b..224b19ef 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -104,8 +104,16 @@ constexpr double UTM_NORTH_FALSE_NORTHING = 0.0;
 constexpr double UTM_SOUTH_FALSE_NORTHING = 10000000.0;
 
 static const std::string INVERSE_OF = "Inverse of ";
-static const char *NULL_GEOCENTRIC_TRANSLATION = "Null geocentric translation";
-static const char *NULL_GEOGRAPHIC_OFFSET = "Null geographic offset";
+static const char *BALLPARK_GEOCENTRIC_TRANSLATION =
+    "Ballpark geocentric translation";
+static const char *BALLPARK_GEOGRAPHIC_OFFSET = "Ballpark geographic offset";
+static const char *BALLPARK_VERTICAL_TRANSFORMATION_PREFIX =
+    " (ballpark vertical transformation";
+static const char *BALLPARK_VERTICAL_TRANSFORMATION =
+    " (ballpark vertical transformation)";
+static const char *BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT =
+    " (ballpark vertical transformation, without ellipsoid height to vertical "
+    "height correction)";
 //! @endcond
 
 //! @cond Doxygen_Suppress
@@ -513,6 +521,7 @@ struct CoordinateOperation::Private {
     crs::CRSPtr interpolationCRS_{};
     util::optional<common::DataEpoch> sourceCoordinateEpoch_{};
     util::optional<common::DataEpoch> targetCoordinateEpoch_{};
+    bool hasBallparkTransformation_ = false;
 
     // do not set this for a ProjectedCRS.definingConversion
     struct CRSStrongRef {
@@ -706,7 +715,7 @@ void CoordinateOperation::setAccuracies(
  * a PROJ pipeline, checking in particular that referenced grids are
  * available.
  */
-bool CoordinateOperation::isPROJInstanciable(
+bool CoordinateOperation::isPROJInstantiable(
     const io::DatabaseContextPtr &databaseContext) const {
     try {
         exportToPROJString(io::PROJStringFormatter::create().get());
@@ -723,6 +732,27 @@ bool CoordinateOperation::isPROJInstanciable(
 
 // ---------------------------------------------------------------------------
 
+/** \brief Return whether a coordinate operation has a "ballpark"
+ * transformation,
+ * that is a very approximate one, due to lack of more accurate transformations.
+ *
+ * Typically a null geographic offset between two horizontal datum, or a
+ * null vertical offset (or limited to unit changes) between two vertical
+ * datum. Errors of several tens to one hundred meters might be expected,
+ * compared to more accurate transformations.
+ */
+bool CoordinateOperation::hasBallparkTransformation() const {
+    return d->hasBallparkTransformation_;
+}
+
+// ---------------------------------------------------------------------------
+
+void CoordinateOperation::setHasBallparkTransformation(bool b) {
+    d->hasBallparkTransformation_ = b;
+}
+
+// ---------------------------------------------------------------------------
+
 //! @cond Doxygen_Suppress
 struct OperationMethod::Private {
     util::optional<std::string> formula_{};
@@ -1523,11 +1553,12 @@ static SingleOperationNNPtr createPROJBased(
     const util::PropertyMap &properties,
     const io::IPROJStringExportableNNPtr &projExportable,
     const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
-    const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies =
-        std::vector<metadata::PositionalAccuracyNNPtr>()) {
+    const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies,
+    bool hasBallparkTransformation) {
     return util::nn_static_pointer_cast<SingleOperation>(
         PROJBasedOperation::create(properties, projExportable, false, sourceCRS,
-                                   targetCRS, accuracies));
+                                   targetCRS, accuracies,
+                                   hasBallparkTransformation));
 }
 //! @endcond
 
@@ -5644,7 +5675,9 @@ void Conversion::_exportToPROJString(
             common::UnitOfMeasure(std::string(), 1.0 / convFactor,
                                   common::UnitOfMeasure::Type::LINEAR)
                 .exportToPROJString();
-        if (!uom.empty()) {
+        if (uom == "m") {
+            // do nothing
+        } else if (!uom.empty()) {
             formatter->addStep("unitconvert");
             formatter->addParam("z_in", uom);
             formatter->addParam("z_out", "m");
@@ -6141,6 +6174,11 @@ TransformationNNPtr Transformation::create(
         accuracies);
     conv->assignSelf(conv);
     conv->setProperties(properties);
+    std::string name;
+    if (properties.getStringValue(common::IdentifiedObject::NAME_KEY, name) &&
+        ci_find(name, "ballpark") != std::string::npos) {
+        conv->setHasBallparkTransformation(true);
+    }
     return conv;
 }
 
@@ -7054,6 +7092,39 @@ TransformationNNPtr Transformation::createVerticalOffset(
 
 // ---------------------------------------------------------------------------
 
+/** \brief Instantiate a transformation based on the Change of Vertical Unit
+ * method.
+ *
+ * This method is defined as [EPSG:1069]
+ * (https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:method:EPSG::1069)
+ *
+ * @param properties See \ref general_properties of the conversion. If the name
+ * is not provided, it is automatically set.
+ * @param sourceCRSIn Source CRS.
+ * @param targetCRSIn Target CRS.
+ * @param factor Conversion factor
+ * @param accuracies Vector of positional accuracy (might be empty).
+ * @return a new Transformation.
+ */
+TransformationNNPtr Transformation::createChangeVerticalUnit(
+    const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
+    const crs::CRSNNPtr &targetCRSIn, const common::Scale &factor,
+    const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies) {
+    return create(
+        properties, sourceCRSIn, targetCRSIn, nullptr,
+        createMethodMapNameEPSGCode(EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT),
+        VectorOfParameters{
+            createOpParamNameEPSGCode(
+                EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR),
+        },
+        VectorOfValues{
+            factor,
+        },
+        accuracies);
+}
+
+// ---------------------------------------------------------------------------
+
 static const char *getCRSQualifierStr(const crs::CRSPtr &crs) {
     auto geod = dynamic_cast<crs::GeodeticCRS *>(crs.get());
     if (geod) {
@@ -7114,10 +7185,10 @@ createPropertiesForInverse(const CoordinateOperation *op, bool derivedFrom,
     // Forge a name for the inverse, either from the forward name, or
     // from the source and target CRS names
     const char *opType;
-    if (starts_with(forwardName, NULL_GEOCENTRIC_TRANSLATION)) {
-        opType = NULL_GEOCENTRIC_TRANSLATION;
-    } else if (starts_with(forwardName, NULL_GEOGRAPHIC_OFFSET)) {
-        opType = NULL_GEOGRAPHIC_OFFSET;
+    if (starts_with(forwardName, BALLPARK_GEOCENTRIC_TRANSLATION)) {
+        opType = BALLPARK_GEOCENTRIC_TRANSLATION;
+    } else if (starts_with(forwardName, BALLPARK_GEOGRAPHIC_OFFSET)) {
+        opType = BALLPARK_GEOGRAPHIC_OFFSET;
     } else if (dynamic_cast<const Transformation *>(op) ||
                starts_with(forwardName, "Transformation from ")) {
         opType = "Transformation";
@@ -7481,6 +7552,17 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
                                  coordinateOperationAccuracies()));
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT) {
+        const double convFactor = parameterValueNumericAsSI(
+            EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR);
+        return d->registerInv(
+            shared_from_this(),
+            createChangeVerticalUnit(
+                createPropertiesForInverse(this, false, false), l_targetCRS,
+                l_sourceCRS, common::Scale(1.0 / convFactor),
+                coordinateOperationAccuracies()));
+    }
+
     return InverseTransformation::create(NN_NO_CHECK(
         util::nn_dynamic_pointer_cast<Transformation>(shared_from_this())));
 }
@@ -8001,6 +8083,46 @@ TransformationNNPtr Transformation::substitutePROJAlternativeGridNames(
         }
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_VERTCON) {
+        auto fileParameter =
+            parameterValue(EPSG_NAME_PARAMETER_VERTICAL_OFFSET_FILE,
+                           EPSG_CODE_PARAMETER_VERTICAL_OFFSET_FILE);
+        if (fileParameter &&
+            fileParameter->type() == ParameterValue::Type::FILENAME) {
+
+            auto filename = fileParameter->valueFile();
+            if (databaseContext->lookForGridAlternative(
+                    filename, projFilename, projGridFormat, inverseDirection)) {
+
+                if (filename == projFilename) {
+                    assert(!inverseDirection);
+                    return self;
+                }
+
+                auto parameters = std::vector<OperationParameterNNPtr>{
+                    createOpParamNameEPSGCode(
+                        EPSG_CODE_PARAMETER_VERTICAL_OFFSET_FILE)};
+                if (inverseDirection) {
+                    return create(createPropertiesForInverse(
+                                      self.as_nullable().get(), true, false),
+                                  targetCRS(), sourceCRS(), nullptr,
+                                  createSimilarPropertiesMethod(method()),
+                                  parameters, {ParameterValue::createFilename(
+                                                  projFilename)},
+                                  coordinateOperationAccuracies())
+                        ->inverseAsTransformation();
+                } else {
+                    return create(
+                        createSimilarPropertiesTransformation(self),
+                        sourceCRS(), targetCRS(), nullptr,
+                        createSimilarPropertiesMethod(method()), parameters,
+                        {ParameterValue::createFilename(projFilename)},
+                        coordinateOperationAccuracies());
+                }
+            }
+        }
+    }
+
     return self;
 }
 // ---------------------------------------------------------------------------
@@ -8016,7 +8138,7 @@ static void ThrowExceptionNotGeodeticGeographic(const char *trfrm_name) {
 // ---------------------------------------------------------------------------
 
 static void setupPROJGeodeticSourceCRS(io::PROJStringFormatter *formatter,
-                                       const crs::CRSNNPtr &crs,
+                                       const crs::CRSNNPtr &crs, bool addPushV3,
                                        const char *trfrm_name) {
     auto sourceCRSGeog = dynamic_cast<const crs::GeographicCRS *>(crs.get());
     if (sourceCRSGeog) {
@@ -8024,6 +8146,11 @@ static void setupPROJGeodeticSourceCRS(io::PROJStringFormatter *formatter,
         sourceCRSGeog->_exportToPROJString(formatter);
         formatter->stopInversion();
 
+        if (addPushV3) {
+            formatter->addStep("push");
+            formatter->addParam("v_3");
+        }
+
         formatter->addStep("cart");
         sourceCRSGeog->ellipsoid()->_exportToPROJString(formatter);
     } else {
@@ -8039,7 +8166,7 @@ static void setupPROJGeodeticSourceCRS(io::PROJStringFormatter *formatter,
 // ---------------------------------------------------------------------------
 
 static void setupPROJGeodeticTargetCRS(io::PROJStringFormatter *formatter,
-                                       const crs::CRSNNPtr &crs,
+                                       const crs::CRSNNPtr &crs, bool addPopV3,
                                        const char *trfrm_name) {
     auto targetCRSGeog = dynamic_cast<const crs::GeographicCRS *>(crs.get());
     if (targetCRSGeog) {
@@ -8047,6 +8174,11 @@ static void setupPROJGeodeticTargetCRS(io::PROJStringFormatter *formatter,
         formatter->setCurrentStepInverted(true);
         targetCRSGeog->ellipsoid()->_exportToPROJString(formatter);
 
+        if (addPopV3) {
+            formatter->addStep("pop");
+            formatter->addParam("v_3");
+        }
+
         targetCRSGeog->_exportToPROJString(formatter);
     } else {
         auto targetCRSGeod = dynamic_cast<const crs::GeodeticCRS *>(crs.get());
@@ -8138,19 +8270,19 @@ void Transformation::_exportToPROJString(
         double z =
             parameterValueNumericAsSI(EPSG_CODE_PARAMETER_Z_AXIS_TRANSLATION);
 
-        if (methodEPSGCode == EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOGRAPHIC_2D ||
-            methodEPSGCode == EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D) {
-            formatter->addStep("push");
-            formatter->addParam("v_3");
-        }
+        bool addPushPopV3 =
+            (methodEPSGCode ==
+                 EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D ||
+             methodEPSGCode ==
+                 EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOGRAPHIC_2D ||
+             methodEPSGCode == EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D ||
+             methodEPSGCode ==
+                 EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOGRAPHIC_2D ||
+             methodEPSGCode ==
+                 EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D);
 
-        setupPROJGeodeticSourceCRS(formatter, sourceCRS(), "Helmert");
+        setupPROJGeodeticSourceCRS(formatter, sourceCRS(), addPushPopV3,
+                                   "Helmert");
 
         formatter->addStep("helmert");
         formatter->addParam("x", x);
@@ -8214,19 +8346,8 @@ void Transformation::_exportToPROJString(
             }
         }
 
-        setupPROJGeodeticTargetCRS(formatter, targetCRS(), "Helmert");
-
-        if (methodEPSGCode == EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOGRAPHIC_2D ||
-            methodEPSGCode == EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D) {
-            formatter->addStep("pop");
-            formatter->addParam("v_3");
-        }
+        setupPROJGeodeticTargetCRS(formatter, targetCRS(), addPushPopV3,
+                                   "Helmert");
 
         return;
     }
@@ -8274,15 +8395,13 @@ void Transformation::_exportToPROJString(
         double pz = parameterValueNumericAsSI(
             EPSG_CODE_PARAMETER_ORDINATE_3_EVAL_POINT);
 
-        if (methodEPSGCode ==
-                EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_PV_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_CF_GEOGRAPHIC_2D) {
-            formatter->addStep("push");
-            formatter->addParam("v_3");
-        }
+        bool addPushPopV3 =
+            (methodEPSGCode ==
+                 EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_PV_GEOGRAPHIC_2D ||
+             methodEPSGCode ==
+                 EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_CF_GEOGRAPHIC_2D);
 
-        setupPROJGeodeticSourceCRS(formatter, sourceCRS(),
+        setupPROJGeodeticSourceCRS(formatter, sourceCRS(), addPushPopV3,
                                    "Molodensky-Badekas");
 
         formatter->addStep("molobadekas");
@@ -8302,17 +8421,9 @@ void Transformation::_exportToPROJString(
             formatter->addParam("convention", "coordinate_frame");
         }
 
-        setupPROJGeodeticTargetCRS(formatter, targetCRS(),
+        setupPROJGeodeticTargetCRS(formatter, targetCRS(), addPushPopV3,
                                    "Molodensky-Badekas");
 
-        if (methodEPSGCode ==
-                EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_PV_GEOGRAPHIC_2D ||
-            methodEPSGCode ==
-                EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_CF_GEOGRAPHIC_2D) {
-            formatter->addStep("pop");
-            formatter->addParam("v_3");
-        }
-
         return;
     }
 
@@ -8795,6 +8906,33 @@ bool SingleOperation::exportToPROJStringGeneric(
                                       "conversion");
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT) {
+        double convFactor = parameterValueNumericAsSI(
+            EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR);
+        auto uom = common::UnitOfMeasure(std::string(), convFactor,
+                                         common::UnitOfMeasure::Type::LINEAR)
+                       .exportToPROJString();
+        auto reverse_uom =
+            common::UnitOfMeasure(std::string(), 1.0 / convFactor,
+                                  common::UnitOfMeasure::Type::LINEAR)
+                .exportToPROJString();
+        if (uom == "m") {
+            // do nothing
+        } else if (!uom.empty()) {
+            formatter->addStep("unitconvert");
+            formatter->addParam("z_in", uom);
+            formatter->addParam("z_out", "m");
+        } else if (!reverse_uom.empty()) {
+            formatter->addStep("unitconvert");
+            formatter->addParam("z_in", "m");
+            formatter->addParam("z_out", reverse_uom);
+        } else {
+            formatter->addStep("affine");
+            formatter->addParam("s33", convFactor);
+        }
+        return true;
+    }
+
     return false;
 }
 
@@ -9080,7 +9218,9 @@ CoordinateOperationNNPtr ConcatenatedOperation::createComputeMetadata(
     }
 
     std::vector<CoordinateOperationNNPtr> flattenOps;
+    bool hasBallparkTransformation = false;
     for (const auto &subOp : operationsIn) {
+        hasBallparkTransformation |= subOp->hasBallparkTransformation();
         auto subOpConcat =
             dynamic_cast<const ConcatenatedOperation *>(subOp.get());
         if (subOpConcat) {
@@ -9120,6 +9260,7 @@ CoordinateOperationNNPtr ConcatenatedOperation::createComputeMetadata(
     }
 
     auto op = create(properties, flattenOps, accuracies);
+    op->setHasBallparkTransformation(hasBallparkTransformation);
     op->d->computedName_ = true;
     return op;
 }
@@ -9627,14 +9768,18 @@ struct PrecomputedOpCharacteristics {
     bool gridsAvailable_ = false;
     bool gridsKnown_ = false;
     size_t stepCount_ = 0;
+    bool isApprox_ = false;
+    bool isNullTransformation_ = false;
 
     PrecomputedOpCharacteristics() = default;
     PrecomputedOpCharacteristics(double area, double accuracy, bool hasGrids,
                                  bool gridsAvailable, bool gridsKnown,
-                                 size_t stepCount)
+                                 size_t stepCount, bool isApprox,
+                                 bool isNullTransformation)
         : area_(area), accuracy_(accuracy), hasGrids_(hasGrids),
           gridsAvailable_(gridsAvailable), gridsKnown_(gridsKnown),
-          stepCount_(stepCount) {}
+          stepCount_(stepCount), isApprox_(isApprox),
+          isNullTransformation_(isNullTransformation) {}
 };
 
 // ---------------------------------------------------------------------------
@@ -9661,6 +9806,22 @@ struct SortFunction {
         // CAUTION: the order of the comparisons is extremely important
         // to get the intended result.
 
+        if (!iterA->second.isApprox_ && iterB->second.isApprox_) {
+            return true;
+        }
+        if (iterA->second.isApprox_ && !iterB->second.isApprox_) {
+            return false;
+        }
+
+        if (!iterA->second.isNullTransformation_ &&
+            iterB->second.isNullTransformation_) {
+            return true;
+        }
+        if (iterA->second.isNullTransformation_ &&
+            !iterB->second.isNullTransformation_) {
+            return false;
+        }
+
         if (iterA->second.hasGrids_ && iterB->second.hasGrids_) {
             // Operations where grids are all available go before other
             if (iterA->second.gridsAvailable_ &&
@@ -9691,6 +9852,17 @@ struct SortFunction {
             return false;
         }
 
+        if (accuracyA < 0 && accuracyB < 0) {
+            // unknown accuracy ? then prefer operations with grids, which
+            // are likely to have best practical accuracy
+            if (iterA->second.hasGrids_ && !iterB->second.hasGrids_) {
+                return true;
+            }
+            if (!iterA->second.hasGrids_ && iterB->second.hasGrids_) {
+                return false;
+            }
+        }
+
         // Operations with larger non-zero area of use go before those with
         // lower one
         const double areaA = iterA->second.area_;
@@ -9722,15 +9894,6 @@ struct SortFunction {
             if (iterA->second.hasGrids_ && !iterB->second.hasGrids_) {
                 return false;
             }
-        } else if (accuracyA < 0 && accuracyB < 0) {
-            // unknown accuracy ? then prefer operations with grids, which
-            // are likely to have best practical accuracy
-            if (iterA->second.hasGrids_ && !iterB->second.hasGrids_) {
-                return true;
-            }
-            if (!iterA->second.hasGrids_ && iterB->second.hasGrids_) {
-                return false;
-            }
         }
 
         // The less intermediate steps, the better
@@ -9882,11 +10045,7 @@ struct FilterResults {
                 bool extentContains =
                     extent->contains(NN_NO_CHECK(areaOfInterest));
                 if (extentContains) {
-                    const auto &name = op->nameStr();
-                    if (name.find(NULL_GEOGRAPHIC_OFFSET) ==
-                            std::string::npos &&
-                        name.find(NULL_GEOCENTRIC_TRANSLATION) ==
-                            std::string::npos) {
+                    if (!op->hasBallparkTransformation()) {
                         hasOpThatContainsAreaOfInterest = true;
                     }
                 }
@@ -9917,11 +10076,7 @@ struct FilterResults {
                     !targetCRSExtent ||
                     extent->contains(NN_NO_CHECK(targetCRSExtent));
                 if (extentContainsSource && extentContainsTarget) {
-                    const auto &name = op->nameStr();
-                    if (name.find(NULL_GEOGRAPHIC_OFFSET) ==
-                            std::string::npos &&
-                        name.find(NULL_GEOCENTRIC_TRANSLATION) ==
-                            std::string::npos) {
+                    if (!op->hasBallparkTransformation()) {
                         hasOpThatContainsAreaOfInterest = true;
                     }
                 }
@@ -10002,13 +10157,7 @@ struct FilterResults {
             bool hasGrids = false;
             bool gridsAvailable = true;
             bool gridsKnown = true;
-            if (context->getAuthorityFactory() &&
-                (gridAvailabilityUse ==
-                     CoordinateOperationContext::GridAvailabilityUse::
-                         USE_FOR_SORTING ||
-                 gridAvailabilityUse ==
-                     CoordinateOperationContext::GridAvailabilityUse::
-                         IGNORE_GRID_AVAILABILITY)) {
+            if (context->getAuthorityFactory()) {
                 const auto gridsNeeded = op->gridsNeeded(
                     context->getAuthorityFactory()->databaseContext());
                 for (const auto &gridDesc : gridsNeeded) {
@@ -10027,9 +10176,17 @@ struct FilterResults {
 
             const auto stepCount = getStepCount(op);
 
+            const bool isApprox =
+                op->nameStr().find(BALLPARK_VERTICAL_TRANSFORMATION_PREFIX) !=
+                std::string::npos;
+            const bool isNullTransformation =
+                op->nameStr().find(BALLPARK_GEOGRAPHIC_OFFSET) !=
+                    std::string::npos ||
+                op->nameStr().find(BALLPARK_GEOCENTRIC_TRANSLATION) !=
+                    std::string::npos;
             map[op.get()] = PrecomputedOpCharacteristics(
                 area, getAccuracy(op), hasGrids, gridsAvailable, gridsKnown,
-                stepCount);
+                stepCount, isApprox, isNullTransformation);
         }
 
         // Sort !
@@ -10041,13 +10198,16 @@ struct FilterResults {
     void removeSyntheticNullTransforms() {
 
         // If we have more than one result, and than the last result is the
-        // default "Null geographic offset" or "Null geocentric translation"
-        // operations we have synthetized, remove it as
+        // default "Ballpark geographic offset" or "Ballpark geocentric
+        // translation"
+        // operations we have synthetized, and that at least one operation
+        // has the desired area of interest, remove it as
         // all previous results are necessarily better
         if (hasOpThatContainsAreaOfInterest && res.size() > 1) {
             const std::string &name = res.back()->nameStr();
-            if (name.find(NULL_GEOGRAPHIC_OFFSET) != std::string::npos ||
-                name.find(NULL_GEOCENTRIC_TRANSLATION) != std::string::npos) {
+            if (name.find(BALLPARK_GEOGRAPHIC_OFFSET) != std::string::npos ||
+                name.find(BALLPARK_GEOCENTRIC_TRANSLATION) !=
+                    std::string::npos) {
                 std::vector<CoordinateOperationNNPtr> resTemp;
                 for (size_t i = 0; i < res.size() - 1; i++) {
                     resTemp.emplace_back(res[i]);
@@ -10460,9 +10620,9 @@ static std::vector<CoordinateOperationNNPtr> findsOpsInRegistryWithIntermediate(
 
 //! @cond Doxygen_Suppress
 static TransformationNNPtr
-createNullGeographicOffset(const crs::CRSNNPtr &sourceCRS,
-                           const crs::CRSNNPtr &targetCRS) {
-    std::string name(NULL_GEOGRAPHIC_OFFSET);
+createBallparkGeographicOffset(const crs::CRSNNPtr &sourceCRS,
+                               const crs::CRSNNPtr &targetCRS) {
+    std::string name(BALLPARK_GEOGRAPHIC_OFFSET);
     name += " from ";
     name += sourceCRS->nameStr();
     name += " to ";
@@ -10653,7 +10813,7 @@ createGeodToGeodPROJBased(const crs::CRSNNPtr &geodSrc,
     auto properties = util::PropertyMap().set(
         common::IdentifiedObject::NAME_KEY,
         buildTransfName(geodSrc->nameStr(), geodDst->nameStr()));
-    return createPROJBased(properties, exportable, geodSrc, geodDst);
+    return createPROJBased(properties, exportable, geodSrc, geodDst, {}, false);
 }
 
 // ---------------------------------------------------------------------------
@@ -10690,7 +10850,9 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased(
     }
 
     return createPROJBased(properties, exportable, sourceCRS, targetCRS,
-                           accuracies);
+                           accuracies,
+                           horizTransform->hasBallparkTransformation() ||
+                               verticalTransform->hasBallparkTransformation());
 }
 
 // ---------------------------------------------------------------------------
@@ -10708,8 +10870,17 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased(
             interpolationGeogCRS);
 
     bool dummy = false;
-    auto ops = std::vector<CoordinateOperationNNPtr>{
-        opSrcCRSToGeogCRS, verticalTransform, opGeogCRStoDstCRS};
+    auto ops = opSrcCRSToGeogCRS->sourceCRS()->_isEquivalentTo(
+                   opSrcCRSToGeogCRS->targetCRS().get())
+                   ? std::vector<CoordinateOperationNNPtr>{verticalTransform,
+                                                           opGeogCRStoDstCRS}
+                   : std::vector<CoordinateOperationNNPtr>{opSrcCRSToGeogCRS,
+                                                           verticalTransform,
+                                                           opGeogCRStoDstCRS};
+    bool hasBallparkTransformation = false;
+    for (const auto &op : ops) {
+        hasBallparkTransformation |= op->hasBallparkTransformation();
+    }
     auto extent = getExtent(ops, true, dummy);
     auto properties = util::PropertyMap();
     properties.set(common::IdentifiedObject::NAME_KEY,
@@ -10728,7 +10899,7 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased(
     }
 
     return createPROJBased(properties, exportable, sourceCRS, targetCRS,
-                           accuracies);
+                           accuracies, hasBallparkTransformation);
 }
 
 //! @endcond
@@ -10784,6 +10955,11 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog(
 
     std::string name(buildTransfName(geogSrc->nameStr(), geogDst->nameStr()));
 
+    const bool sameDatum =
+        geogSrc->datum() != nullptr && geogDst->datum() != nullptr &&
+        geogSrc->datum()->_isEquivalentTo(
+            geogDst->datum().get(), util::IComparable::Criterion::EQUIVALENT);
+
     // Do they differ by vertical units ?
     if (vconvSrc != vconvDst &&
         geogSrc->ellipsoid()->_isEquivalentTo(
@@ -10799,18 +10975,19 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog(
                                         name),
                 common::Scale(factor));
             conv->setCRSs(sourceCRS, targetCRS, nullptr);
+            conv->setHasBallparkTransformation(!sameDatum);
             res.push_back(conv);
             return res;
         } else {
-            res.emplace_back(createGeodToGeodPROJBased(sourceCRS, targetCRS));
+            auto op = createGeodToGeodPROJBased(sourceCRS, targetCRS);
+            op->setHasBallparkTransformation(!sameDatum);
+            res.emplace_back(op);
             return res;
         }
     }
 
     // Do the CRS differ only by their axis order ?
-    if (geogSrc->datum() != nullptr && geogDst->datum() != nullptr &&
-        geogSrc->datum()->_isEquivalentTo(
-            geogDst->datum().get(), util::IComparable::Criterion::EQUIVALENT) &&
+    if (sameDatum &&
         !srcCS->_isEquivalentTo(dstCS.get(),
                                 util::IComparable::Criterion::EQUIVALENT)) {
         auto srcOrder = srcCS->axisOrder();
@@ -10871,7 +11048,8 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog(
                          metadata::Extent::WORLD),
                 datum, dstCS));
 
-        steps.emplace_back(createNullGeographicOffset(sourceCRS, interm_crs));
+        steps.emplace_back(
+            createBallparkGeographicOffset(sourceCRS, interm_crs));
 
         steps.emplace_back(Transformation::createLongitudeRotation(
             util::PropertyMap()
@@ -10906,15 +11084,17 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog(
                          metadata::Extent::WORLD),
                 sourceCRS, interm_crs, offset_pm));
             steps.emplace_back(
-                createNullGeographicOffset(interm_crs, targetCRS));
+                createBallparkGeographicOffset(interm_crs, targetCRS));
         } else {
             steps.emplace_back(
-                createNullGeographicOffset(sourceCRS, targetCRS));
+                createBallparkGeographicOffset(sourceCRS, targetCRS));
         }
     }
 
-    res.emplace_back(ConcatenatedOperation::createComputeMetadata(
-        steps, !allowEmptyIntersection));
+    auto op = ConcatenatedOperation::createComputeMetadata(
+        steps, !allowEmptyIntersection);
+    op->setHasBallparkTransformation(!sameDatum);
+    res.emplace_back(op);
     return res;
 }
 
@@ -10966,8 +11146,8 @@ findCandidateGeodCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
 
 static bool isNullTransformation(const std::string &name) {
 
-    return starts_with(name, NULL_GEOCENTRIC_TRANSLATION) ||
-           starts_with(name, NULL_GEOGRAPHIC_OFFSET);
+    return starts_with(name, BALLPARK_GEOCENTRIC_TRANSLATION) ||
+           starts_with(name, BALLPARK_GEOGRAPHIC_OFFSET);
 }
 
 // ---------------------------------------------------------------------------
@@ -11085,9 +11265,9 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
 // ---------------------------------------------------------------------------
 
 static CoordinateOperationNNPtr
-createNullGeocentricTranslation(const crs::CRSNNPtr &sourceCRS,
-                                const crs::CRSNNPtr &targetCRS) {
-    std::string name(NULL_GEOCENTRIC_TRANSLATION);
+createBallparkGeocentricTranslation(const crs::CRSNNPtr &sourceCRS,
+                                    const crs::CRSNNPtr &targetCRS) {
+    std::string name(BALLPARK_GEOCENTRIC_TRANSLATION);
     name += " from ";
     name += sourceCRS->nameStr();
     name += " to ";
@@ -11309,7 +11489,7 @@ CoordinateOperationFactory::Private::createOperations(
                             util::nn_dynamic_pointer_cast<cs::CartesianCS>(
                                 geodSrc->coordinateSystem()))));
                 auto opFirst =
-                    createNullGeocentricTranslation(sourceCRS, interm_crs);
+                    createBallparkGeocentricTranslation(sourceCRS, interm_crs);
                 auto opSecond =
                     createGeographicGeocentric(interm_crs, targetCRS);
                 res.emplace_back(ConcatenatedOperation::createComputeMetadata(
@@ -11324,7 +11504,7 @@ CoordinateOperationFactory::Private::createOperations(
 
         if (isSrcGeocentric && isTargetGeocentric) {
             res.emplace_back(
-                createNullGeocentricTranslation(sourceCRS, targetCRS));
+                createBallparkGeocentricTranslation(sourceCRS, targetCRS));
             return res;
         }
 
@@ -11526,30 +11706,35 @@ CoordinateOperationFactory::Private::createOperations(
     if (vertSrc && vertDst) {
         const auto &srcDatum = vertSrc->datum();
         const auto &dstDatum = vertDst->datum();
-        if (srcDatum && dstDatum &&
-            srcDatum->_isEquivalentTo(
-                dstDatum.get(), util::IComparable::Criterion::EQUIVALENT)) {
-            const double convSrc = vertSrc->coordinateSystem()
-                                       ->axisList()[0]
-                                       ->unit()
-                                       .conversionToSI();
-            const double convDst = vertDst->coordinateSystem()
-                                       ->axisList()[0]
-                                       ->unit()
-                                       .conversionToSI();
-            if (convSrc != convDst) {
-                const double factor = convSrc / convDst;
-                auto conv = Conversion::createChangeVerticalUnit(
-                    util::PropertyMap().set(
-                        common::IdentifiedObject::NAME_KEY,
-                        buildTransfName(sourceCRS->nameStr(),
-                                        targetCRS->nameStr())),
-                    common::Scale(factor));
-                conv->setCRSs(sourceCRS, targetCRS, nullptr);
-                res.push_back(conv);
-                return res;
-            }
+        const bool equivalentVDatum =
+            (srcDatum && dstDatum &&
+             srcDatum->_isEquivalentTo(
+                 dstDatum.get(), util::IComparable::Criterion::EQUIVALENT));
+
+        const double convSrc =
+            vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI();
+        const double convDst =
+            vertDst->coordinateSystem()->axisList()[0]->unit().conversionToSI();
+
+        const double factor = convSrc / convDst;
+        auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
+        if (!equivalentVDatum) {
+            name += BALLPARK_VERTICAL_TRANSFORMATION;
+            auto conv = Transformation::createChangeVerticalUnit(
+                util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                        name),
+                sourceCRS, targetCRS, common::Scale(factor), {});
+            conv->setHasBallparkTransformation(true);
+            res.push_back(conv);
+        } else if (convSrc != convDst) {
+            auto conv = Conversion::createChangeVerticalUnit(
+                util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                        name),
+                common::Scale(factor));
+            conv->setCRSs(sourceCRS, targetCRS, nullptr);
+            res.push_back(conv);
         }
+        return res;
     }
 
     // A bit odd case as we are comparing apples to oranges, but in case
@@ -11562,17 +11747,17 @@ CoordinateOperationFactory::Private::createOperations(
         if (geogAxis.size() == 3) {
             convDst = geogAxis[2]->unit().conversionToSI();
         }
-        if (convSrc != convDst) {
-            const double factor = convSrc / convDst;
-            auto conv = Conversion::createChangeVerticalUnit(
-                util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
-                                        buildTransfName(sourceCRS->nameStr(),
-                                                        targetCRS->nameStr())),
-                common::Scale(factor));
-            conv->setCRSs(sourceCRS, targetCRS, nullptr);
-            res.push_back(conv);
-            return res;
-        }
+
+        const double factor = convSrc / convDst;
+        auto conv = Transformation::createChangeVerticalUnit(
+            util::PropertyMap().set(
+                common::IdentifiedObject::NAME_KEY,
+                buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) +
+                    BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT),
+            sourceCRS, targetCRS, common::Scale(factor), {});
+        conv->setHasBallparkTransformation(true);
+        res.push_back(conv);
+        return res;
     }
 
     // reverse of previous case
@@ -11877,7 +12062,8 @@ PROJBasedOperationNNPtr PROJBasedOperation::create(
     const util::PropertyMap &properties,
     const io::IPROJStringExportableNNPtr &projExportable, bool inverse,
     const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
-    const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies) {
+    const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies,
+    bool hasBallparkTransformation) {
 
     auto formatter = io::PROJStringFormatter::create();
     if (inverse) {
@@ -11903,6 +12089,7 @@ PROJBasedOperationNNPtr PROJBasedOperation::create(
     op->setAccuracies(accuracies);
     op->projStringExportable_ = projExportable.as_nullable();
     op->inverse_ = inverse;
+    op->setHasBallparkTransformation(hasBallparkTransformation);
     return op;
 }
 
@@ -11916,7 +12103,7 @@ CoordinateOperationNNPtr PROJBasedOperation::inverse() const {
                 createPropertiesForInverse(this, false, false),
                 NN_NO_CHECK(projStringExportable_), !inverse_,
                 NN_NO_CHECK(targetCRS()), NN_NO_CHECK(sourceCRS()),
-                coordinateOperationAccuracies()));
+                coordinateOperationAccuracies(), hasBallparkTransformation()));
     }
 
     auto formatter = io::PROJStringFormatter::create();
@@ -11929,11 +12116,11 @@ CoordinateOperationNNPtr PROJBasedOperation::inverse() const {
     }
     formatter->stopInversion();
 
-    return util::nn_static_pointer_cast<CoordinateOperation>(
-        PROJBasedOperation::create(
-            createPropertiesForInverse(this, false, false),
-            formatter->toString(), targetCRS(), sourceCRS(),
-            coordinateOperationAccuracies()));
+    auto op = PROJBasedOperation::create(
+        createPropertiesForInverse(this, false, false), formatter->toString(),
+        targetCRS(), sourceCRS(), coordinateOperationAccuracies());
+    op->setHasBallparkTransformation(hasBallparkTransformation());
+    return util::nn_static_pointer_cast<CoordinateOperation>(op);
 }
 
 // ---------------------------------------------------------------------------
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 01a588e3..aabb15a1 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -453,7 +453,7 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
                 auto transf =
                     util::nn_dynamic_pointer_cast<operation::Transformation>(
                         op);
-                if (transf && !starts_with(transf->nameStr(), "Null geo")) {
+                if (transf && !starts_with(transf->nameStr(), "Ballpark geo")) {
                     try {
                         transf->getTOWGS84Parameters();
                     } catch (const std::exception &) {
@@ -486,7 +486,7 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
                                     operation::Transformation>(subops[1]);
                                 if (transf &&
                                     !starts_with(transf->nameStr(),
-                                                 "Null geo")) {
+                                                 "Ballpark geo")) {
                                     try {
                                         transf->getTOWGS84Parameters();
                                     } catch (const std::exception &) {
@@ -3855,7 +3855,7 @@ BoundCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                 for (const auto &op : ops) {
                     std::string opTransfPROJString;
                     bool opTransfPROJStringValid = false;
-                    if (op->nameStr().find("Null geographic") == 0) {
+                    if (op->nameStr().find("Ballpark geographic") == 0) {
                         if (isTOWGS84Compatible()) {
                             auto params =
                                 transformation()->getTOWGS84Parameters();
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index 81abcdf1..bafa1b5d 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -3873,7 +3873,7 @@ AuthorityFactory::getDescriptionText(const std::string &code) const {
 /** \brief Return a list of information on CRS objects
  *
  * This is functionnaly equivalent of listing the codes from an authority,
- * instanciating
+ * instantiating
  * a CRS object for each of them and getting the information from this CRS
  * object, but this implementation has much less overhead.
  *
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index e14239b0..a1608464 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -5102,71 +5102,6 @@ const std::string &PROJStringFormatter::toString() const {
                 }
             }
 
-            // axisswap order=2,1, pop/push v_3, axisswap order=2,1 -> can
-            // suppress axisswap
-            if (i + 1 < d->steps_.size() && prevStep.name == "axisswap" &&
-                (curStep.name == "push" || curStep.name == "pop") &&
-                prevStepParamCount == 1 &&
-                prevStep.paramValues[0].equals("order", "2,1") &&
-                curStepParamCount == 1 && curStep.paramValues[0].key == "v_3") {
-                auto iterNext = iterCur;
-                ++iterNext;
-                auto &nextStep = *iterNext;
-                if (nextStep.name == "axisswap" &&
-                    nextStep.paramValues.size() == 1 &&
-                    nextStep.paramValues[0].equals("order", "2,1")) {
-                    d->steps_.erase(iterPrev);
-                    d->steps_.erase(iterNext);
-                    changeDone = true;
-                    break;
-                }
-            }
-
-            // push v_3, axisswap order=2,1, pop v_3 -> can suppress push/pop
-            if (i + 1 < d->steps_.size() && prevStep.name == "push" &&
-                prevStepParamCount == 1 &&
-                prevStep.paramValues[0].key == "v_3" &&
-                curStep.name == "axisswap" && curStepParamCount == 1 &&
-                curStep.paramValues[0].equals("order", "2,1")) {
-                auto iterNext = iterCur;
-                ++iterNext;
-                auto &nextStep = *iterNext;
-                if (nextStep.name == "pop" &&
-                    nextStep.paramValues.size() == 1 &&
-                    nextStep.paramValues[0].key == "v_3") {
-                    d->steps_.erase(iterPrev);
-                    d->steps_.erase(iterNext);
-                    changeDone = true;
-                    break;
-                }
-            }
-
-            // unitconvert xy_in=A xy_out=B, pop/push v_3, unitconvert xy_in=B
-            // xy_out=A -> can suppress unitconvert
-            if (i + 1 < d->steps_.size() && prevStep.name == "unitconvert" &&
-                (curStep.name == "push" || curStep.name == "pop") &&
-                prevStepParamCount == 2 &&
-                prevStep.paramValues[0].key == "xy_in" &&
-                prevStep.paramValues[1].key == "xy_out" &&
-                curStepParamCount == 1 && curStep.paramValues[0].key == "v_3") {
-                auto iterNext = iterCur;
-                ++iterNext;
-                auto &nextStep = *iterNext;
-                if (nextStep.name == "unitconvert" &&
-                    nextStep.paramValues.size() == 2 &&
-                    nextStep.paramValues[0].key == "xy_in" &&
-                    nextStep.paramValues[1].key == "xy_out" &&
-                    nextStep.paramValues[0].value ==
-                        prevStep.paramValues[1].value &&
-                    nextStep.paramValues[1].value ==
-                        prevStep.paramValues[0].value) {
-                    d->steps_.erase(iterPrev);
-                    d->steps_.erase(iterNext);
-                    changeDone = true;
-                    break;
-                }
-            }
-
             // for practical purposes WGS84 and GRS80 ellipsoids are
             // equivalents (cartesian transform between both lead to differences
             // of the order of 1e-14 deg..).
diff --git a/src/proj.h b/src/proj.h
index af22c341..1f71e115 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -999,7 +999,10 @@ int PROJ_DLL proj_coordoperation_get_method_info(PJ_CONTEXT *ctx,
                                                  const char **out_method_auth_name,
                                                  const char **out_method_code);
 
-int PROJ_DLL proj_coordoperation_is_instanciable(PJ_CONTEXT *ctx,
+int PROJ_DLL proj_coordoperation_is_instantiable(PJ_CONTEXT *ctx,
+                                                 const PJ *coordoperation);
+
+int PROJ_DLL proj_coordoperation_has_ballpark_transformation(PJ_CONTEXT *ctx,
                                                  const PJ *coordoperation);
 
 int PROJ_DLL proj_coordoperation_get_param_count(PJ_CONTEXT *ctx,
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 14b69492..448b65c8 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -203,7 +203,7 @@ PJ_COORD PROJ_DLL pj_approx_3D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD co
 /* Grid functionality */
 int             proj_vgrid_init(PJ *P, const char *grids);
 int             proj_hgrid_init(PJ *P, const char *grids);
-double          proj_vgrid_value(PJ *P, PJ_LP lp);
+double          proj_vgrid_value(PJ *P, PJ_LP lp, double vmultiplier);
 PJ_LP           proj_hgrid_value(PJ *P, PJ_LP lp);
 PJ_LP           proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction);
 
diff --git a/src/proj_symbol_rename.h b/src/proj_symbol_rename.h
index 0ce342fd..b6e887ca 100644
--- a/src/proj_symbol_rename.h
+++ b/src/proj_symbol_rename.h
@@ -115,7 +115,7 @@
 #define proj_coordoperation_get_param_count internal_proj_coordoperation_get_param_count
 #define proj_coordoperation_get_param_index internal_proj_coordoperation_get_param_index
 #define proj_coordoperation_get_towgs84_values internal_proj_coordoperation_get_towgs84_values
-#define proj_coordoperation_is_instanciable internal_proj_coordoperation_is_instanciable
+#define proj_coordoperation_is_instantiable internal_proj_coordoperation_is_instantiable
 #define proj_create internal_proj_create
 #define proj_create_argv internal_proj_create_argv
 #define proj_create_cartesian_2D_cs internal_proj_create_cartesian_2D_cs
diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp
index c28e1489..0e0d641c 100644
--- a/src/transformations/deformation.cpp
+++ b/src/transformations/deformation.cpp
@@ -92,7 +92,7 @@ static PJ_XYZ get_grid_shift(PJ* P, PJ_XYZ cartesian) {
 
     /* look up correction values in grids */
     shift.lp    = proj_hgrid_value(P, geodetic.lp);
-    shift.enu.u = proj_vgrid_value(P, geodetic.lp);
+    shift.enu.u = proj_vgrid_value(P, geodetic.lp, 1.0);
 
     if (proj_errno(P) == PJD_ERR_GRID_AREA)
         proj_log_debug(P, "deformation: coordinate (%.3f, %.3f) outside deformation model",
diff --git a/src/transformations/vgridshift.cpp b/src/transformations/vgridshift.cpp
index fda38ec3..4cd48fb6 100644
--- a/src/transformations/vgridshift.cpp
+++ b/src/transformations/vgridshift.cpp
@@ -26,7 +26,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
     if (P->vgridlist_geoid != nullptr) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.xyz.z += Q->forward_multiplier * proj_vgrid_value(P, point.lp);
+        point.xyz.z += proj_vgrid_value(P, point.lp, Q->forward_multiplier);
     }
 
     return point.xyz;
@@ -41,7 +41,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
     if (P->vgridlist_geoid != nullptr) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.xyz.z -= Q->forward_multiplier * proj_vgrid_value(P, point.lp);
+        point.xyz.z -= proj_vgrid_value(P, point.lp, Q->forward_multiplier);
     }
 
     return point.lpz;