diff --git a/alg/gdaltransformer.cpp b/alg/gdaltransformer.cpp
index 72d0c6884264..6c8487915163 100644
--- a/alg/gdaltransformer.cpp
+++ b/alg/gdaltransformer.cpp
@@ -2149,6 +2149,12 @@ GDALCreateGenImgProjTransformer2( GDALDatasetH hSrcDS, GDALDatasetH hDstDS,
             InsertCenterLong( hSrcDS, &oSrcSRS, aosOptions );
         }
 
+        if( CPLFetchBool(papszOptions, "PROMOTE_TO_3D", false) )
+        {
+            oSrcSRS.PromoteTo3D(nullptr);
+            oDstSRS.PromoteTo3D(nullptr);
+        }
+
         if( !(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&
               dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0) )
         {
diff --git a/apps/gdalwarp_lib.cpp b/apps/gdalwarp_lib.cpp
index fac5aa06a1c5..94f9ac3a7b47 100644
--- a/apps/gdalwarp_lib.cpp
+++ b/apps/gdalwarp_lib.cpp
@@ -559,27 +559,16 @@ GDALWarpAppOptions* GDALWarpAppOptionsClone(const GDALWarpAppOptions *psOptionsI
 
 #ifdef USE_PROJ_BASED_VERTICAL_SHIFT_METHOD
 
-/************************************************************************/
-/*                      ApplyVerticalShift()                            */
-/************************************************************************/
-
-static bool ApplyVerticalShift( GDALDatasetH hWrkSrcDS,
-                                const GDALWarpAppOptions* psOptions,
-                                GDALWarpOptions *psWO )
+static bool MustApplyVerticalShift( GDALDatasetH hWrkSrcDS,
+                                    const GDALWarpAppOptions* psOptions,
+                                    OGRSpatialReference& oSRSSrc,
+                                    OGRSpatialReference& oSRSDst,
+                                    bool& bSrcHasVertAxis,
+                                    bool& bDstHasVertAxis )
 {
-    bool bApplyVShift = false;
-
-    if( psOptions->bVShift )
-    {
-        bApplyVShift = true;
-        psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions,
-                                                 "APPLY_VERTICAL_SHIFT",
-                                                 "YES");
-    }
+    bool bApplyVShift = psOptions->bVShift;
 
     // Check if we must do vertical shift grid transform
-    OGRSpatialReference oSRSSrc;
-    OGRSpatialReference oSRSDst;
     const char* pszSrcWKT = CSLFetchNameValue(psOptions->papszTO, "SRC_SRS");
     if( pszSrcWKT )
         oSRSSrc.SetFromUserInput( pszSrcWKT );
@@ -594,14 +583,46 @@ static bool ApplyVerticalShift( GDALDatasetH hWrkSrcDS,
     if( pszDstWKT )
         oSRSDst.SetFromUserInput( pszDstWKT );
 
-    const bool bSrcHasVertAxis =
+    bSrcHasVertAxis =
         oSRSSrc.IsCompound() ||
         ((oSRSSrc.IsProjected() || oSRSSrc.IsGeographic()) && oSRSSrc.GetAxesCount() == 3);
 
-    const bool bDstHasVertAxis =
+    bDstHasVertAxis =
         oSRSDst.IsCompound() ||
         ((oSRSDst.IsProjected() || oSRSDst.IsGeographic()) && oSRSDst.GetAxesCount() == 3);
 
+    if( (GDALGetRasterCount(hWrkSrcDS) == 1 || psOptions->bVShift) &&
+        (bSrcHasVertAxis || bDstHasVertAxis) )
+    {
+        bApplyVShift = true;
+    }
+    return bApplyVShift;
+}
+
+/************************************************************************/
+/*                      ApplyVerticalShift()                            */
+/************************************************************************/
+
+static bool ApplyVerticalShift( GDALDatasetH hWrkSrcDS,
+                                const GDALWarpAppOptions* psOptions,
+                                GDALWarpOptions *psWO )
+{
+    if( psOptions->bVShift )
+    {
+        psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions,
+                                                 "APPLY_VERTICAL_SHIFT",
+                                                 "YES");
+    }
+
+    OGRSpatialReference oSRSSrc;
+    OGRSpatialReference oSRSDst;
+    bool bSrcHasVertAxis = false;
+    bool bDstHasVertAxis = false;
+    bool bApplyVShift = MustApplyVerticalShift( hWrkSrcDS, psOptions,
+                                                oSRSSrc, oSRSDst,
+                                                bSrcHasVertAxis,
+                                                bDstHasVertAxis );
+
     if( (GDALGetRasterCount(hWrkSrcDS) == 1 || psOptions->bVShift) &&
         (bSrcHasVertAxis || bDstHasVertAxis) )
     {
@@ -2228,6 +2249,22 @@ GDALDatasetH GDALWarpDirect( const char *pszDest, GDALDatasetH hDstDS,
         psOptions->papszTO = CSLSetNameValue(psOptions->papszTO,
                                              "STRIP_VERT_CS", "YES");
     }
+    else if( nSrcCount > 0 )
+    {
+        bool bSrcHasVertAxis = false;
+        bool bDstHasVertAxis = false;
+        OGRSpatialReference oSRSSrc;
+        OGRSpatialReference oSRSDst;
+        bool bApplyVShift = MustApplyVerticalShift( pahSrcDS[0], psOptions,
+                                                    oSRSSrc, oSRSDst,
+                                                    bSrcHasVertAxis,
+                                                    bDstHasVertAxis );
+        if( bApplyVShift )
+        {
+            psOptions->papszTO = CSLSetNameValue(psOptions->papszTO,
+                                                 "PROMOTE_TO_3D", "YES");
+        }
+    }
 #else
     psOptions->papszTO = CSLSetNameValue(psOptions->papszTO,
                                          "STRIP_VERT_CS", "YES");
@@ -2715,7 +2752,7 @@ GDALDatasetH GDALWarpDirect( const char *pszDest, GDALDatasetH hDstDS,
 #ifdef USE_PROJ_BASED_VERTICAL_SHIFT_METHOD
         if( !psOptions->bNoVShift )
         {
-            // Can modify both psWO->papszWarpOptions
+            // Can modify psWO->papszWarpOptions
             if( ApplyVerticalShift( hWrkSrcDS, psOptions, psWO ) )
             {
                 bUseApproxTransformer = false;
diff --git a/autotest/alg/applyverticalshiftgrid.py b/autotest/alg/applyverticalshiftgrid.py
index 2c15fb78762e..b95731af19a8 100755
--- a/autotest/alg/applyverticalshiftgrid.py
+++ b/autotest/alg/applyverticalshiftgrid.py
@@ -295,7 +295,8 @@ def test_applyverticalshiftgrid_6():
     grid_ds = None
 
     ds = gdal.Warp('', '../gcore/data/byte.tif', format='MEM',
-                   dstSRS='+proj=utm +zone=11 +datum=NAD27 +geoidgrids=./tmp/applyverticalshiftgrid_6.gtx +vunits=m +no_defs')
+                   srcSRS='EPSG:32611',
+                   dstSRS='+proj=utm +zone=11 +datum=WGS84 +geoidgrids=./tmp/applyverticalshiftgrid_6.gtx +vunits=m +no_defs')
     cs = ds.GetRasterBand(1).Checksum()
     assert cs == 4783
 
@@ -314,7 +315,8 @@ def test_applyverticalshiftgrid_7():
     grid_ds = None
 
     ds = gdal.Warp('', '../gcore/data/byte.tif', format='MEM',
-                   dstSRS='+proj=utm +zone=11 +datum=NAD27 +geoidgrids=./tmp/applyverticalshiftgrid_7.gtx +vunits=m +no_defs')
+                   srcSRS='EPSG:32611',
+                   dstSRS='+proj=utm +zone=11 +datum=WGS84 +geoidgrids=./tmp/applyverticalshiftgrid_7.gtx +vunits=m +no_defs')
     cs = ds.GetRasterBand(1).Checksum()
     assert cs == 4783