QGIS API Documentation 3.40.0-Bratislava (b56115d8743)
Loading...
Searching...
No Matches
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 "qgs3dutils.h"
19#include "qgslogger.h"
20#include "qgsrasterlayer.h"
22#include "qgsgdalutils.h"
23
24
26{
28
29 // the whole world is projected to a square:
30 // X going from 180 W to 180 E
31 // Y going from ~85 N to ~85 S (=atan(sinh(pi)) ... to get a square)
32 QgsCoordinateTransform ct( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3857" ) ), transformContext );
34 const QgsPointXY topLeftLonLat( -180, 180.0 / M_PI * std::atan( std::sinh( M_PI ) ) );
35 const QgsPointXY bottomRightLonLat( 180, 180.0 / M_PI * std::atan( std::sinh( -M_PI ) ) );
36 const QgsPointXY topLeft = ct.transform( topLeftLonLat );
37 const QgsPointXY bottomRight = ct.transform( bottomRightLonLat );
38 mXSpan = ( bottomRight.x() - topLeft.x() );
39}
40
42
44{
45 // using terrain tiles stored on AWS and listed within Registry of Open Data on AWS
46 // see https://registry.opendata.aws/terrain-tiles/
47 //
48 // tiles are generated using a variety of sources (SRTM, ETOPO1 and more detailed data for some countries)
49 // for more details and attribution see https://github.com/tilezen/joerd/blob/master/docs/data-sources.md
50
51 DataSource ds;
52 ds.uri = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png";
53 ds.zMin = 0;
54 ds.zMax = 15;
55 return ds;
56}
57
59{
60 mDataSource = ds;
61 const QString uri = QString( "type=xyz&url=%1&zmin=%2&zmax=%3" ).arg( mDataSource.uri ).arg( mDataSource.zMin ).arg( mDataSource.zMax );
62 mOnlineDtm.reset( new QgsRasterLayer( uri, "terrarium", "wms" ) );
63}
64
65
66void QgsTerrainDownloader::adjustExtentAndResolution( double mupp, const QgsRectangle &extentOrig, QgsRectangle &extent, int &res )
67{
68 const double xMin = floor( extentOrig.xMinimum() / mupp ) * mupp;
69 const double xMax = ceil( extentOrig.xMaximum() / mupp ) * mupp;
70
71 const double yMin = floor( extentOrig.yMinimum() / mupp ) * mupp;
72 const double yMax = ceil( extentOrig.yMaximum() / mupp ) * mupp;
73
74 extent = QgsRectangle( xMin, yMin, xMax, yMax );
75 res = round( ( xMax - xMin ) / mupp );
76}
77
78
79double QgsTerrainDownloader::findBestTileResolution( double requestedMupp ) const
80{
81 int zoom = 0;
82 for ( ; zoom <= 15; ++zoom )
83 {
84 const double tileMupp = mXSpan / ( 256 * ( 1 << zoom ) );
85 if ( tileMupp <= requestedMupp )
86 break;
87 }
88
89 if ( zoom > 15 ) zoom = 15;
90 const double finalMupp = mXSpan / ( 256 * ( 1 << zoom ) );
91 return finalMupp;
92}
93
94
95void QgsTerrainDownloader::tileImageToHeightMap( const QImage &img, QByteArray &heightMap )
96{
97 // for description of the "terrarium" format:
98 // https://github.com/tilezen/joerd/blob/master/docs/formats.md
99
100 // assuming ARGB premultiplied but with alpha 255
101 const QRgb *rgb = reinterpret_cast<const QRgb *>( img.constBits() );
102 const int count = img.width() * img.height();
103 heightMap.resize( sizeof( float ) * count );
104 float *hData = reinterpret_cast<float *>( heightMap.data() );
105 for ( int i = 0; i < count; ++i )
106 {
107 const QRgb c = rgb[i];
108 if ( qAlpha( c ) == 255 )
109 {
110 const float h = qRed( c ) * 256 + qGreen( c ) + qBlue( c ) / 256.f - 32768;
111 *hData++ = h;
112 }
113 else
114 {
115 *hData++ = std::numeric_limits<float>::quiet_NaN();
116 }
117 }
118}
119
120
121QByteArray QgsTerrainDownloader::getHeightMap( const QgsRectangle &extentOrig, int res, const QgsCoordinateReferenceSystem &destCrs, const QgsCoordinateTransformContext &context, QString tmpFilenameImg, QString tmpFilenameTif )
122{
123 if ( !mOnlineDtm || !mOnlineDtm->isValid() )
124 {
125 QgsDebugError( "missing a valid data source" );
126 return QByteArray();
127 }
128
129 QgsRectangle extentTr = Qgs3DUtils::tryReprojectExtent2D( extentOrig, destCrs, mOnlineDtm->crs(), context );
130 const double requestedMupp = extentTr.width() / res;
131 const double finalMupp = findBestTileResolution( requestedMupp );
132
133 // adjust extent to match native resolution of terrain tiles
134
135 QgsRectangle extent;
136 const int resOrig = res;
137 adjustExtentAndResolution( finalMupp, extentTr, extent, res );
138
139 // request tile
140
141 QgsRasterBlock *b = mOnlineDtm->dataProvider()->block( 1, extent, res, res );
142 const QImage img = b->image();
143 delete b;
144 if ( !tmpFilenameImg.isEmpty() )
145 img.save( tmpFilenameImg );
146
147 // convert to height data
148
149 QByteArray heightMap;
150 tileImageToHeightMap( img, heightMap );
151
152 // prepare source/destination datasets for resampling
153
154 const gdal::dataset_unique_ptr hSrcDS( QgsGdalUtils::createSingleBandMemoryDataset( GDT_Float32, extent, res, res, mOnlineDtm->crs() ) );
156 if ( !tmpFilenameTif.isEmpty() )
157 hDstDS = QgsGdalUtils::createSingleBandTiffDataset( tmpFilenameTif, GDT_Float32, extentOrig, resOrig, resOrig, destCrs );
158 else
159 hDstDS = QgsGdalUtils::createSingleBandMemoryDataset( GDT_Float32, extentOrig, resOrig, resOrig, destCrs );
160
161 if ( !hSrcDS || !hDstDS )
162 {
163 QgsDebugError( "failed to create GDAL dataset for heightmap" );
164 return QByteArray();
165 }
166
167 const CPLErr err = GDALRasterIO( GDALGetRasterBand( hSrcDS.get(), 1 ), GF_Write, 0, 0, res, res, heightMap.data(), res, res, GDT_Float32, 0, 0 );
168 if ( err != CE_None )
169 {
170 QgsDebugError( "failed to write heightmap data to GDAL dataset" );
171 return QByteArray();
172 }
173
174 // resample to the desired extent + resolution
175 QgsGdalUtils::resampleSingleBandRaster( hSrcDS.get(), hDstDS.get(), GRA_Bilinear,
176 context.calculateCoordinateOperation( mOnlineDtm->crs(), destCrs ).toUtf8().constData() );
177
178 QByteArray heightMapOut;
179 heightMapOut.resize( resOrig * resOrig * sizeof( float ) );
180 char *data = heightMapOut.data();
181
182 // read the data back
183
184 const CPLErr err2 = GDALRasterIO( GDALGetRasterBand( hDstDS.get(), 1 ), GF_Read, 0, 0, resOrig, resOrig, data, resOrig, resOrig, GDT_Float32, 0, 0 );
185 if ( err2 != CE_None )
186 {
187 QgsDebugError( "failed to read heightmap data from GDAL dataset" );
188 return QByteArray();
189 }
190
191 return heightMapOut;
192}
static QgsRectangle tryReprojectExtent2D(const QgsRectangle &extent, const QgsCoordinateReferenceSystem &crs1, const QgsCoordinateReferenceSystem &crs2, const QgsCoordinateTransformContext &context)
Reprojects extent from crs1 to crs2 coordinate reference system with context context.
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.
QgsPointXY transform(const QgsPointXY &point, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward) const
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:60
double x
Definition qgspointxy.h:63
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.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
double width() const
Returns the width of the rectangle.
double xMaximum() const
Returns the x maximum value (right side of rectangle).
double yMaximum() const
Returns the y maximum value (top side of rectangle).
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.
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 QgsDebugError(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.