QGIS API Documentation 3.28.0-Firenze (ed3ad0430f)
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
17
18#include "qgslogger.h"
19#include "qgsrasterlayer.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
65void 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
78double QgsTerrainDownloader::findBestTileResolution( double requestedMupp ) const
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
94void 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
120QByteArray 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:140
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.