[3d] Reproject DEM for terrain if the DEM has different CRS (fixes #1…

…7514)

src/3d/qgs3dmapsettings.cpp

@@ -84,7 +84,9 @@ void Qgs3DMapSettings::readXml( const QDomElement &elem, const QgsReadWriteConte
   QString terrainGenType = elemTerrainGenerator.attribute( "type" );
   if ( terrainGenType == "dem" )
   {
-    mTerrainGenerator.reset( new QgsDemTerrainGenerator );
+    QgsDemTerrainGenerator *demTerrainGenerator = new QgsDemTerrainGenerator;
+    demTerrainGenerator->setCrs( mCrs );
+    mTerrainGenerator.reset( demTerrainGenerator );
   }
   else if ( terrainGenType == "quantized-mesh" )
   {

src/3d/terrain/qgsdemterraingenerator.cpp

@@ -35,9 +35,16 @@ QgsRasterLayer *QgsDemTerrainGenerator::layer() const
   return qobject_cast<QgsRasterLayer *>( mLayer.layer.data() );
 }
 
+void QgsDemTerrainGenerator::setCrs( const QgsCoordinateReferenceSystem &crs )
+{
+  mCrs = crs;
+  updateGenerator();
+}
+
 QgsTerrainGenerator *QgsDemTerrainGenerator::clone() const
 {
   QgsDemTerrainGenerator *cloned = new QgsDemTerrainGenerator;
+  cloned->mCrs = mCrs;
   cloned->mLayer = mLayer;
   cloned->mResolution = mResolution;
   cloned->mSkirtHeight = mSkirtHeight;
@@ -66,13 +73,17 @@ void QgsDemTerrainGenerator::writeXml( QDomElement &elem ) const
   elem.setAttribute( "layer", mLayer.layerId );
   elem.setAttribute( "resolution", mResolution );
   elem.setAttribute( "skirt-height", mSkirtHeight );
+
+  // crs is not read/written - it should be the same as destination crs of the map
 }
 
 void QgsDemTerrainGenerator::readXml( const QDomElement &elem )
 {
   mLayer = QgsMapLayerRef( elem.attribute( "layer" ) );
   mResolution = elem.attribute( "resolution" ).toInt();
   mSkirtHeight = elem.attribute( "skirt-height" ).toFloat();
+
+  // crs is not read/written - it should be the same as destination crs of the map
 }
 
 void QgsDemTerrainGenerator::resolveReferences( const QgsProject &project )
@@ -91,7 +102,11 @@ void QgsDemTerrainGenerator::updateGenerator()
   QgsRasterLayer *dem = layer();
   if ( dem )
   {
-    mTerrainTilingScheme = QgsTilingScheme( dem->extent(), dem->crs() );
+    QgsRectangle te = dem->extent();
+    QgsCoordinateTransform terrainToMapTransform( dem->crs(), mCrs );
+    te = terrainToMapTransform.transformBoundingBox( te );
+
+    mTerrainTilingScheme = QgsTilingScheme( te, mCrs );
     delete mHeightMapGenerator;
     mHeightMapGenerator = new QgsDemHeightMapGenerator( dem, mTerrainTilingScheme, mResolution );
   }

src/3d/terrain/qgsdemterraingenerator.h

@@ -45,6 +45,9 @@ class _3D_EXPORT QgsDemTerrainGenerator : public QgsTerrainGenerator
     //! Returns raster layer with elevation model to be used for terrain generation
     QgsRasterLayer *layer() const;
 
+    //! Sets CRS of the terrain
+    void setCrs( const QgsCoordinateReferenceSystem &crs );
+
     //! Sets resolution of the generator (how many elevation samples on one side of a terrain tile)
     void setResolution( int resolution ) { mResolution = resolution; updateGenerator(); }
     //! Returns resolution of the generator (how many elevation samples on one side of a terrain tile)
@@ -73,6 +76,7 @@ class _3D_EXPORT QgsDemTerrainGenerator : public QgsTerrainGenerator
 
     QgsDemHeightMapGenerator *mHeightMapGenerator = nullptr;
 
+    QgsCoordinateReferenceSystem mCrs;
     //! source layer for heights
     QgsMapLayerRef mLayer;
     //! how many vertices to place on one side of the tile

src/3d/terrain/qgsdemterraintileloader_p.cpp

@@ -113,6 +113,7 @@ void QgsDemTerrainTileLoader::onHeightMapReady( int jobId, const QByteArray &hei
 // ---------------------
 
 #include <qgsrasterlayer.h>
+#include <qgsrasterprojector.h>
 #include <QtConcurrent/QtConcurrentRun>
 #include <QFutureWatcher>
 
@@ -132,13 +133,22 @@ QgsDemHeightMapGenerator::~QgsDemHeightMapGenerator()
 
 #include <QElapsedTimer>
 
-static QByteArray _readDtmData( QgsRasterDataProvider *provider, const QgsRectangle &extent, int res )
+static QByteArray _readDtmData( QgsRasterDataProvider *provider, const QgsRectangle &extent, int res, const QgsCoordinateReferenceSystem &destCrs )
 {
   QElapsedTimer t;
   t.start();
 
   // TODO: use feedback object? (but GDAL currently does not support cancelation anyway)
-  QgsRasterBlock *block = provider->block( 1, extent, res, res );
+  QgsRasterInterface *input = provider;
+  std::unique_ptr<QgsRasterProjector> projector;
+  if ( provider->crs() != destCrs )
+  {
+    projector.reset( new QgsRasterProjector );
+    projector->setCrs( provider->crs(), destCrs );
+    projector->setInput( provider );
+    input = projector.get();
+  }
+  QgsRasterBlock *block = input->block( 1, extent, res, res );
 
   QByteArray data;
   if ( block )
@@ -168,7 +178,7 @@ int QgsDemHeightMapGenerator::render( int x, int y, int z )
   jd.extent = extent;
   jd.timer.start();
   // make a clone of the data provider so it is safe to use in worker thread
-  jd.future = QtConcurrent::run( _readDtmData, mClonedProvider, extent, mResolution );
+  jd.future = QtConcurrent::run( _readDtmData, mClonedProvider, extent, mResolution, mTilingScheme.crs() );
 
   QFutureWatcher<QByteArray> *fw = new QFutureWatcher<QByteArray>;
   fw->setFuture( jd.future );

src/app/3d/qgs3dmapconfigwidget.cpp

@@ -90,6 +90,7 @@ void Qgs3DMapConfigWidget::apply()
     if ( tGenNeedsUpdate )
     {
       QgsDemTerrainGenerator *demTerrainGen = new QgsDemTerrainGenerator;
+      demTerrainGen->setCrs( mMap->crs() );
       demTerrainGen->setLayer( demLayer );
       demTerrainGen->setResolution( spinTerrainResolution->value() );
       demTerrainGen->setSkirtHeight( spinTerrainSkirtHeight->value() );
@@ -139,6 +140,7 @@ void Qgs3DMapConfigWidget::updateMaxZoomLevel()
   if ( demLayer )
   {
     QgsDemTerrainGenerator *demTerrainGen = new QgsDemTerrainGenerator;
+    demTerrainGen->setCrs( mMap->crs() );
     demTerrainGen->setLayer( demLayer );
     demTerrainGen->setResolution( spinTerrainResolution->value() );
     tGen.reset( demTerrainGen );