diff --git a/src/3d/terrain/qgsdemterraintileloader_p.cpp b/src/3d/terrain/qgsdemterraintileloader_p.cpp
index b4cdd37b7dbe..87bce73a8e05 100644
--- a/src/3d/terrain/qgsdemterraintileloader_p.cpp
+++ b/src/3d/terrain/qgsdemterraintileloader_p.cpp
@@ -157,6 +157,7 @@ QgsDemHeightMapGenerator::QgsDemHeightMapGenerator( QgsRasterLayer *dtm, const Q
   , mResolution( resolution )
   , mLastJobId( 0 )
   , mDownloader( dtm ? nullptr : new QgsTerrainDownloader( transformContext ) )
+  , mTransformContext( transformContext )
 {
 }
 
@@ -205,9 +206,9 @@ static QByteArray _readDtmData( QgsRasterDataProvider *provider, const QgsRectan
   return data;
 }
 
-static QByteArray _readOnlineDtm( QgsTerrainDownloader *downloader, const QgsRectangle &extent, int res, const QgsCoordinateReferenceSystem &destCrs )
+static QByteArray _readOnlineDtm( QgsTerrainDownloader *downloader, const QgsRectangle &extent, int res, const QgsCoordinateReferenceSystem &destCrs, const QgsCoordinateTransformContext &context )
 {
-  return downloader->getHeightMap( extent, res, destCrs );
+  return downloader->getHeightMap( extent, res, destCrs, context );
 }
 
 int QgsDemHeightMapGenerator::render( const QgsChunkNodeId &nodeId )
@@ -231,7 +232,7 @@ int QgsDemHeightMapGenerator::render( const QgsChunkNodeId &nodeId )
   if ( mDtm )
     jd.future = QtConcurrent::run( _readDtmData, mClonedProvider, extent, mResolution, mTilingScheme.crs() );
   else
-    jd.future = QtConcurrent::run( _readOnlineDtm, mDownloader.get(), extent, mResolution, mTilingScheme.crs() );
+    jd.future = QtConcurrent::run( _readOnlineDtm, mDownloader.get(), extent, mResolution, mTilingScheme.crs(), mTransformContext );
 
   QFutureWatcher<QByteArray> *fw = new QFutureWatcher<QByteArray>( nullptr );
   fw->setFuture( jd.future );
diff --git a/src/3d/terrain/qgsdemterraintileloader_p.h b/src/3d/terrain/qgsdemterraintileloader_p.h
index 75ce55537a6c..d309b7818613 100644
--- a/src/3d/terrain/qgsdemterraintileloader_p.h
+++ b/src/3d/terrain/qgsdemterraintileloader_p.h
@@ -35,6 +35,7 @@
 #include <QMutex>
 
 #include "qgschunknode_p.h"
+#include "qgscoordinatetransformcontext.h"
 #include "qgsrectangle.h"
 #include "qgsterraintileloader_p.h"
 #include "qgstilingscheme.h"
@@ -138,6 +139,8 @@ class QgsDemHeightMapGenerator : public QObject
     mutable QMutex mLazyLoadDtmCoarseDataMutex;
     //! used for height queries
     QByteArray mDtmCoarseData;
+
+    QgsCoordinateTransformContext mTransformContext;
 };
 
 /// @endcond
diff --git a/src/3d/terrain/qgsterraindownloader.cpp b/src/3d/terrain/qgsterraindownloader.cpp
index fcd54e3400a8..5fa23656baf1 100644
--- a/src/3d/terrain/qgsterraindownloader.cpp
+++ b/src/3d/terrain/qgsterraindownloader.cpp
@@ -177,8 +177,8 @@ QByteArray QgsTerrainDownloader::getHeightMap( const QgsRectangle &extentOrig, i
   }
 
   // resample to the desired extent + resolution
-
-  QgsGdalUtils::resampleSingleBandRaster( hSrcDS.get(), hDstDS.get(), GRA_Bilinear );
+  QgsGdalUtils::resampleSingleBandRaster( hSrcDS.get(), hDstDS.get(), GRA_Bilinear,
+                                          context.calculateCoordinateOperation( mOnlineDtm->crs(), destCrs ).toUtf8().constData() );
 
   QByteArray heightMapOut;
   heightMapOut.resize( resOrig * resOrig * sizeof( float ) );
diff --git a/src/core/qgsgdalutils.cpp b/src/core/qgsgdalutils.cpp
index 79cd17880fd2..33f5c61bc84e 100644
--- a/src/core/qgsgdalutils.cpp
+++ b/src/core/qgsgdalutils.cpp
@@ -145,7 +145,7 @@ gdal::dataset_unique_ptr QgsGdalUtils::imageToMemoryDataset( const QImage &image
   return hSrcDS;
 }
 
-bool QgsGdalUtils::resampleSingleBandRaster( GDALDatasetH hSrcDS, GDALDatasetH hDstDS, GDALResampleAlg resampleAlg )
+bool QgsGdalUtils::resampleSingleBandRaster( GDALDatasetH hSrcDS, GDALDatasetH hDstDS, GDALResampleAlg resampleAlg, const char *pszCoordinateOperation )
 {
   gdal::warp_options_unique_ptr psWarpOptions( GDALCreateWarpOptions() );
   psWarpOptions->hSrcDS = hSrcDS;
@@ -160,10 +160,12 @@ bool QgsGdalUtils::resampleSingleBandRaster( GDALDatasetH hSrcDS, GDALDatasetH h
   psWarpOptions->eResampleAlg = resampleAlg;
 
   // Establish reprojection transformer.
-  psWarpOptions->pTransformerArg =
-    GDALCreateGenImgProjTransformer( hSrcDS, GDALGetProjectionRef( hSrcDS ),
-                                     hDstDS, GDALGetProjectionRef( hDstDS ),
-                                     FALSE, 0.0, 1 );
+  char **papszOptions = nullptr;
+  if ( pszCoordinateOperation != nullptr )
+    papszOptions = CSLSetNameValue( papszOptions, "COORDINATE_OPERATION", pszCoordinateOperation );
+  psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer2( hSrcDS, hDstDS, papszOptions );
+  CSLDestroy( papszOptions );
+
   if ( ! psWarpOptions->pTransformerArg )
   {
     return false;
diff --git a/src/core/qgsgdalutils.h b/src/core/qgsgdalutils.h
index 82f518806651..4ad637cf275d 100644
--- a/src/core/qgsgdalutils.h
+++ b/src/core/qgsgdalutils.h
@@ -66,7 +66,7 @@ class CORE_EXPORT QgsGdalUtils
      * \returns TRUE on success
      * \since QGIS 3.8
      */
-    static bool resampleSingleBandRaster( GDALDatasetH hSrcDS, GDALDatasetH hDstDS, GDALResampleAlg resampleAlg );
+    static bool resampleSingleBandRaster( GDALDatasetH hSrcDS, GDALDatasetH hDstDS, GDALResampleAlg resampleAlg, const char *pszCoordinateOperation );
 
     /**
      * Resamples a QImage \a image using GDAL resampler.
diff --git a/src/providers/postgres/raster/qgspostgresrasterprovider.cpp b/src/providers/postgres/raster/qgspostgresrasterprovider.cpp
index dca3c2f39317..b9c7c19289b0 100644
--- a/src/providers/postgres/raster/qgspostgresrasterprovider.cpp
+++ b/src/providers/postgres/raster/qgspostgresrasterprovider.cpp
@@ -492,7 +492,7 @@ bool QgsPostgresRasterProvider::readBlock( int bandNo, const QgsRectangle &viewE
     }
 
     // Resample the raster to the final bounds and resolution
-    if ( ! QgsGdalUtils::resampleSingleBandRaster( tmpDS.get(), dstDS.get(), GDALResampleAlg::GRA_NearestNeighbour ) )
+    if ( ! QgsGdalUtils::resampleSingleBandRaster( tmpDS.get(), dstDS.get(), GDALResampleAlg::GRA_NearestNeighbour, nullptr ) )
     {
       const QString lastError = QString::fromUtf8( CPLGetLastErrorMsg() ) ;
       QgsMessageLog::logMessage( tr( "Unable to resample and transform destination raster for tiles from %1: %2" )
diff --git a/tests/src/core/testqgsgdalutils.cpp b/tests/src/core/testqgsgdalutils.cpp
index b4d319ae2e22..c3032698ca6c 100644
--- a/tests/src/core/testqgsgdalutils.cpp
+++ b/tests/src/core/testqgsgdalutils.cpp
@@ -194,7 +194,7 @@ void TestQgsGdalUtils::testResampleSingleBandRaster()
   gdal::dataset_unique_ptr dstDS = QgsGdalUtils::createSingleBandTiffDataset( outputFilename, GDT_Float32, outputExtent, 2, 2, QgsCoordinateReferenceSystem( "EPSG:4326" ) );
   QVERIFY( dstDS );
 
-  QgsGdalUtils::resampleSingleBandRaster( srcDS.get(), dstDS.get(), GRA_NearestNeighbour );
+  QgsGdalUtils::resampleSingleBandRaster( srcDS.get(), dstDS.get(), GRA_NearestNeighbour, nullptr );
   dstDS.reset();
 
   std::unique_ptr<QgsRasterLayer> layer( new QgsRasterLayer( outputFilename, "test", "gdal" ) );