QGIS API Documentation  3.23.0-Master (eb871beae0)
qgsterraindownloader.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsterraindownloader.cpp
3  --------------------------------------
4  Date : March 2019
5  Copyright : (C) 2019 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 
16 #include "qgsterraindownloader.h"
17 
18 #include "qgslogger.h"
19 #include "qgsrasterlayer.h"
20 #include "qgscoordinatetransform.h"
21 #include "qgsgdalutils.h"
22 
23 
25 {
27 
28  // the whole world is projected to a square:
29  // X going from 180 W to 180 E
30  // Y going from ~85 N to ~85 S (=atan(sinh(pi)) ... to get a square)
31  const QgsCoordinateTransform ct( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3857" ) ), transformContext );
32  const QgsPointXY topLeftLonLat( -180, 180.0 / M_PI * std::atan( std::sinh( M_PI ) ) );
33  const QgsPointXY bottomRightLonLat( 180, 180.0 / M_PI * std::atan( std::sinh( -M_PI ) ) );
34  const QgsPointXY topLeft = ct.transform( topLeftLonLat );
35  const QgsPointXY bottomRight = ct.transform( bottomRightLonLat );
36  mXSpan = ( bottomRight.x() - topLeft.x() );
37 }
38 
40 
42 {
43  // using terrain tiles stored on AWS and listed within Registry of Open Data on AWS
44  // see https://registry.opendata.aws/terrain-tiles/
45  //
46  // tiles are generated using a variety of sources (SRTM, ETOPO1 and more detailed data for some countries)
47  // for more details and attribution see https://github.com/tilezen/joerd/blob/master/docs/data-sources.md
48 
49  DataSource ds;
50  ds.uri = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png";
51  ds.zMin = 0;
52  ds.zMax = 15;
53  return ds;
54 }
55 
57 {
58  mDataSource = ds;
59  const QString uri = QString( "type=xyz&url=%1&zmin=%2&zmax=%3" ).arg( mDataSource.uri ).arg( mDataSource.zMin ).arg( mDataSource.zMax );
60  mOnlineDtm.reset( new QgsRasterLayer( uri, "terrarium", "wms" ) );
61 }
62 
63 
64 void QgsTerrainDownloader::adjustExtentAndResolution( double mupp, const QgsRectangle &extentOrig, QgsRectangle &extent, int &res )
65 {
66  const double xMin = floor( extentOrig.xMinimum() / mupp ) * mupp;
67  const double xMax = ceil( extentOrig.xMaximum() / mupp ) * mupp;
68 
69  const double yMin = floor( extentOrig.yMinimum() / mupp ) * mupp;
70  const double yMax = ceil( extentOrig.yMaximum() / mupp ) * mupp;
71 
72  extent = QgsRectangle( xMin, yMin, xMax, yMax );
73  res = round( ( xMax - xMin ) / mupp );
74 }
75 
76 
77 double QgsTerrainDownloader::findBestTileResolution( double requestedMupp )
78 {
79  int zoom = 0;
80  for ( ; zoom <= 15; ++zoom )
81  {
82  const double tileMupp = mXSpan / ( 256 * ( 1 << zoom ) );
83  if ( tileMupp <= requestedMupp )
84  break;
85  }
86 
87  if ( zoom > 15 ) zoom = 15;
88  const double finalMupp = mXSpan / ( 256 * ( 1 << zoom ) );
89  return finalMupp;
90 }
91 
92 
93 void QgsTerrainDownloader::tileImageToHeightMap( const QImage &img, QByteArray &heightMap )
94 {
95  // for description of the "terrarium" format:
96  // https://github.com/tilezen/joerd/blob/master/docs/formats.md
97 
98  // assuming ARGB premultiplied but with alpha 255
99  const QRgb *rgb = reinterpret_cast<const QRgb *>( img.constBits() );
100  const int count = img.width() * img.height();
101  heightMap.resize( sizeof( float ) * count );
102  float *hData = reinterpret_cast<float *>( heightMap.data() );
103  for ( int i = 0; i < count; ++i )
104  {
105  const QRgb c = rgb[i];
106  if ( qAlpha( c ) == 255 )
107  {
108  const float h = qRed( c ) * 256 + qGreen( c ) + qBlue( c ) / 256.f - 32768;
109  *hData++ = h;
110  }
111  else
112  {
113  *hData++ = std::numeric_limits<float>::quiet_NaN();
114  }
115  }
116 }
117 
118 
119 QByteArray QgsTerrainDownloader::getHeightMap( const QgsRectangle &extentOrig, int res, const QgsCoordinateReferenceSystem &destCrs, const QgsCoordinateTransformContext &context, QString tmpFilenameImg, QString tmpFilenameTif )
120 {
121  if ( !mOnlineDtm || !mOnlineDtm->isValid() )
122  {
123  QgsDebugMsg( "missing a valid data source" );
124  return QByteArray();
125  }
126 
127  QgsRectangle extentTr = extentOrig;
128  if ( destCrs != mOnlineDtm->crs() )
129  {
130  // if in different CRS - need to reproject extent and resolution
131  const QgsCoordinateTransform ct( destCrs, mOnlineDtm->crs(), context );
132  extentTr = ct.transformBoundingBox( extentOrig );
133  }
134 
135  const double requestedMupp = extentTr.width() / res;
136  const double finalMupp = findBestTileResolution( requestedMupp );
137 
138  // adjust extent to match native resolution of terrain tiles
139 
140  QgsRectangle extent;
141  const int resOrig = res;
142  adjustExtentAndResolution( finalMupp, extentTr, extent, res );
143 
144  // request tile
145 
146  QgsRasterBlock *b = mOnlineDtm->dataProvider()->block( 1, extent, res, res );
147  const QImage img = b->image();
148  delete b;
149  if ( !tmpFilenameImg.isEmpty() )
150  img.save( tmpFilenameImg );
151 
152  // convert to height data
153 
154  QByteArray heightMap;
155  tileImageToHeightMap( img, heightMap );
156 
157  // prepare source/destination datasets for resampling
158 
159  const gdal::dataset_unique_ptr hSrcDS( QgsGdalUtils::createSingleBandMemoryDataset( GDT_Float32, extent, res, res, mOnlineDtm->crs() ) );
161  if ( !tmpFilenameTif.isEmpty() )
162  hDstDS = QgsGdalUtils::createSingleBandTiffDataset( tmpFilenameTif, GDT_Float32, extentOrig, resOrig, resOrig, destCrs );
163  else
164  hDstDS = QgsGdalUtils::createSingleBandMemoryDataset( GDT_Float32, extentOrig, resOrig, resOrig, destCrs );
165 
166  if ( !hSrcDS || !hDstDS )
167  {
168  QgsDebugMsg( "failed to create GDAL dataset for heightmap" );
169  return QByteArray();
170  }
171 
172  const CPLErr err = GDALRasterIO( GDALGetRasterBand( hSrcDS.get(), 1 ), GF_Write, 0, 0, res, res, heightMap.data(), res, res, GDT_Float32, 0, 0 );
173  if ( err != CE_None )
174  {
175  QgsDebugMsg( "failed to write heightmap data to GDAL dataset" );
176  return QByteArray();
177  }
178 
179  // resample to the desired extent + resolution
180  QgsGdalUtils::resampleSingleBandRaster( hSrcDS.get(), hDstDS.get(), GRA_Bilinear,
181  context.calculateCoordinateOperation( mOnlineDtm->crs(), destCrs ).toUtf8().constData() );
182 
183  QByteArray heightMapOut;
184  heightMapOut.resize( resOrig * resOrig * sizeof( float ) );
185  char *data = heightMapOut.data();
186 
187  // read the data back
188 
189  const CPLErr err2 = GDALRasterIO( GDALGetRasterBand( hDstDS.get(), 1 ), GF_Read, 0, 0, resOrig, resOrig, data, resOrig, resOrig, GDT_Float32, 0, 0 );
190  if ( err2 != CE_None )
191  {
192  QgsDebugMsg( "failed to read heightmap data from GDAL dataset" );
193  return QByteArray();
194  }
195 
196  return heightMapOut;
197 }
This class represents a coordinate reference system (CRS).
Contains information about the context in which a coordinate transform is executed.
QString calculateCoordinateOperation(const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination) const
Returns the Proj coordinate operation string to use when transforming from the specified source CRS t...
Class for doing transforms between two map coordinate systems.
QgsRectangle transformBoundingBox(const QgsRectangle &rectangle, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward, bool handle180Crossover=false) const SIP_THROW(QgsCsException)
Transforms a rectangle from the source CRS to the destination CRS.
QgsPointXY transform(const QgsPointXY &point, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward) const SIP_THROW(QgsCsException)
Transform the point from the source CRS to the destination CRS.
static bool resampleSingleBandRaster(GDALDatasetH hSrcDS, GDALDatasetH hDstDS, GDALResampleAlg resampleAlg, const char *pszCoordinateOperation)
Resamples a single band raster to the destination dataset with different resolution (and possibly wit...
static gdal::dataset_unique_ptr createSingleBandTiffDataset(const QString &filename, GDALDataType dataType, const QgsRectangle &extent, int width, int height, const QgsCoordinateReferenceSystem &crs)
Creates a new single band TIFF dataset with given parameters.
static gdal::dataset_unique_ptr createSingleBandMemoryDataset(GDALDataType dataType, const QgsRectangle &extent, int width, int height, const QgsCoordinateReferenceSystem &crs)
Creates a new single band memory dataset with given parameters.
A class to represent a 2D point.
Definition: qgspointxy.h:59
Q_GADGET double x
Definition: qgspointxy.h:62
Raster data container.
QImage image() const
Returns an image containing the block data, if the block's data type is color.
Represents a raster layer.
A rectangle specified with double values.
Definition: qgsrectangle.h:42
double yMaximum() const SIP_HOLDGIL
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:193
double xMaximum() const SIP_HOLDGIL
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:183
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
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:223
QgsTerrainDownloader(const QgsCoordinateTransformContext &transformContext)
Constructs a QgsTerrainDownloader object.
static DataSource defaultDataSource()
Returns the data source used by default.
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,...
void setDataSource(const DataSource &ds)
Configures data source to be used for download of terrain tiles.
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
Definition: qgsogrutils.h:138
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
Definition of data source for terrain tiles (assuming "terrarium" data encoding with usual XYZ tiling...
QString uri
HTTP(S) template for XYZ tiles requests (e.g. http://example.com/{z}/{x}/{y}.png)
int zMin
Minimum zoom level (Z) with valid data.
int zMax
Maximum zoom level (Z) with valid data.