Merge pull request #1279 from rouault/vertcrs_transform_improvements

Vertcrs transform improvements

24 files changed, 986 insertions, 530 deletions

diff --git a/data/sql/grid_alternatives.sql b/data/sql/grid_alternatives.sql
index a45c43f1..99f20d0e 100644
--- a/data/sql/grid_alternatives.sql
+++ b/data/sql/grid_alternatives.sql
@@ -142,6 +142,53 @@ INSERT INTO grid_alternatives(original_grid_name,
                               'proj-datumgrid-europe',
                               NULL, NULL, NULL, NULL);
 
+-- Continental USA VERTCON: NGVD (19)29 height to NAVD (19)88 height
+
+INSERT INTO grid_alternatives(original_grid_name,
+                              proj_grid_name,
+                              proj_grid_format,
+                              proj_method,
+                              inverse_direction,
+                              package_name,
+                              url, direct_download, open_license, directory)
+                      VALUES ('vertconw.94',
+                              'vertconw.gtx',
+                              'GTX',
+                              'vgridshift',
+                              0,
+                              'proj-datumgrid-north-america',
+                              NULL, NULL, NULL, NULL);
+
+INSERT INTO grid_alternatives(original_grid_name,
+                              proj_grid_name,
+                              proj_grid_format,
+                              proj_method,
+                              inverse_direction,
+                              package_name,
+                              url, direct_download, open_license, directory)
+                      VALUES ('vertconc.94',
+                              'vertconc.gtx',
+                              'GTX',
+                              'vgridshift',
+                              0,
+                              'proj-datumgrid-north-america',
+                              NULL, NULL, NULL, NULL);
+
+INSERT INTO grid_alternatives(original_grid_name,
+                              proj_grid_name,
+                              proj_grid_format,
+                              proj_method,
+                              inverse_direction,
+                              package_name,
+                              url, direct_download, open_license, directory)
+                      VALUES ('vertcone.94',
+                              'vertcone.gtx',
+                              'GTX',
+                              'vgridshift',
+                              0,
+                              'proj-datumgrid-north-america',
+                              NULL, NULL, NULL, NULL);
+
 -- EGM models
 
 INSERT INTO grid_alternatives(original_grid_name,
diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp
index 92b655f9..fd737dcc 100644
--- a/include/proj/coordinateoperation.hpp
+++ b/include/proj/coordinateoperation.hpp
@@ -146,7 +146,9 @@ class PROJ_GCC_DLL CoordinateOperation : public common::ObjectUsage,
     gridsNeeded(const io::DatabaseContextPtr &databaseContext) const = 0;
 
     PROJ_DLL bool
-    isPROJInstanciable(const io::DatabaseContextPtr &databaseContext) const;
+    isPROJInstantiable(const io::DatabaseContextPtr &databaseContext) const;
+
+    PROJ_DLL bool hasBallparkTransformation() const;
 
   protected:
     PROJ_INTERNAL CoordinateOperation();
@@ -167,6 +169,7 @@ class PROJ_GCC_DLL CoordinateOperation : public common::ObjectUsage,
     PROJ_INTERNAL
     void setAccuracies(
         const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies);
+    PROJ_INTERNAL void setHasBallparkTransformation(bool b);
 
   private:
     PROJ_OPAQUE_PRIVATE_DATA
@@ -1493,6 +1496,11 @@ class PROJ_GCC_DLL Transformation : public SingleOperation {
     PROJ_DLL TransformationNNPtr substitutePROJAlternativeGridNames(
         io::DatabaseContextNNPtr databaseContext) const;
 
+    PROJ_DLL static TransformationNNPtr createChangeVerticalUnit(
+        const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
+        const crs::CRSNNPtr &targetCRSIn, const common::Scale &factor,
+        const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies);
+
     PROJ_PRIVATE :
         //! @cond Doxygen_Suppress
         PROJ_INTERNAL const std::string &
diff --git a/include/proj/internal/coordinateoperation_internal.hpp b/include/proj/internal/coordinateoperation_internal.hpp
index 8428b8bf..65b786b2 100644
--- a/include/proj/internal/coordinateoperation_internal.hpp
+++ b/include/proj/internal/coordinateoperation_internal.hpp
@@ -246,7 +246,8 @@ class PROJBasedOperation : public SingleOperation {
     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 hasRoughTransformation);
 
     std::set<GridDescription>
     gridsNeeded(const io::DatabaseContextPtr &databaseContext) const override;
diff --git a/scripts/reformat_cpp.sh b/scripts/reformat_cpp.sh
index b3140ebe..cc10b53e 100755
--- a/scripts/reformat_cpp.sh
+++ b/scripts/reformat_cpp.sh
@@ -15,7 +15,7 @@ esac
 
 TOPDIR="$SCRIPT_DIR/.."
 
-for i in "$TOPDIR"/include/proj/*.hpp "$TOPDIR"/include/proj/internal/*.hpp "$TOPDIR"/src/iso19111/*.cpp "$TOPDIR"/test/unit/*.cpp; do
+for i in "$TOPDIR"/include/proj/*.hpp "$TOPDIR"/include/proj/internal/*.hpp "$TOPDIR"/src/iso19111/*.cpp "$TOPDIR"/test/unit/*.cpp "$TOPDIR"/src/apps/projinfo.cpp; do
     if ! echo "$i" | grep -q "lru_cache.hpp"; then
         "$SCRIPT_DIR"/reformat.sh "$i";
     fi
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;
diff --git a/test/cli/testprojinfo b/test/cli/testprojinfo
index 111c071e..24a1fdd5 100755
--- a/test/cli/testprojinfo
+++ b/test/cli/testprojinfo
@@ -53,6 +53,10 @@ echo "Testing projinfo -s EPSG:4326 -t EPSG:32631" >> ${OUT}
 $EXE -s EPSG:4326 -t EPSG:32631 >>${OUT}
 echo "" >>${OUT}
 
+echo "Testing projinfo -s NAD27 -t NAD83" >> ${OUT}
+$EXE -s NAD27 -t NAD83 >>${OUT}
+echo "" >>${OUT}
+
 echo "Testing projinfo -s NAD27 -t NAD83 --grid-check none --spatial-test intersects --summary" >> ${OUT}
 $EXE -s NAD27 -t NAD83 --grid-check none --spatial-test intersects --summary >>${OUT}
 echo "" >>${OUT}
diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist
index ebc59c40..674e9631 100644
--- a/test/cli/testprojinfo_out.dist
+++ b/test/cli/testprojinfo_out.dist
@@ -136,6 +136,10 @@ GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.25722
 
 
 Testing projinfo -s EPSG:4326 -t EPSG:32631
+Candidate operations found: 1
+-------------------------------------
+Operation n°1:
+
 EPSG:16031, UTM zone 31N, 0 m, World - N hemisphere - 0°E to 6°E
 
 PROJ string:
@@ -162,17 +166,82 @@ CONVERSION["UTM zone 31N",
         ID["EPSG",8807]],
     ID["EPSG",16031]]
 
+Testing projinfo -s NAD27 -t NAD83
+Candidate operations found: 1
+Note: using '--spatial-test intersects' would bring more results (7)
+-------------------------------------
+Operation n°1:
+
+unknown id, Ballpark geographic offset from NAD27 to NAD83, unknown accuracy, World, has ballpark transformation
+
+PROJ string:
+
+
+WKT2_2018 string:
+COORDINATEOPERATION["Ballpark geographic offset from NAD27 to NAD83",
+    SOURCECRS[
+        GEOGCRS["NAD27",
+            DATUM["North American Datum 1927",
+                ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
+                    LENGTHUNIT["metre",1]]],
+            PRIMEM["Greenwich",0,
+                ANGLEUNIT["degree",0.0174532925199433]],
+            CS[ellipsoidal,2],
+                AXIS["geodetic latitude (Lat)",north,
+                    ORDER[1],
+                    ANGLEUNIT["degree",0.0174532925199433]],
+                AXIS["geodetic longitude (Lon)",east,
+                    ORDER[2],
+                    ANGLEUNIT["degree",0.0174532925199433]],
+            USAGE[
+                SCOPE["unknown"],
+                AREA["North America - NAD27"],
+                BBOX[7.15,167.65,83.17,-47.74]],
+            ID["EPSG",4267]]],
+    TARGETCRS[
+        GEOGCRS["NAD83",
+            DATUM["North American Datum 1983",
+                ELLIPSOID["GRS 1980",6378137,298.257222101,
+                    LENGTHUNIT["metre",1]]],
+            PRIMEM["Greenwich",0,
+                ANGLEUNIT["degree",0.0174532925199433]],
+            CS[ellipsoidal,2],
+                AXIS["geodetic latitude (Lat)",north,
+                    ORDER[1],
+                    ANGLEUNIT["degree",0.0174532925199433]],
+                AXIS["geodetic longitude (Lon)",east,
+                    ORDER[2],
+                    ANGLEUNIT["degree",0.0174532925199433]],
+            USAGE[
+                SCOPE["unknown"],
+                AREA["North America - NAD83"],
+                BBOX[14.92,167.65,86.46,-47.74]],
+            ID["EPSG",4269]]],
+    METHOD["Geographic2D offsets",
+        ID["EPSG",9619]],
+    PARAMETER["Latitude offset",0,
+        ANGLEUNIT["degree",0.0174532925199433],
+        ID["EPSG",8601]],
+    PARAMETER["Longitude offset",0,
+        ANGLEUNIT["degree",0.0174532925199433],
+        ID["EPSG",8602]],
+    USAGE[
+        SCOPE["unknown"],
+        AREA["World"],
+        BBOX[-90,-180,90,180]]]
+
 Testing projinfo -s NAD27 -t NAD83 --grid-check none --spatial-test intersects --summary
 Candidate operations found: 7
 DERIVED_FROM(EPSG):1312, NAD27 to NAD83 (3), 1.0 m, Canada
 DERIVED_FROM(EPSG):1313, NAD27 to NAD83 (4), 1.5 m, Canada - NAD27
 DERIVED_FROM(EPSG):1241, NAD27 to NAD83 (1), 0.15 m, USA - CONUS including EEZ
 DERIVED_FROM(EPSG):1243, NAD27 to NAD83 (2), 0.5 m, USA - Alaska including EEZ
-unknown id, Null geographic offset from NAD27 to NAD83, unknown accuracy, World
 EPSG:1462, NAD27 to NAD83 (5), 1.0 m, Canada - Quebec
 EPSG:1573, NAD27 to NAD83 (6), 1.5 m, Canada - Quebec
+unknown id, Ballpark geographic offset from NAD27 to NAD83, unknown accuracy, World, has ballpark transformation
 
 Testing projinfo -s NAD27 -t NAD83 --grid-check none --spatial-test intersects
+Candidate operations found: 7
 -------------------------------------
 Operation n°1:
 
@@ -366,13 +435,13 @@ COORDINATEOPERATION["NAD27 to NAD83 (2)",
 -------------------------------------
 Operation n°5:
 
-unknown id, Null geographic offset from NAD27 to NAD83, unknown accuracy, World
+EPSG:1462, NAD27 to NAD83 (5), 1.0 m, Canada - Quebec
 
 PROJ string:
-
++proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=GS2783v1.QUE +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 WKT2_2018 string:
-COORDINATEOPERATION["Null geographic offset from NAD27 to NAD83",
+COORDINATEOPERATION["NAD27 to NAD83 (5)",
     SOURCECRS[
         GEOGCRS["NAD27",
             DATUM["North American Datum 1927",
@@ -386,12 +455,7 @@ COORDINATEOPERATION["Null geographic offset from NAD27 to NAD83",
                     ANGLEUNIT["degree",0.0174532925199433]],
                 AXIS["geodetic longitude (Lon)",east,
                     ORDER[2],
-                    ANGLEUNIT["degree",0.0174532925199433]],
-            USAGE[
-                SCOPE["unknown"],
-                AREA["North America - NAD27"],
-                BBOX[7.15,167.65,83.17,-47.74]],
-            ID["EPSG",4267]]],
+                    ANGLEUNIT["degree",0.0174532925199433]]]],
     TARGETCRS[
         GEOGCRS["NAD83",
             DATUM["North American Datum 1983",
@@ -405,35 +469,27 @@ COORDINATEOPERATION["Null geographic offset from NAD27 to NAD83",
                     ANGLEUNIT["degree",0.0174532925199433]],
                 AXIS["geodetic longitude (Lon)",east,
                     ORDER[2],
-                    ANGLEUNIT["degree",0.0174532925199433]],
-            USAGE[
-                SCOPE["unknown"],
-                AREA["North America - NAD83"],
-                BBOX[14.92,167.65,86.46,-47.74]],
-            ID["EPSG",4269]]],
-    METHOD["Geographic2D offsets",
-        ID["EPSG",9619]],
-    PARAMETER["Latitude offset",0,
-        ANGLEUNIT["degree",0.0174532925199433],
-        ID["EPSG",8601]],
-    PARAMETER["Longitude offset",0,
-        ANGLEUNIT["degree",0.0174532925199433],
-        ID["EPSG",8602]],
+                    ANGLEUNIT["degree",0.0174532925199433]]]],
+    METHOD["NTv1",
+        ID["EPSG",9614]],
+    PARAMETERFILE["Latitude and longitude difference file","GS2783v1.QUE"],
+    OPERATIONACCURACY[1.0],
     USAGE[
         SCOPE["unknown"],
-        AREA["World"],
-        BBOX[-90,-180,90,180]]]
+        AREA["Canada - Quebec"],
+        BBOX[44.99,-79.85,62.62,-57.1]],
+    ID["EPSG",1462]]
 
 -------------------------------------
 Operation n°6:
 
-EPSG:1462, NAD27 to NAD83 (5), 1.0 m, Canada - Quebec
+EPSG:1573, NAD27 to NAD83 (6), 1.5 m, Canada - Quebec
 
 PROJ string:
-+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=GS2783v1.QUE +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
++proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=QUE27-83.gsb +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 WKT2_2018 string:
-COORDINATEOPERATION["NAD27 to NAD83 (5)",
+COORDINATEOPERATION["NAD27 to NAD83 (6)",
     SOURCECRS[
         GEOGCRS["NAD27",
             DATUM["North American Datum 1927",
@@ -462,26 +518,26 @@ COORDINATEOPERATION["NAD27 to NAD83 (5)",
                 AXIS["geodetic longitude (Lon)",east,
                     ORDER[2],
                     ANGLEUNIT["degree",0.0174532925199433]]]],
-    METHOD["NTv1",
-        ID["EPSG",9614]],
-    PARAMETERFILE["Latitude and longitude difference file","GS2783v1.QUE"],
-    OPERATIONACCURACY[1.0],
+    METHOD["NTv2",
+        ID["EPSG",9615]],
+    PARAMETERFILE["Latitude and longitude difference file","QUE27-83.gsb"],
+    OPERATIONACCURACY[1.5],
     USAGE[
         SCOPE["unknown"],
         AREA["Canada - Quebec"],
         BBOX[44.99,-79.85,62.62,-57.1]],
-    ID["EPSG",1462]]
+    ID["EPSG",1573]]
 
 -------------------------------------
 Operation n°7:
 
-EPSG:1573, NAD27 to NAD83 (6), 1.5 m, Canada - Quebec
+unknown id, Ballpark geographic offset from NAD27 to NAD83, unknown accuracy, World, has ballpark transformation
 
 PROJ string:
-+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=QUE27-83.gsb +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
+
 
 WKT2_2018 string:
-COORDINATEOPERATION["NAD27 to NAD83 (6)",
+COORDINATEOPERATION["Ballpark geographic offset from NAD27 to NAD83",
     SOURCECRS[
         GEOGCRS["NAD27",
             DATUM["North American Datum 1927",
@@ -495,7 +551,12 @@ COORDINATEOPERATION["NAD27 to NAD83 (6)",
                     ANGLEUNIT["degree",0.0174532925199433]],
                 AXIS["geodetic longitude (Lon)",east,
                     ORDER[2],
-                    ANGLEUNIT["degree",0.0174532925199433]]]],
+                    ANGLEUNIT["degree",0.0174532925199433]],
+            USAGE[
+                SCOPE["unknown"],
+                AREA["North America - NAD27"],
+                BBOX[7.15,167.65,83.17,-47.74]],
+            ID["EPSG",4267]]],
     TARGETCRS[
         GEOGCRS["NAD83",
             DATUM["North American Datum 1983",
@@ -509,16 +570,24 @@ COORDINATEOPERATION["NAD27 to NAD83 (6)",
                     ANGLEUNIT["degree",0.0174532925199433]],
                 AXIS["geodetic longitude (Lon)",east,
                     ORDER[2],
-                    ANGLEUNIT["degree",0.0174532925199433]]]],
-    METHOD["NTv2",
-        ID["EPSG",9615]],
-    PARAMETERFILE["Latitude and longitude difference file","QUE27-83.gsb"],
-    OPERATIONACCURACY[1.5],
+                    ANGLEUNIT["degree",0.0174532925199433]],
+            USAGE[
+                SCOPE["unknown"],
+                AREA["North America - NAD83"],
+                BBOX[14.92,167.65,86.46,-47.74]],
+            ID["EPSG",4269]]],
+    METHOD["Geographic2D offsets",
+        ID["EPSG",9619]],
+    PARAMETER["Latitude offset",0,
+        ANGLEUNIT["degree",0.0174532925199433],
+        ID["EPSG",8601]],
+    PARAMETER["Longitude offset",0,
+        ANGLEUNIT["degree",0.0174532925199433],
+        ID["EPSG",8602]],
     USAGE[
         SCOPE["unknown"],
-        AREA["Canada - Quebec"],
-        BBOX[44.99,-79.85,62.62,-57.1]],
-    ID["EPSG",1573]]
+        AREA["World"],
+        BBOX[-90,-180,90,180]]]
 
 Testing projinfo -s EPSG:4230 -t EPSG:4258 --bbox 8,54.51,15.24,57.8 --summary
 Candidate operations found: 1
diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp
index ad637786..3f99b1b3 100644
--- a/test/unit/gie_self_tests.cpp
+++ b/test/unit/gie_self_tests.cpp
@@ -421,9 +421,11 @@ TEST(gie, info_functions) {
     /* check a few key characteristics of the Mercator projection */
     EXPECT_NEAR(factors.angular_distortion, 0.0, 1e-7)
         << factors.angular_distortion; /* angular distortion should be 0 */
+
+    /* Meridian/parallel angle should be 90 deg */
     EXPECT_NEAR(factors.meridian_parallel_angle, M_PI_2, 1e-7)
-        << factors.meridian_parallel_angle; /* Meridian/parallel angle should be
-                                               90 deg */
+        << factors.meridian_parallel_angle;
+
     EXPECT_EQ(factors.meridian_convergence,
               0.0); /* meridian convergence should be 0 */
 
@@ -701,13 +703,16 @@ TEST(gie, proj_create_crs_to_crs_PULKOVO42_ETRS89) {
     EXPECT_NEAR(c.xy.x, 44.999701238, 1e-9);
     EXPECT_NEAR(c.xy.y, 24.998474948, 1e-9);
     EXPECT_EQ(std::string(proj_pj_info(P).definition),
-              "proj=pipeline step proj=push v_3 step proj=axisswap order=2,1 "
-              "step proj=unitconvert xy_in=deg xy_out=rad step proj=cart "
+              "proj=pipeline step proj=axisswap order=2,1 "
+              "step proj=unitconvert xy_in=deg xy_out=rad "
+              "step proj=push v_3 "
+              "step proj=cart "
               "ellps=krass step proj=helmert x=2.3287 y=-147.0425 z=-92.0802 "
               "rx=0.3092483 ry=-0.32482185 rz=-0.49729934 s=5.68906266 "
-              "convention=coordinate_frame step inv proj=cart ellps=GRS80 step "
-              "proj=unitconvert xy_in=rad xy_out=deg step proj=axisswap "
-              "order=2,1 step proj=pop v_3");
+              "convention=coordinate_frame step inv proj=cart ellps=GRS80 "
+              "step proj=pop v_3 "
+              "step proj=unitconvert xy_in=rad xy_out=deg step proj=axisswap "
+              "order=2,1");
 
     c = proj_trans(P, PJ_INV, c);
     EXPECT_NEAR(c.xy.x, 45, 1e-8);
@@ -730,12 +735,15 @@ TEST(gie, proj_create_crs_to_crs_PULKOVO42_ETRS89) {
     EXPECT_NEAR(c.xy.x, 51.999714150, 1e-9);
     EXPECT_NEAR(c.xy.y, 19.998187811, 1e-9);
     EXPECT_EQ(std::string(proj_pj_info(P).definition),
-              "proj=pipeline step proj=push v_3 step proj=axisswap order=2,1 "
-              "step proj=unitconvert xy_in=deg xy_out=rad step proj=cart "
+              "proj=pipeline step proj=axisswap order=2,1 "
+              "step proj=unitconvert xy_in=deg xy_out=rad "
+              "step proj=push v_3 "
+              "step proj=cart "
               "ellps=krass step proj=helmert x=33.4 y=-146.6 z=-76.3 rx=-0.359 "
               "ry=-0.053 rz=0.844 s=-0.84 convention=position_vector step inv "
-              "proj=cart ellps=GRS80 step proj=unitconvert xy_in=rad "
-              "xy_out=deg step proj=axisswap order=2,1 step proj=pop v_3");
+              "proj=cart ellps=GRS80 step proj=pop v_3 "
+              "step proj=unitconvert xy_in=rad "
+              "xy_out=deg step proj=axisswap order=2,1");
 
     proj_destroy(P);
 }
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 6205a9b8..486ab0c7 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -838,14 +838,16 @@ TEST_F(CApi, proj_create_from_database) {
         ASSERT_NE(info.definition, nullptr);
         EXPECT_EQ(
             info.definition,
-            std::string("proj=pipeline step proj=push v_3 step proj=axisswap "
+            std::string("proj=pipeline step proj=axisswap "
                         "order=2,1 step proj=unitconvert xy_in=deg xy_out=rad "
+                        "step proj=push v_3 "
                         "step proj=cart ellps=bessel step proj=helmert "
                         "x=601.705 y=84.263 z=485.227 rx=-4.7354 ry=-1.3145 "
-                        "rz=-5.393 s=-2.3887 convention=coordinate_frame step "
-                        "inv proj=cart ellps=GRS80 step proj=unitconvert "
-                        "xy_in=rad xy_out=deg step proj=axisswap order=2,1 "
-                        "step proj=pop v_3"));
+                        "rz=-5.393 s=-2.3887 convention=coordinate_frame "
+                        "step inv proj=cart ellps=GRS80 "
+                        "step proj=pop v_3 "
+                        "step proj=unitconvert xy_in=rad xy_out=deg "
+                        "step proj=axisswap order=2,1"));
         EXPECT_EQ(info.accuracy, 1);
     }
 }
@@ -1300,13 +1302,13 @@ TEST_F(CApi, proj_coordoperation_get_grid_used) {
 
 // ---------------------------------------------------------------------------
 
-TEST_F(CApi, proj_coordoperation_is_instanciable) {
+TEST_F(CApi, proj_coordoperation_is_instantiable) {
     auto op = proj_create_from_database(m_ctxt, "EPSG", "1671",
                                         PJ_CATEGORY_COORDINATE_OPERATION, true,
                                         nullptr);
     ASSERT_NE(op, nullptr);
     ObjectKeeper keeper(op);
-    EXPECT_EQ(proj_coordoperation_is_instanciable(m_ctxt, op), 1);
+    EXPECT_EQ(proj_coordoperation_is_instantiable(m_ctxt, op), 1);
 }
 
 // ---------------------------------------------------------------------------
@@ -1343,6 +1345,7 @@ TEST_F(CApi, proj_create_operations) {
     auto op = proj_list_get(m_ctxt, res, 0);
     ASSERT_NE(op, nullptr);
     ObjectKeeper keeper_op(op);
+    EXPECT_FALSE(proj_coordoperation_has_ballpark_transformation(m_ctxt, op));
 
     EXPECT_EQ(proj_get_name(op), std::string("NAD27 to NAD83 (3)"));
 }
@@ -1399,8 +1402,9 @@ TEST_F(CApi, proj_create_operations_with_pivot) {
         ASSERT_NE(op, nullptr);
         ObjectKeeper keeper_op(op);
 
-        EXPECT_EQ(proj_get_name(op),
-                  std::string("Null geographic offset from WGS 84 to JGD2011"));
+        EXPECT_EQ(
+            proj_get_name(op),
+            std::string("Ballpark geographic offset from WGS 84 to JGD2011"));
     }
 
     // Restrict pivot to Tokyo CRS
@@ -1421,7 +1425,7 @@ TEST_F(CApi, proj_create_operations_with_pivot) {
         ASSERT_NE(res, nullptr);
         ObjListKeeper keeper_res(res);
         EXPECT_EQ(proj_list_get_count(res), 7);
-        auto op = proj_list_get(m_ctxt, res, 1);
+        auto op = proj_list_get(m_ctxt, res, 0);
         ASSERT_NE(op, nullptr);
         ObjectKeeper keeper_op(op);
 
@@ -1451,7 +1455,7 @@ TEST_F(CApi, proj_create_operations_with_pivot) {
         ASSERT_NE(res, nullptr);
         ObjListKeeper keeper_res(res);
         // includes results from ESRI
-        EXPECT_EQ(proj_list_get_count(res), 5);
+        EXPECT_EQ(proj_list_get_count(res), 4);
         auto op = proj_list_get(m_ctxt, res, 0);
         ASSERT_NE(op, nullptr);
         ObjectKeeper keeper_op(op);
diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp
index 80de017f..944e0ebe 100644
--- a/test/unit/test_factory.cpp
+++ b/test/unit/test_factory.cpp
@@ -660,13 +660,13 @@ TEST(factory, AuthorityFactory_createCoordinateOperation_helmert_3) {
                  NoSuchAuthorityCodeException);
     auto op = factory->createCoordinateOperation("1113", false);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-              "+inv +proj=longlat +a=6378249.145 +rf=293.4663077 +step "
-              "+proj=cart +a=6378249.145 +rf=293.4663077 +step +proj=helmert "
-              "+x=-143 +y=-90 +z=-294 +step +inv +proj=cart +ellps=WGS84 +step "
-              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
-              "+order=2,1 +step +proj=pop +v_3");
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +inv "
+              "+proj=longlat +a=6378249.145 +rf=293.4663077 +step +proj=push "
+              "+v_3 +step +proj=cart +a=6378249.145 +rf=293.4663077 +step "
+              "+proj=helmert +x=-143 +y=-90 +z=-294 +step +inv +proj=cart "
+              "+ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert "
+              "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
 }
 
 // ---------------------------------------------------------------------------
@@ -675,13 +675,13 @@ TEST(factory, AuthorityFactory_createCoordinateOperation_helmert_7_CF) {
     auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
     auto op = factory->createCoordinateOperation("7676", false);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-              "+proj=cart +ellps=bessel +step +proj=helmert +x=577.88891 "
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=bessel +step +proj=helmert +x=577.88891 "
               "+y=165.22205 +z=391.18289 +rx=-4.9145 +ry=0.94729 +rz=13.05098 "
               "+s=7.78664 +convention=coordinate_frame +step +inv +proj=cart "
-              "+ellps=WGS84 +step +proj=unitconvert +xy_in=rad +xy_out=deg "
-              "+step +proj=axisswap +order=2,1 +step +proj=pop +v_3");
+              "+ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert "
+              "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
 }
 
 // ---------------------------------------------------------------------------
@@ -852,15 +852,15 @@ TEST(factory,
     EXPECT_TRUE(so->validateParameters().empty());
 
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-              "+proj=cart +ellps=bessel +step +proj=molobadekas +x=593.032 "
-              "+y=26 +z=478.741 +rx=0.409394387439237 +ry=-0.359705195614311 "
-              "+rz=1.86849100035057 +s=4.0772 +px=3903453.148 +py=368135.313 "
-              "+pz=5012970.306 +convention=coordinate_frame +step +inv "
-              "+proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad "
-              "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop "
-              "+v_3");
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=bessel +step +proj=molobadekas "
+              "+x=593.032 +y=26 +z=478.741 +rx=0.409394387439237 "
+              "+ry=-0.359705195614311 +rz=1.86849100035057 +s=4.0772 "
+              "+px=3903453.148 +py=368135.313 +pz=5012970.306 "
+              "+convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 "
+              "+step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad "
+              "+xy_out=deg +step +proj=axisswap +order=2,1");
 }
 
 // ---------------------------------------------------------------------------
@@ -2071,12 +2071,12 @@ TEST_F(FactoryWithTmpDatabase, AuthorityFactory_wkt_based_transformation) {
     ASSERT_EQ(res.size(), 1U);
     EXPECT_EQ(res[0]->nameStr(), "My WKT string based op");
     EXPECT_EQ(res[0]->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-              "+proj=cart +ellps=WGS84 +step +proj=helmert +x=1 +y=2 +z=3 "
-              "+step +inv +proj=cart +ellps=WGS84 +step +proj=unitconvert "
-              "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 +step "
-              "+proj=pop +v_3");
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=WGS84 +step +proj=helmert +x=1 +y=2 "
+              "+z=3 +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
+              "+proj=axisswap +order=2,1");
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 90deb661..112b46e3 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -565,12 +565,12 @@ TEST(operation, transformation_createGeocentricTranslations) {
     EXPECT_EQ(inv_transf_as_transf->getTOWGS84Parameters(), expected_inv);
 
     EXPECT_EQ(transf->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-              "+proj=cart +ellps=GRS80 +step +proj=helmert +x=1 +y=2 +z=3 "
-              "+step +inv +proj=cart +ellps=WGS84 +step +proj=unitconvert "
-              "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 +step "
-              "+proj=pop +v_3");
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=GRS80 +step +proj=helmert +x=1 +y=2 "
+              "+z=3 +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
+              "+proj=axisswap +order=2,1");
 }
 
 // ---------------------------------------------------------------------------
@@ -673,13 +673,13 @@ TEST(operation, transformation_createPositionVector) {
     EXPECT_EQ(transf->getTOWGS84Parameters(), expected);
 
     EXPECT_EQ(transf->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-              "+proj=cart +ellps=GRS80 +step +proj=helmert +x=1 +y=2 +z=3 "
-              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
-              "+proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
-              "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop "
-              "+v_3");
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=GRS80 +step +proj=helmert +x=1 +y=2 "
+              "+z=3 +rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step "
+              "+inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step "
+              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
+              "+order=2,1");
 
     auto inv_transf = transf->inverse();
     ASSERT_EQ(inv_transf->coordinateOperationAccuracies().size(), 1U);
@@ -696,12 +696,12 @@ TEST(operation, transformation_createPositionVector) {
 #else
     EXPECT_EQ(
         inv_transf->exportToPROJString(PROJStringFormatter::create().get()),
-        "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap +order=2,1 "
-        "+step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-        "+ellps=WGS84 +step +inv +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 "
-        "+rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
-        "+ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
-        "+proj=axisswap +order=2,1 +step +proj=pop +v_3");
+        "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+        "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step "
+        "+proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=1 +y=2 +z=3 +rx=4 "
+        "+ry=5 +rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
+        "+ellps=GRS80 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad "
+        "+xy_out=deg +step +proj=axisswap +order=2,1");
 
     // In WKT, use approximate formula
     auto wkt = inv_transf->exportToWKT(WKTFormatter::create().get());
@@ -742,13 +742,13 @@ TEST(operation, transformation_createCoordinateFrameRotation) {
     EXPECT_EQ(params, expected);
 
     EXPECT_EQ(transf->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-              "+proj=cart +ellps=GRS80 +step +proj=helmert +x=1 +y=2 +z=3 "
-              "+rx=-4 +ry=-5 +rz=-6 +s=7 +convention=coordinate_frame +step "
-              "+inv +proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
-              "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop "
-              "+v_3");
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=GRS80 +step +proj=helmert +x=1 +y=2 "
+              "+z=3 +rx=-4 +ry=-5 +rz=-6 +s=7 +convention=coordinate_frame "
+              "+step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step "
+              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
+              "+order=2,1");
 
     auto inv_transf = transf->inverse();
     ASSERT_EQ(inv_transf->coordinateOperationAccuracies().size(), 0U);
@@ -765,12 +765,12 @@ TEST(operation, transformation_createCoordinateFrameRotation) {
 #else
     EXPECT_EQ(
         inv_transf->exportToPROJString(PROJStringFormatter::create().get()),
-        "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap +order=2,1 "
-        "+step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-        "+ellps=WGS84 +step +inv +proj=helmert +x=1 +y=2 +z=3 +rx=-4 +ry=-5 "
-        "+rz=-6 +s=7 +convention=coordinate_frame +step +inv +proj=cart "
-        "+ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
-        "+proj=axisswap +order=2,1 +step +proj=pop +v_3");
+        "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+        "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step "
+        "+proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=1 +y=2 +z=3 "
+        "+rx=-4 +ry=-5 +rz=-6 +s=7 +convention=coordinate_frame +step +inv "
+        "+proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=unitconvert "
+        "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
 
     // In WKT, use approximate formula
     auto wkt = inv_transf->exportToWKT(WKTFormatter::create().get());
@@ -4199,17 +4199,18 @@ TEST(operation, geogCRS_to_geogCRS_context_default) {
         EXPECT_EQ(list[0]->getEPSGCode(), 15994); // Romania - 3m
         EXPECT_EQ(list[1]->getEPSGCode(), 1644);  // Poland - 1m
         EXPECT_EQ(list[2]->nameStr(),
-                  "Null geographic offset from Pulkovo 1942(58) to ETRS89");
+                  "Ballpark geographic offset from Pulkovo 1942(58) to ETRS89");
 
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-            "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-            "+proj=cart +ellps=krass +step +proj=helmert +x=2.3287 "
+            "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+            "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+            "+step +proj=cart +ellps=krass +step +proj=helmert +x=2.3287 "
             "+y=-147.0425 +z=-92.0802 +rx=0.3092483 +ry=-0.32482185 "
             "+rz=-0.49729934 +s=5.68906266 +convention=coordinate_frame +step "
-            "+inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad "
-            "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop +v_3");
+            "+inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step "
+            "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
+            "+order=2,1");
     }
 
     // Reverse case
@@ -4224,13 +4225,14 @@ TEST(operation, geogCRS_to_geogCRS_context_default) {
 
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-            "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-            "+proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=2.3287 "
+            "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+            "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+            "+step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=2.3287 "
             "+y=-147.0425 +z=-92.0802 +rx=0.3092483 +ry=-0.32482185 "
             "+rz=-0.49729934 +s=5.68906266 +convention=coordinate_frame +step "
-            "+inv +proj=cart +ellps=krass +step +proj=unitconvert +xy_in=rad "
-            "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop +v_3");
+            "+inv +proj=cart +ellps=krass +step +proj=pop +v_3 +step "
+            "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
+            "+order=2,1");
     }
 }
 
@@ -4362,26 +4364,20 @@ TEST(operation, geogCRS_to_geogCRS_context_inverse_needed) {
             authFactory->createCoordinateReferenceSystem("4275"), // NTF
             authFactory->createCoordinateReferenceSystem("4258"), // ETRS89
             ctxt);
-        ASSERT_EQ(list.size(), 3U);
+        ASSERT_EQ(list.size(), 2U);
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-            "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-            "+proj=cart +ellps=clrk80ign +step +proj=helmert +x=-168 +y=-60 "
-            "+z=320 +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert "
-            "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 +step "
-            "+proj=pop +v_3");
+            "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+            "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+            "+step +proj=cart +ellps=clrk80ign +step +proj=helmert +x=-168 "
+            "+y=-60 +z=320 +step +inv +proj=cart +ellps=GRS80 +step +proj=pop "
+            "+v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
+            "+proj=axisswap +order=2,1");
         EXPECT_EQ(list[1]->exportToPROJString(
                       PROJStringFormatter::create(
                           PROJStringFormatter::Convention::PROJ_5,
                           authFactory->databaseContext())
                           .get()),
-                  "");
-        EXPECT_EQ(list[2]->exportToPROJString(
-                      PROJStringFormatter::create(
-                          PROJStringFormatter::Convention::PROJ_5,
-                          authFactory->databaseContext())
-                          .get()),
                   "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
                   "+proj=unitconvert +xy_in=deg +xy_out=rad +step "
                   "+proj=hgridshift +grids=ntf_r93.gsb +step +proj=unitconvert "
@@ -4511,7 +4507,8 @@ TEST(operation, geogCRS_to_geogCRS_noop) {
     auto op = CoordinateOperationFactory::create()->createOperation(
         GeographicCRS::EPSG_4326, GeographicCRS::EPSG_4326);
     ASSERT_TRUE(op != nullptr);
-    EXPECT_EQ(op->nameStr(), "Null geographic offset from WGS 84 to WGS 84");
+    EXPECT_EQ(op->nameStr(),
+              "Ballpark geographic offset from WGS 84 to WGS 84");
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), "");
     EXPECT_EQ(op->inverse()->nameStr(), op->nameStr());
 }
@@ -4772,9 +4769,10 @@ TEST(operation, geocentricCRS_to_geogCRS_different_datum) {
     auto op = CoordinateOperationFactory::create()->createOperation(
         createGeocentricDatumWGS84(), GeographicCRS::EPSG_4269);
     ASSERT_TRUE(op != nullptr);
-    EXPECT_EQ(op->nameStr(), "Null geocentric translation from WGS 84 to NAD83 "
-                             "(geocentric) + Conversion from NAD83 "
-                             "(geocentric) to NAD83");
+    EXPECT_EQ(op->nameStr(),
+              "Ballpark geocentric translation from WGS 84 to NAD83 "
+              "(geocentric) + Conversion from NAD83 "
+              "(geocentric) to NAD83");
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +inv +proj=cart +ellps=GRS80 +step "
               "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
@@ -4789,7 +4787,7 @@ TEST(operation, geogCRS_to_geocentricCRS_different_datum) {
         GeographicCRS::EPSG_4269, createGeocentricDatumWGS84());
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->nameStr(), "Conversion from NAD83 to NAD83 (geocentric) + "
-                             "Null geocentric translation from NAD83 "
+                             "Ballpark geocentric translation from NAD83 "
                              "(geocentric) to WGS 84");
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
@@ -4805,7 +4803,7 @@ TEST(operation, geocentricCRS_to_geocentricCRS_noop) {
         createGeocentricDatumWGS84(), createGeocentricDatumWGS84());
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->nameStr(),
-              "Null geocentric translation from WGS 84 to WGS 84");
+              "Ballpark geocentric translation from WGS 84 to WGS 84");
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), "");
     EXPECT_EQ(op->inverse()->nameStr(), op->nameStr());
 }
@@ -5324,14 +5322,14 @@ TEST(operation, boundCRS_of_geogCRS_to_geogCRS) {
         boundCRS, GeographicCRS::EPSG_4326);
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=grad +xy_out=rad "
-              "+step +inv +proj=longlat +ellps=clrk80ign +pm=paris +step "
-              "+proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 +z=3 "
-              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
-              "+proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
-              "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop "
-              "+v_3");
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=grad +xy_out=rad +step +inv "
+              "+proj=longlat +ellps=clrk80ign +pm=paris +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 "
+              "+z=3 +rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step "
+              "+inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step "
+              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
+              "+order=2,1");
 }
 
 // ---------------------------------------------------------------------------
@@ -5345,13 +5343,13 @@ TEST(operation, boundCRS_of_geogCRS_to_geogCRS_with_area) {
         boundCRS, authFactory->createCoordinateReferenceSystem("4326"));
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-              "+proj=cart +ellps=clrk66 +step +proj=helmert +x=1 +y=2 +z=3 "
-              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
-              "+proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
-              "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop "
-              "+v_3");
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=clrk66 +step +proj=helmert +x=1 +y=2 "
+              "+z=3 +rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step "
+              "+inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step "
+              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
+              "+order=2,1");
 }
 
 // ---------------------------------------------------------------------------
@@ -5387,8 +5385,8 @@ TEST(operation, createOperation_boundCRS_identified_by_datum) {
         NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest));
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=unitconvert "
-              "+xy_in=deg +xy_out=rad +step +proj=cart +ellps=WGS84 +step "
+              "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step "
               "+proj=helmert +x=263 +y=-6 +z=-431 +step +inv +proj=cart "
               "+ellps=clrk80ign +step +proj=pop +v_3 +step +proj=utm +zone=32 "
               "+ellps=clrk80ign");
@@ -5461,13 +5459,11 @@ TEST(operation, boundCRS_of_projCRS_to_geogCRS) {
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +inv +proj=utm +zone=31 +ellps=clrk80ign "
-              "+pm=paris +step +proj=longlat +ellps=clrk80ign +pm=paris +step "
-              "+proj=push +v_3 +step +inv +proj=longlat +ellps=clrk80ign "
-              "+pm=paris +step +proj=cart +ellps=clrk80ign +step +proj=helmert "
-              "+x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
-              "+convention=position_vector +step +inv +proj=cart +ellps=WGS84 "
-              "+step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
-              "+proj=axisswap +order=2,1 +step +proj=pop +v_3");
+              "+pm=paris +step +proj=push +v_3 +step +proj=cart "
+              "+ellps=clrk80ign +step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 "
+              "+rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
+              "+ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert "
+              "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
 }
 
 // ---------------------------------------------------------------------------
@@ -5483,13 +5479,13 @@ TEST(operation, boundCRS_of_geogCRS_to_projCRS) {
         CoordinateOperationFactory::create()->createOperation(boundCRS, utm31);
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=grad +xy_out=rad "
-              "+step +inv +proj=longlat +ellps=clrk80ign +pm=paris +step "
-              "+proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 +z=3 "
-              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
-              "+proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=utm "
-              "+zone=31 +ellps=WGS84");
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=grad +xy_out=rad +step +inv "
+              "+proj=longlat +ellps=clrk80ign +pm=paris +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 "
+              "+z=3 +rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step "
+              "+inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step "
+              "+proj=utm +zone=31 +ellps=WGS84");
 }
 
 // ---------------------------------------------------------------------------
@@ -5501,14 +5497,14 @@ TEST(operation, geogCRS_to_boundCRS_of_geogCRS) {
         GeographicCRS::EPSG_4326, boundCRS);
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-              "+proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=1 +y=2 +z=3 "
-              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
-              "+proj=cart +ellps=clrk80ign +step +proj=longlat "
-              "+ellps=clrk80ign +pm=paris +step +proj=unitconvert +xy_in=rad "
-              "+xy_out=grad +step +proj=axisswap +order=2,1 +step +proj=pop "
-              "+v_3");
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=1 "
+              "+y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector "
+              "+step +inv +proj=cart +ellps=clrk80ign +step +proj=pop +v_3 "
+              "+step +proj=longlat +ellps=clrk80ign +pm=paris +step "
+              "+proj=unitconvert +xy_in=rad +xy_out=grad +step +proj=axisswap "
+              "+order=2,1");
 }
 
 // ---------------------------------------------------------------------------
@@ -5531,14 +5527,12 @@ TEST(operation, boundCRS_to_boundCRS) {
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +inv +proj=utm +zone=31 +ellps=clrk80ign "
-              "+pm=paris +step +proj=longlat +ellps=clrk80ign +pm=paris +step "
-              "+proj=push +v_3 +step +inv +proj=longlat +ellps=clrk80ign "
-              "+pm=paris +step +proj=cart +ellps=clrk80ign +step +proj=helmert "
-              "+x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
-              "+convention=position_vector +step +inv +proj=helmert +x=8 +y=9 "
-              "+z=10 +rx=11 +ry=12 +rz=13 +s=14 +convention=position_vector "
-              "+step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step "
-              "+proj=utm +zone=32 +ellps=GRS80");
+              "+pm=paris +step +proj=push +v_3 +step +proj=cart "
+              "+ellps=clrk80ign +step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 "
+              "+rz=6 +s=7 +convention=position_vector +step +inv +proj=helmert "
+              "+x=8 +y=9 +z=10 +rx=11 +ry=12 +rz=13 +s=14 "
+              "+convention=position_vector +step +inv +proj=cart +ellps=GRS80 "
+              "+step +proj=pop +v_3 +step +proj=utm +zone=32 +ellps=GRS80");
 }
 
 // ---------------------------------------------------------------------------
@@ -5552,12 +5546,12 @@ TEST(operation, boundCRS_to_boundCRS_noop_for_TOWGS84) {
                                                                     boundCRS2);
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=grad +xy_out=rad "
-              "+step +inv +proj=longlat +ellps=clrk80ign +pm=paris +step "
-              "+proj=cart +ellps=clrk80ign +step +inv +proj=cart +ellps=GRS80 "
-              "+step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
-              "+proj=axisswap +order=2,1 +step +proj=pop +v_3");
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=grad +xy_out=rad +step +inv "
+              "+proj=longlat +ellps=clrk80ign +pm=paris +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=clrk80ign +step +inv +proj=cart "
+              "+ellps=GRS80 +step +proj=pop +v_3 +step +proj=unitconvert "
+              "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
 }
 
 // ---------------------------------------------------------------------------
@@ -5955,13 +5949,13 @@ TEST(operation, compoundCRS_with_boundGeogCRS_to_geogCRS) {
         compound, GeographicCRS::EPSG_4979);
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-              "+proj=cart +ellps=WGS84 +step +proj=helmert +x=1 +y=2 +z=3 "
-              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
-              "+proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
-              "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop "
-              "+v_3");
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=WGS84 +step +proj=helmert +x=1 +y=2 "
+              "+z=3 +rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step "
+              "+inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step "
+              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
+              "+order=2,1");
 }
 
 // ---------------------------------------------------------------------------
@@ -5979,12 +5973,12 @@ TEST(operation, compoundCRS_with_boundGeogCRS_and_boundVerticalCRS_to_geogCRS) {
     // Not completely sure the order of horizontal and vertical operations
     // makes sense
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
-              "+order=2,1 +step +proj=unitconvert +xy_in=grad +xy_out=rad "
-              "+step +inv +proj=longlat +ellps=clrk80ign +pm=paris +step "
-              "+proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 +z=3 "
-              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
-              "+proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step "
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=grad +xy_out=rad +step +inv "
+              "+proj=longlat +ellps=clrk80ign +pm=paris +step +proj=push +v_3 "
+              "+step +proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 "
+              "+z=3 +rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step "
+              "+inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step "
               "+proj=vgridshift +grids=egm08_25.gtx +multiplier=1 +step "
               "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
               "+order=2,1");
@@ -6018,14 +6012,12 @@ TEST(operation, compoundCRS_with_boundProjCRS_and_boundVerticalCRS_to_geogCRS) {
     // makes sense
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +inv +proj=utm +zone=31 +ellps=clrk80ign "
-              "+pm=paris +step +proj=longlat +ellps=clrk80ign +pm=paris +step "
-              "+proj=push +v_3 +step +inv +proj=longlat +ellps=clrk80ign "
-              "+pm=paris +step +proj=cart +ellps=clrk80ign +step +proj=helmert "
-              "+x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
-              "+convention=position_vector +step +inv +proj=cart +ellps=WGS84 "
-              "+step +proj=pop +v_3 +step +proj=vgridshift +grids=egm08_25.gtx "
-              "+multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg "
-              "+step +proj=axisswap +order=2,1");
+              "+pm=paris +step +proj=push +v_3 +step +proj=cart "
+              "+ellps=clrk80ign +step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 "
+              "+rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
+              "+ellps=WGS84 +step +proj=pop +v_3 +step +proj=vgridshift "
+              "+grids=egm08_25.gtx +multiplier=1 +step +proj=unitconvert "
+              "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
 
     auto opInverse = CoordinateOperationFactory::create()->createOperation(
         GeographicCRS::EPSG_4979, compound);
@@ -6125,6 +6117,63 @@ TEST(operation, compoundCRS_to_compoundCRS_with_vertical_transform) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, compoundCRS_to_compoundCRS_context) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        // NAD27 + NGVD29 height (ftUS)
+        authFactory->createCoordinateReferenceSystem("7406"),
+        // NAD83(NSRS2007) + NAVD88 height
+        authFactory->createCoordinateReferenceSystem("5500"), ctxt);
+    // 152 or 155 depending if the VERTCON grids are there
+    ASSERT_GE(list.size(), 152U);
+    EXPECT_FALSE(list[0]->hasBallparkTransformation());
+    EXPECT_EQ(list[0]->nameStr(), "NGVD29 height (ftUS) to NAVD88 height (3) + "
+                                  "NAD27 to WGS 84 (79) + Inverse of "
+                                  "NAD83(NSRS2007) to WGS 84 (1)");
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=rad +z_out=m "
+              "+step +proj=vgridshift +grids=vertcone.gtx +multiplier=0.001 "
+              "+step +proj=hgridshift +grids=conus +step "
+              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
+              "+order=2,1");
+
+    bool foundApprox = false;
+    for (size_t i = 0; i < list.size(); i++) {
+        auto projString =
+            list[i]->exportToPROJString(PROJStringFormatter::create().get());
+        EXPECT_TRUE(
+            projString.find("+proj=pipeline +step +proj=axisswap +order=2,1 "
+                            "+step +proj=unitconvert +xy_in=deg +z_in=us-ft "
+                            "+xy_out=rad +z_out=m") == 0)
+            << list[i]->nameStr();
+        if (list[i]->nameStr().find("Transformation from NGVD29 height (ftUS) "
+                                    "to NAVD88 height (ballpark vertical "
+                                    "transformation)") == 0) {
+            EXPECT_TRUE(list[i]->hasBallparkTransformation());
+            EXPECT_EQ(list[i]->nameStr(),
+                      "Transformation from NGVD29 height (ftUS) to NAVD88 "
+                      "height (ballpark vertical transformation) + NAD27 to "
+                      "WGS 84 (79) + Inverse of NAD83(NSRS2007) to WGS 84 (1)");
+            EXPECT_EQ(projString,
+                      "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+                      "+proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=rad "
+                      "+z_out=m +step +proj=hgridshift +grids=conus "
+                      "+step +proj=unitconvert +xy_in=rad "
+                      "+xy_out=deg +step +proj=axisswap +order=2,1");
+            foundApprox = true;
+            break;
+        }
+    }
+    EXPECT_TRUE(foundApprox);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, vertCRS_to_vertCRS) {
 
     auto vertcrs_m_obj = PROJStringParser().createFromPROJString("+vunits=m");
@@ -6174,6 +6223,25 @@ TEST(operation, vertCRS_to_vertCRS) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, vertCRS_to_vertCRS_context) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        // NGVD29 height (m)
+        authFactory->createCoordinateReferenceSystem("7968"),
+        // NAVD88 height (1)
+        authFactory->createCoordinateReferenceSystem("5703"), ctxt);
+    ASSERT_EQ(list.size(), 3U);
+    EXPECT_EQ(list[0]->nameStr(), "NGVD29 height (m) to NAVD88 height (3)");
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=vgridshift +grids=vertcone.gtx +multiplier=0.001");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, compoundCRS_to_geogCRS_3D) {
 
     auto compoundcrs_ft_obj = PROJStringParser().createFromPROJString(
@@ -6190,6 +6258,7 @@ TEST(operation, compoundCRS_to_geogCRS_3D) {
         auto op = CoordinateOperationFactory::create()->createOperation(
             NN_CHECK_ASSERT(compoundcrs_ft), NN_CHECK_ASSERT(geogcrs_m));
         ASSERT_TRUE(op != nullptr);
+        EXPECT_TRUE(op->hasBallparkTransformation());
         EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
                   "+proj=pipeline +step +inv +proj=merc +lon_0=0 +k=1 +x_0=0 "
                   "+y_0=0 +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
@@ -6200,6 +6269,7 @@ TEST(operation, compoundCRS_to_geogCRS_3D) {
         auto op = CoordinateOperationFactory::create()->createOperation(
             NN_CHECK_ASSERT(geogcrs_m), NN_CHECK_ASSERT(compoundcrs_ft));
         ASSERT_TRUE(op != nullptr);
+        EXPECT_TRUE(op->hasBallparkTransformation());
         EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
                   "+proj=pipeline +step +proj=unitconvert +xy_in=deg +z_in=m "
                   "+xy_out=rad +z_out=ft +step +proj=merc +lon_0=0 +k=1 +x_0=0 "
@@ -6209,6 +6279,64 @@ TEST(operation, compoundCRS_to_geogCRS_3D) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, compoundCRS_to_geogCRS_3D_context) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    // CompoundCRS to Geog3DCRS, with vertical unit change, but without
+    // ellipsoid height <--> vertical height correction
+    {
+        auto ctxt =
+            CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            authFactory->createCoordinateReferenceSystem(
+                "7406"), // NAD27 + NGVD29 height (ftUS)
+            authFactory->createCoordinateReferenceSystem("4979"), // WGS 84
+            ctxt);
+        ASSERT_GE(list.size(), 1U);
+        EXPECT_TRUE(list[0]->hasBallparkTransformation());
+        EXPECT_EQ(list[0]->nameStr(),
+                  "NAD27 to WGS 84 (79) + Transformation from NGVD29 height "
+                  "(ftUS) to WGS 84 (ballpark vertical transformation, without "
+                  "ellipsoid height to vertical height correction)");
+        EXPECT_EQ(list[0]->exportToPROJString(
+                      PROJStringFormatter::create(
+                          PROJStringFormatter::Convention::PROJ_5,
+                          authFactory->databaseContext())
+                          .get()),
+                  "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+                  "+proj=unitconvert +xy_in=deg +xy_out=rad +step "
+                  "+proj=hgridshift +grids=conus +step +proj=unitconvert "
+                  "+xy_in=rad +z_in=us-ft +xy_out=deg +z_out=m +step "
+                  "+proj=axisswap +order=2,1");
+    }
+
+    // CompoundCRS to Geog3DCRS, with same vertical unit, but without
+    // ellipsoid height <--> vertical height correction
+    {
+        auto ctxt =
+            CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            authFactory->createCoordinateReferenceSystem(
+                "5500"), // NAD83(NSRS2007) + NAVD88 height
+            authFactory->createCoordinateReferenceSystem("4979"), // WGS 84
+            ctxt);
+        ASSERT_GE(list.size(), 1U);
+        EXPECT_TRUE(list[0]->hasBallparkTransformation());
+        EXPECT_EQ(list[0]->nameStr(),
+                  "NAD83(NSRS2007) to WGS 84 (1) + Transformation from NAVD88 "
+                  "height to WGS 84 (ballpark vertical transformation, without "
+                  "ellipsoid height to vertical height correction)");
+        EXPECT_EQ(list[0]->exportToPROJString(
+                      PROJStringFormatter::create(
+                          PROJStringFormatter::Convention::PROJ_5,
+                          authFactory->databaseContext())
+                          .get()),
+                  "");
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, IGNF_LAMB1_TO_EPSG_4326) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), std::string());
@@ -6223,6 +6351,7 @@ TEST(operation, IGNF_LAMB1_TO_EPSG_4326) {
         ctxt);
     ASSERT_EQ(list.size(), 2U);
 
+    EXPECT_FALSE(list[0]->hasBallparkTransformation());
     EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +inv +proj=lcc +lat_1=49.5 +lat_0=49.5 "
               "+lon_0=0 +k_0=0.99987734 +x_0=600000 +y_0=200000 "
@@ -6230,14 +6359,15 @@ TEST(operation, IGNF_LAMB1_TO_EPSG_4326) {
               "+grids=ntf_r93.gsb +step +proj=unitconvert +xy_in=rad "
               "+xy_out=deg +step +proj=axisswap +order=2,1");
 
+    EXPECT_FALSE(list[1]->hasBallparkTransformation());
     EXPECT_EQ(list[1]->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +inv +proj=lcc +lat_1=49.5 +lat_0=49.5 "
               "+lon_0=0 +k_0=0.99987734 +x_0=600000 +y_0=200000 "
               "+ellps=clrk80ign +pm=paris +step +proj=push +v_3 +step "
               "+proj=cart +ellps=clrk80ign +step +proj=helmert +x=-168 +y=-60 "
-              "+z=320 +step +inv +proj=cart +ellps=WGS84 +step "
-              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
-              "+order=2,1 +step +proj=pop +v_3");
+              "+z=320 +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
+              "+proj=axisswap +order=2,1");
 
     auto list2 = CoordinateOperationFactory::create()->createOperations(
         AuthorityFactory::create(DatabaseContext::create(), "EPSG")
@@ -6264,14 +6394,14 @@ TEST(operation, IGNF_LAMB1_TO_EPSG_4326) {
 
 // ---------------------------------------------------------------------------
 
-TEST(operation, isPROJInstanciable) {
+TEST(operation, isPROJInstantiable) {
 
     {
         auto transformation = Transformation::createGeocentricTranslations(
             PropertyMap(), GeographicCRS::EPSG_4269, GeographicCRS::EPSG_4326,
             1.0, 2.0, 3.0, {});
         EXPECT_TRUE(
-            transformation->isPROJInstanciable(DatabaseContext::create()));
+            transformation->isPROJInstantiable(DatabaseContext::create()));
     }
 
     // Missing grid
@@ -6280,7 +6410,7 @@ TEST(operation, isPROJInstanciable) {
             PropertyMap(), GeographicCRS::EPSG_4807, GeographicCRS::EPSG_4326,
             "foo.gsb", std::vector<PositionalAccuracyNNPtr>());
         EXPECT_FALSE(
-            transformation->isPROJInstanciable(DatabaseContext::create()));
+            transformation->isPROJInstantiable(DatabaseContext::create()));
     }
 
     // Unsupported method
@@ -6292,7 +6422,7 @@ TEST(operation, isPROJInstanciable) {
             std::vector<GeneralParameterValueNNPtr>{},
             std::vector<PositionalAccuracyNNPtr>{});
         EXPECT_FALSE(
-            transformation->isPROJInstanciable(DatabaseContext::create()));
+            transformation->isPROJInstantiable(DatabaseContext::create()));
     }
 }