1 /***************************************************************************
2 testqgsgdalutils.cpp
3 --------------------
4 begin : September 2018
5 copyright : (C) 2018 Even Rouault
6 email : even.rouault at spatialys.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 #include "qgstest.h"
16 #include <QObject>
17 #include <QString>
18 #include <QStringList>
19 #include <QSettings>
20
21 #include <gdal.h>
22
23 #include "qgsgdalutils.h"
24 #include "qgsapplication.h"
25 #include "qgsrasterlayer.h"
26
27
28 class TestQgsGdalUtils: public QObject
29 {
30 Q_OBJECT
31
32 private slots:
33 void initTestCase();// will be called before the first testfunction is executed.
34 void cleanupTestCase();// will be called after the last testfunction was executed.
35 void init();// will be called before each testfunction is executed.
36 void cleanup();// will be called after every testfunction.
37 void supportsRasterCreate();
38 void testCreateSingleBandMemoryDataset();
39 void testCreateMultiBandMemoryDataset();
40 void testCreateSingleBandTiffDataset();
41 void testResampleSingleBandRaster();
42 void testImageToDataset();
43 void testResampleImageToImage();
44
45 private:
46
47 double identify( GDALDatasetH dataset, int band, int px, int py );
48 };
49
initTestCase()50 void TestQgsGdalUtils::initTestCase()
51 {
52 QgsApplication::init();
53 QgsApplication::initQgis();
54 GDALAllRegister();
55 }
56
cleanupTestCase()57 void TestQgsGdalUtils::cleanupTestCase()
58 {
59 QgsApplication::exitQgis();
60 }
61
init()62 void TestQgsGdalUtils::init()
63 {
64
65 }
66
cleanup()67 void TestQgsGdalUtils::cleanup()
68 {
69
70 }
71
supportsRasterCreate()72 void TestQgsGdalUtils::supportsRasterCreate()
73 {
74 QVERIFY( QgsGdalUtils::supportsRasterCreate( GDALGetDriverByName( "GTiff" ) ) );
75 QVERIFY( QgsGdalUtils::supportsRasterCreate( GDALGetDriverByName( "GPKG" ) ) );
76
77 // special case
78 QVERIFY( !QgsGdalUtils::supportsRasterCreate( GDALGetDriverByName( "SQLite" ) ) );
79
80 // create-only
81 QVERIFY( !QgsGdalUtils::supportsRasterCreate( GDALGetDriverByName( "DTED" ) ) );
82
83 // vector-only
84 QVERIFY( !QgsGdalUtils::supportsRasterCreate( GDALGetDriverByName( "ESRI Shapefile" ) ) );
85 }
86
87 #if PROJ_VERSION_MAJOR>=6
88 #define EPSG_4326_WKT \
89 "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]]," \
90 "AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433," \
91 "AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"
92 #else
93 #define EPSG_4326_WKT \
94 "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]]," \
95 "AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433," \
96 "AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]"
97 #endif
98
testCreateSingleBandMemoryDataset()99 void TestQgsGdalUtils::testCreateSingleBandMemoryDataset()
100 {
101 gdal::dataset_unique_ptr ds1 = QgsGdalUtils::createSingleBandMemoryDataset( GDT_Float32, QgsRectangle( 1, 1, 21, 11 ), 40, 20, QgsCoordinateReferenceSystem( "EPSG:4326" ) );
102 QVERIFY( ds1 );
103
104 QCOMPARE( GDALGetRasterCount( ds1.get() ), 1 );
105 QCOMPARE( GDALGetRasterXSize( ds1.get() ), 40 );
106 QCOMPARE( GDALGetRasterYSize( ds1.get() ), 20 );
107
108 #if PROJ_VERSION_MAJOR>=6
109 QCOMPARE( GDALGetProjectionRef( ds1.get() ), R"""(GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]])""" );
110 #else
111 QCOMPARE( GDALGetProjectionRef( ds1.get() ), EPSG_4326_WKT );
112 #endif
113 double geoTransform[6];
114 double geoTransformExpected[] = { 1, 0.5, 0, 11, 0, -0.5 };
115 QCOMPARE( GDALGetGeoTransform( ds1.get(), geoTransform ), CE_None );
116 QVERIFY( memcmp( geoTransform, geoTransformExpected, sizeof( double ) * 6 ) == 0 );
117
118 QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( ds1.get(), 1 ) ), GDT_Float32 );
119 }
120
testCreateMultiBandMemoryDataset()121 void TestQgsGdalUtils::testCreateMultiBandMemoryDataset()
122 {
123 gdal::dataset_unique_ptr ds1 = QgsGdalUtils::createMultiBandMemoryDataset( GDT_Float32, 4, QgsRectangle( 1, 1, 21, 11 ), 40, 20, QgsCoordinateReferenceSystem( "EPSG:4326" ) );
124 QVERIFY( ds1 );
125
126 QCOMPARE( GDALGetRasterCount( ds1.get() ), 4 );
127 QCOMPARE( GDALGetRasterXSize( ds1.get() ), 40 );
128 QCOMPARE( GDALGetRasterYSize( ds1.get() ), 20 );
129
130 #if PROJ_VERSION_MAJOR>=6
131 QCOMPARE( GDALGetProjectionRef( ds1.get() ), R"""(GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]])""" );
132 #else
133 QCOMPARE( GDALGetProjectionRef( ds1.get() ), EPSG_4326_WKT );
134 #endif
135 double geoTransform[6];
136 double geoTransformExpected[] = { 1, 0.5, 0, 11, 0, -0.5 };
137 QCOMPARE( GDALGetGeoTransform( ds1.get(), geoTransform ), CE_None );
138 QVERIFY( memcmp( geoTransform, geoTransformExpected, sizeof( double ) * 6 ) == 0 );
139
140 QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( ds1.get(), 1 ) ), GDT_Float32 );
141 QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( ds1.get(), 2 ) ), GDT_Float32 );
142 QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( ds1.get(), 3 ) ), GDT_Float32 );
143 QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( ds1.get(), 4 ) ), GDT_Float32 );
144 }
145
testCreateSingleBandTiffDataset()146 void TestQgsGdalUtils::testCreateSingleBandTiffDataset()
147 {
148 QString filename = QDir::tempPath() + "/qgis_test_single_band_raster.tif";
149 QFile::remove( filename );
150 QVERIFY( !QFile::exists( filename ) );
151
152 gdal::dataset_unique_ptr ds1 = QgsGdalUtils::createSingleBandTiffDataset( filename, GDT_Float32, QgsRectangle( 1, 1, 21, 11 ), 40, 20, QgsCoordinateReferenceSystem( "EPSG:4326" ) );
153 QVERIFY( ds1 );
154
155 QCOMPARE( GDALGetRasterCount( ds1.get() ), 1 );
156 QCOMPARE( GDALGetRasterXSize( ds1.get() ), 40 );
157 QCOMPARE( GDALGetRasterYSize( ds1.get() ), 20 );
158
159 #if PROJ_VERSION_MAJOR>=6
160 QCOMPARE( GDALGetProjectionRef( ds1.get() ), R"""(GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]])""" );
161 #else
162 QCOMPARE( GDALGetProjectionRef( ds1.get() ), EPSG_4326_WKT );
163 #endif
164
165 double geoTransform[6];
166 double geoTransformExpected[] = { 1, 0.5, 0, 11, 0, -0.5 };
167 QCOMPARE( GDALGetGeoTransform( ds1.get(), geoTransform ), CE_None );
168 QVERIFY( memcmp( geoTransform, geoTransformExpected, sizeof( double ) * 6 ) == 0 );
169
170 QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( ds1.get(), 1 ) ), GDT_Float32 );
171
172 ds1.reset(); // makes sure the file is fully written
173
174 QVERIFY( QFile::exists( filename ) );
175
176 std::unique_ptr<QgsRasterLayer> layer( new QgsRasterLayer( filename, "test", "gdal" ) );
177 QVERIFY( layer->isValid() );
178 QCOMPARE( layer->extent(), QgsRectangle( 1, 1, 21, 11 ) );
179 QCOMPARE( layer->width(), 40 );
180 QCOMPARE( layer->height(), 20 );
181
182 layer.reset(); // let's clean up before removing the file
183 QFile::remove( filename );
184 }
185
testResampleSingleBandRaster()186 void TestQgsGdalUtils::testResampleSingleBandRaster()
187 {
188 QString inputFilename = QString( TEST_DATA_DIR ) + "/float1-16.tif";
189 gdal::dataset_unique_ptr srcDS( GDALOpen( inputFilename.toUtf8().constData(), GA_ReadOnly ) );
190 QVERIFY( srcDS );
191
192 QString outputFilename = QDir::tempPath() + "/qgis_test_float1-16_resampled.tif";
193 QgsRectangle outputExtent( 106.25, -6.75, 106.55, -6.45 );
194 gdal::dataset_unique_ptr dstDS = QgsGdalUtils::createSingleBandTiffDataset( outputFilename, GDT_Float32, outputExtent, 2, 2, QgsCoordinateReferenceSystem( "EPSG:4326" ) );
195 QVERIFY( dstDS );
196
197 QgsGdalUtils::resampleSingleBandRaster( srcDS.get(), dstDS.get(), GRA_NearestNeighbour, nullptr );
198 dstDS.reset();
199
200 std::unique_ptr<QgsRasterLayer> layer( new QgsRasterLayer( outputFilename, "test", "gdal" ) );
201 QVERIFY( layer );
202 std::unique_ptr<QgsRasterBlock> block( layer->dataProvider()->block( 1, outputExtent, 2, 2 ) );
203 QVERIFY( block );
204 QCOMPARE( block->value( 0, 0 ), 6. );
205 QCOMPARE( block->value( 1, 1 ), 11. );
206
207 layer.reset();
208 QFile::remove( outputFilename );
209 }
210
testImageToDataset()211 void TestQgsGdalUtils::testImageToDataset()
212 {
213 QString inputFilename = QString( TEST_DATA_DIR ) + "/rgb256x256.png";
214 QImage src = QImage( inputFilename );
215 src = src.convertToFormat( QImage::Format_ARGB32 );
216 QVERIFY( !src.isNull() );
217
218 gdal::dataset_unique_ptr dstDS = QgsGdalUtils::imageToMemoryDataset( QImage() );
219 QVERIFY( !dstDS );
220
221 dstDS = QgsGdalUtils::imageToMemoryDataset( src );
222 QVERIFY( dstDS );
223
224 QCOMPARE( GDALGetRasterCount( dstDS.get() ), 4 );
225 QCOMPARE( GDALGetRasterXSize( dstDS.get() ), 256 );
226 QCOMPARE( GDALGetRasterYSize( dstDS.get() ), 256 );
227
228 QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( dstDS.get(), 1 ) ), GDT_Byte );
229 QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( dstDS.get(), 2 ) ), GDT_Byte );
230 QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( dstDS.get(), 3 ) ), GDT_Byte );
231 QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( dstDS.get(), 4 ) ), GDT_Byte );
232
233 QCOMPARE( identify( dstDS.get(), 1, 50, 50 ), 255.0 );
234 QCOMPARE( identify( dstDS.get(), 1, 200, 50 ), 255.0 );
235 QCOMPARE( identify( dstDS.get(), 1, 50, 200 ), 0.0 );
236 QCOMPARE( identify( dstDS.get(), 1, 200, 200 ), 0.0 );
237
238 QCOMPARE( identify( dstDS.get(), 2, 50, 50 ), 255.0 );
239 QCOMPARE( identify( dstDS.get(), 2, 200, 50 ), 0.0 );
240 QCOMPARE( identify( dstDS.get(), 2, 50, 200 ), 255.0 );
241 QCOMPARE( identify( dstDS.get(), 2, 200, 200 ), 0.0 );
242
243 QCOMPARE( identify( dstDS.get(), 3, 50, 50 ), 0.0 );
244 QCOMPARE( identify( dstDS.get(), 3, 200, 50 ), 0.0 );
245 QCOMPARE( identify( dstDS.get(), 3, 50, 200 ), 0.0 );
246 QCOMPARE( identify( dstDS.get(), 3, 200, 200 ), 255.0 );
247
248 QCOMPARE( identify( dstDS.get(), 4, 50, 50 ), 255.0 );
249 QCOMPARE( identify( dstDS.get(), 4, 200, 50 ), 255.0 );
250 QCOMPARE( identify( dstDS.get(), 4, 50, 200 ), 255.0 );
251 QCOMPARE( identify( dstDS.get(), 4, 200, 200 ), 255.0 );
252 }
253
testResampleImageToImage()254 void TestQgsGdalUtils::testResampleImageToImage()
255 {
256 QString inputFilename = QString( TEST_DATA_DIR ) + "/rgb256x256.png";
257 QImage src = QImage( inputFilename );
258 src = src.convertToFormat( QImage::Format_ARGB32 );
259 QVERIFY( !src.isNull() );
260
261 QImage res = QgsGdalUtils::resampleImage( QImage(), QSize( 50, 50 ), GRIORA_NearestNeighbour );
262 QVERIFY( res.isNull() );
263
264 res = QgsGdalUtils::resampleImage( src, QSize( 50, 50 ), GRIORA_NearestNeighbour );
265 QVERIFY( !res.isNull() );
266 QCOMPARE( res.width(), 50 );
267 QCOMPARE( res.height(), 50 );
268
269 QCOMPARE( qRed( res.pixel( 10, 10 ) ), 255 );
270 QCOMPARE( qGreen( res.pixel( 10, 10 ) ), 255 );
271 QCOMPARE( qBlue( res.pixel( 10, 10 ) ), 0 );
272 QCOMPARE( qAlpha( res.pixel( 10, 10 ) ), 255 );
273
274 QCOMPARE( qRed( res.pixel( 40, 10 ) ), 255 );
275 QCOMPARE( qGreen( res.pixel( 40, 10 ) ), 0 );
276 QCOMPARE( qBlue( res.pixel( 40, 10 ) ), 0 );
277 QCOMPARE( qAlpha( res.pixel( 40, 10 ) ), 255 );
278
279 QCOMPARE( qRed( res.pixel( 10, 40 ) ), 0 );
280 QCOMPARE( qGreen( res.pixel( 10, 40 ) ), 255 );
281 QCOMPARE( qBlue( res.pixel( 10, 40 ) ), 0 );
282 QCOMPARE( qAlpha( res.pixel( 10, 40 ) ), 255 );
283
284 QCOMPARE( qRed( res.pixel( 40, 40 ) ), 0 );
285 QCOMPARE( qGreen( res.pixel( 40, 40 ) ), 0 );
286 QCOMPARE( qBlue( res.pixel( 40, 40 ) ), 255 );
287 QCOMPARE( qAlpha( res.pixel( 40, 40 ) ), 255 );
288 }
289
identify(GDALDatasetH dataset,int band,int px,int py)290 double TestQgsGdalUtils::identify( GDALDatasetH dataset, int band, int px, int py )
291 {
292 GDALRasterBandH hBand = GDALGetRasterBand( dataset, band );
293
294 float *pafScanline = ( float * ) CPLMalloc( sizeof( float ) );
295 CPLErr err = GDALRasterIO( hBand, GF_Read, px, py, 1, 1,
296 pafScanline, 1, 1, GDT_Float32, 0, 0 );
297 double value = err == CE_None ? pafScanline[0] : std::numeric_limits<double>::quiet_NaN();
298 CPLFree( pafScanline );
299
300 return value;
301 }
302
303 QGSTEST_MAIN( TestQgsGdalUtils )
304 #include "testqgsgdalutils.moc"
305