QGIS API Documentation 3.29.0-Master (19d7edcfed)
qgsdemterraintileloader_p.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsdemterraintileloader_p.cpp
3 --------------------------------------
4 Date : July 2017
5 Copyright : (C) 2017 by Martin Dobias
6 Email : wonder dot sk at gmail dot com
7 ***************************************************************************
8 * *
9 * This program is free software; you can redistribute it and/or modify *
10 * it under the terms of the GNU General Public License as published by *
11 * the Free Software Foundation; either version 2 of the License, or *
12 * (at your option) any later version. *
13 * *
14 ***************************************************************************/
15
17
18#include "qgs3dmapsettings.h"
19#include "qgschunknode_p.h"
22#include "qgseventtracing.h"
24#include "qgsterrainentity_p.h"
27#include "qgsterraingenerator.h"
28
29#include <Qt3DRender/QGeometryRenderer>
30#include <Qt3DCore/QTransform>
31#include <QMutexLocker>
32
34
35static void _heightMapMinMax( const QByteArray &heightMap, float &zMin, float &zMax )
36{
37 const float *zBits = ( const float * ) heightMap.constData();
38 int zCount = heightMap.count() / sizeof( float );
39 bool first = true;
40
41 zMin = zMax = std::numeric_limits<float>::quiet_NaN();
42 for ( int i = 0; i < zCount; ++i )
43 {
44 float z = zBits[i];
45 if ( std::isnan( z ) )
46 continue;
47 if ( first )
48 {
49 zMin = zMax = z;
50 first = false;
51 }
52 zMin = std::min( zMin, z );
53 zMax = std::max( zMax, z );
54 }
55}
56
57
58QgsDemTerrainTileLoader::QgsDemTerrainTileLoader( QgsTerrainEntity *terrain, QgsChunkNode *node, QgsTerrainGenerator *terrainGenerator )
59 : QgsTerrainTileLoader( terrain, node )
60 , mResolution( 0 )
61{
62
63 QgsDemHeightMapGenerator *heightMapGenerator = nullptr;
64 if ( terrainGenerator->type() == QgsTerrainGenerator::Dem )
65 {
66 QgsDemTerrainGenerator *generator = static_cast<QgsDemTerrainGenerator *>( terrainGenerator );
67 heightMapGenerator = generator->heightMapGenerator();
68 mSkirtHeight = generator->skirtHeight();
69 }
70 else if ( terrainGenerator->type() == QgsTerrainGenerator::Online )
71 {
72 QgsOnlineTerrainGenerator *generator = static_cast<QgsOnlineTerrainGenerator *>( terrainGenerator );
73 heightMapGenerator = generator->heightMapGenerator();
74 mSkirtHeight = generator->skirtHeight();
75 }
76 else
77 Q_ASSERT( false );
78
79 // get heightmap asynchronously
80 connect( heightMapGenerator, &QgsDemHeightMapGenerator::heightMapReady, this, &QgsDemTerrainTileLoader::onHeightMapReady );
81 mHeightMapJobId = heightMapGenerator->render( node->tileId() );
82 mResolution = heightMapGenerator->resolution();
83}
84
85Qt3DCore::QEntity *QgsDemTerrainTileLoader::createEntity( Qt3DCore::QEntity *parent )
86{
87 float zMin, zMax;
88 _heightMapMinMax( mHeightMap, zMin, zMax );
89
90 if ( std::isnan( zMin ) || std::isnan( zMax ) )
91 {
92 // no data available for this tile
93 return nullptr;
94 }
95
96 const Qgs3DMapSettings &map = terrain()->map3D();
97 QgsChunkNodeId nodeId = mNode->tileId();
98 QgsRectangle extent = map.terrainGenerator()->tilingScheme().tileToExtent( nodeId );
99 double x0 = extent.xMinimum() - map.origin().x();
100 double y0 = extent.yMinimum() - map.origin().y();
101 double side = extent.width();
102 double half = side / 2;
103
104
105 QgsTerrainTileEntity *entity = new QgsTerrainTileEntity( nodeId );
106
107 // create geometry renderer
108
109 Qt3DRender::QGeometryRenderer *mesh = new Qt3DRender::QGeometryRenderer;
110 mesh->setGeometry( new DemTerrainTileGeometry( mResolution, side, map.terrainVerticalScale(), mSkirtHeight, mHeightMap, mesh ) );
111 entity->addComponent( mesh ); // takes ownership if the component has no parent
112
113 // create material
114
115 createTextureComponent( entity, map.isTerrainShadingEnabled(), map.terrainShadingMaterial(), !map.layers().empty() );
116
117 // create transform
118
119 Qt3DCore::QTransform *transform = nullptr;
120 transform = new Qt3DCore::QTransform();
121 entity->addComponent( transform );
122
123 transform->setScale( side );
124 transform->setTranslation( QVector3D( x0 + half, 0, - ( y0 + half ) ) );
125
126 mNode->setExactBbox( QgsAABB( x0, zMin * map.terrainVerticalScale(), -y0, x0 + side, zMax * map.terrainVerticalScale(), -( y0 + side ) ) );
127 mNode->updateParentBoundingBoxesRecursively();
128
129 entity->setParent( parent );
130 return entity;
131}
132
133void QgsDemTerrainTileLoader::onHeightMapReady( int jobId, const QByteArray &heightMap )
134{
135 if ( mHeightMapJobId == jobId )
136 {
137 this->mHeightMap = heightMap;
138 mHeightMapJobId = -1;
139
140 // continue loading - texture
141 loadTexture();
142 }
143}
144
145
146// ---------------------
147
148#include <qgsrasterlayer.h>
149#include <qgsrasterprojector.h>
150#include <QtConcurrent/QtConcurrentRun>
151#include <QFutureWatcher>
152#include "qgsterraindownloader.h"
153
154QgsDemHeightMapGenerator::QgsDemHeightMapGenerator( QgsRasterLayer *dtm, const QgsTilingScheme &tilingScheme, int resolution, const QgsCoordinateTransformContext &transformContext )
155 : mDtmExtent( dtm ? dtm->extent() : QgsRectangle() )
156 , mClonedProvider( dtm ? qgis::down_cast<QgsRasterDataProvider *>( dtm->dataProvider()->clone() ) : nullptr )
157 , mTilingScheme( tilingScheme )
158 , mResolution( resolution )
159 , mLastJobId( 0 )
160 , mDownloader( dtm ? nullptr : new QgsTerrainDownloader( transformContext ) )
161 , mTransformContext( transformContext )
162{
163}
164
165QgsDemHeightMapGenerator::~QgsDemHeightMapGenerator()
166{
167 delete mClonedProvider;
168}
169
170
171static QByteArray _readDtmData( QgsRasterDataProvider *provider, const QgsRectangle &extent, int res, const QgsCoordinateReferenceSystem &destCrs, const QgsRectangle &clippingExtent )
172{
173 provider->moveToThread( QThread::currentThread() );
174
175 QgsEventTracing::ScopedEvent e( QStringLiteral( "3D" ), QStringLiteral( "DEM" ) );
176
177 // TODO: use feedback object? (but GDAL currently does not support cancellation anyway)
178 QgsRasterInterface *input = provider;
179 std::unique_ptr<QgsRasterProjector> projector;
180 if ( provider->crs() != destCrs )
181 {
182 projector.reset( new QgsRasterProjector );
183 projector->setCrs( provider->crs(), destCrs, provider->transformContext() );
184 projector->setInput( provider );
185 input = projector.get();
186 }
187 std::unique_ptr< QgsRasterBlock > block( input->block( 1, extent, res, res ) );
188
189 QByteArray data;
190 if ( block )
191 {
192 block->convert( Qgis::DataType::Float32 ); // currently we expect just floats
193
194 // set noData outside our clippingExtent
195 const QRect subRect = QgsRasterBlock::subRect( extent, block->width(), block->height(), clippingExtent );
196 if ( !block->hasNoDataValue() )
197 {
198 // QgsRasterBlock::setIsNoDataExcept() misbehaves without a defined no data value
199 // see https://github.com/qgis/QGIS/issues/51285
200 block->setNoDataValue( std::numeric_limits<float>::lowest() );
201 }
202 block->setIsNoDataExcept( subRect );
203
204 data = block->data();
205 data.detach(); // this should make a deep copy
206
207 if ( block->hasNoData() )
208 {
209 // turn all no-data values into NaN in the output array
210 float *floatData = reinterpret_cast<float *>( data.data() );
211 Q_ASSERT( data.count() % sizeof( float ) == 0 );
212 int count = data.count() / sizeof( float );
213 for ( int i = 0; i < count; ++i )
214 {
215 if ( block->isNoData( i ) )
216 floatData[i] = std::numeric_limits<float>::quiet_NaN();
217 }
218 }
219 }
220
221 delete provider;
222 return data;
223}
224
225static QByteArray _readOnlineDtm( QgsTerrainDownloader *downloader, const QgsRectangle &extent, int res, const QgsCoordinateReferenceSystem &destCrs, const QgsCoordinateTransformContext &context )
226{
227 return downloader->getHeightMap( extent, res, destCrs, context );
228}
229
230int QgsDemHeightMapGenerator::render( const QgsChunkNodeId &nodeId )
231{
232 QgsEventTracing::addEvent( QgsEventTracing::AsyncBegin, QStringLiteral( "3D" ), QStringLiteral( "DEM" ), nodeId.text() );
233
234 // extend the rect by half-pixel on each side? to get the values in "corners"
235 QgsRectangle extent = mTilingScheme.tileToExtent( nodeId );
236 float mapUnitsPerPixel = extent.width() / mResolution;
237 extent.grow( mapUnitsPerPixel / 2 );
238 // but make sure not to go beyond the root tile's full extent (returns invalid values)
239 QgsRectangle rootTileExtent = mTilingScheme.tileToExtent( 0, 0, 0 );
240 extent = extent.intersect( rootTileExtent );
241
242 JobData jd;
243 jd.jobId = ++mLastJobId;
244 jd.tileId = nodeId;
245 jd.extent = extent;
246 jd.timer.start();
247 QFutureWatcher<QByteArray> *fw = new QFutureWatcher<QByteArray>( nullptr );
248 connect( fw, &QFutureWatcher<QByteArray>::finished, this, &QgsDemHeightMapGenerator::onFutureFinished );
249 connect( fw, &QFutureWatcher<QByteArray>::finished, fw, &QObject::deleteLater );
250 if ( mClonedProvider )
251 {
252 // make a clone of the data provider so it is safe to use in worker thread
253 std::unique_ptr< QgsRasterDataProvider > clonedProviderClone( mClonedProvider->clone() );
254 clonedProviderClone->moveToThread( nullptr );
255 jd.future = QtConcurrent::run( _readDtmData, clonedProviderClone.release(), extent, mResolution, mTilingScheme.crs(), mTilingScheme.fullExtent() );
256 }
257 else
258 {
259 jd.future = QtConcurrent::run( _readOnlineDtm, mDownloader.get(), extent, mResolution, mTilingScheme.crs(), mTransformContext );
260 }
261
262 fw->setFuture( jd.future );
263
264 mJobs.insert( fw, jd );
265
266 return jd.jobId;
267}
268
269void QgsDemHeightMapGenerator::waitForFinished()
270{
271 for ( auto it = mJobs.keyBegin(); it != mJobs.keyEnd(); it++ )
272 {
273 QFutureWatcher<QByteArray> *fw = *it;
274 disconnect( fw, &QFutureWatcher<QByteArray>::finished, this, &QgsDemHeightMapGenerator::onFutureFinished );
275 disconnect( fw, &QFutureWatcher<QByteArray>::finished, fw, &QObject::deleteLater );
276 }
277 QVector<QFutureWatcher<QByteArray>*> toBeDeleted;
278 for ( auto it = mJobs.keyBegin(); it != mJobs.keyEnd(); it++ )
279 {
280 QFutureWatcher<QByteArray> *fw = *it;
281 fw->waitForFinished();
282 JobData jobData = mJobs.value( fw );
283 toBeDeleted.push_back( fw );
284
285 QByteArray data = jobData.future.result();
286 emit heightMapReady( jobData.jobId, data );
287 }
288
289 for ( QFutureWatcher<QByteArray> *fw : toBeDeleted )
290 {
291 mJobs.remove( fw );
292 fw->deleteLater();
293 }
294}
295
296void QgsDemHeightMapGenerator::lazyLoadDtmCoarseData( int res, const QgsRectangle &rect )
297{
298 QMutexLocker locker( &mLazyLoadDtmCoarseDataMutex );
299 if ( mDtmCoarseData.isEmpty() )
300 {
301 std::unique_ptr< QgsRasterBlock > block( mClonedProvider->block( 1, rect, res, res ) );
302 block->convert( Qgis::DataType::Float32 );
303 mDtmCoarseData = block->data();
304 mDtmCoarseData.detach(); // make a deep copy
305 }
306}
307
308float QgsDemHeightMapGenerator::heightAt( double x, double y )
309{
310 if ( !mClonedProvider )
311 return 0; // TODO: calculate heights for online DTM
312
313 // TODO: this is quite a primitive implementation: better to use heightmaps currently in use
314 int res = 1024;
315 lazyLoadDtmCoarseData( res, mDtmExtent );
316
317 int cellX = ( int )( ( x - mDtmExtent.xMinimum() ) / mDtmExtent.width() * res + .5f );
318 int cellY = ( int )( ( mDtmExtent.yMaximum() - y ) / mDtmExtent.height() * res + .5f );
319 cellX = std::clamp( cellX, 0, res - 1 );
320 cellY = std::clamp( cellY, 0, res - 1 );
321
322 const float *data = ( const float * ) mDtmCoarseData.constData();
323 return data[cellX + cellY * res];
324}
325
326void QgsDemHeightMapGenerator::onFutureFinished()
327{
328 QFutureWatcher<QByteArray> *fw = static_cast<QFutureWatcher<QByteArray>*>( sender() );
329 Q_ASSERT( fw );
330 Q_ASSERT( mJobs.contains( fw ) );
331 JobData jobData = mJobs.value( fw );
332
333 mJobs.remove( fw );
334 fw->deleteLater();
335
336 QgsEventTracing::addEvent( QgsEventTracing::AsyncEnd, QStringLiteral( "3D" ), QStringLiteral( "DEM" ), jobData.tileId.text() );
337
338 QByteArray data = jobData.future.result();
339 emit heightMapReady( jobData.jobId, data );
340}
341
@ Float32
Thirty two bit floating point (float)
double terrainVerticalScale() const
Returns vertical scale (exaggeration) of terrain.
QgsTerrainGenerator * terrainGenerator() const
Returns the terrain generator.
bool isTerrainShadingEnabled() const
Returns whether terrain shading is enabled.
QgsPhongMaterialSettings terrainShadingMaterial() const
Returns terrain shading material.
QList< QgsMapLayer * > layers() const
Returns the list of 3D map layers to be rendered in the scene.
QgsVector3D origin() const
Returns coordinates in map CRS at which 3D scene has origin (0,0,0)
3
Definition: qgsaabb.h:34
This class represents a coordinate reference system (CRS).
Contains information about the context in which a coordinate transform is executed.
QgsCoordinateTransformContext transformContext() const
Returns data provider coordinate transform context.
virtual QgsCoordinateReferenceSystem crs() const =0
Returns the coordinate system for the data source.
QgsDemHeightMapGenerator * heightMapGenerator()
Returns height map generator object - takes care of extraction of elevations from the layer)
float skirtHeight() const
Returns skirt height (in world units). Skirts at the edges of terrain tiles help hide cracks between ...
QgsDemHeightMapGenerator * heightMapGenerator()
Returns height map generator object - takes care of extraction of elevations from the layer)
float skirtHeight() const
Returns skirt height (in world units). Skirts at the edges of terrain tiles help hide cracks between ...
static QRect subRect(const QgsRectangle &extent, int width, int height, const QgsRectangle &subExtent)
For extent and width, height find rectangle covered by subextent.
Base class for raster data providers.
bool setInput(QgsRasterInterface *input) override
Set input.
Base class for processing filters like renderers, reprojector, resampler etc.
virtual QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr)=0
Read block of data using given extent and size.
Represents a raster layer.
QgsRasterProjector implements approximate projection support for it calculates grid of points in sour...
A rectangle specified with double values.
Definition: qgsrectangle.h:42
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:188
double yMinimum() const SIP_HOLDGIL
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:198
void grow(double delta)
Grows the rectangle in place by the specified amount.
Definition: qgsrectangle.h:296
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:223
QgsRectangle intersect(const QgsRectangle &rect) const
Returns the intersection with the given rectangle.
Definition: qgsrectangle.h:333
QByteArray getHeightMap(const QgsRectangle &extentOrig, int res, const QgsCoordinateReferenceSystem &destCrs, const QgsCoordinateTransformContext &context=QgsCoordinateTransformContext(), QString tmpFilenameImg=QString(), QString tmpFilenameTif=QString())
For given extent and resolution (number of pixels for width/height) in specified CRS,...
@ Dem
Terrain is built from raster layer with digital elevation model.
@ Online
Terrain is built from downloaded tiles with digital elevation model.
virtual Type type() const =0
What texture generator implementation is this.
const QgsTilingScheme & tilingScheme() const
Returns tiling scheme of the terrain.
QgsRectangle tileToExtent(int x, int y, int z) const
Returns map coordinates of the extent of a tile.
double y() const
Returns Y coordinate.
Definition: qgsvector3d.h:51
double x() const
Returns X coordinate.
Definition: qgsvector3d.h:49