data/Makefile.am

@@ -64,7 +64,7 @@ proj.db: $(DATAPATH)/sql/*.sql
 	 for x in $(SQL_ORDERED_LIST); do \
 		export SQL_EXPANDED_LIST="$${SQL_EXPANDED_LIST} $(DATAPATH)/$$x"; \
 	 done; \
-	 if test "x$(PROJ_DB_CACHE_DIR)" != "x" -a -f "$(PROJ_DB_CACHE_DIR)/proj.db" -a -f "$(PROJ_DB_CACHE_DIR)/proj.db.sql.md5" ; then \
+	 if test "x$(PROJ_DB_CACHE_DIR)" != "x" -a -x "$$(command -v md5sum)" -a -f "$(PROJ_DB_CACHE_DIR)/proj.db" -a -f "$(PROJ_DB_CACHE_DIR)/proj.db.sql.md5" ; then \
 	    cat $${SQL_EXPANDED_LIST} | md5sum | diff - "$(PROJ_DB_CACHE_DIR)/proj.db.sql.md5" > /dev/null \
 	    && (echo "Reusing cached proj.db"; cp "$(PROJ_DB_CACHE_DIR)/proj.db" proj.db); \
 	 fi; \
@@ -84,7 +84,7 @@ proj.db: $(DATAPATH)/sql/*.sql
 		$(RM) proj.db; \
 		exit 1; \
 	 else \
-	   if test "x$(PROJ_DB_CACHE_DIR)" != "x" ; then \
+	   if test "x$(PROJ_DB_CACHE_DIR)" != "x" -a -x "$$(command -v md5sum)" ; then \
 		mkdir -p "$(PROJ_DB_CACHE_DIR)"; \
 		cat $${SQL_EXPANDED_LIST} | md5sum > "$(PROJ_DB_CACHE_DIR)/proj.db.sql.md5"; \
 		cp proj.db "$(PROJ_DB_CACHE_DIR)"; \

data/Makefile.in

@@ -576,7 +576,7 @@ proj.db: $(DATAPATH)/sql/*.sql
 	 for x in $(SQL_ORDERED_LIST); do \
 		export SQL_EXPANDED_LIST="$${SQL_EXPANDED_LIST} $(DATAPATH)/$$x"; \
 	 done; \
-	 if test "x$(PROJ_DB_CACHE_DIR)" != "x" -a -f "$(PROJ_DB_CACHE_DIR)/proj.db" -a -f "$(PROJ_DB_CACHE_DIR)/proj.db.sql.md5" ; then \
+	 if test "x$(PROJ_DB_CACHE_DIR)" != "x" -a -x "$$(command -v md5sum)" -a -f "$(PROJ_DB_CACHE_DIR)/proj.db" -a -f "$(PROJ_DB_CACHE_DIR)/proj.db.sql.md5" ; then \
 	    cat $${SQL_EXPANDED_LIST} | md5sum | diff - "$(PROJ_DB_CACHE_DIR)/proj.db.sql.md5" > /dev/null \
 	    && (echo "Reusing cached proj.db"; cp "$(PROJ_DB_CACHE_DIR)/proj.db" proj.db); \
 	 fi; \
@@ -596,7 +596,7 @@ proj.db: $(DATAPATH)/sql/*.sql
 		$(RM) proj.db; \
 		exit 1; \
 	 else \
-	   if test "x$(PROJ_DB_CACHE_DIR)" != "x" ; then \
+	   if test "x$(PROJ_DB_CACHE_DIR)" != "x" -a -x "$$(command -v md5sum)" ; then \
 		mkdir -p "$(PROJ_DB_CACHE_DIR)"; \
 		cat $${SQL_EXPANDED_LIST} | md5sum > "$(PROJ_DB_CACHE_DIR)/proj.db.sql.md5"; \
 		cp proj.db "$(PROJ_DB_CACHE_DIR)"; \

data/sql/customizations.sql

@@ -19,6 +19,25 @@ DELETE FROM "supersession" WHERE superseded_table_name = 'grid_transformation' A
                                  replacement_code = '8885' AND
                                  source = 'EPSG';
 
+-- ('EPSG','7001','ETRS89 to NAP height (1)') lacks an interpolationCRS with Amersfoort / EPSG:4289
+-- See https://salsa.debian.org/debian-gis-team/proj-rdnap/blob/debian/2008-8/Use%20of%20RDTRANS2008%20and%20NAPTRANS2008.pdf
+-- "The naptrans2008 VDatum-grid is referenced to the Bessel-1841 ellipsoid"
+CREATE TABLE dummy(foo);
+CREATE TRIGGER check_grid_transformation_epsg_7001
+BEFORE INSERT ON dummy
+FOR EACH ROW BEGIN
+    SELECT RAISE(ABORT, 'grid_transformation EPSG:7001 entry is not ETRS89 to NAP height (1)')
+        WHERE NOT EXISTS(SELECT 1 FROM grid_transformation WHERE auth_name = 'EPSG' AND code = '7001' AND name = 'ETRS89 to NAP height (1)');
+    SELECT RAISE(ABORT, 'grid_transformation EPSG:7001 entry has already an interpolationCRS')
+        WHERE EXISTS(SELECT 1 FROM grid_transformation WHERE auth_name = 'EPSG' AND code = '7001' AND interpolation_crs_auth_name IS NOT NULL);
+END;
+INSERT INTO dummy DEFAULT VALUES;
+DROP TRIGGER check_grid_transformation_epsg_7001;
+DROP TABLE dummy;
+UPDATE grid_transformation SET interpolation_crs_auth_name = 'EPSG',
+                               interpolation_crs_code = '4289'
+                           WHERE auth_name = 'EPSG' AND code = '7001';
+
 -- Define the allowed authorities, and their precedence, when researching a
 -- coordinate operation
 

data/sql/grid_alternatives.sql

@@ -798,3 +798,41 @@ INSERT INTO grid_alternatives(original_grid_name,
                               0,
                               'proj-datumgrid-oceania',
                               NULL, NULL, NULL, NULL);
+
+-- Netherlands / RDNAP (non-free grids)
+
+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 ('naptrans2008.gtx',
+                              'naptrans2008.gtx',
+                              'GTX',
+                              'vgridshift',
+                              0,
+                              NULL, -- package name
+                              'https://salsa.debian.org/debian-gis-team/proj-rdnap/raw/upstream/2008/naptrans2008.gtx',
+                              1, -- direct download
+                              0, -- non-freely licensed. See https://salsa.debian.org/debian-gis-team/proj-rdnap/raw/master/debian/copyright
+                              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 ('rdtrans2008.gsb',
+                              'rdtrans2008.gsb',
+                              'NTv2',
+                              'hgridshift',
+                              0,
+                              NULL, -- package name
+                              'https://salsa.debian.org/debian-gis-team/proj-rdnap/raw/upstream/2008/rdtrans2008.gsb',
+                              1, -- direct download
+                              0, -- non-freely licensed. See https://salsa.debian.org/debian-gis-team/proj-rdnap/raw/master/debian/copyright
+                              NULL);

debian/changelog

+proj (6.1.0~rc2-1~exp1) experimental; urgency=medium
+
+  * New upstream release candidate.
+  * Update 6.1.0~rc1 symbols for other architectures.
+  * Drop patches, applied/included upstream.
+
+ -- Bas Couwenberg <sebastic@debian.org>  Wed, 08 May 2019 06:20:12 +0200
+
 proj (6.1.0~rc1-1~exp2) experimental; urgency=medium
 
   * Add patch from PR #1454 to fix proj-rdnap test failures.

debian/libproj15.symbols

-# SymbolsHelper-Confirmed: 6.1.0~rc1 amd64 armel armhf ia64 mips mipsel powerpc powerpcspe
+# SymbolsHelper-Confirmed: 6.1.0~rc1 amd64 arm64 armel armhf hppa hurd-i386 i386 ia64 kfreebsd-amd64 kfreebsd-i386 m68k mips mips64el mipsel powerpc powerpcspe ppc64 ppc64el riscv64 s390x sh4 sparc64 x32
 libproj.so.15 #PACKAGE# #MINVER#
 * Build-Depends-Package: libproj-dev
  _Z10pj_ell_setP9projCtx_tP8ARG_listPdS3_@Base 6.0.0
@@ -540,7 +540,7 @@ libproj.so.15 #PACKAGE# #MINVER#
  _ZN5osgeo4proj6common9DataEpochD2Ev@Base 6.0.0
  _ZN5osgeo4proj8internal10replaceAllERKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEES9_S9_@Base 6.0.0
  _ZN5osgeo4proj8internal13c_locale_stodERKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEE@Base 6.0.0
- (arch=amd64)_ZN5osgeo4proj8internal5splitERKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEES9_@Base 6.1.0~rc1
+ _ZN5osgeo4proj8internal5splitERKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEES9_@Base 6.1.0~rc1
  _ZN5osgeo4proj8internal5splitERKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEEc@Base 6.0.0
  _ZN5osgeo4proj8internal7ci_findERKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEEPKc@Base 6.0.0
  _ZN5osgeo4proj8internal7tolowerERKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEE@Base 6.0.0
@@ -696,7 +696,7 @@ libproj.so.15 #PACKAGE# #MINVER#
  _ZN5osgeo4proj9operation14Transformation33createTimeDependentPositionVectorERKNS0_4util11PropertyMapERKN7dropbox6oxygen2nnISt10shared_ptrINS0_3crs3CRSEEEESG_dddddddddddddddRKSt6vectorINS9_ISA_INS0_8metadata18PositionalAccuracyEEEESaISL_EE@Base 6.0.0
  _ZN5osgeo4proj9operation14Transformation35createGeographic2DWithHeightOffsetsERKNS0_4util11PropertyMapERKN7dropbox6oxygen2nnISt10shared_ptrINS0_3crs3CRSEEEESG_RKNS0_6common5AngleESK_RKNSH_6LengthERKSt6vectorINS9_ISA_INS0_8metadata18PositionalAccuracyEEEESaISS_EE@Base 6.0.0
  (arch=!amd64)_ZN5osgeo4proj9operation14Transformation40createGravityRelatedHeightToGeographic3DERKNS0_4util11PropertyMapERKN7dropbox6oxygen2nnISt10shared_ptrINS0_3crs3CRSEEEESG_RKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEERKSt6vectorINS9_ISA_INS0_8metadata18PositionalAccuracyEEEESaIST_EE@Base 6.1.0~rc1
- (arch=amd64)_ZN5osgeo4proj9operation14Transformation40createGravityRelatedHeightToGeographic3DERKNS0_4util11PropertyMapERKN7dropbox6oxygen2nnISt10shared_ptrINS0_3crs3CRSEEEESG_RKSD_RKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEERKSt6vectorINS9_ISA_INS0_8metadata18PositionalAccuracyEEEESaISV_EE@Base 6.1.0~rc1
+ _ZN5osgeo4proj9operation14Transformation40createGravityRelatedHeightToGeographic3DERKNS0_4util11PropertyMapERKN7dropbox6oxygen2nnISt10shared_ptrINS0_3crs3CRSEEEESG_RKSD_RKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEERKSt6vectorINS9_ISA_INS0_8metadata18PositionalAccuracyEEEESaISV_EE@Base 6.1.0~rc1
  _ZN5osgeo4proj9operation14Transformation42createTimeDependentCoordinateFrameRotationERKNS0_4util11PropertyMapERKN7dropbox6oxygen2nnISt10shared_ptrINS0_3crs3CRSEEEESG_dddddddddddddddRKSt6vectorINS9_ISA_INS0_8metadata18PositionalAccuracyEEEESaISL_EE@Base 6.0.0
  _ZN5osgeo4proj9operation14Transformation6createERKNS0_4util11PropertyMapERKN7dropbox6oxygen2nnISt10shared_ptrINS0_3crs3CRSEEEESG_RKSD_RKNS9_ISA_INS1_15OperationMethodEEEERKSt6vectorINS9_ISA_INS1_21GeneralParameterValueEEEESaISR_EERKSO_INS9_ISA_INS0_8metadata18PositionalAccuracyEEEESaISZ_EE@Base 6.0.0
  _ZN5osgeo4proj9operation14Transformation6createERKNS0_4util11PropertyMapERKN7dropbox6oxygen2nnISt10shared_ptrINS0_3crs3CRSEEEESG_RKSD_S6_RKSt6vectorINS9_ISA_INS1_18OperationParameterEEEESaISM_EERKSJ_INS9_ISA_INS1_14ParameterValueEEEESaIST_EERKSJ_INS9_ISA_INS0_8metadata18PositionalAccuracyEEEESaIS11_EE@Base 6.0.0
@@ -1077,8 +1077,8 @@ libproj.so.15 #PACKAGE# #MINVER#
  (optional=templinst)_ZNSt12__shared_ptrIN5osgeo4proj5datum22GeodeticReferenceFrameELN9__gnu_cxx12_Lock_policyE2EEC2ERKS6_@Base 6.1.0~rc1
  (optional=templinst|arch=armel)_ZNSt12__shared_ptrIN5osgeo4proj6common16IdentifiedObjectELN9__gnu_cxx12_Lock_policyE1EEC1INS1_9operation15OperationMethodEvEERKS_IT_LS5_1EE@Base 6.0.0
  (optional=templinst|arch=armel)_ZNSt12__shared_ptrIN5osgeo4proj6common16IdentifiedObjectELN9__gnu_cxx12_Lock_policyE1EEC2INS1_9operation15OperationMethodEvEERKS_IT_LS5_1EE@Base 6.0.0
- (optional=templinst|arch=!armel)_ZNSt12__shared_ptrIN5osgeo4proj6common16IdentifiedObjectELN9__gnu_cxx12_Lock_policyE2EEC1INS1_9operation15OperationMethodEvEERKS_IT_LS5_2EE@Base 6.0.0
- (optional=templinst|arch=!armel)_ZNSt12__shared_ptrIN5osgeo4proj6common16IdentifiedObjectELN9__gnu_cxx12_Lock_policyE2EEC2INS1_9operation15OperationMethodEvEERKS_IT_LS5_2EE@Base 6.0.0
+ (optional=templinst|arch=!armel !riscv64)_ZNSt12__shared_ptrIN5osgeo4proj6common16IdentifiedObjectELN9__gnu_cxx12_Lock_policyE2EEC1INS1_9operation15OperationMethodEvEERKS_IT_LS5_2EE@Base 6.0.0
+ (optional=templinst|arch=!armel !riscv64)_ZNSt12__shared_ptrIN5osgeo4proj6common16IdentifiedObjectELN9__gnu_cxx12_Lock_policyE2EEC2INS1_9operation15OperationMethodEvEERKS_IT_LS5_2EE@Base 6.0.0
  (optional=templinst|arch=armel)_ZNSt12__shared_ptrIN5osgeo4proj8metadata6ExtentELN9__gnu_cxx12_Lock_policyE1EEC1ERKS6_@Base 6.0.0
  (optional=templinst|arch=armel)_ZNSt12__shared_ptrIN5osgeo4proj8metadata6ExtentELN9__gnu_cxx12_Lock_policyE1EEC2ERKS6_@Base 6.0.0
  (optional=templinst|arch=!armel !riscv64)_ZNSt12__shared_ptrIN5osgeo4proj8metadata6ExtentELN9__gnu_cxx12_Lock_policyE2EEC1ERKS6_@Base 6.0.0

debian/patches/series

-spelling-errors.patch
-pr1454-createOperations-fix-case-of-ETRS89-3D-to-proj-string-with-nadgrids-and-geoidgrids.patch

debian/patches/spelling-errors.patch

-Description: Fix spelling errors.
- * unknow -> unknown
-Author: Bas Couwenberg <sebastic@debian.org>
-Forwarded: https://github.com/OSGeo/proj.4/pull/1452
-Applied-Upstream: https://github.com/OSGeo/proj.4/commit/a16cb05010115690ecca8938d95ca544f2ad65c6
-
---- a/src/strerrno.cpp
-+++ b/src/strerrno.cpp
-@@ -70,7 +70,7 @@ pj_err_list[] = {
-     "argument not numerical or out of range",                          /* -58 */
-     "inconsistent unit type between input and output",                 /* -59 */
-     "arguments are mutually exclusive",                                /* -60 */
--    "generic error of unknow origin",                                  /* -61 */
-+    "generic error of unknown origin",                                 /* -61 */
- 
-     /* When adding error messages, remember to update ID defines in
-        projects.h, and transient_error array in pj_transform                  */

include/proj/coordinateoperation.hpp

@@ -1476,7 +1476,8 @@ class PROJ_GCC_DLL Transformation : public SingleOperation {
     PROJ_DLL static TransformationNNPtr
     createGravityRelatedHeightToGeographic3D(
         const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
-        const crs::CRSNNPtr &targetCRSIn, const std::string &filename,
+        const crs::CRSNNPtr &targetCRSIn, const crs::CRSPtr &interpolationCRSIn,
+        const std::string &filename,
         const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies);
 
     PROJ_DLL static TransformationNNPtr createVERTCON(

include/proj/internal/internal.hpp

@@ -140,6 +140,9 @@ std::string toupper(const std::string &osStr);
 PROJ_FOR_TEST std::vector<std::string> split(const std::string &osStr,
                                              char separator);
 
+PROJ_FOR_TEST std::vector<std::string> split(const std::string &osStr,
+                                             const std::string &separator);
+
 bool ci_equal(const char *a, const char *b) noexcept;
 
 #ifdef SUPPORT_DELETED_FUNCTION

include/proj/io.hpp

@@ -427,7 +427,8 @@ class PROJ_GCC_DLL PROJStringFormatter {
     PROJ_INTERNAL void setOmitProjLongLatIfPossible(bool omit);
     PROJ_INTERNAL bool omitProjLongLatIfPossible() const;
 
-    PROJ_INTERNAL void setOmitZUnitConversion(bool omit);
+    PROJ_INTERNAL void pushOmitZUnitConversion();
+    PROJ_INTERNAL void popOmitZUnitConversion();
     PROJ_INTERNAL bool omitZUnitConversion() const;
 
     PROJ_INTERNAL void setDropEarlyBindingsTerms(bool drop);

src/iso19111/coordinateoperation.cpp

@@ -6958,11 +6958,11 @@ TransformationNNPtr Transformation::createNTv2(
 static TransformationNNPtr _createGravityRelatedHeightToGeographic3D(
     const util::PropertyMap &properties, bool inverse,
     const crs::CRSNNPtr &sourceCRSIn, const crs::CRSNNPtr &targetCRSIn,
-    const std::string &filename,
+    const crs::CRSPtr &interpolationCRSIn, const std::string &filename,
     const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies) {
     return Transformation::create(
-        properties, sourceCRSIn, targetCRSIn, nullptr,
+        properties, sourceCRSIn, targetCRSIn, interpolationCRSIn,
         util::PropertyMap().set(
             common::IdentifiedObject::NAME_KEY,
             inverse ? INVERSE_OF + PROJ_WKT2_NAME_METHOD_HEIGHT_TO_GEOG3D
@@ -6981,17 +6981,20 @@ static TransformationNNPtr _createGravityRelatedHeightToGeographic3D(
  * At minimum the name should be defined.
  * @param sourceCRSIn Source CRS.
  * @param targetCRSIn Target CRS.
+ * @param interpolationCRSIn Interpolation CRS. (might be null)
  * @param filename GRID filename.
  * @param accuracies Vector of positional accuracy (might be empty).
  * @return new Transformation.
  */
 TransformationNNPtr Transformation::createGravityRelatedHeightToGeographic3D(
     const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
-    const crs::CRSNNPtr &targetCRSIn, const std::string &filename,
+    const crs::CRSNNPtr &targetCRSIn, const crs::CRSPtr &interpolationCRSIn,
+    const std::string &filename,
     const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies) {
     return _createGravityRelatedHeightToGeographic3D(
-        properties, false, sourceCRSIn, targetCRSIn, filename, accuracies);
+        properties, false, sourceCRSIn, targetCRSIn, interpolationCRSIn,
+        filename, accuracies);
 }
 // ---------------------------------------------------------------------------
@@ -7302,8 +7305,20 @@ createPropertiesForInverse(const CoordinateOperation *op, bool derivedFrom,
     auto targetCRS = op->targetCRS();
     std::string name;
     if (!forwardName.empty()) {
-        if (starts_with(forwardName, INVERSE_OF)) {
-            name = forwardName.substr(INVERSE_OF.size());
+        if (starts_with(forwardName, INVERSE_OF) ||
+            forwardName.find(" + ") != std::string::npos) {
+            auto tokens = split(forwardName, " + ");
+            for (size_t i = tokens.size(); i > 0;) {
+                i--;
+                if (!name.empty()) {
+                    name += " + ";
+                }
+                if (starts_with(tokens[i], INVERSE_OF)) {
+                    name += tokens[i].substr(INVERSE_OF.size());
+                } else {
+                    name += INVERSE_OF + tokens[i];
+                }
+            }
         } else if (!sourceCRS || !targetCRS ||
                    forwardName != buildOpName(opType, sourceCRS, targetCRS)) {
             name = INVERSE_OF + forwardName;
@@ -8195,13 +8210,14 @@ TransformationNNPtr Transformation::substitutePROJAlternativeGridNames(
                 return createGravityRelatedHeightToGeographic3D(
                            createPropertiesForInverse(self.as_nullable().get(),
                                                       true, false),
-                           targetCRS(), sourceCRS(), projFilename,
-                           coordinateOperationAccuracies())
+                           targetCRS(), sourceCRS(), interpolationCRS(),
+                           projFilename, coordinateOperationAccuracies())
                     ->inverseAsTransformation();
             } else {
                 return createGravityRelatedHeightToGeographic3D(
                     createSimilarPropertiesTransformation(self), sourceCRS(),
-                    targetCRS(), projFilename, coordinateOperationAccuracies());
+                    targetCRS(), interpolationCRS(), projFilename,
+                    coordinateOperationAccuracies());
             }
         }
     }
@@ -10972,19 +10988,19 @@ struct MyPROJStringExportableHorizVertical final
         // cppcheck-suppress functionStatic
         _exportToPROJString(io::PROJStringFormatter *formatter) const override {
-        formatter->setOmitZUnitConversion(true);
+        formatter->pushOmitZUnitConversion();
         horizTransform->_exportToPROJString(formatter);
         formatter->startInversion();
         geogDst->addAngularUnitConvertAndAxisSwap(formatter);
         formatter->stopInversion();
-        formatter->setOmitZUnitConversion(false);
+        formatter->popOmitZUnitConversion();
         verticalTransform->_exportToPROJString(formatter);
-        formatter->setOmitZUnitConversion(true);
+        formatter->pushOmitZUnitConversion();
         geogDst->addAngularUnitConvertAndAxisSwap(formatter);
-        formatter->setOmitZUnitConversion(false);
+        formatter->popOmitZUnitConversion();
     }
 };
@@ -11016,7 +11032,7 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final
         // cppcheck-suppress functionStatic
         _exportToPROJString(io::PROJStringFormatter *formatter) const override {
-        formatter->setOmitZUnitConversion(true);
+        formatter->pushOmitZUnitConversion();
         opSrcCRSToGeogCRS->_exportToPROJString(formatter);
@@ -11024,17 +11040,17 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final
         interpolationGeogCRS->addAngularUnitConvertAndAxisSwap(formatter);
         formatter->stopInversion();
-        formatter->setOmitZUnitConversion(false);
+        formatter->popOmitZUnitConversion();
         verticalTransform->_exportToPROJString(formatter);
-        formatter->setOmitZUnitConversion(true);
+        formatter->pushOmitZUnitConversion();
         interpolationGeogCRS->addAngularUnitConvertAndAxisSwap(formatter);
         opGeogCRStoDstCRS->_exportToPROJString(formatter);
-        formatter->setOmitZUnitConversion(false);
+        formatter->popOmitZUnitConversion();
     }
 };
@@ -12135,6 +12151,37 @@ CoordinateOperationFactory::Private::createOperations(
             }
         }
+        auto vertCRSOfBaseOfBoundSrc =
+            dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get());
+        if (vertCRSOfBaseOfBoundSrc && hubSrcGeog &&
+            hubSrcGeog->coordinateSystem()->axisList().size() == 3 &&
+            geogDst->coordinateSystem()->axisList().size() == 3) {
+            auto opsFirst = createOperations(sourceCRS, hubSrc, context);
+            auto opsSecond = createOperations(hubSrc, targetCRS, context);
+            if (!opsFirst.empty() && !opsSecond.empty()) {
+                for (const auto &opFirst : opsFirst) {
+                    for (const auto &opLast : opsSecond) {
+                        // Exclude artificial transformations from the hub
+                        // to the target CRS
+                        if (!opLast->hasBallparkTransformation()) {
+                            try {
+                                res.emplace_back(
+                                    ConcatenatedOperation::
+                                        createComputeMetadata(
+                                            {opFirst, opLast},
+                                            !allowEmptyIntersection));
+                            } catch (
+                                const InvalidOperationEmptyIntersection &) {
+                            }
+                        }
+                    }
+                }
+                if (!res.empty()) {
+                    return res;
+                }
+            }
+        }
+
         return createOperations(boundSrc->baseCRS(), targetCRS, context);
     }
@@ -12349,7 +12396,8 @@ CoordinateOperationFactory::Private::createOperations(
         const auto &componentsSrc = compoundSrc->componentReferenceSystems();
         if (!componentsSrc.empty()) {
             std::vector<CoordinateOperationNNPtr> horizTransforms;
-            if (componentsSrc[0]->extractGeographicCRS()) {
+            auto srcGeogCRS = componentsSrc[0]->extractGeographicCRS();
+            if (srcGeogCRS) {
                 horizTransforms =
                     createOperations(componentsSrc[0], targetCRS, context);
             }
@@ -12363,6 +12411,55 @@ CoordinateOperationFactory::Private::createOperations(
                 for (const auto &horizTransform : horizTransforms) {
                     for (const auto &verticalTransform : verticalTransforms) {
+                        crs::GeographicCRSPtr interpolationGeogCRS;
+                        auto transformationVerticalTransform =
+                            dynamic_cast<const Transformation *>(
+                                verticalTransform.get());
+                        if (transformationVerticalTransform) {
+                            auto interpTransformCRS =
+                                transformationVerticalTransform
+                                    ->interpolationCRS();
+                            if (interpTransformCRS) {
+                                auto nn_interpTransformCRS =
+                                    NN_NO_CHECK(interpTransformCRS);
+                                if (dynamic_cast<const crs::GeographicCRS *>(
+                                        nn_interpTransformCRS.get())) {
+                                    interpolationGeogCRS =
+                                        util::nn_dynamic_pointer_cast<
+                                            crs::GeographicCRS>(
+                                            nn_interpTransformCRS);
+                                }
+                            }
+                        }
+                        bool done = false;
+                        if (interpolationGeogCRS &&
+                            (interpolationGeogCRS->_isEquivalentTo(
+                                 srcGeogCRS.get(),
+                                 util::IComparable::Criterion::EQUIVALENT) ||
+                             interpolationGeogCRS->is2DPartOf3D(
+                                 NN_NO_CHECK(srcGeogCRS.get())))) {
+                            auto srcToInterp = createOperations(
+                                componentsSrc[0],
+                                NN_NO_CHECK(interpolationGeogCRS), context);
+                            auto interpToCompoundHoriz = createOperations(
+                                NN_NO_CHECK(interpolationGeogCRS),
+                                componentsSrc[0], context);
+                            if (!srcToInterp.empty() &&
+                                !interpToCompoundHoriz.empty()) {
+                                auto op = createHorizVerticalHorizPROJBased(
+                                    sourceCRS, componentsSrc[0],
+                                    srcToInterp.front(), verticalTransform,
+                                    interpToCompoundHoriz.front(),
+                                    interpolationGeogCRS);
+                                done = true;
+                                res.emplace_back(
+                                    ConcatenatedOperation::
+                                        createComputeMetadata(
+                                            {op, horizTransform},
+                                            !allowEmptyIntersection));
+                            }
+                        }
+                        if (!done) {
                             auto op = createHorizVerticalPROJBased(
                                 sourceCRS, targetCRS, horizTransform,
                                 verticalTransform);
@@ -12370,6 +12467,7 @@ CoordinateOperationFactory::Private::createOperations(
                             res.emplace_back(op);
                         }
                     }
+                }
                 return res;
             } else {
                 return horizTransforms;

src/iso19111/internal.cpp

@@ -298,6 +298,21 @@ std::vector<std::string> split(const std::string &str, char separator) {
 
 // ---------------------------------------------------------------------------
 
+std::vector<std::string> split(const std::string &str,
+                               const std::string &separator) {
+    std::vector<std::string> res;
+    size_t lastPos = 0;
+    size_t newPos = 0;
+    while ((newPos = str.find(separator, lastPos)) != std::string::npos) {
+        res.push_back(str.substr(lastPos, newPos - lastPos));
+        lastPos = newPos + separator.size();
+    }
+    res.push_back(str.substr(lastPos));
+    return res;
+}
+
+// ---------------------------------------------------------------------------
+
 #ifdef _WIN32
 
 // For some reason, sqlite3_snprintf() in the sqlite3 builds used on AppVeyor

src/iso19111/io.cpp

@@ -3818,7 +3818,7 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {
                     Transformation::createGravityRelatedHeightToGeographic3D(
                         PropertyMap().set(IdentifiedObject::NAME_KEY,
                                           transformationName),
-                        crs, GeographicCRS::EPSG_4979,
+                        crs, GeographicCRS::EPSG_4979, nullptr,
                         stripQuotes(extensionChildren[1]),
                         std::vector<PositionalAccuracyNNPtr>());
                 return nn_static_pointer_cast<CRS>(BoundCRS::create(
@@ -4922,7 +4922,7 @@ struct PROJStringFormatter::Private {
     };
     std::vector<InversionStackElt> inversionStack_{InversionStackElt()};
     bool omitProjLongLatIfPossible_ = false;
-    bool omitZUnitConversion_ = false;
+    std::vector<bool> omitZUnitConversion_{false};
     DatabaseContextPtr dbContext_{};
     bool useApproxTMerc_ = false;
     bool addNoDefs_ = true;
@@ -5939,15 +5939,21 @@ bool PROJStringFormatter::omitProjLongLatIfPossible() const {
 
 // ---------------------------------------------------------------------------
 
-void PROJStringFormatter::setOmitZUnitConversion(bool omit) {
-    assert(d->omitZUnitConversion_ ^ omit);
-    d->omitZUnitConversion_ = omit;
+void PROJStringFormatter::pushOmitZUnitConversion() {
+    d->omitZUnitConversion_.push_back(true);
+}
+
+// ---------------------------------------------------------------------------
+
+void PROJStringFormatter::popOmitZUnitConversion() {
+    assert(d->omitZUnitConversion_.size() > 1);
+    d->omitZUnitConversion_.pop_back();
 }
 
 // ---------------------------------------------------------------------------
 
 bool PROJStringFormatter::omitZUnitConversion() const {
-    return d->omitZUnitConversion_;
+    return d->omitZUnitConversion_.back();
 }
 
 // ---------------------------------------------------------------------------
@@ -6995,7 +7001,7 @@ PROJStringParser::Private::buildBoundOrCompoundCRSIfNeeded(int iStep,
             Transformation::createGravityRelatedHeightToGeographic3D(
                 PropertyMap().set(IdentifiedObject::NAME_KEY,
                                   "unknown to WGS84 ellipsoidal height"),
-                crs, GeographicCRS::EPSG_4979, geoidgrids,
+                crs, GeographicCRS::EPSG_4979, nullptr, geoidgrids,
                 std::vector<PositionalAccuracyNNPtr>());
         auto boundvcrs =
             BoundCRS::create(vcrs, GeographicCRS::EPSG_4979, transformation);

src/strerrno.cpp

@@ -70,7 +70,7 @@ pj_err_list[] = {
     "argument not numerical or out of range",                          /* -58 */
     "inconsistent unit type between input and output",                 /* -59 */
     "arguments are mutually exclusive",                                /* -60 */
-    "generic error of unknow origin",                                  /* -61 */
+    "generic error of unknown origin",                                 /* -61 */
 
     /* When adding error messages, remember to update ID defines in
        projects.h, and transient_error array in pj_transform                  */

test/unit/test_operation.cpp

@@ -5903,6 +5903,35 @@ TEST(operation, boundCRS_with_basecrs_with_extent_to_geogCRS) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, ETRS89_3D_to_proj_string_with_geoidgrids_nadgrids) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    // ETRS89 3D
+    auto src = authFactory->createCoordinateReferenceSystem("4937");
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 "
+        "+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel "
+        "+nadgrids=rdtrans2008.gsb +geoidgrids=naptrans2008.gtx +units=m "
+        "+type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        src, NN_NO_CHECK(dst), ctxt);
+    ASSERT_EQ(list.size(), 1U);
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +inv +proj=vgridshift +grids=naptrans2008.gtx "
+              "+multiplier=1 "
+              "+step +inv +proj=hgridshift +grids=rdtrans2008.gsb "
+              "+step +proj=sterea +lat_0=52.1561605555556 "
+              "+lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 "
+              "+y_0=463000 +ellps=bessel");
+}
+
+// ---------------------------------------------------------------------------
+
 static VerticalCRSNNPtr createVerticalCRS() {
     PropertyMap propertiesVDatum;
     propertiesVDatum.set(Identifier::CODESPACE_KEY, "EPSG")
@@ -5941,8 +5970,8 @@ static BoundCRSNNPtr createBoundVerticalCRS() {
     auto vertCRS = createVerticalCRS();
     auto transformation =
         Transformation::createGravityRelatedHeightToGeographic3D(
-            PropertyMap(), vertCRS, GeographicCRS::EPSG_4979, "egm08_25.gtx",
-            std::vector<PositionalAccuracyNNPtr>());
+            PropertyMap(), vertCRS, GeographicCRS::EPSG_4979, nullptr,
+            "egm08_25.gtx", std::vector<PositionalAccuracyNNPtr>());
     return BoundCRS::create(vertCRS, GeographicCRS::EPSG_4979, transformation);
 }
 
@@ -6727,6 +6756,22 @@ TEST(operation, compoundCRS_from_WKT2_no_id_to_geogCRS_3D_context) {
     auto list =
         CoordinateOperationFactory::create()->createOperations(src, dst, ctxt);
     ASSERT_GE(list.size(), 1U);
+
+    {
+        // Important here is vgridshift before hgridshift
+        auto op_proj =
+            list[0]->exportToPROJString(PROJStringFormatter::create().get());
+        EXPECT_EQ(
+            op_proj,
+            "+proj=pipeline +step +inv +proj=sterea +lat_0=52.1561605555556 "
+            "+lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 "
+            "+ellps=bessel "
+            "+step +proj=vgridshift +grids=naptrans2008.gtx +multiplier=1 "
+            "+step +proj=hgridshift +grids=rdtrans2008.gsb "
+            "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
+            "+step +proj=axisswap +order=2,1");
+    }
+
     auto wkt2 =
         "COMPOUNDCRS[\"unknown\",\n"
         "  PROJCRS[\"unknown\",\n"
@@ -6759,8 +6804,12 @@ TEST(operation, compoundCRS_from_WKT2_no_id_to_geogCRS_3D_context) {
     for (size_t i = 0; i < list.size(); i++) {
         const auto &op = list[i];
         const auto &op2 = list2[i];
-        EXPECT_TRUE(
-            op->isEquivalentTo(op2.get(), IComparable::Criterion::EQUIVALENT));
+        auto op_proj =
+            op->exportToPROJString(PROJStringFormatter::create().get());
+        auto op2_proj =
+            op2->exportToPROJString(PROJStringFormatter::create().get());
+        EXPECT_EQ(op_proj, op2_proj) << "op=" << op->nameStr()
+                                     << " op2=" << op2->nameStr();
     }
 }