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