QGIS API Documentation  3.8.0-Zanzibar (11aff65)
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"
20 #include "qgsdemterraingenerator.h"
23 #include "qgsterrainentity_p.h"
25 #include "qgsterraintileentity_p.h"
26 
27 #include <Qt3DRender/QGeometryRenderer>
28 
30 
31 static void _heightMapMinMax( const QByteArray &heightMap, float &zMin, float &zMax )
32 {
33  const float *zBits = ( const float * ) heightMap.constData();
34  int zCount = heightMap.count() / sizeof( float );
35  bool first = true;
36 
37  zMin = zMax = std::numeric_limits<float>::quiet_NaN();
38  for ( int i = 0; i < zCount; ++i )
39  {
40  float z = zBits[i];
41  if ( std::isnan( z ) )
42  continue;
43  if ( first )
44  {
45  zMin = zMax = z;
46  first = false;
47  }
48  zMin = std::min( zMin, z );
49  zMax = std::max( zMax, z );
50  }
51 }
52 
53 
54 QgsDemTerrainTileLoader::QgsDemTerrainTileLoader( QgsTerrainEntity *terrain, QgsChunkNode *node )
55  : QgsTerrainTileLoader( terrain, node )
56  , mResolution( 0 )
57 {
58 
59  const Qgs3DMapSettings &map = terrain->map3D();
60 
61  QgsDemHeightMapGenerator *heightMapGenerator = nullptr;
63  {
64  QgsDemTerrainGenerator *generator = static_cast<QgsDemTerrainGenerator *>( map.terrainGenerator() );
65  heightMapGenerator = generator->heightMapGenerator();
66  mSkirtHeight = generator->skirtHeight();
67  }
68  else if ( map.terrainGenerator()->type() == QgsTerrainGenerator::Online )
69  {
70  QgsOnlineTerrainGenerator *generator = static_cast<QgsOnlineTerrainGenerator *>( map.terrainGenerator() );
71  heightMapGenerator = generator->heightMapGenerator();
72  mSkirtHeight = generator->skirtHeight();
73  }
74  else
75  Q_ASSERT( false );
76 
77  // get heightmap asynchronously
78  connect( heightMapGenerator, &QgsDemHeightMapGenerator::heightMapReady, this, &QgsDemTerrainTileLoader::onHeightMapReady );
79  mHeightMapJobId = heightMapGenerator->render( node->tileX(), node->tileY(), node->tileZ() );
80  mResolution = heightMapGenerator->resolution();
81 }
82 
83 Qt3DCore::QEntity *QgsDemTerrainTileLoader::createEntity( Qt3DCore::QEntity *parent )
84 {
85  float zMin, zMax;
86  _heightMapMinMax( mHeightMap, zMin, zMax );
87 
88  if ( std::isnan( zMin ) || std::isnan( zMax ) )
89  {
90  // no data available for this tile
91  return nullptr;
92  }
93 
94  const Qgs3DMapSettings &map = terrain()->map3D();
95  QgsRectangle extent = map.terrainGenerator()->tilingScheme().tileToExtent( mNode->tileX(), mNode->tileY(), mNode->tileZ() ); //node->extent;
96  double x0 = extent.xMinimum() - map.origin().x();
97  double y0 = extent.yMinimum() - map.origin().y();
98  double side = extent.width();
99  double half = side / 2;
100 
101 
102  QgsTerrainTileEntity *entity = new QgsTerrainTileEntity;
103 
104  // create geometry renderer
105 
106  Qt3DRender::QGeometryRenderer *mesh = new Qt3DRender::QGeometryRenderer;
107  mesh->setGeometry( new DemTerrainTileGeometry( mResolution, side, map.terrainVerticalScale(), mSkirtHeight, mHeightMap, mesh ) );
108  entity->addComponent( mesh ); // takes ownership if the component has no parent
109 
110  // create material
111 
112  createTextureComponent( entity, map.isTerrainShadingEnabled(), map.terrainShadingMaterial() );
113 
114  // create transform
115 
116  Qt3DCore::QTransform *transform = nullptr;
117  transform = new Qt3DCore::QTransform();
118  entity->addComponent( transform );
119 
120  transform->setScale( side );
121  transform->setTranslation( QVector3D( x0 + half, 0, - ( y0 + half ) ) );
122 
123  mNode->setExactBbox( QgsAABB( x0, zMin * map.terrainVerticalScale(), -y0, x0 + side, zMax * map.terrainVerticalScale(), -( y0 + side ) ) );
124 
125  entity->setEnabled( false );
126  entity->setParent( parent );
127  return entity;
128 }
129 
130 void QgsDemTerrainTileLoader::onHeightMapReady( int jobId, const QByteArray &heightMap )
131 {
132  if ( mHeightMapJobId == jobId )
133  {
134  this->mHeightMap = heightMap;
135  mHeightMapJobId = -1;
136 
137  // continue loading - texture
138  loadTexture();
139  }
140 }
141 
142 
143 // ---------------------
144 
145 #include <qgsrasterlayer.h>
146 #include <qgsrasterprojector.h>
147 #include <QtConcurrent/QtConcurrentRun>
148 #include <QFutureWatcher>
149 #include "qgsterraindownloader.h"
150 
151 QgsDemHeightMapGenerator::QgsDemHeightMapGenerator( QgsRasterLayer *dtm, const QgsTilingScheme &tilingScheme, int resolution, const QgsCoordinateTransformContext &transformContext )
152  : mDtm( dtm )
153  , mClonedProvider( dtm ? ( QgsRasterDataProvider * )dtm->dataProvider()->clone() : nullptr )
154  , mTilingScheme( tilingScheme )
155  , mResolution( resolution )
156  , mLastJobId( 0 )
157  , mDownloader( dtm ? nullptr : new QgsTerrainDownloader( transformContext ) )
158 {
159 }
160 
161 QgsDemHeightMapGenerator::~QgsDemHeightMapGenerator()
162 {
163  delete mClonedProvider;
164 }
165 
166 #include <QElapsedTimer>
167 
168 static QByteArray _readDtmData( QgsRasterDataProvider *provider, const QgsRectangle &extent, int res, const QgsCoordinateReferenceSystem &destCrs )
169 {
170  QElapsedTimer t;
171  t.start();
172 
173  // TODO: use feedback object? (but GDAL currently does not support cancellation anyway)
174  QgsRasterInterface *input = provider;
175  std::unique_ptr<QgsRasterProjector> projector;
176  if ( provider->crs() != destCrs )
177  {
178  projector.reset( new QgsRasterProjector );
179  projector->setCrs( provider->crs(), destCrs, provider->transformContext() );
180  projector->setInput( provider );
181  input = projector.get();
182  }
183  std::unique_ptr< QgsRasterBlock > block( input->block( 1, extent, res, res ) );
184 
185  QByteArray data;
186  if ( block )
187  {
188  block->convert( Qgis::Float32 ); // currently we expect just floats
189  data = block->data();
190  data.detach(); // this should make a deep copy
191 
192  if ( block->hasNoData() )
193  {
194  // turn all no-data values into NaN in the output array
195  float *floatData = reinterpret_cast<float *>( data.data() );
196  Q_ASSERT( data.count() % sizeof( float ) == 0 );
197  int count = data.count() / sizeof( float );
198  for ( int i = 0; i < count; ++i )
199  {
200  if ( block->isNoData( i ) )
201  floatData[i] = std::numeric_limits<float>::quiet_NaN();
202  }
203  }
204  }
205  return data;
206 }
207 
208 
209 static QByteArray _readOnlineDtm( QgsTerrainDownloader *downloader, const QgsRectangle &extent, int res, const QgsCoordinateReferenceSystem &destCrs )
210 {
211  return downloader->getHeightMap( extent, res, destCrs );
212 }
213 
214 int QgsDemHeightMapGenerator::render( int x, int y, int z )
215 {
216  Q_ASSERT( mJobs.isEmpty() ); // should be always just one active job...
217 
218  // extend the rect by half-pixel on each side? to get the values in "corners"
219  QgsRectangle extent = mTilingScheme.tileToExtent( x, y, z );
220  float mapUnitsPerPixel = extent.width() / mResolution;
221  extent.grow( mapUnitsPerPixel / 2 );
222  // but make sure not to go beyond the full extent (returns invalid values)
223  QgsRectangle fullExtent = mTilingScheme.tileToExtent( 0, 0, 0 );
224  extent = extent.intersect( fullExtent );
225 
226  JobData jd;
227  jd.jobId = ++mLastJobId;
228  jd.extent = extent;
229  jd.timer.start();
230  // make a clone of the data provider so it is safe to use in worker thread
231  if ( mDtm )
232  jd.future = QtConcurrent::run( _readDtmData, mClonedProvider, extent, mResolution, mTilingScheme.crs() );
233  else
234  jd.future = QtConcurrent::run( _readOnlineDtm, mDownloader.get(), extent, mResolution, mTilingScheme.crs() );
235 
236  QFutureWatcher<QByteArray> *fw = new QFutureWatcher<QByteArray>( nullptr );
237  fw->setFuture( jd.future );
238  connect( fw, &QFutureWatcher<QByteArray>::finished, this, &QgsDemHeightMapGenerator::onFutureFinished );
239 
240  mJobs.insert( fw, jd );
241 
242  return jd.jobId;
243 }
244 
245 QByteArray QgsDemHeightMapGenerator::renderSynchronously( int x, int y, int z )
246 {
247  // extend the rect by half-pixel on each side? to get the values in "corners"
248  QgsRectangle extent = mTilingScheme.tileToExtent( x, y, z );
249  float mapUnitsPerPixel = extent.width() / mResolution;
250  extent.grow( mapUnitsPerPixel / 2 );
251  // but make sure not to go beyond the full extent (returns invalid values)
252  QgsRectangle fullExtent = mTilingScheme.tileToExtent( 0, 0, 0 );
253  extent = extent.intersect( fullExtent );
254 
255  std::unique_ptr< QgsRasterBlock > block( mDtm->dataProvider()->block( 1, extent, mResolution, mResolution ) );
256 
257  QByteArray data;
258  if ( block )
259  {
260  block->convert( Qgis::Float32 ); // currently we expect just floats
261  data = block->data();
262  data.detach(); // this should make a deep copy
263  }
264 
265  return data;
266 }
267 
268 float QgsDemHeightMapGenerator::heightAt( double x, double y )
269 {
270  if ( !mDtm )
271  return 0; // TODO: calculate heights for online DTM
272 
273  // TODO: this is quite a primitive implementation: better to use heightmaps currently in use
274  int res = 1024;
275  QgsRectangle rect = mDtm->extent();
276  if ( mDtmCoarseData.isEmpty() )
277  {
278  std::unique_ptr< QgsRasterBlock > block( mDtm->dataProvider()->block( 1, rect, res, res ) );
279  block->convert( Qgis::Float32 );
280  mDtmCoarseData = block->data();
281  mDtmCoarseData.detach(); // make a deep copy
282  }
283 
284  int cellX = ( int )( ( x - rect.xMinimum() ) / rect.width() * res + .5f );
285  int cellY = ( int )( ( rect.yMaximum() - y ) / rect.height() * res + .5f );
286  cellX = qBound( 0, cellX, res - 1 );
287  cellY = qBound( 0, cellY, res - 1 );
288 
289  const float *data = ( const float * ) mDtmCoarseData.constData();
290  return data[cellX + cellY * res];
291 }
292 
293 void QgsDemHeightMapGenerator::onFutureFinished()
294 {
295  QFutureWatcher<QByteArray> *fw = static_cast<QFutureWatcher<QByteArray>*>( sender() );
296  Q_ASSERT( fw );
297  Q_ASSERT( mJobs.contains( fw ) );
298  JobData jobData = mJobs.value( fw );
299 
300  mJobs.remove( fw );
301  fw->deleteLater();
302 
303  QByteArray data = jobData.future.result();
304  emit heightMapReady( jobData.jobId, data );
305 }
306 
3 Axis-aligned bounding box - in world coords.
Definition: qgsaabb.h:30
A rectangle specified with double values.
Definition: qgsrectangle.h:41
bool setInput(QgsRasterInterface *input) override
Set input.
QgsTerrainGenerator * clone() const override
Makes a copy of the current instance.
bool isTerrainShadingEnabled() const
Returns whether terrain shading is enabled.
QgsDemHeightMapGenerator * heightMapGenerator()
Returns height map generator object - takes care of extraction of elevations from the layer) ...
QgsPhongMaterialSettings terrainShadingMaterial() const
Returns terrain shading material.
This class provides qgis with the ability to render raster datasets onto the mapcanvas.
3 Implementation of terrain generator that uses online resources to download heightmaps.
virtual Type type() const =0
What texture generator implementation is this.
int resolution() const
Returns resolution of the generator (how many elevation samples on one side of a terrain tile) ...
Thirty two bit floating point (float)
Definition: qgis.h:87
float skirtHeight() const
Returns skirt height (in world units). Skirts at the edges of terrain tiles help hide cracks between ...
double y() const
Returns Y coordinate.
Definition: qgsvector3d.h:51
3 Definition of the world
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, download necessary tile images (if not cached already) and produce height map out of them (byte array of res*res float values)
QgsVector3D origin() const
Returns coordinates in map CRS at which 3D scene has origin (0,0,0)
Terrain is built from raster layer with digital elevation model.
void grow(double delta)
Grows the rectangle in place by the specified amount.
Definition: qgsrectangle.h:275
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.
virtual QgsCoordinateReferenceSystem crs() const =0
Returns the coordinate system for the data source.
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:202
QgsRectangle intersect(const QgsRectangle &rect) const
Returns the intersection with the given rectangle.
Definition: qgsrectangle.h:312
float skirtHeight() const
Returns skirt height (in world units). Skirts at the edges of terrain tiles help hide cracks between ...
QgsRectangle extent() const override
extent of the terrain in terrain&#39;s CRS
const QgsTilingScheme & tilingScheme() const
Returns tiling scheme of the terrain.
Contains information about the context in which a coordinate transform is executed.
Base class for processing filters like renderers, reprojector, resampler etc.
3 Takes care of downloading terrain data from a publicly available data source.
double terrainVerticalScale() const
Returns vertical scale (exaggeration) of terrain.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:177
QgsRasterProjector implements approximate projection support for it calculates grid of points in sour...
QgsCoordinateTransformContext transformContext() const
Returns data provider coordinate transform context.
QgsDemHeightMapGenerator * heightMapGenerator()
Returns height map generator object - takes care of extraction of elevations from the layer) ...
3 Implementation of terrain generator that uses a raster layer with DEM to build terrain.
This class represents a coordinate reference system (CRS).
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:167
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:172
Terrain is built from downloaded tiles with digital elevation model.
3 The class encapsulates tiling scheme (just like with WMTS / TMS / XYZ layers).
QgsRectangle tileToExtent(int x, int y, int z) const
Returns map coordinates of the extent of a tile.
double x() const
Returns X coordinate.
Definition: qgsvector3d.h:49
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:209
Base class for raster data providers.
QgsTerrainGenerator * terrainGenerator() const
Returns terrain generator. It takes care of producing terrain tiles from the input data...