QGIS API Documentation  3.22.4-Białowieża (ce8e65e95e)
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  QgsCoordinateTransform ct( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3857" ) ), transformContext );
33  const QgsPointXY topLeftLonLat( -180, 180.0 / M_PI * std::atan( std::sinh( M_PI ) ) );
34  const QgsPointXY bottomRightLonLat( 180, 180.0 / M_PI * std::atan( std::sinh( -M_PI ) ) );
35  const QgsPointXY topLeft = ct.transform( topLeftLonLat );
36  const QgsPointXY bottomRight = ct.transform( bottomRightLonLat );
37  mXSpan = ( bottomRight.x() - topLeft.x() );
38 }
39 
41 
43 {
44  // using terrain tiles stored on AWS and listed within Registry of Open Data on AWS
45  // see https://registry.opendata.aws/terrain-tiles/
46  //
47  // tiles are generated using a variety of sources (SRTM, ETOPO1 and more detailed data for some countries)
48  // for more details and attribution see https://github.com/tilezen/joerd/blob/master/docs/data-sources.md
49 
50  DataSource ds;
51  ds.uri = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png";
52  ds.zMin = 0;
53  ds.zMax = 15;
54  return ds;
55 }
56 
58 {
59  mDataSource = ds;
60  const QString uri = QString( "type=xyz&url=%1&zmin=%2&zmax=%3" ).arg( mDataSource.uri ).arg( mDataSource.zMin ).arg( mDataSource.zMax );
61  mOnlineDtm.reset( new QgsRasterLayer( uri, "terrarium", "wms" ) );
62 }
63 
64 
65 void QgsTerrainDownloader::adjustExtentAndResolution( double mupp, const QgsRectangle &extentOrig, QgsRectangle &extent, int &res )
66 {
67  const double xMin = floor( extentOrig.xMinimum() / mupp ) * mupp;
68  const double xMax = ceil( extentOrig.xMaximum() / mupp ) * mupp;
69 
70  const double yMin = floor( extentOrig.yMinimum() / mupp ) * mupp;
71  const double yMax = ceil( extentOrig.yMaximum() / mupp ) * mupp;
72 
73  extent = QgsRectangle( xMin, yMin, xMax, yMax );
74  res = round( ( xMax - xMin ) / mupp );
75 }
76 
77 
78 double QgsTerrainDownloader::findBestTileResolution( double requestedMupp )
79 {
80  int zoom = 0;
81  for ( ; zoom <= 15; ++zoom )
82  {
83  const double tileMupp = mXSpan / ( 256 * ( 1 << zoom ) );
84  if ( tileMupp <= requestedMupp )
85  break;
86  }
87 
88  if ( zoom > 15 ) zoom = 15;
89  const double finalMupp = mXSpan / ( 256 * ( 1 << zoom ) );
90  return finalMupp;
91 }
92 
93 
94 void QgsTerrainDownloader::tileImageToHeightMap( const QImage &img, QByteArray &heightMap )
95 {
96  // for description of the "terrarium" format:
97  // https://github.com/tilezen/joerd/blob/master/docs/formats.md
98 
99  // assuming ARGB premultiplied but with alpha 255
100  const QRgb *rgb = reinterpret_cast<const QRgb *>( img.constBits() );
101  const int count = img.width() * img.height();
102  heightMap.resize( sizeof( float ) * count );
103  float *hData = reinterpret_cast<float *>( heightMap.data() );
104  for ( int i = 0; i < count; ++i )
105  {
106  const QRgb c = rgb[i];
107  if ( qAlpha( c ) == 255 )
108  {
109  const float h = qRed( c ) * 256 + qGreen( c ) + qBlue( c ) / 256.f - 32768;
110  *hData++ = h;
111  }
112  else
113  {
114  *hData++ = std::numeric_limits<float>::quiet_NaN();
115  }
116  }
117 }
118 
119 
120 QByteArray QgsTerrainDownloader::getHeightMap( const QgsRectangle &extentOrig, int res, const QgsCoordinateReferenceSystem &destCrs, const QgsCoordinateTransformContext &context, QString tmpFilenameImg, QString tmpFilenameTif )
121 {
122  if ( !mOnlineDtm || !mOnlineDtm->isValid() )
123  {
124  QgsDebugMsg( "missing a valid data source" );
125  return QByteArray();
126  }
127 
128  QgsRectangle extentTr = extentOrig;
129  if ( destCrs != mOnlineDtm->crs() )
130  {
131  // if in different CRS - need to reproject extent and resolution
132  QgsCoordinateTransform ct( destCrs, mOnlineDtm->crs(), context );
134  extentTr = ct.transformBoundingBox( extentOrig );
135  }
136 
137  const double requestedMupp = extentTr.width() / res;
138  const double finalMupp = findBestTileResolution( requestedMupp );
139 
140  // adjust extent to match native resolution of terrain tiles
141 
142  QgsRectangle extent;
143  const int resOrig = res;
144  adjustExtentAndResolution( finalMupp, extentTr, extent, res );
145 
146  // request tile
147 
148  QgsRasterBlock *b = mOnlineDtm->dataProvider()->block( 1, extent, res, res );
149  const QImage img = b->image();
150  delete b;
151  if ( !tmpFilenameImg.isEmpty() )
152  img.save( tmpFilenameImg );
153 
154  // convert to height data
155 
156  QByteArray heightMap;
157  tileImageToHeightMap( img, heightMap );
158 
159  // prepare source/destination datasets for resampling
160 
161  const gdal::dataset_unique_ptr hSrcDS( QgsGdalUtils::createSingleBandMemoryDataset( GDT_Float32, extent, res, res, mOnlineDtm->crs() ) );
163  if ( !tmpFilenameTif.isEmpty() )
164  hDstDS = QgsGdalUtils::createSingleBandTiffDataset( tmpFilenameTif, GDT_Float32, extentOrig, resOrig, resOrig, destCrs );
165  else
166  hDstDS = QgsGdalUtils::createSingleBandMemoryDataset( GDT_Float32, extentOrig, resOrig, resOrig, destCrs );
167 
168  if ( !hSrcDS || !hDstDS )
169  {
170  QgsDebugMsg( "failed to create GDAL dataset for heightmap" );
171  return QByteArray();
172  }
173 
174  const CPLErr err = GDALRasterIO( GDALGetRasterBand( hSrcDS.get(), 1 ), GF_Write, 0, 0, res, res, heightMap.data(), res, res, GDT_Float32, 0, 0 );
175  if ( err != CE_None )
176  {
177  QgsDebugMsg( "failed to write heightmap data to GDAL dataset" );
178  return QByteArray();
179  }
180 
181  // resample to the desired extent + resolution
182  QgsGdalUtils::resampleSingleBandRaster( hSrcDS.get(), hDstDS.get(), GRA_Bilinear,
183  context.calculateCoordinateOperation( mOnlineDtm->crs(), destCrs ).toUtf8().constData() );
184 
185  QByteArray heightMapOut;
186  heightMapOut.resize( resOrig * resOrig * sizeof( float ) );
187  char *data = heightMapOut.data();
188 
189  // read the data back
190 
191  const CPLErr err2 = GDALRasterIO( GDALGetRasterBand( hDstDS.get(), 1 ), GF_Read, 0, 0, resOrig, resOrig, data, resOrig, resOrig, GDT_Float32, 0, 0 );
192  if ( err2 != CE_None )
193  {
194  QgsDebugMsg( "failed to read heightmap data from GDAL dataset" );
195  return QByteArray();
196  }
197 
198  return heightMapOut;
199 }
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.
void setBallparkTransformsAreAppropriate(bool appropriate)
Sets whether approximate "ballpark" results are appropriate for this coordinate transform.
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.