Skip to content

Commit

Permalink
Add tests, fix issues from review, customizable terrain downloader
Browse files Browse the repository at this point in the history
  • Loading branch information
@wonder-sk
wonder-sk committed Mar 19, 2019
1 parent 407952f commit 3032a65
Show file tree
Hide file tree
Showing 5 changed files with 154 additions and 17 deletions.
10 changes: 4 additions & 6 deletions src/3d/terrain/qgsonlineterraingenerator.cpp
View file
Expand Up @@ -18,10 +18,9 @@
#include "qgsdemterraintileloader_p.h"


QgsOnlineTerrainGenerator::~QgsOnlineTerrainGenerator()
{
delete mHeightMapGenerator;
}
QgsOnlineTerrainGenerator::QgsOnlineTerrainGenerator() = default;

QgsOnlineTerrainGenerator::~QgsOnlineTerrainGenerator() = default;

QgsChunkLoader *QgsOnlineTerrainGenerator::createChunkLoader( QgsChunkNode *node ) const
{
Expand Down Expand Up @@ -114,6 +113,5 @@ void QgsOnlineTerrainGenerator::updateGenerator()
mTerrainTilingScheme = QgsTilingScheme( mExtent, mCrs );
}

delete mHeightMapGenerator;
mHeightMapGenerator = new QgsDemHeightMapGenerator( nullptr, mTerrainTilingScheme, mResolution );
mHeightMapGenerator.reset( new QgsDemHeightMapGenerator( nullptr, mTerrainTilingScheme, mResolution ) );
}
6 changes: 3 additions & 3 deletions src/3d/terrain/qgsonlineterraingenerator.h
View file
Expand Up @@ -33,7 +33,7 @@ class _3D_EXPORT QgsOnlineTerrainGenerator : public QgsTerrainGenerator
{
public:
//! Constructor for QgsOnlineTerrainGenerator
QgsOnlineTerrainGenerator() = default;
QgsOnlineTerrainGenerator();
~QgsOnlineTerrainGenerator() override;

//! Sets extent of the terrain
Expand All @@ -55,7 +55,7 @@ class _3D_EXPORT QgsOnlineTerrainGenerator : public QgsTerrainGenerator
float skirtHeight() const { return mSkirtHeight; }

//! Returns height map generator object - takes care of extraction of elevations from the layer)
QgsDemHeightMapGenerator *heightMapGenerator() { return mHeightMapGenerator; }
QgsDemHeightMapGenerator *heightMapGenerator() { return mHeightMapGenerator.get(); }

QgsTerrainGenerator *clone() const override SIP_FACTORY;
Type type() const override;
Expand All @@ -80,7 +80,7 @@ class _3D_EXPORT QgsOnlineTerrainGenerator : public QgsTerrainGenerator
//! height of the "skirts" at the edges of tiles to hide cracks between adjacent cracks
float mSkirtHeight = 10.f;

QgsDemHeightMapGenerator *mHeightMapGenerator = nullptr;
std::unique_ptr<QgsDemHeightMapGenerator> mHeightMapGenerator;
};

#endif // QGSONLINETERRAINGENERATOR_H
42 changes: 36 additions & 6 deletions src/3d/terrain/qgsterraindownloader.cpp
View file
Expand Up @@ -23,8 +23,7 @@

QgsTerrainDownloader::QgsTerrainDownloader()
{
QString uri = "type=xyz&url=http://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png&zmax=15&zmin=0";
onlineDtm.reset( new QgsRasterLayer( uri, "terrarium", "wms" ) );
setDataSource( defaultDataSource() );

// the whole world is projected to a square:
// X going from 180 W to 180 E
Expand All @@ -41,6 +40,28 @@ QgsTerrainDownloader::QgsTerrainDownloader()

QgsTerrainDownloader::~QgsTerrainDownloader() = default;

QgsTerrainDownloader::DataSource QgsTerrainDownloader::defaultDataSource()
{
// using terrain tiles stored on AWS and listed within Registry of Open Data on AWS
// see https://registry.opendata.aws/terrain-tiles/
//
// tiles are generated using a variety of sources (SRTM, ETOPO1 and more detailed data for some countries)
// for more details and attribution see https://github.com/tilezen/joerd/blob/master/docs/data-sources.md

DataSource ds;
ds.uri = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png";
ds.zMin = 0;
ds.zMax = 15;
return ds;
}

void QgsTerrainDownloader::setDataSource( const QgsTerrainDownloader::DataSource &ds )
{
mDataSource = ds;
QString uri = QString( "type=xyz&url=%1&zmin=%2&zmax=%3" ).arg( mDataSource.uri ).arg( mDataSource.zMin ).arg( mDataSource.zMax );
mOnlineDtm.reset( new QgsRasterLayer( uri, "terrarium", "wms" ) );
}


void QgsTerrainDownloader::adjustExtentAndResolution( double mupp, const QgsRectangle &extentOrig, QgsRectangle &extent, int &res )
{
Expand Down Expand Up @@ -73,6 +94,9 @@ double QgsTerrainDownloader::findBestTileResolution( double requestedMupp )

void QgsTerrainDownloader::tileImageToHeightMap( const QImage &img, QByteArray &heightMap )
{
// for description of the "terrarium" format:
// https://github.com/tilezen/joerd/blob/master/docs/formats.md

// assuming ARGB premultiplied but with alpha 255
const QRgb *rgb = reinterpret_cast<const QRgb *>( img.constBits() );
int count = img.width() * img.height();
Expand All @@ -96,11 +120,17 @@ void QgsTerrainDownloader::tileImageToHeightMap( const QImage &img, QByteArray &

QByteArray QgsTerrainDownloader::getHeightMap( const QgsRectangle &extentOrig, int res, const QgsCoordinateReferenceSystem &destCrs, const QgsCoordinateTransformContext &context, QString tmpFilenameImg, QString tmpFilenameTif )
{
if ( !mOnlineDtm || !mOnlineDtm->isValid() )
{
QgsDebugMsg( "missing a valid data source" );
return QByteArray();
}

QgsRectangle extentTr = extentOrig;
if ( destCrs != onlineDtm->crs() )
if ( destCrs != mOnlineDtm->crs() )
{
// if in different CRS - need to reproject extent and resolution
QgsCoordinateTransform ct( destCrs, onlineDtm->crs(), context );
QgsCoordinateTransform ct( destCrs, mOnlineDtm->crs(), context );
extentTr = ct.transformBoundingBox( extentOrig );
}

Expand All @@ -115,7 +145,7 @@ QByteArray QgsTerrainDownloader::getHeightMap( const QgsRectangle &extentOrig, i

// request tile

QgsRasterBlock *b = onlineDtm->dataProvider()->block( 1, extent, res, res );
QgsRasterBlock *b = mOnlineDtm->dataProvider()->block( 1, extent, res, res );
QImage img = b->image();
delete b;
if ( !tmpFilenameImg.isEmpty() )
Expand All @@ -128,7 +158,7 @@ QByteArray QgsTerrainDownloader::getHeightMap( const QgsRectangle &extentOrig, i

// prepare source/destination datasets for resampling

gdal::dataset_unique_ptr hSrcDS( QgsGdalUtils::createSingleBandMemoryDataset( GDT_Float32, extent, res, res, onlineDtm->crs() ) );
gdal::dataset_unique_ptr hSrcDS( QgsGdalUtils::createSingleBandMemoryDataset( GDT_Float32, extent, res, res, mOnlineDtm->crs() ) );
gdal::dataset_unique_ptr hDstDS;
if ( !tmpFilenameTif.isEmpty() )
hDstDS = QgsGdalUtils::createSingleBandTiffDataset( tmpFilenameTif, GDT_Float32, extentOrig, resOrig, resOrig, destCrs );
Expand Down
25 changes: 23 additions & 2 deletions src/3d/terrain/qgsterraindownloader.h
View file
Expand Up @@ -32,7 +32,10 @@ class QgsRasterLayer;
* \ingroup 3d
* Takes care of downloading terrain data from a publicly available data source.
*
* Currently using terrain tiles in GeoTIFF format hosted on AWS.
* Currently using terrain tiles in Terrarium format hosted on AWS. More info:
* - data format: https://github.com/tilezen/joerd/blob/master/docs/formats.md
* - data sources: https://github.com/tilezen/joerd/blob/master/docs/data-sources.md
* - hosting: https://registry.opendata.aws/terrain-tiles/
*
* \since QGIS 3.8
*/
Expand All @@ -43,6 +46,23 @@ class _3D_EXPORT QgsTerrainDownloader
QgsTerrainDownloader();
~QgsTerrainDownloader();

//! Definition of data source for terrain tiles (assuming "terrarium" data encoding with usual XYZ tiling scheme)
typedef struct
{
QString uri; //!< HTTP(S) template for XYZ tiles requests (e.g. http://example.com/{z}/{x}/{y}.png)
int zMin = 0; //!< Minimum zoom level (Z) with valid data
int zMax = 0; //!< Maximum zoom level (Z) with valid data
} DataSource;

//! Returns the data source used by default
static DataSource defaultDataSource();

//! Configures data source to be used for download of terrain tiles
void setDataSource( const DataSource &ds );

//! Returns currently configured data source
DataSource dataSource() const { return mDataSource; }

/**
* 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)
Expand All @@ -69,7 +89,8 @@ class _3D_EXPORT QgsTerrainDownloader
static void tileImageToHeightMap( const QImage &img, QByteArray &heightMap );

private:
std::unique_ptr<QgsRasterLayer> onlineDtm;
DataSource mDataSource;
std::unique_ptr<QgsRasterLayer> mOnlineDtm;
double mXSpan = 0; //!< Width of the tile at zoom level 0 in map units
};

Expand Down
88 changes: 88 additions & 0 deletions tests/src/core/testqgsgdalutils.cpp
View file
Expand Up @@ -22,6 +22,8 @@

#include "qgsgdalutils.h"
#include "qgsapplication.h"
#include "qgsrasterlayer.h"


class TestQgsGdalUtils: public QObject
{
Expand All @@ -33,6 +35,9 @@ class TestQgsGdalUtils: public QObject
void init();// will be called before each testfunction is executed.
void cleanup();// will be called after every testfunction.
void supportsRasterCreate();
void testCreateSingleBandMemoryDataset();
void testCreateSingleBandTiffDataset();
void testResampleSingleBandRaster();

private:
};
Expand Down Expand Up @@ -74,5 +79,88 @@ void TestQgsGdalUtils::supportsRasterCreate()
QVERIFY( !QgsGdalUtils::supportsRasterCreate( GDALGetDriverByName( "ESRI Shapefile" ) ) );
}

#define EPSG_4326_WKT \
"GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]]," \
"AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433," \
"AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]"

void TestQgsGdalUtils::testCreateSingleBandMemoryDataset()
{
gdal::dataset_unique_ptr ds1 = QgsGdalUtils::createSingleBandMemoryDataset( GDT_Float32, QgsRectangle( 1, 1, 21, 11 ), 40, 20, QgsCoordinateReferenceSystem( "EPSG:4326" ) );
QVERIFY( ds1 );

QCOMPARE( GDALGetRasterCount( ds1.get() ), 1 );
QCOMPARE( GDALGetRasterXSize( ds1.get() ), 40 );
QCOMPARE( GDALGetRasterYSize( ds1.get() ), 20 );

QCOMPARE( GDALGetProjectionRef( ds1.get() ), EPSG_4326_WKT );
double geoTransform[6];
double geoTransformExpected[] = { 1, 0.5, 0, 11, 0, -0.5 };
QCOMPARE( GDALGetGeoTransform( ds1.get(), geoTransform ), CE_None );
QVERIFY( memcmp( geoTransform, geoTransformExpected, sizeof( double ) * 6 ) == 0 );

QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( ds1.get(), 1 ) ), GDT_Float32 );
}

void TestQgsGdalUtils::testCreateSingleBandTiffDataset()
{
QString filename = QDir::tempPath() + "/qgis_test_single_band_raster.tif";
QFile::remove( filename );
QVERIFY( !QFile::exists( filename ) );

gdal::dataset_unique_ptr ds1 = QgsGdalUtils::createSingleBandTiffDataset( filename, GDT_Float32, QgsRectangle( 1, 1, 21, 11 ), 40, 20, QgsCoordinateReferenceSystem( "EPSG:4326" ) );
QVERIFY( ds1 );

QCOMPARE( GDALGetRasterCount( ds1.get() ), 1 );
QCOMPARE( GDALGetRasterXSize( ds1.get() ), 40 );
QCOMPARE( GDALGetRasterYSize( ds1.get() ), 20 );

QCOMPARE( GDALGetProjectionRef( ds1.get() ), EPSG_4326_WKT );
double geoTransform[6];
double geoTransformExpected[] = { 1, 0.5, 0, 11, 0, -0.5 };
QCOMPARE( GDALGetGeoTransform( ds1.get(), geoTransform ), CE_None );
QVERIFY( memcmp( geoTransform, geoTransformExpected, sizeof( double ) * 6 ) == 0 );

QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( ds1.get(), 1 ) ), GDT_Float32 );

ds1.reset(); // makes sure the file is fully written

QVERIFY( QFile::exists( filename ) );

std::unique_ptr<QgsRasterLayer> layer( new QgsRasterLayer( filename, "test", "gdal" ) );
QVERIFY( layer->isValid() );
QCOMPARE( layer->extent(), QgsRectangle( 1, 1, 21, 11 ) );
QCOMPARE( layer->width(), 40 );
QCOMPARE( layer->height(), 20 );

layer.reset(); // let's clean up before removing the file
QFile::remove( filename );
}

void TestQgsGdalUtils::testResampleSingleBandRaster()
{
QString inputFilename = QString( TEST_DATA_DIR ) + "/float1-16.tif";
gdal::dataset_unique_ptr srcDS( GDALOpen( inputFilename.toUtf8().constData(), GA_ReadOnly ) );
QVERIFY( srcDS );

QString outputFilename = QDir::tempPath() + "/qgis_test_float1-16_resampled.tif";
QgsRectangle outputExtent( 106.25, -6.75, 106.55, -6.45 );
gdal::dataset_unique_ptr dstDS = QgsGdalUtils::createSingleBandTiffDataset( outputFilename, GDT_Float32, outputExtent, 2, 2, QgsCoordinateReferenceSystem( "EPSG:4326" ) );
QVERIFY( dstDS );

QgsGdalUtils::resampleSingleBandRaster( srcDS.get(), dstDS.get(), GRA_NearestNeighbour );
dstDS.reset();

std::unique_ptr<QgsRasterLayer> layer( new QgsRasterLayer( outputFilename, "test", "gdal" ) );
QVERIFY( layer );
std::unique_ptr<QgsRasterBlock> block( layer->dataProvider()->block( 1, outputExtent, 2, 2 ) );
QVERIFY( block );
QCOMPARE( block->value( 0, 0 ), 6. );
QCOMPARE( block->value( 1, 1 ), 11. );

layer.reset();
QFile::remove( outputFilename );
}

QGSTEST_MAIN( TestQgsGdalUtils )
#include "testqgsgdalutils.moc"

0 comments on commit 3032a65

Please sign in to comment.