Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsgdalprovider.cpp - QGIS Data provider for
3 : : GDAL rasters
4 : : -------------------
5 : : begin : November, 2010
6 : : copyright : (C) 2010 by Radim Blazek
7 : : email : radim dot blazek at gmail dot com
8 : : ***************************************************************************/
9 : :
10 : : /***************************************************************************
11 : : * *
12 : : * This program is free software; you can redistribute it and/or modify *
13 : : * it under the terms of the GNU General Public License as published by *
14 : : * the Free Software Foundation; either version 2 of the License, or *
15 : : * (at your option) any later version. *
16 : : * *
17 : : ***************************************************************************/
18 : :
19 : : #include "qgsgdalprovider.h"
20 : : ///@cond PRIVATE
21 : :
22 : : #include "qgslogger.h"
23 : : #include "qgsgdalproviderbase.h"
24 : : #include "qgsgdalprovider.h"
25 : : #include "qgsconfig.h"
26 : :
27 : : #include "qgsgdalutils.h"
28 : : #include "qgsapplication.h"
29 : : #include "qgsauthmanager.h"
30 : : #include "qgscoordinatetransform.h"
31 : : #include "qgsdataitem.h"
32 : : #include "qgsdataitemprovider.h"
33 : : #include "qgsdatasourceuri.h"
34 : : #include "qgsgdaldataitems.h"
35 : : #include "qgshtmlutils.h"
36 : : #include "qgsmessagelog.h"
37 : : #include "qgsrectangle.h"
38 : : #include "qgscoordinatereferencesystem.h"
39 : : #include "qgsrasterbandstats.h"
40 : : #include "qgsrasteridentifyresult.h"
41 : : #include "qgsrasterlayer.h"
42 : : #include "qgsrasterpyramid.h"
43 : : #include "qgspointxy.h"
44 : : #include "qgssettings.h"
45 : : #include "qgsogrutils.h"
46 : : #include "qgsruntimeprofiler.h"
47 : :
48 : : #include <QImage>
49 : : #include <QColor>
50 : : #include <QProcess>
51 : : #include <QMessageBox>
52 : : #include <QDir>
53 : : #include <QFileInfo>
54 : : #include <QFile>
55 : : #include <QHash>
56 : : #include <QTime>
57 : : #include <QTextDocument>
58 : : #include <QDebug>
59 : : #include <QRegularExpression>
60 : :
61 : : #include <gdalwarper.h>
62 : : #include <gdal.h>
63 : : #include <ogr_srs_api.h>
64 : : #include <cpl_conv.h>
65 : : #include <cpl_string.h>
66 : :
67 : : #define ERRMSG(message) QGS_ERROR_MESSAGE(message,"GDAL provider")
68 : : #define ERR(message) QgsError(message,"GDAL provider")
69 : :
70 : : #define PROVIDER_KEY QStringLiteral( "gdal" )
71 : : #define PROVIDER_DESCRIPTION QStringLiteral( "GDAL data provider" )
72 : :
73 : : // To avoid potential races when destroying related instances ("main" and clones)
74 : : #if QT_VERSION < QT_VERSION_CHECK(5, 14, 0)
75 : : Q_GLOBAL_STATIC_WITH_ARGS( QMutex, sGdalProviderMutex, ( QMutex::Recursive ) )
76 : : #else
77 : 0 : Q_GLOBAL_STATIC( QRecursiveMutex, sGdalProviderMutex )
78 : : #endif
79 : :
80 : 5 : QHash< QgsGdalProvider *, QVector<QgsGdalProvider::DatasetPair> > QgsGdalProvider::mgDatasetCache;
81 : :
82 : : int QgsGdalProvider::mgDatasetCacheSize = 0;
83 : :
84 : : // Number of cached datasets from which we will try to do eviction when a
85 : : // provider has 2 or more cached datasets
86 : : const int MIN_THRESHOLD_FOR_CACHE_CLEANUP = 10;
87 : :
88 : : // Maximum number of cached datasets
89 : : // We try to keep at least 1 cached dataset per parent provider between
90 : : // MIN_THRESHOLD_FOR_CACHE_CLEANUP and MAX_CACHE_SIZE. But we don't want to
91 : : // maintain more than MAX_CACHE_SIZE datasets opened to avoid running short of
92 : : // file descriptors.
93 : : const int MAX_CACHE_SIZE = 50;
94 : :
95 : 0 : struct QgsGdalProgress
96 : : {
97 : : int type;
98 : 0 : QgsGdalProvider *provider = nullptr;
99 : 0 : QgsRasterBlockFeedback *feedback = nullptr;
100 : : };
101 : : //
102 : : // global callback function
103 : : //
104 : 0 : int CPL_STDCALL progressCallback( double dfComplete,
105 : : const char *pszMessage,
106 : : void *pProgressArg )
107 : : {
108 : : Q_UNUSED( pszMessage )
109 : :
110 : : static double sDfLastComplete = -1.0;
111 : :
112 : 0 : QgsGdalProgress *prog = static_cast<QgsGdalProgress *>( pProgressArg );
113 : :
114 : 0 : if ( sDfLastComplete > dfComplete )
115 : : {
116 : 0 : if ( sDfLastComplete >= 1.0 )
117 : 0 : sDfLastComplete = -1.0;
118 : : else
119 : 0 : sDfLastComplete = dfComplete;
120 : 0 : }
121 : :
122 : 0 : if ( std::floor( sDfLastComplete * 10 ) != std::floor( dfComplete * 10 ) )
123 : : {
124 : 0 : if ( prog->feedback )
125 : 0 : prog->feedback->setProgress( dfComplete * 100 );
126 : 0 : }
127 : 0 : sDfLastComplete = dfComplete;
128 : :
129 : 0 : if ( prog->feedback && prog->feedback->isCanceled() )
130 : 0 : return false;
131 : :
132 : 0 : return true;
133 : 0 : }
134 : :
135 : 0 : QgsGdalProvider::QgsGdalProvider( const QString &uri, const QgsError &error )
136 : 0 : : QgsRasterDataProvider( uri, QgsDataProvider::ProviderOptions() )
137 : 0 : , mpRefCounter( new QAtomicInt( 1 ) )
138 : 0 : , mpLightRefCounter( new QAtomicInt( 1 ) )
139 : 0 : , mUpdate( false )
140 : 0 : {
141 : 0 : mGeoTransform[0] = 0;
142 : 0 : mGeoTransform[1] = 1;
143 : 0 : mGeoTransform[2] = 0;
144 : 0 : mGeoTransform[3] = 0;
145 : 0 : mGeoTransform[4] = 0;
146 : 0 : mGeoTransform[5] = -1;
147 : 0 : setError( error );
148 : 0 : }
149 : :
150 : 0 : QgsGdalProvider::QgsGdalProvider( const QString &uri, const ProviderOptions &options, bool update, GDALDatasetH dataset )
151 : 0 : : QgsRasterDataProvider( uri, options )
152 : 0 : , mpRefCounter( new QAtomicInt( 1 ) )
153 : : #if QT_VERSION < QT_VERSION_CHECK(5, 14, 0)
154 : : , mpMutex( new QMutex( QMutex::Recursive ) )
155 : : #else
156 : 0 : , mpMutex( new QRecursiveMutex() )
157 : : #endif
158 : 0 : , mpParent( new QgsGdalProvider * ( this ) )
159 : 0 : , mpLightRefCounter( new QAtomicInt( 1 ) )
160 : 0 : , mUpdate( update )
161 : 0 : {
162 : 0 : mGeoTransform[0] = 0;
163 : 0 : mGeoTransform[1] = 1;
164 : 0 : mGeoTransform[2] = 0;
165 : 0 : mGeoTransform[3] = 0;
166 : 0 : mGeoTransform[4] = 0;
167 : 0 : mGeoTransform[5] = -1;
168 : :
169 : 0 : QgsDebugMsgLevel( "constructing with uri '" + uri + "'.", 2 );
170 : :
171 : 0 : QgsGdalProviderBase::registerGdalDrivers();
172 : :
173 : : #ifndef QT_NO_NETWORKPROXY
174 : 0 : QgsGdalUtils::setupProxy();
175 : : #endif
176 : :
177 : 0 : std::unique_ptr< QgsScopedRuntimeProfile > profile;
178 : 0 : if ( QgsApplication::profiler()->groupIsActive( QStringLiteral( "projectload" ) ) )
179 : 0 : profile = std::make_unique< QgsScopedRuntimeProfile >( tr( "Open data source" ), QStringLiteral( "projectload" ) );
180 : :
181 : 0 : if ( !CPLGetConfigOption( "AAIGRID_DATATYPE", nullptr ) )
182 : : {
183 : : // GDAL tends to open AAIGrid as Float32 which results in lost precision
184 : : // and confusing values shown to users, force Float64
185 : 0 : CPLSetConfigOption( "AAIGRID_DATATYPE", "Float64" );
186 : 0 : }
187 : :
188 : : #if !(GDAL_VERSION_MAJOR > 2 || (GDAL_VERSION_MAJOR == 2 && GDAL_VERSION_MINOR >= 3))
189 : : if ( !CPLGetConfigOption( "VRT_SHARED_SOURCE", nullptr ) )
190 : : {
191 : : // GDAL < 2.3 has issues with use of VRT in multi-threaded
192 : : // scenarios. See https://github.com/qgis/QGIS/issues/24413 /
193 : : // https://trac.osgeo.org/gdal/ticket/6939
194 : : CPLSetConfigOption( "VRT_SHARED_SOURCE", "NO" );
195 : : }
196 : : #endif
197 : :
198 : : // To get buildSupportedRasterFileFilter the provider is called with empty uri
199 : 0 : if ( uri.isEmpty() )
200 : : {
201 : 0 : return;
202 : : }
203 : :
204 : 0 : mGdalDataset = nullptr;
205 : 0 : if ( dataset )
206 : : {
207 : 0 : mGdalBaseDataset = dataset;
208 : 0 : initBaseDataset();
209 : 0 : }
210 : : else
211 : : {
212 : 0 : ( void )initIfNeeded();
213 : : }
214 : 0 : }
215 : :
216 : 0 : QgsGdalProvider::QgsGdalProvider( const QgsGdalProvider &other )
217 : 0 : : QgsRasterDataProvider( other.dataSourceUri(), QgsDataProvider::ProviderOptions() )
218 : 0 : , mUpdate( false )
219 : 0 : {
220 : 0 : mDriverName = other.mDriverName;
221 : :
222 : : // The JP2OPENJPEG driver might consume too much memory on large datasets
223 : : // so make sure to really use a single one.
224 : : // The PostGISRaster driver internally uses a per-thread connection cache.
225 : : // This can lead to crashes if two datasets created by the same thread are used at the same time.
226 : 0 : bool forceUseSameDataset = ( mDriverName.toUpper() == QLatin1String( "JP2OPENJPEG" ) ||
227 : 0 : mDriverName == QLatin1String( "PostGISRaster" ) ||
228 : 0 : CSLTestBoolean( CPLGetConfigOption( "QGIS_GDAL_FORCE_USE_SAME_DATASET", "FALSE" ) ) );
229 : :
230 : 0 : if ( forceUseSameDataset )
231 : : {
232 : 0 : ++ ( *other.mpRefCounter );
233 : : // cppcheck-suppress copyCtorPointerCopying
234 : 0 : mpRefCounter = other.mpRefCounter;
235 : 0 : mpMutex = other.mpMutex;
236 : 0 : mpLightRefCounter = new QAtomicInt( 1 );
237 : 0 : mHasInit = other.mHasInit;
238 : 0 : mValid = other.mValid;
239 : 0 : mGdalBaseDataset = other.mGdalBaseDataset;
240 : 0 : mGdalDataset = other.mGdalDataset;
241 : 0 : }
242 : : else
243 : : {
244 : :
245 : 0 : ++ ( *other.mpLightRefCounter );
246 : 0 :
247 : 0 : mpRefCounter = new QAtomicInt( 1 );
248 : 0 : mpLightRefCounter = other.mpLightRefCounter;
249 : 0 : #if QT_VERSION < QT_VERSION_CHECK(5, 14, 0)
250 : : mpMutex = new QMutex( QMutex::Recursive );
251 : : #else
252 : 0 : mpMutex = new QRecursiveMutex();
253 : : #endif
254 : :
255 : 0 : mpParent = other.mpParent;
256 : :
257 : 0 : if ( getCachedGdalHandles( const_cast<QgsGdalProvider *>( &other ), mGdalBaseDataset, mGdalDataset ) )
258 : : {
259 : 0 : QgsDebugMsgLevel( QStringLiteral( "recycling already opened dataset" ), 5 );
260 : 0 : mHasInit = true;
261 : 0 : mValid = other.mValid;
262 : 0 : }
263 : 0 : else
264 : : {
265 : 0 : QgsDebugMsgLevel( QStringLiteral( "will need to open new dataset" ), 5 );
266 : 0 : mHasInit = false;
267 : 0 : mValid = false;
268 : : }
269 : :
270 : : }
271 : :
272 : 0 : mHasPyramids = other.mHasPyramids;
273 : 0 : mGdalDataType = other.mGdalDataType;
274 : 0 : mExtent = other.mExtent;
275 : 0 : mWidth = other.mWidth;
276 : 0 : mHeight = other.mHeight;
277 : 0 : mXBlockSize = other.mXBlockSize;
278 : 0 : mYBlockSize = other.mYBlockSize;
279 : 0 : memcpy( mGeoTransform, other.mGeoTransform, sizeof( mGeoTransform ) );
280 : 0 : mCrs = other.mCrs;
281 : 0 : mPyramidList = other.mPyramidList;
282 : 0 : mSubLayers = other.mSubLayers;
283 : 0 : mMaskBandExposedAsAlpha = other.mMaskBandExposedAsAlpha;
284 : 0 : mBandCount = other.mBandCount;
285 : 0 : copyBaseSettings( other );
286 : 0 : }
287 : :
288 : 0 : QString QgsGdalProvider::dataSourceUri( bool expandAuthConfig ) const
289 : : {
290 : 0 : if ( expandAuthConfig && QgsDataProvider::dataSourceUri( ).contains( QLatin1String( "authcfg" ) ) )
291 : 0 : {
292 : 0 : QString uri( QgsDataProvider::dataSourceUri() );
293 : : // Check for authcfg
294 : 0 : QRegularExpression authcfgRe( R"raw(authcfg='?([^'\s]+)'?)raw" );
295 : 0 : QRegularExpressionMatch match;
296 : 0 : if ( uri.contains( authcfgRe, &match ) )
297 : : {
298 : 0 : uri = uri.replace( match.captured( 0 ), QString() );
299 : 0 : QString configId( match.captured( 1 ) );
300 : 0 : QStringList connectionItems;
301 : 0 : connectionItems << uri;
302 : 0 : if ( QgsApplication::authManager()->updateDataSourceUriItems( connectionItems, configId, QStringLiteral( "gdal" ) ) )
303 : : {
304 : 0 : uri = connectionItems.first( );
305 : 0 : }
306 : 0 : }
307 : 0 : return uri;
308 : 0 : }
309 : : else
310 : : {
311 : 0 : return QgsDataProvider::dataSourceUri();
312 : : }
313 : 0 : }
314 : :
315 : 0 : QgsGdalProvider *QgsGdalProvider::clone() const
316 : : {
317 : 0 : return new QgsGdalProvider( *this );
318 : 0 : }
319 : :
320 : 0 : bool QgsGdalProvider::getCachedGdalHandles( QgsGdalProvider *provider,
321 : : GDALDatasetH &gdalBaseDataset,
322 : : GDALDatasetH &gdalDataset )
323 : : {
324 : 0 : QMutexLocker locker( sGdalProviderMutex() );
325 : :
326 : 0 : auto iter = mgDatasetCache.find( provider );
327 : 0 : if ( iter == mgDatasetCache.end() )
328 : : {
329 : 0 : return false;
330 : : }
331 : :
332 : 0 : if ( !iter.value().isEmpty() )
333 : : {
334 : 0 : DatasetPair pair = iter.value().takeFirst();
335 : 0 : mgDatasetCacheSize --;
336 : 0 : gdalBaseDataset = pair.mGdalBaseDataset;
337 : 0 : gdalDataset = pair.mGdalDataset;
338 : 0 : return true;
339 : : }
340 : 0 : return false;
341 : 0 : }
342 : :
343 : 0 : bool QgsGdalProvider::cacheGdalHandlesForLaterReuse( QgsGdalProvider *provider,
344 : 0 : GDALDatasetH gdalBaseDataset,
345 : : GDALDatasetH gdalDataset )
346 : : {
347 : 0 : QMutexLocker locker( sGdalProviderMutex() );
348 : :
349 : : // If the cache size goes above the soft limit, try to do evict a cached
350 : : // dataset for the provider that has the most cached entries
351 : 0 : if ( mgDatasetCacheSize >= MIN_THRESHOLD_FOR_CACHE_CLEANUP )
352 : 0 : {
353 : 0 : auto iter = mgDatasetCache.find( provider );
354 : 0 : if ( iter == mgDatasetCache.end() || iter.value().isEmpty() )
355 : : {
356 : 0 : QgsGdalProvider *candidateProvider = nullptr;
357 : 0 : int nLargestCountOfCachedDatasets = 0;
358 : 0 : for ( iter = mgDatasetCache.begin(); iter != mgDatasetCache.end(); ++iter )
359 : : {
360 : 0 : if ( iter.value().size() > nLargestCountOfCachedDatasets )
361 : : {
362 : 0 : candidateProvider = iter.key();
363 : 0 : nLargestCountOfCachedDatasets = iter.value().size();
364 : 0 : }
365 : 0 : }
366 : :
367 : : Q_ASSERT( candidateProvider );
368 : : Q_ASSERT( !mgDatasetCache[ candidateProvider ].isEmpty() );
369 : :
370 : : // If the candidate is ourselves, then do nothing
371 : 0 : if ( candidateProvider == provider )
372 : 0 : return false;
373 : :
374 : : // If the candidate provider has at least 2 cached datasets, then
375 : : // we can evict one.
376 : : // In the case where providers have at most one cached dataset, then
377 : : // evict one arbitrarily
378 : 0 : if ( nLargestCountOfCachedDatasets >= 2 ||
379 : 0 : mgDatasetCacheSize >= MAX_CACHE_SIZE )
380 : : {
381 : 0 : mgDatasetCacheSize --;
382 : 0 : DatasetPair pair = mgDatasetCache[ candidateProvider ].takeLast();
383 : 0 : if ( pair.mGdalBaseDataset != pair.mGdalDataset )
384 : : {
385 : 0 : GDALDereferenceDataset( pair.mGdalBaseDataset );
386 : 0 : }
387 : 0 : if ( pair.mGdalDataset )
388 : : {
389 : 0 : GDALClose( pair.mGdalDataset );
390 : 0 : }
391 : 0 : }
392 : 0 : }
393 : : else
394 : : {
395 : 0 : return false;
396 : : }
397 : 0 : }
398 : :
399 : : // Add handles to the cache
400 : 0 : auto iter = mgDatasetCache.find( provider );
401 : 0 : if ( iter == mgDatasetCache.end() )
402 : : {
403 : 0 : mgDatasetCache[provider] = QVector<DatasetPair>();
404 : 0 : iter = mgDatasetCache.find( provider );
405 : 0 : }
406 : :
407 : 0 : mgDatasetCacheSize ++;
408 : : DatasetPair pair;
409 : 0 : pair.mGdalBaseDataset = gdalBaseDataset;
410 : 0 : pair.mGdalDataset = gdalDataset;
411 : 0 : iter.value().push_back( pair );
412 : :
413 : 0 : return true;
414 : 0 : }
415 : :
416 : 0 : void QgsGdalProvider::closeCachedGdalHandlesFor( QgsGdalProvider *provider )
417 : : {
418 : 0 : QMutexLocker locker( sGdalProviderMutex() );
419 : 0 : auto iter = mgDatasetCache.find( provider );
420 : 0 : if ( iter != mgDatasetCache.end() )
421 : : {
422 : 0 : while ( !iter.value().isEmpty() )
423 : : {
424 : 0 : mgDatasetCacheSize --;
425 : 0 : DatasetPair pair = iter.value().takeLast();
426 : 0 : if ( pair.mGdalBaseDataset != pair.mGdalDataset )
427 : : {
428 : 0 : GDALDereferenceDataset( pair.mGdalBaseDataset );
429 : 0 : }
430 : 0 : if ( pair.mGdalDataset )
431 : : {
432 : 0 : GDALClose( pair.mGdalDataset );
433 : 0 : }
434 : : }
435 : 0 : mgDatasetCache.erase( iter );
436 : 0 : }
437 : 0 : }
438 : :
439 : :
440 : 0 : QgsGdalProvider::~QgsGdalProvider()
441 : 0 : {
442 : 0 : QMutexLocker locker( sGdalProviderMutex() );
443 : :
444 : 0 : if ( mGdalTransformerArg )
445 : 0 : GDALDestroyTransformer( mGdalTransformerArg );
446 : :
447 : 0 : int lightRefCounter = -- ( *mpLightRefCounter );
448 : 0 : int refCounter = -- ( *mpRefCounter );
449 : 0 : if ( refCounter == 0 )
450 : : {
451 : 0 : if ( mpParent && *mpParent && *mpParent != this && mGdalBaseDataset &&
452 : 0 : cacheGdalHandlesForLaterReuse( *mpParent, mGdalBaseDataset, mGdalDataset ) )
453 : : {
454 : : // do nothing
455 : 0 : }
456 : : else
457 : : {
458 : 0 : if ( mGdalBaseDataset != mGdalDataset )
459 : : {
460 : 0 : GDALDereferenceDataset( mGdalBaseDataset );
461 : 0 : }
462 : 0 : if ( mGdalDataset )
463 : : {
464 : : // Check if already a PAM (persistent auxiliary metadata) file exists
465 : 0 : QString pamFile = dataSourceUri( true ) + QLatin1String( ".aux.xml" );
466 : 0 : bool pamFileAlreadyExists = QFileInfo::exists( pamFile );
467 : :
468 : 0 : GDALClose( mGdalDataset );
469 : :
470 : : // If GDAL created a PAM file right now by using estimated metadata, delete it right away
471 : 0 : if ( !mStatisticsAreReliable && !pamFileAlreadyExists && QFileInfo::exists( pamFile ) )
472 : 0 : QFile( pamFile ).remove();
473 : 0 : }
474 : :
475 : 0 : if ( mpParent && *mpParent == this )
476 : : {
477 : 0 : *mpParent = nullptr;
478 : 0 : closeCachedGdalHandlesFor( this );
479 : 0 : }
480 : : }
481 : 0 : delete mpMutex;
482 : 0 : delete mpRefCounter;
483 : 0 : if ( lightRefCounter == 0 )
484 : : {
485 : 0 : delete mpLightRefCounter;
486 : 0 : delete mpParent;
487 : 0 : }
488 : 0 : }
489 : 0 : }
490 : :
491 : :
492 : 0 : void QgsGdalProvider::closeDataset()
493 : : {
494 : 0 : if ( !mValid )
495 : : {
496 : 0 : return;
497 : : }
498 : 0 : mValid = false;
499 : :
500 : 0 : if ( mGdalBaseDataset != mGdalDataset )
501 : : {
502 : 0 : GDALDereferenceDataset( mGdalBaseDataset );
503 : 0 : }
504 : 0 : mGdalBaseDataset = nullptr;
505 : :
506 : 0 : GDALClose( mGdalDataset );
507 : 0 : mGdalDataset = nullptr;
508 : :
509 : 0 : closeCachedGdalHandlesFor( this );
510 : 0 : }
511 : :
512 : 0 : void QgsGdalProvider::reloadProviderData()
513 : : {
514 : 0 : QMutexLocker locker( mpMutex );
515 : 0 : closeDataset();
516 : :
517 : 0 : mHasInit = false;
518 : 0 : ( void )initIfNeeded();
519 : 0 : }
520 : :
521 : 0 : QString QgsGdalProvider::htmlMetadata()
522 : : {
523 : 0 : QMutexLocker locker( mpMutex );
524 : 0 : if ( !initIfNeeded() )
525 : 0 : return QString();
526 : :
527 : 0 : GDALDriverH hDriver = GDALGetDriverByName( mDriverName.toLocal8Bit().constData() );
528 : 0 : if ( !hDriver )
529 : 0 : return QString();
530 : :
531 : 0 : QString myMetadata;
532 : :
533 : : // GDAL Driver description
534 : 0 : myMetadata += QStringLiteral( "<tr><td class=\"highlight\">" ) + tr( "GDAL Driver Description" ) + QStringLiteral( "</td><td>" ) + mDriverName + QStringLiteral( "</td></tr>\n" );
535 : :
536 : : // GDAL Driver Metadata
537 : 0 : myMetadata += QStringLiteral( "<tr><td class=\"highlight\">" ) + tr( "GDAL Driver Metadata" ) + QStringLiteral( "</td><td>" ) + QString( GDALGetMetadataItem( hDriver, GDAL_DMD_LONGNAME, nullptr ) ) + QStringLiteral( "</td></tr>\n" );
538 : :
539 : : // Dataset description
540 : 0 : myMetadata += QStringLiteral( "<tr><td class=\"highlight\">" ) + tr( "Dataset Description" ) + QStringLiteral( "</td><td>" ) + QString::fromUtf8( GDALGetDescription( mGdalDataset ) ) + QStringLiteral( "</td></tr>\n" );
541 : :
542 : : // compression
543 : 0 : QString compression = QString( GDALGetMetadataItem( mGdalDataset, "COMPRESSION", "IMAGE_STRUCTURE" ) );
544 : 0 : myMetadata += QStringLiteral( "<tr><td class=\"highlight\">" ) + tr( "Compression" ) + QStringLiteral( "</td><td>" ) + compression + QStringLiteral( "</td></tr>\n" );
545 : :
546 : : // Band details
547 : 0 : for ( int i = 1; i <= GDALGetRasterCount( mGdalDataset ); ++i )
548 : : {
549 : 0 : GDALRasterBandH gdalBand = GDALGetRasterBand( mGdalDataset, i );
550 : 0 : char **GDALmetadata = GDALGetMetadata( gdalBand, nullptr );
551 : 0 : myMetadata += QStringLiteral( "<tr><td class=\"highlight\">" ) + tr( "Band %1" ).arg( i ) + QStringLiteral( "</td><td>" );
552 : 0 : if ( GDALmetadata )
553 : : {
554 : 0 : QStringList metadata = QgsOgrUtils::cStringListToQStringList( GDALmetadata );
555 : 0 : myMetadata += QgsHtmlUtils::buildBulletList( metadata );
556 : 0 : }
557 : :
558 : 0 : char **GDALcategories = GDALGetRasterCategoryNames( gdalBand );
559 : :
560 : 0 : if ( GDALcategories )
561 : : {
562 : 0 : QStringList categories = QgsOgrUtils::cStringListToQStringList( GDALcategories );
563 : 0 : myMetadata += QgsHtmlUtils::buildBulletList( categories );
564 : 0 : }
565 : 0 : myMetadata += QLatin1String( "</td></tr>" );
566 : 0 : }
567 : :
568 : : // More information
569 : 0 : myMetadata += QStringLiteral( "<tr><td class=\"highlight\">" ) + tr( "More information" ) + QStringLiteral( "</td><td>\n" );
570 : :
571 : 0 : if ( mMaskBandExposedAsAlpha )
572 : : {
573 : 0 : myMetadata += tr( "Mask band (exposed as alpha band)" ) + QStringLiteral( "<br />\n" );
574 : 0 : }
575 : :
576 : 0 : char **GDALmetadata = GDALGetMetadata( mGdalDataset, nullptr );
577 : 0 : if ( GDALmetadata )
578 : : {
579 : 0 : QStringList metadata = QgsOgrUtils::cStringListToQStringList( GDALmetadata );
580 : 0 : myMetadata += QgsHtmlUtils::buildBulletList( metadata );
581 : 0 : }
582 : :
583 : : //just use the first band
584 : 0 : if ( GDALGetRasterCount( mGdalDataset ) > 0 )
585 : : {
586 : 0 : GDALRasterBandH myGdalBand = GDALGetRasterBand( mGdalDataset, 1 );
587 : 0 : if ( GDALGetOverviewCount( myGdalBand ) > 0 )
588 : : {
589 : : int myOverviewInt;
590 : 0 : for ( myOverviewInt = 0; myOverviewInt < GDALGetOverviewCount( myGdalBand ); myOverviewInt++ )
591 : : {
592 : : GDALRasterBandH myOverview;
593 : 0 : myOverview = GDALGetOverview( myGdalBand, myOverviewInt );
594 : 0 : QStringList metadata;
595 : 0 : metadata.append( QStringLiteral( "X : " ) + QString::number( GDALGetRasterBandXSize( myOverview ) ) );
596 : 0 : metadata.append( QStringLiteral( "Y : " ) + QString::number( GDALGetRasterBandYSize( myOverview ) ) );
597 : 0 : myMetadata += QgsHtmlUtils::buildBulletList( metadata );
598 : 0 : }
599 : 0 : }
600 : 0 : }
601 : :
602 : : // End more information
603 : 0 : myMetadata += QLatin1String( "</td></tr>\n" );
604 : :
605 : : // Dimensions
606 : 0 : myMetadata += QStringLiteral( "<tr><td class=\"highlight\">" ) + tr( "Dimensions" ) + QStringLiteral( "</td><td>" );
607 : 0 : myMetadata += tr( "X: %1 Y: %2 Bands: %3" )
608 : 0 : .arg( GDALGetRasterXSize( mGdalDataset ) )
609 : 0 : .arg( GDALGetRasterYSize( mGdalDataset ) )
610 : 0 : .arg( GDALGetRasterCount( mGdalDataset ) );
611 : 0 : myMetadata += QLatin1String( "</td></tr>\n" );
612 : :
613 : 0 : if ( GDALGetGeoTransform( mGdalDataset, mGeoTransform ) != CE_None )
614 : : {
615 : : // if the raster does not have a valid transform we need to use
616 : : // a pixel size of (1,-1), but GDAL returns (1,1)
617 : 0 : mGeoTransform[5] = -1;
618 : 0 : }
619 : : else
620 : : {
621 : : // Origin
622 : 0 : myMetadata += QStringLiteral( "<tr><td class=\"highlight\">" ) + tr( "Origin" ) + QStringLiteral( "</td><td>" ) + QString::number( mGeoTransform[0] ) + QStringLiteral( "," ) + QString::number( mGeoTransform[3] ) + QStringLiteral( "</td></tr>\n" );
623 : :
624 : : // Pixel size
625 : 0 : myMetadata += QStringLiteral( "<tr><td class=\"highlight\">" ) + tr( "Pixel Size" ) + QStringLiteral( "</td><td>" ) + QString::number( mGeoTransform[1], 'g', 19 ) + QStringLiteral( "," ) + QString::number( mGeoTransform[5], 'g', 19 ) + QStringLiteral( "</td></tr>\n" );
626 : : }
627 : :
628 : 0 : return myMetadata;
629 : 0 : }
630 : :
631 : 0 : QgsRasterBlock *QgsGdalProvider::block( int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback )
632 : : {
633 : 0 : std::unique_ptr< QgsRasterBlock > block = std::make_unique< QgsRasterBlock >( dataType( bandNo ), width, height );
634 : 0 : if ( !initIfNeeded() )
635 : 0 : return block.release();
636 : 0 : if ( sourceHasNoDataValue( bandNo ) && useSourceNoDataValue( bandNo ) )
637 : : {
638 : 0 : block->setNoDataValue( sourceNoDataValue( bandNo ) );
639 : 0 : }
640 : :
641 : 0 : if ( block->isEmpty() )
642 : : {
643 : 0 : return block.release();
644 : : }
645 : :
646 : 0 : if ( !mExtent.intersects( extent ) )
647 : : {
648 : : // the requested extent is completely outside of the raster's extent - nothing to do
649 : 0 : block->setIsNoData();
650 : 0 : return block.release();
651 : : }
652 : :
653 : 0 : if ( !mExtent.contains( extent ) )
654 : : {
655 : 0 : QRect subRect = QgsRasterBlock::subRect( extent, width, height, mExtent );
656 : 0 : block->setIsNoDataExcept( subRect );
657 : 0 : }
658 : 0 : if ( !readBlock( bandNo, extent, width, height, block->bits(), feedback ) )
659 : : {
660 : 0 : block->setError( { tr( "Error occurred while reading block." ), QStringLiteral( "GDAL" ) } );
661 : 0 : block->setIsNoData();
662 : 0 : block->setValid( false );
663 : 0 : return block.release();
664 : : }
665 : : // apply scale and offset
666 : : Q_ASSERT( block ); // to make cppcheck happy
667 : 0 : block->applyScaleOffset( bandScale( bandNo ), bandOffset( bandNo ) );
668 : 0 : block->applyNoDataValues( userNoDataValues( bandNo ) );
669 : 0 : return block.release();
670 : 0 : }
671 : :
672 : 0 : bool QgsGdalProvider::readBlock( int bandNo, int xBlock, int yBlock, void *data )
673 : : {
674 : 0 : QMutexLocker locker( mpMutex );
675 : 0 : if ( !initIfNeeded() )
676 : 0 : return false;
677 : :
678 : : // TODO!!!: Check data alignment!!! May it happen that nearest value which
679 : : // is not nearest is assigned to an output cell???
680 : :
681 : 0 : GDALRasterBandH myGdalBand = getBand( bandNo );
682 : : //GDALReadBlock( myGdalBand, xBlock, yBlock, block );
683 : :
684 : : // We have to read with correct data type consistent with other readBlock functions
685 : 0 : int xOff = xBlock * mXBlockSize;
686 : 0 : int yOff = yBlock * mYBlockSize;
687 : 0 : CPLErr err = gdalRasterIO( myGdalBand, GF_Read, xOff, yOff, mXBlockSize, mYBlockSize, data, mXBlockSize, mYBlockSize, ( GDALDataType ) mGdalDataType.at( bandNo - 1 ), 0, 0 );
688 : 0 : if ( err != CPLE_None )
689 : : {
690 : 0 : QgsLogger::warning( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
691 : 0 : return false;
692 : : }
693 : :
694 : 0 : return true;
695 : 0 : }
696 : :
697 : 0 : bool QgsGdalProvider::canDoResampling(
698 : : int bandNo,
699 : : const QgsRectangle &reqExtent,
700 : : int bufferWidthPix,
701 : : int bufferHeightPix )
702 : : {
703 : 0 : QMutexLocker locker( mpMutex );
704 : 0 : if ( !initIfNeeded() )
705 : 0 : return false;
706 : :
707 : 0 : GDALRasterBandH gdalBand = getBand( bandNo );
708 : 0 : if ( GDALGetRasterColorTable( gdalBand ) )
709 : 0 : return false;
710 : :
711 : 0 : const double reqXRes = reqExtent.width() / bufferWidthPix;
712 : 0 : const double reqYRes = reqExtent.height() / bufferHeightPix;
713 : 0 : const double srcXRes = mGeoTransform[1];
714 : 0 : const double srcYRes = fabs( mGeoTransform[5] );
715 : : // Compute the resampling factor :
716 : : // < 1 means upsampling / zoom-in
717 : : // > 1 means downsampling / zoom-out
718 : 0 : const double resamplingFactor = std::max( reqXRes / srcXRes, reqYRes / srcYRes );
719 : :
720 : 0 : if ( resamplingFactor < 1 )
721 : : {
722 : : // upsampling
723 : 0 : return mZoomedInResamplingMethod != QgsRasterDataProvider::ResamplingMethod::Nearest;
724 : : }
725 : :
726 : 0 : if ( resamplingFactor < 1.1 )
727 : : {
728 : : // very close to nominal resolution ==> check compatibility of zoom-in or zoom-out resampler with what GDAL can do
729 : 0 : return mZoomedInResamplingMethod != QgsRasterDataProvider::ResamplingMethod::Nearest ||
730 : 0 : mZoomedOutResamplingMethod != QgsRasterDataProvider::ResamplingMethod::Nearest;
731 : : }
732 : :
733 : : // if no zoom out resampling, exit now
734 : 0 : if ( mZoomedOutResamplingMethod == QgsRasterDataProvider::ResamplingMethod::Nearest )
735 : : {
736 : 0 : return false;
737 : : }
738 : :
739 : : // if the resampling factor is not too large, we can do the downsampling
740 : : // from the full resolution with acceptable performance
741 : 0 : if ( resamplingFactor <= ( mMaxOversampling + .1 ) )
742 : : {
743 : 0 : return true;
744 : : }
745 : :
746 : : // otherwise, we have to check that we have an overview band that satisfies
747 : : // that criterion
748 : 0 : const int nOvrCount = GDALGetOverviewCount( gdalBand );
749 : 0 : for ( int i = 0; i < nOvrCount; i++ )
750 : : {
751 : 0 : GDALRasterBandH hOvrBand = GDALGetOverview( gdalBand, i );
752 : 0 : const double ovrResamplingFactor = xSize() / static_cast<double>( GDALGetRasterBandXSize( hOvrBand ) );
753 : 0 : if ( resamplingFactor <= ( mMaxOversampling + .1 ) * ovrResamplingFactor )
754 : : {
755 : 0 : return true;
756 : : }
757 : 0 : }
758 : :
759 : : // too much zoomed out
760 : 0 : return false;
761 : 0 : }
762 : :
763 : 0 : static GDALRIOResampleAlg getGDALResamplingAlg( QgsGdalProvider::ResamplingMethod method )
764 : : {
765 : 0 : GDALRIOResampleAlg eResampleAlg = GRIORA_NearestNeighbour;
766 : 0 : switch ( method )
767 : : {
768 : : case QgsGdalProvider::ResamplingMethod::Nearest:
769 : 0 : eResampleAlg = GRIORA_NearestNeighbour;
770 : 0 : break;
771 : :
772 : : case QgsGdalProvider::ResamplingMethod::Bilinear:
773 : 0 : eResampleAlg = GRIORA_Bilinear;
774 : 0 : break;
775 : :
776 : : case QgsGdalProvider::ResamplingMethod::Cubic:
777 : 0 : eResampleAlg = GRIORA_Cubic;
778 : 0 : break;
779 : :
780 : : case QgsRasterDataProvider::ResamplingMethod::CubicSpline:
781 : 0 : eResampleAlg = GRIORA_CubicSpline;
782 : 0 : break;
783 : :
784 : : case QgsRasterDataProvider::ResamplingMethod::Lanczos:
785 : 0 : eResampleAlg = GRIORA_Lanczos;
786 : 0 : break;
787 : :
788 : : case QgsRasterDataProvider::ResamplingMethod::Average:
789 : 0 : eResampleAlg = GRIORA_Average;
790 : 0 : break;
791 : :
792 : : case QgsRasterDataProvider::ResamplingMethod::Mode:
793 : 0 : eResampleAlg = GRIORA_Mode;
794 : 0 : break;
795 : :
796 : : case QgsRasterDataProvider::ResamplingMethod::Gauss:
797 : 0 : eResampleAlg = GRIORA_Gauss;
798 : 0 : break;
799 : : }
800 : :
801 : 0 : return eResampleAlg;
802 : : }
803 : :
804 : 0 : bool QgsGdalProvider::readBlock( int bandNo, QgsRectangle const &reqExtent, int bufferWidthPix, int bufferHeightPix, void *data, QgsRasterBlockFeedback *feedback )
805 : : {
806 : 0 : QMutexLocker locker( mpMutex );
807 : 0 : if ( !initIfNeeded() )
808 : 0 : return false;
809 : :
810 : 0 : QgsDebugMsgLevel( "bufferWidthPix = " + QString::number( bufferWidthPix ), 5 );
811 : 0 : QgsDebugMsgLevel( "bufferHeightPix = " + QString::number( bufferHeightPix ), 5 );
812 : 0 : QgsDebugMsgLevel( "reqExtent: " + reqExtent.toString(), 5 );
813 : :
814 : 0 : for ( int i = 0; i < 6; i++ )
815 : : {
816 : 0 : QgsDebugMsgLevel( QStringLiteral( "transform : %1" ).arg( mGeoTransform[i] ), 5 );
817 : 0 : }
818 : :
819 : 0 : const size_t dataSize = static_cast<size_t>( dataTypeSize( bandNo ) );
820 : :
821 : 0 : QgsRectangle intersectExtent = reqExtent.intersect( mExtent );
822 : 0 : if ( intersectExtent.isEmpty() )
823 : : {
824 : 0 : QgsDebugMsgLevel( QStringLiteral( "draw request outside view extent." ), 2 );
825 : 0 : return false;
826 : : }
827 : 0 : QgsDebugMsgLevel( "extent: " + mExtent.toString(), 5 );
828 : 0 : QgsDebugMsgLevel( "intersectExtent: " + intersectExtent.toString(), 5 );
829 : :
830 : 0 : const double reqXRes = reqExtent.width() / bufferWidthPix;
831 : 0 : const double reqYRes = reqExtent.height() / bufferHeightPix;
832 : :
833 : 0 : const double srcXRes = mGeoTransform[1];
834 : 0 : const double srcYRes = mGeoTransform[5]; // may be negative?
835 : 0 : QgsDebugMsgLevel( QStringLiteral( "reqXRes = %1 reqYRes = %2 srcXRes = %3 srcYRes = %4" ).arg( reqXRes ).arg( reqYRes ).arg( srcXRes ).arg( srcYRes ), 5 );
836 : 0 : const double resamplingFactor = std::max( reqXRes / srcXRes, reqYRes / srcYRes );
837 : :
838 : 0 : GDALRasterBandH gdalBand = getBand( bandNo );
839 : 0 : const GDALDataType type = static_cast<GDALDataType>( mGdalDataType.at( bandNo - 1 ) );
840 : :
841 : : // Find top, bottom rows and left, right column the raster extent covers
842 : : // These are limits in target grid space
843 : 0 : QRect subRect = QgsRasterBlock::subRect( reqExtent, bufferWidthPix, bufferHeightPix, intersectExtent );
844 : 0 : const int tgtTopOri = subRect.top();
845 : 0 : const int tgtBottomOri = subRect.bottom();
846 : 0 : const int tgtLeftOri = subRect.left();
847 : 0 : const int tgtRightOri = subRect.right();
848 : :
849 : : // Calculate rows/cols limits in source raster grid space
850 : 0 : int srcLeft = 0;
851 : 0 : int srcTop = 0;
852 : 0 : int srcBottom = ySize() - 1;
853 : 0 : int srcRight = xSize() - 1;
854 : :
855 : : // Get necessary src extent aligned to src resolution
856 : 0 : if ( mExtent.xMinimum() < intersectExtent.xMinimum() )
857 : : {
858 : 0 : srcLeft = static_cast<int>( std::floor( ( intersectExtent.xMinimum() - mExtent.xMinimum() ) / srcXRes ) );
859 : 0 : }
860 : 0 : if ( mExtent.xMaximum() > intersectExtent.xMaximum() )
861 : : {
862 : : // Clamp to raster width to avoid potential rounding errors (see GH #34435)
863 : 0 : srcRight = std::min( mWidth - 1, static_cast<int>( std::floor( ( intersectExtent.xMaximum() - mExtent.xMinimum() ) / srcXRes ) ) );
864 : 0 : }
865 : :
866 : : // GDAL states that mGeoTransform[3] is top, may it also be bottom and mGeoTransform[5] positive?
867 : 0 : if ( mExtent.yMaximum() > intersectExtent.yMaximum() )
868 : : {
869 : 0 : srcTop = static_cast<int>( std::floor( -1. * ( mExtent.yMaximum() - intersectExtent.yMaximum() ) / srcYRes ) );
870 : 0 : }
871 : 0 : if ( mExtent.yMinimum() < intersectExtent.yMinimum() )
872 : : {
873 : : // Clamp to raster height to avoid potential rounding errors (see GH #34435)
874 : 0 : srcBottom = std::min( mHeight - 1, static_cast<int>( std::floor( -1. * ( mExtent.yMaximum() - intersectExtent.yMinimum() ) / srcYRes ) ) );
875 : 0 : }
876 : :
877 : 0 : const int srcWidth = srcRight - srcLeft + 1;
878 : 0 : const int srcHeight = srcBottom - srcTop + 1;
879 : :
880 : : // Use GDAL resampling if asked and possible
881 : 0 : if ( mProviderResamplingEnabled &&
882 : 0 : canDoResampling( bandNo, reqExtent, bufferWidthPix, bufferHeightPix ) )
883 : : {
884 : 0 : int tgtTop = tgtTopOri;
885 : 0 : int tgtBottom = tgtBottomOri;
886 : 0 : int tgtLeft = tgtLeftOri;
887 : 0 : int tgtRight = tgtRightOri;
888 : :
889 : : // Deal with situations that requires adjusting tgt coordinates.
890 : : // This is due rounding in QgsRasterBlock::subRect() and
891 : : // when the difference between the raster extent and the
892 : : // request extent is less than half of the request resolution.
893 : : // This is to avoid to send dfXOff, dfYOff, dfXSize, dfYSize values that
894 : : // would be outside of the raster extent, as GDAL (at least currently)
895 : : // rejects them
896 : 0 : if ( reqExtent.xMinimum() + tgtLeft * reqXRes < mExtent.xMinimum() )
897 : : {
898 : 0 : if ( GDALRasterIO( gdalBand, GF_Read, 0, srcTop, 1, srcHeight,
899 : 0 : static_cast<char *>( data ) + tgtTopOri * bufferWidthPix * dataSize,
900 : 0 : 1, tgtBottomOri - tgtTopOri + 1, type,
901 : 0 : dataSize, dataSize * bufferWidthPix ) != CE_None )
902 : : {
903 : 0 : return false;
904 : : }
905 : 0 : tgtLeft ++;
906 : 0 : }
907 : 0 : if ( reqExtent.yMaximum() - tgtTop * reqYRes > mExtent.yMaximum() )
908 : : {
909 : 0 : if ( GDALRasterIO( gdalBand, GF_Read, srcLeft, 0, srcWidth, 1,
910 : 0 : static_cast<char *>( data ) + tgtLeftOri * dataSize,
911 : 0 : tgtRightOri - tgtLeftOri + 1, 1, type,
912 : 0 : dataSize, dataSize * bufferWidthPix ) != CE_None )
913 : : {
914 : 0 : return false;
915 : : }
916 : 0 : tgtTop ++;
917 : 0 : }
918 : 0 : if ( reqExtent.xMinimum() + ( tgtRight + 1 ) * reqXRes > mExtent.xMaximum() )
919 : : {
920 : 0 : if ( GDALRasterIO( gdalBand, GF_Read, xSize() - 1, srcTop, 1, srcHeight,
921 : 0 : static_cast<char *>( data ) + ( tgtTopOri * bufferWidthPix + tgtRightOri ) * dataSize,
922 : 0 : 1, tgtBottomOri - tgtTopOri + 1, type,
923 : 0 : dataSize, dataSize * bufferWidthPix ) != CE_None )
924 : : {
925 : 0 : return false;
926 : : }
927 : 0 : tgtRight --;
928 : 0 : }
929 : 0 : if ( reqExtent.yMaximum() - ( tgtBottom + 1 ) * reqYRes < mExtent.yMinimum() )
930 : : {
931 : 0 : if ( GDALRasterIO( gdalBand, GF_Read, srcLeft, ySize() - 1, srcWidth, 1,
932 : 0 : static_cast<char *>( data ) + ( tgtBottomOri * bufferWidthPix + tgtLeftOri ) * dataSize,
933 : 0 : tgtRightOri - tgtLeftOri + 1, 1, type,
934 : 0 : dataSize, dataSize * bufferWidthPix ) != CE_None )
935 : : {
936 : 0 : return false;
937 : : }
938 : 0 : tgtBottom --;
939 : 0 : }
940 : :
941 : 0 : int tgtWidth = tgtRight - tgtLeft + 1;
942 : 0 : int tgtHeight = tgtBottom - tgtTop + 1;
943 : 0 : if ( tgtWidth > 0 && tgtHeight > 0 )
944 : : {
945 : 0 : QgsDebugMsgLevel( QStringLiteral( "using GDAL resampling path" ), 5 );
946 : : GDALRasterIOExtraArg sExtraArg;
947 : 0 : INIT_RASTERIO_EXTRA_ARG( sExtraArg );
948 : :
949 : : ResamplingMethod method;
950 : 0 : if ( resamplingFactor < 1 )
951 : : {
952 : 0 : method = mZoomedInResamplingMethod;
953 : 0 : }
954 : 0 : else if ( resamplingFactor < 1.1 )
955 : : {
956 : : // very close to nominal resolution ==> use either zoomed out resampler or zoomed in resampler
957 : 0 : if ( mZoomedOutResamplingMethod != ResamplingMethod::Nearest )
958 : 0 : method = mZoomedOutResamplingMethod;
959 : : else
960 : 0 : method = mZoomedInResamplingMethod;
961 : 0 : }
962 : : else
963 : : {
964 : 0 : method = mZoomedOutResamplingMethod;
965 : : }
966 : 0 : sExtraArg.eResampleAlg = getGDALResamplingAlg( method );
967 : :
968 : 0 : if ( mMaskBandExposedAsAlpha &&
969 : 0 : bandNo == GDALGetRasterCount( mGdalDataset ) + 1 &&
970 : 0 : sExtraArg.eResampleAlg != GRIORA_NearestNeighbour &&
971 : 0 : sExtraArg.eResampleAlg != GRIORA_Bilinear )
972 : : {
973 : : // As time of writing in GDAL up to 3.1, there's a difference of behavior
974 : : // when using non-nearest resampling on mask bands.
975 : : // With bilinear, values are returned in [0,255] range
976 : : // whereas with cubic (and other convolution-based methods), there are in [0,1] range.
977 : : // As we want [0,255] range, fallback to Bilinear for the mask band
978 : 0 : sExtraArg.eResampleAlg = GRIORA_Bilinear;
979 : 0 : }
980 : :
981 : 0 : sExtraArg.bFloatingPointWindowValidity = true;
982 : 0 : sExtraArg.dfXOff = ( reqExtent.xMinimum() + tgtLeft * reqXRes - mExtent.xMinimum() ) / srcXRes;
983 : 0 : sExtraArg.dfYOff = ( mExtent.yMaximum() - ( reqExtent.yMaximum() - tgtTop * reqYRes ) ) / -srcYRes;
984 : 0 : sExtraArg.dfXSize = tgtWidth * reqXRes / srcXRes;
985 : 0 : sExtraArg.dfYSize = tgtHeight * reqYRes / -srcYRes;
986 : 0 : return GDALRasterIOEx( gdalBand, GF_Read,
987 : 0 : static_cast<int>( std::floor( sExtraArg.dfXOff ) ),
988 : 0 : static_cast<int>( std::floor( sExtraArg.dfYOff ) ),
989 : 0 : std::max( 1, static_cast<int>( std::floor( sExtraArg.dfXSize ) ) ),
990 : 0 : std::max( 1, static_cast<int>( std::floor( sExtraArg.dfYSize ) ) ),
991 : 0 : static_cast<char *>( data ) +
992 : 0 : ( tgtTop * bufferWidthPix + tgtLeft ) * dataSize,
993 : 0 : tgtWidth,
994 : 0 : tgtHeight,
995 : 0 : type,
996 : 0 : dataSize,
997 : 0 : dataSize * bufferWidthPix,
998 : 0 : &sExtraArg ) == CE_None;
999 : : }
1000 : 0 : }
1001 : : // Provider resampling was asked but we cannot do it in a performant way
1002 : : // (too much downsampling compared to the allowed maximum resampling factor),
1003 : : // so fallback to something replicating QgsRasterResampleFilter behavior
1004 : 0 : else if ( mProviderResamplingEnabled &&
1005 : 0 : mZoomedOutResamplingMethod != QgsRasterDataProvider::ResamplingMethod::Nearest &&
1006 : 0 : resamplingFactor > 1 )
1007 : : {
1008 : : // Do the resampling in two steps:
1009 : : // - downsample with nearest neighbour down to tgtWidth * mMaxOversampling, tgtHeight * mMaxOversampling
1010 : : // - then downsample with mZoomedOutResamplingMethod down to tgtWidth, tgtHeight
1011 : 0 : const int tgtWidth = tgtRightOri - tgtLeftOri + 1;
1012 : 0 : const int tgtHeight = tgtBottomOri - tgtTopOri + 1;
1013 : :
1014 : 0 : const int tmpWidth = static_cast<int>( tgtWidth * mMaxOversampling + 0.5 );
1015 : 0 : const int tmpHeight = static_cast<int>( tgtHeight * mMaxOversampling + 0.5 );
1016 : :
1017 : : // Allocate temporary block
1018 : 0 : size_t bufferSize = dataSize * static_cast<size_t>( tmpWidth ) * static_cast<size_t>( tmpHeight );
1019 : : #ifdef Q_PROCESSOR_X86_32
1020 : : // Safety check for 32 bit systems
1021 : : qint64 _buffer_size = dataSize * static_cast<qint64>( tmpWidth ) * static_cast<qint64>( tmpHeight );
1022 : : if ( _buffer_size != static_cast<qint64>( bufferSize ) )
1023 : : {
1024 : : QgsDebugMsg( QStringLiteral( "Integer overflow calculating buffer size on a 32 bit system." ) );
1025 : : return false;
1026 : : }
1027 : : #endif
1028 : 0 : char *tmpBlock = static_cast<char *>( qgsMalloc( bufferSize ) );
1029 : 0 : if ( ! tmpBlock )
1030 : : {
1031 : 0 : QgsDebugMsgLevel( QStringLiteral( "Couldn't allocate temporary buffer of %1 bytes" ).arg( dataSize * tmpWidth * tmpHeight ), 5 );
1032 : 0 : return false;
1033 : : }
1034 : 0 : CPLErrorReset();
1035 : :
1036 : :
1037 : 0 : CPLErr err = gdalRasterIO( gdalBand, GF_Read,
1038 : 0 : srcLeft, srcTop, srcWidth, srcHeight,
1039 : 0 : static_cast<void *>( tmpBlock ),
1040 : 0 : tmpWidth, tmpHeight, type,
1041 : 0 : 0, 0, feedback );
1042 : :
1043 : 0 : if ( err != CPLE_None )
1044 : : {
1045 : 0 : const QString lastError = QString::fromUtf8( CPLGetLastErrorMsg() ) ;
1046 : 0 : if ( feedback )
1047 : 0 : feedback->appendError( lastError );
1048 : :
1049 : 0 : QgsLogger::warning( "RasterIO error: " + lastError );
1050 : 0 : qgsFree( tmpBlock );
1051 : 0 : return false;
1052 : 0 : }
1053 : :
1054 : 0 : GDALDriverH hDriverMem = GDALGetDriverByName( "MEM" );
1055 : 0 : if ( !hDriverMem )
1056 : : {
1057 : 0 : qgsFree( tmpBlock );
1058 : 0 : return false;
1059 : : }
1060 : 0 : gdal::dataset_unique_ptr hSrcDS( GDALCreate(
1061 : 0 : hDriverMem, "", tmpWidth, tmpHeight, 0, GDALGetRasterDataType( gdalBand ), nullptr ) );
1062 : :
1063 : 0 : char **papszOptions = QgsGdalUtils::papszFromStringList( QStringList()
1064 : 0 : << QStringLiteral( "PIXELOFFSET=%1" ).arg( dataSize )
1065 : 0 : << QStringLiteral( "LINEOFFSET=%1" ).arg( dataSize * tmpWidth )
1066 : 0 : << QStringLiteral( "DATAPOINTER=%1" ).arg( reinterpret_cast< qulonglong >( tmpBlock ) ) );
1067 : 0 : GDALAddBand( hSrcDS.get(), GDALGetRasterDataType( gdalBand ), papszOptions );
1068 : 0 : CSLDestroy( papszOptions );
1069 : :
1070 : : GDALRasterIOExtraArg sExtraArg;
1071 : 0 : INIT_RASTERIO_EXTRA_ARG( sExtraArg );
1072 : :
1073 : 0 : sExtraArg.eResampleAlg = getGDALResamplingAlg( mZoomedOutResamplingMethod );
1074 : :
1075 : 0 : CPLErr eErr = GDALRasterIOEx( GDALGetRasterBand( hSrcDS.get(), 1 ),
1076 : : GF_Read,
1077 : 0 : 0, 0, tmpWidth, tmpHeight,
1078 : 0 : static_cast<char *>( data ) +
1079 : 0 : ( tgtTopOri * bufferWidthPix + tgtLeftOri ) * dataSize,
1080 : 0 : tgtWidth,
1081 : 0 : tgtHeight,
1082 : 0 : type,
1083 : 0 : dataSize,
1084 : 0 : dataSize * bufferWidthPix,
1085 : : &sExtraArg );
1086 : :
1087 : 0 : qgsFree( tmpBlock );
1088 : :
1089 : 0 : return eErr == CE_None;
1090 : 0 : }
1091 : :
1092 : 0 : const int tgtTop = tgtTopOri;
1093 : 0 : const int tgtBottom = tgtBottomOri;
1094 : 0 : const int tgtLeft = tgtLeftOri;
1095 : 0 : const int tgtRight = tgtRightOri;
1096 : 0 : QgsDebugMsgLevel( QStringLiteral( "tgtTop = %1 tgtBottom = %2 tgtLeft = %3 tgtRight = %4" ).arg( tgtTop ).arg( tgtBottom ).arg( tgtLeft ).arg( tgtRight ), 5 );
1097 : : // target size in pizels
1098 : 0 : const int tgtWidth = tgtRight - tgtLeft + 1;
1099 : 0 : const int tgtHeight = tgtBottom - tgtTop + 1;
1100 : :
1101 : 0 : QgsDebugMsgLevel( QStringLiteral( "srcTop = %1 srcBottom = %2 srcWidth = %3 srcHeight = %4" ).arg( srcTop ).arg( srcBottom ).arg( srcWidth ).arg( srcHeight ), 5 );
1102 : :
1103 : : // Determine the dimensions of the buffer into which we will ask GDAL to write
1104 : : // pixels.
1105 : : // In downsampling scenarios, we will use the request resolution to compute the dimension
1106 : : // In upsampling (or exactly at raster resolution) scenarios, we will use the raster resolution to compute the dimension
1107 : 0 : int tmpWidth = srcWidth;
1108 : 0 : int tmpHeight = srcHeight;
1109 : :
1110 : 0 : if ( reqXRes > srcXRes )
1111 : : {
1112 : : // downsampling
1113 : 0 : tmpWidth = static_cast<int>( std::round( srcWidth * srcXRes / reqXRes ) );
1114 : 0 : }
1115 : 0 : if ( reqYRes > std::fabs( srcYRes ) )
1116 : : {
1117 : : // downsampling
1118 : 0 : tmpHeight = static_cast<int>( std::round( -1.*srcHeight * srcYRes / reqYRes ) );
1119 : 0 : }
1120 : :
1121 : 0 : const double tmpXMin = mExtent.xMinimum() + srcLeft * srcXRes;
1122 : 0 : const double tmpYMax = mExtent.yMaximum() + srcTop * srcYRes;
1123 : 0 : QgsDebugMsgLevel( QStringLiteral( "tmpXMin = %1 tmpYMax = %2 tmpWidth = %3 tmpHeight = %4" ).arg( tmpXMin ).arg( tmpYMax ).arg( tmpWidth ).arg( tmpHeight ), 5 );
1124 : :
1125 : : // Allocate temporary block
1126 : 0 : size_t bufferSize = dataSize * static_cast<size_t>( tmpWidth ) * static_cast<size_t>( tmpHeight );
1127 : : #ifdef Q_PROCESSOR_X86_32
1128 : : // Safety check for 32 bit systems
1129 : : qint64 _buffer_size = dataSize * static_cast<qint64>( tmpWidth ) * static_cast<qint64>( tmpHeight );
1130 : : if ( _buffer_size != static_cast<qint64>( bufferSize ) )
1131 : : {
1132 : : QgsDebugMsg( QStringLiteral( "Integer overflow calculating buffer size on a 32 bit system." ) );
1133 : : return false;
1134 : : }
1135 : : #endif
1136 : 0 : char *tmpBlock = static_cast<char *>( qgsMalloc( bufferSize ) );
1137 : 0 : if ( ! tmpBlock )
1138 : : {
1139 : 0 : QgsDebugMsgLevel( QStringLiteral( "Couldn't allocate temporary buffer of %1 bytes" ).arg( dataSize * tmpWidth * tmpHeight ), 5 );
1140 : 0 : return false;
1141 : : }
1142 : 0 : CPLErrorReset();
1143 : :
1144 : 0 : CPLErr err = gdalRasterIO( gdalBand, GF_Read,
1145 : 0 : srcLeft, srcTop, srcWidth, srcHeight,
1146 : 0 : static_cast<void *>( tmpBlock ),
1147 : 0 : tmpWidth, tmpHeight, type,
1148 : 0 : 0, 0, feedback );
1149 : :
1150 : 0 : if ( err != CPLE_None )
1151 : : {
1152 : 0 : const QString lastError = QString::fromUtf8( CPLGetLastErrorMsg() ) ;
1153 : 0 : if ( feedback )
1154 : 0 : feedback->appendError( lastError );
1155 : :
1156 : 0 : QgsLogger::warning( "RasterIO error: " + lastError );
1157 : 0 : qgsFree( tmpBlock );
1158 : 0 : return false;
1159 : 0 : }
1160 : :
1161 : 0 : const double tmpXRes = srcWidth * srcXRes / tmpWidth;
1162 : 0 : const double tmpYRes = srcHeight * srcYRes / tmpHeight; // negative
1163 : :
1164 : 0 : double y = intersectExtent.yMaximum() - 0.5 * reqYRes;
1165 : 0 : for ( int row = 0; row < tgtHeight; row++ )
1166 : : {
1167 : 0 : int tmpRow = static_cast<int>( std::floor( -1. * ( tmpYMax - y ) / tmpYRes ) );
1168 : 0 : tmpRow = std::min( tmpRow, tmpHeight - 1 );
1169 : :
1170 : 0 : char *srcRowBlock = tmpBlock + dataSize * tmpRow * tmpWidth;
1171 : 0 : char *dstRowBlock = ( char * )data + dataSize * ( tgtTop + row ) * bufferWidthPix;
1172 : :
1173 : 0 : double x = ( intersectExtent.xMinimum() + 0.5 * reqXRes - tmpXMin ) / tmpXRes; // cell center
1174 : 0 : double increment = reqXRes / tmpXRes;
1175 : :
1176 : 0 : char *dst = dstRowBlock + dataSize * tgtLeft;
1177 : 0 : char *src = srcRowBlock;
1178 : 0 : int tmpCol = 0;
1179 : 0 : int lastCol = 0;
1180 : 0 : for ( int col = 0; col < tgtWidth; ++col )
1181 : : {
1182 : : // std::floor() is quite slow! Use just cast to int.
1183 : 0 : tmpCol = static_cast<int>( x );
1184 : 0 : tmpCol = std::min( tmpCol, tmpWidth - 1 );
1185 : 0 : if ( tmpCol > lastCol )
1186 : : {
1187 : 0 : src += ( tmpCol - lastCol ) * dataSize;
1188 : 0 : lastCol = tmpCol;
1189 : 0 : }
1190 : 0 : memcpy( dst, src, dataSize );
1191 : 0 : dst += dataSize;
1192 : 0 : x += increment;
1193 : 0 : }
1194 : 0 : y -= reqYRes;
1195 : 0 : }
1196 : :
1197 : 0 : qgsFree( tmpBlock );
1198 : 0 : return true;
1199 : 0 : }
1200 : :
1201 : : /**
1202 : : * \param bandNumber the number of the band for which you want a color table
1203 : : * \param list a pointer the object that will hold the color table
1204 : : * \return TRUE if a color table was able to be read, FALSE otherwise
1205 : : */
1206 : 0 : QList<QgsColorRampShader::ColorRampItem> QgsGdalProvider::colorTable( int bandNumber )const
1207 : : {
1208 : 0 : QMutexLocker locker( mpMutex );
1209 : 0 : if ( !const_cast<QgsGdalProvider *>( this )->initIfNeeded() )
1210 : 0 : return QList<QgsColorRampShader::ColorRampItem>();
1211 : 0 : return QgsGdalProviderBase::colorTable( mGdalDataset, bandNumber );
1212 : 0 : }
1213 : :
1214 : 0 : QgsCoordinateReferenceSystem QgsGdalProvider::crs() const
1215 : : {
1216 : 0 : return mCrs;
1217 : : }
1218 : :
1219 : 0 : QgsRectangle QgsGdalProvider::extent() const
1220 : : {
1221 : : //TODO
1222 : : //mExtent = QgsGdal::extent( mGisdbase, mLocation, mMapset, mMapName, QgsGdal::Raster );
1223 : 0 : return mExtent;
1224 : : }
1225 : :
1226 : : // this is only called once when statistics are calculated
1227 : : // TODO
1228 : 0 : int QgsGdalProvider::xBlockSize() const
1229 : : {
1230 : 0 : return mXBlockSize;
1231 : : }
1232 : 0 : int QgsGdalProvider::yBlockSize() const
1233 : : {
1234 : 0 : return mYBlockSize;
1235 : : }
1236 : :
1237 : 0 : int QgsGdalProvider::xSize() const { return mWidth; }
1238 : 0 : int QgsGdalProvider::ySize() const { return mHeight; }
1239 : :
1240 : 0 : QString QgsGdalProvider::generateBandName( int bandNumber ) const
1241 : : {
1242 : 0 : QMutexLocker locker( mpMutex );
1243 : 0 : if ( !const_cast<QgsGdalProvider *>( this )->initIfNeeded() )
1244 : 0 : return QString();
1245 : :
1246 : 0 : if ( mDriverName == QLatin1String( "netCDF" ) || mDriverName == QLatin1String( "GTiff" ) )
1247 : : {
1248 : 0 : char **GDALmetadata = GDALGetMetadata( mGdalDataset, nullptr );
1249 : 0 : if ( GDALmetadata )
1250 : : {
1251 : 0 : QStringList metadata = QgsOgrUtils::cStringListToQStringList( GDALmetadata );
1252 : 0 : QStringList dimExtraValues;
1253 : 0 : QMap<QString, QString> unitsMap;
1254 : 0 : for ( QStringList::const_iterator i = metadata.constBegin(); i != metadata.constEnd(); ++i )
1255 : : {
1256 : 0 : QString val( *i );
1257 : 0 : if ( !val.startsWith( QLatin1String( "NETCDF_DIM_EXTRA" ) ) && !val.startsWith( QLatin1String( "GTIFF_DIM_EXTRA" ) ) && !val.contains( QLatin1String( "#units=" ) ) )
1258 : 0 : continue;
1259 : 0 : QStringList values = val.split( '=' );
1260 : 0 : val = values.at( 1 );
1261 : 0 : if ( values.at( 0 ) == QLatin1String( "NETCDF_DIM_EXTRA" ) || values.at( 0 ) == QLatin1String( "GTIFF_DIM_EXTRA" ) )
1262 : : {
1263 : 0 : dimExtraValues = val.replace( '{', QString() ).replace( '}', QString() ).split( ',' );
1264 : : //http://qt-project.org/doc/qt-4.8/qregexp.html#capturedTexts
1265 : 0 : }
1266 : : else
1267 : : {
1268 : 0 : unitsMap[values.at( 0 ).split( '#' ).at( 0 )] = val;
1269 : : }
1270 : 0 : }
1271 : 0 : if ( !dimExtraValues.isEmpty() )
1272 : : {
1273 : 0 : QStringList bandNameValues;
1274 : 0 : GDALRasterBandH gdalBand = GDALGetRasterBand( mGdalDataset, bandNumber );
1275 : 0 : GDALmetadata = GDALGetMetadata( gdalBand, nullptr );
1276 : 0 : if ( GDALmetadata )
1277 : : {
1278 : 0 : metadata = QgsOgrUtils::cStringListToQStringList( GDALmetadata );
1279 : 0 : for ( QStringList::const_iterator i = metadata.constBegin(); i != metadata.constEnd(); ++i )
1280 : : {
1281 : 0 : QString val( *i );
1282 : 0 : if ( !val.startsWith( QLatin1String( "NETCDF_DIM_" ) ) && !val.startsWith( QLatin1String( "GTIFF_DIM_" ) ) )
1283 : 0 : continue;
1284 : 0 : QStringList values = val.split( '=' );
1285 : 0 : for ( QStringList::const_iterator j = dimExtraValues.constBegin(); j != dimExtraValues.constEnd(); ++j )
1286 : : {
1287 : 0 : QString dim = ( *j );
1288 : 0 : if ( values.at( 0 ) != "NETCDF_DIM_" + dim && values.at( 0 ) != "GTIFF_DIM_" + dim )
1289 : 0 : continue;
1290 : 0 : if ( unitsMap.contains( dim ) && !unitsMap[dim].isEmpty() && unitsMap[dim] != QLatin1String( "none" ) )
1291 : 0 : bandNameValues.append( dim + '=' + values.at( 1 ) + " (" + unitsMap[dim] + ')' );
1292 : : else
1293 : 0 : bandNameValues.append( dim + '=' + values.at( 1 ) );
1294 : 0 : }
1295 : 0 : }
1296 : 0 : }
1297 : 0 : if ( !bandNameValues.isEmpty() )
1298 : : {
1299 : 0 : return tr( "Band" ) + QStringLiteral( " %1: %2" ).arg( bandNumber, 1 + ( int ) std::log10( ( float ) bandCount() ), 10, QChar( '0' ) ).arg( bandNameValues.join( QLatin1String( " / " ) ) );
1300 : : }
1301 : 0 : }
1302 : 0 : }
1303 : 0 : }
1304 : 0 : QString generatedBandName = QgsRasterDataProvider::generateBandName( bandNumber );
1305 : 0 : GDALRasterBandH myGdalBand = getBand( bandNumber );
1306 : 0 : if ( ! myGdalBand )
1307 : : {
1308 : 0 : QgsLogger::warning( QStringLiteral( "Band %1 does not exist." ).arg( bandNumber ) );
1309 : 0 : return QString();
1310 : : }
1311 : :
1312 : 0 : QString gdalBandName( GDALGetDescription( myGdalBand ) );
1313 : 0 : if ( !gdalBandName.isEmpty() )
1314 : : {
1315 : 0 : return generatedBandName + QStringLiteral( ": " ) + gdalBandName;
1316 : : }
1317 : 0 : return generatedBandName;
1318 : 0 : }
1319 : :
1320 : 0 : QgsRasterIdentifyResult QgsGdalProvider::identify( const QgsPointXY &point, QgsRaster::IdentifyFormat format, const QgsRectangle &boundingBox, int width, int height, int /*dpi*/ )
1321 : : {
1322 : 0 : QMutexLocker locker( mpMutex );
1323 : 0 : if ( !initIfNeeded() )
1324 : 0 : return QgsRasterIdentifyResult( ERR( tr( "Cannot read data" ) ) );
1325 : :
1326 : 0 : QgsDebugMsgLevel( QStringLiteral( "thePoint = %1 %2" ).arg( point.x(), 0, 'g', 10 ).arg( point.y(), 0, 'g', 10 ), 3 );
1327 : :
1328 : 0 : QMap<int, QVariant> results;
1329 : :
1330 : 0 : if ( format != QgsRaster::IdentifyFormatValue )
1331 : : {
1332 : 0 : return QgsRasterIdentifyResult( ERR( tr( "Format not supported" ) ) );
1333 : : }
1334 : :
1335 : 0 : if ( !extent().contains( point ) )
1336 : : {
1337 : : // Outside the raster
1338 : 0 : for ( int bandNo = 1; bandNo <= bandCount(); bandNo++ )
1339 : : {
1340 : 0 : results.insert( bandNo, QVariant() ); // null QVariant represents no data
1341 : 0 : }
1342 : 0 : return QgsRasterIdentifyResult( QgsRaster::IdentifyFormatValue, results );
1343 : : }
1344 : :
1345 : 0 : QgsRectangle finalExtent = boundingBox;
1346 : 0 : if ( finalExtent.isEmpty() )
1347 : 0 : finalExtent = extent();
1348 : :
1349 : 0 : QgsDebugMsgLevel( QStringLiteral( "myExtent = %1 " ).arg( finalExtent.toString() ), 3 );
1350 : :
1351 : 0 : if ( width == 0 )
1352 : 0 : width = xSize();
1353 : 0 : if ( height == 0 )
1354 : 0 : height = ySize();
1355 : :
1356 : 0 : QgsDebugMsgLevel( QStringLiteral( "theWidth = %1 height = %2" ).arg( width ).arg( height ), 3 );
1357 : :
1358 : : // Calculate the row / column where the point falls
1359 : 0 : double xres = ( finalExtent.width() ) / width;
1360 : 0 : double yres = ( finalExtent.height() ) / height;
1361 : :
1362 : : // Offset, not the cell index -> std::floor
1363 : 0 : int col = static_cast< int >( std::floor( ( point.x() - finalExtent.xMinimum() ) / xres ) );
1364 : 0 : int row = static_cast< int >( std::floor( ( finalExtent.yMaximum() - point.y() ) / yres ) );
1365 : :
1366 : 0 : QgsDebugMsgLevel( QStringLiteral( "row = %1 col = %2" ).arg( row ).arg( col ), 3 );
1367 : :
1368 : 0 : int r = 0;
1369 : 0 : int c = 0;
1370 : 0 : int w = 1;
1371 : 0 : int h = 1;
1372 : :
1373 : 0 : double xMin = finalExtent.xMinimum() + col * xres;
1374 : 0 : double xMax = xMin + xres * w;
1375 : 0 : double yMax = finalExtent.yMaximum() - row * yres;
1376 : 0 : double yMin = yMax - yres * h;
1377 : 0 : QgsRectangle pixelExtent( xMin, yMin, xMax, yMax );
1378 : :
1379 : 0 : for ( int i = 1; i <= bandCount(); i++ )
1380 : : {
1381 : 0 : std::unique_ptr< QgsRasterBlock > bandBlock( block( i, pixelExtent, w, h ) );
1382 : :
1383 : 0 : if ( !bandBlock )
1384 : : {
1385 : 0 : return QgsRasterIdentifyResult( ERR( tr( "Cannot read data" ) ) );
1386 : : }
1387 : :
1388 : 0 : double value = bandBlock->value( r, c );
1389 : :
1390 : 0 : if ( ( sourceHasNoDataValue( i ) && useSourceNoDataValue( i ) &&
1391 : 0 : ( std::isnan( value ) || qgsDoubleNear( value, sourceNoDataValue( i ) ) ) ) ||
1392 : 0 : ( QgsRasterRange::contains( value, userNoDataValues( i ) ) ) )
1393 : : {
1394 : 0 : results.insert( i, QVariant() ); // null QVariant represents no data
1395 : 0 : }
1396 : : else
1397 : : {
1398 : 0 : if ( sourceDataType( i ) == Qgis::Float32 )
1399 : : {
1400 : : // Insert a float QVariant so that QgsMapToolIdentify::identifyRasterLayer()
1401 : : // can print a string without an excessive precision
1402 : 0 : results.insert( i, static_cast<float>( value ) );
1403 : 0 : }
1404 : : else
1405 : 0 : results.insert( i, value );
1406 : : }
1407 : 0 : }
1408 : 0 : return QgsRasterIdentifyResult( QgsRaster::IdentifyFormatValue, results );
1409 : 0 : }
1410 : :
1411 : 0 : bool QgsGdalProvider::worldToPixel( double x, double y, int &col, int &row ) const
1412 : : {
1413 : : /*
1414 : : * This set of equations solves
1415 : : * Xgeo = GT(0) + XpixelGT(1) + YlineGT(2)
1416 : : * Ygeo = GT(3) + XpixelGT(4) + YlineGT(5)
1417 : : * when Xgeo, Ygeo are given
1418 : : */
1419 : 0 : double div = ( mGeoTransform[2] * mGeoTransform[4] - mGeoTransform[1] * mGeoTransform[5] );
1420 : 0 : if ( div < 2 * std::numeric_limits<double>::epsilon() )
1421 : 0 : return false;
1422 : 0 : double doubleCol = -( mGeoTransform[2] * ( mGeoTransform[3] - y ) + mGeoTransform[5] * ( x - mGeoTransform[0] ) ) / div;
1423 : 0 : double doubleRow = ( mGeoTransform[1] * ( mGeoTransform[3] - y ) + mGeoTransform[4] * ( x - mGeoTransform[0] ) ) / div;
1424 : : // note: we truncate here, not round, otherwise values will be 0.5 pixels off
1425 : 0 : col = static_cast< int >( doubleCol );
1426 : 0 : row = static_cast< int >( doubleRow );
1427 : 0 : return true;
1428 : 0 : }
1429 : :
1430 : 0 : double QgsGdalProvider::sample( const QgsPointXY &point, int band, bool *ok, const QgsRectangle &, int, int, int )
1431 : : {
1432 : 0 : if ( ok )
1433 : 0 : *ok = false;
1434 : :
1435 : 0 : QMutexLocker locker( mpMutex );
1436 : 0 : if ( !initIfNeeded() )
1437 : 0 : return std::numeric_limits<double>::quiet_NaN();
1438 : :
1439 : 0 : if ( !extent().contains( point ) )
1440 : : {
1441 : : // Outside the raster
1442 : 0 : return std::numeric_limits<double>::quiet_NaN();
1443 : : }
1444 : :
1445 : 0 : GDALRasterBandH hBand = GDALGetRasterBand( mGdalDataset, band );
1446 : 0 : if ( !hBand )
1447 : 0 : return std::numeric_limits<double>::quiet_NaN();
1448 : :
1449 : : int row;
1450 : : int col;
1451 : 0 : if ( !worldToPixel( point.x(), point.y(), col, row ) )
1452 : 0 : return std::numeric_limits<double>::quiet_NaN();
1453 : :
1454 : 0 : float value = 0;
1455 : 0 : CPLErr err = GDALRasterIO( hBand, GF_Read, col, row, 1, 1,
1456 : 0 : &value, 1, 1, GDT_Float32, 0, 0 );
1457 : 0 : if ( err != CE_None )
1458 : 0 : return std::numeric_limits<double>::quiet_NaN();
1459 : :
1460 : 0 : if ( ( sourceHasNoDataValue( band ) && useSourceNoDataValue( band ) &&
1461 : 0 : ( std::isnan( value ) || qgsDoubleNear( static_cast< double >( value ), sourceNoDataValue( band ) ) ) ) ||
1462 : 0 : ( QgsRasterRange::contains( static_cast< double >( value ), userNoDataValues( band ) ) ) )
1463 : : {
1464 : 0 : return std::numeric_limits<double>::quiet_NaN();
1465 : : }
1466 : :
1467 : 0 : if ( ok )
1468 : 0 : *ok = true;
1469 : :
1470 : 0 : return static_cast< double >( value ) * bandScale( band ) + bandOffset( band );
1471 : 0 : }
1472 : :
1473 : 0 : int QgsGdalProvider::capabilities() const
1474 : : {
1475 : 0 : QMutexLocker locker( mpMutex );
1476 : 0 : if ( !const_cast<QgsGdalProvider *>( this )->initIfNeeded() )
1477 : 0 : return 0;
1478 : :
1479 : 0 : int capability = QgsRasterDataProvider::Identify
1480 : : | QgsRasterDataProvider::IdentifyValue
1481 : : | QgsRasterDataProvider::Size
1482 : : | QgsRasterDataProvider::BuildPyramids
1483 : : | QgsRasterDataProvider::Create
1484 : : | QgsRasterDataProvider::Remove
1485 : : | QgsRasterDataProvider::Prefetch;
1486 : 0 : if ( mDriverName != QLatin1String( "WMS" ) )
1487 : : {
1488 : 0 : capability |= QgsRasterDataProvider::Size;
1489 : 0 : }
1490 : 0 : return capability;
1491 : 0 : }
1492 : :
1493 : 0 : Qgis::DataType QgsGdalProvider::sourceDataType( int bandNo ) const
1494 : : {
1495 : 0 : QMutexLocker locker( mpMutex );
1496 : 0 : if ( !const_cast<QgsGdalProvider *>( this )->initIfNeeded() )
1497 : 0 : return dataTypeFromGdal( GDT_Byte );
1498 : :
1499 : 0 : if ( mMaskBandExposedAsAlpha && bandNo == GDALGetRasterCount( mGdalDataset ) + 1 )
1500 : 0 : return dataTypeFromGdal( GDT_Byte );
1501 : :
1502 : 0 : GDALRasterBandH myGdalBand = GDALGetRasterBand( mGdalDataset, bandNo );
1503 : 0 : GDALDataType myGdalDataType = GDALGetRasterDataType( myGdalBand );
1504 : 0 : Qgis::DataType myDataType = dataTypeFromGdal( myGdalDataType );
1505 : :
1506 : : // define if the band has scale and offset to apply
1507 : 0 : double myScale = bandScale( bandNo );
1508 : 0 : double myOffset = bandOffset( bandNo );
1509 : 0 : if ( myScale != 1.0 || myOffset != 0.0 )
1510 : : {
1511 : : // if the band has scale or offset to apply change dataType
1512 : 0 : switch ( myDataType )
1513 : : {
1514 : : case Qgis::UnknownDataType:
1515 : : case Qgis::ARGB32:
1516 : : case Qgis::ARGB32_Premultiplied:
1517 : 0 : return myDataType;
1518 : : case Qgis::Byte:
1519 : : case Qgis::UInt16:
1520 : : case Qgis::Int16:
1521 : : case Qgis::UInt32:
1522 : : case Qgis::Int32:
1523 : : case Qgis::Float32:
1524 : : case Qgis::CInt16:
1525 : 0 : myDataType = Qgis::Float32;
1526 : 0 : break;
1527 : : case Qgis::Float64:
1528 : : case Qgis::CInt32:
1529 : : case Qgis::CFloat32:
1530 : 0 : myDataType = Qgis::Float64;
1531 : 0 : break;
1532 : : case Qgis::CFloat64:
1533 : 0 : return myDataType;
1534 : : }
1535 : 0 : }
1536 : 0 : return myDataType;
1537 : 0 : }
1538 : :
1539 : 0 : Qgis::DataType QgsGdalProvider::dataType( int bandNo ) const
1540 : : {
1541 : 0 : if ( mMaskBandExposedAsAlpha && bandNo == mBandCount )
1542 : 0 : return dataTypeFromGdal( GDT_Byte );
1543 : :
1544 : 0 : if ( bandNo <= 0 || bandNo > mGdalDataType.count() ) return Qgis::UnknownDataType;
1545 : :
1546 : 0 : return dataTypeFromGdal( mGdalDataType[bandNo - 1] );
1547 : 0 : }
1548 : :
1549 : 0 : double QgsGdalProvider::bandScale( int bandNo ) const
1550 : : {
1551 : 0 : QMutexLocker locker( mpMutex );
1552 : 0 : if ( !const_cast<QgsGdalProvider *>( this )->initIfNeeded() )
1553 : 0 : return 1.0;
1554 : :
1555 : 0 : GDALRasterBandH myGdalBand = getBand( bandNo );
1556 : : int bGotScale;
1557 : 0 : double myScale = GDALGetRasterScale( myGdalBand, &bGotScale );
1558 : :
1559 : : // if scale==0, ignore both scale and offset
1560 : 0 : if ( bGotScale && !qgsDoubleNear( myScale, 0.0 ) )
1561 : 0 : return myScale;
1562 : : else
1563 : 0 : return 1.0;
1564 : 0 : }
1565 : :
1566 : 0 : double QgsGdalProvider::bandOffset( int bandNo ) const
1567 : : {
1568 : 0 : QMutexLocker locker( mpMutex );
1569 : 0 : if ( !const_cast<QgsGdalProvider *>( this )->initIfNeeded() )
1570 : 0 : return 0.0;
1571 : :
1572 : 0 : GDALRasterBandH myGdalBand = getBand( bandNo );
1573 : :
1574 : : // if scale==0, ignore both scale and offset
1575 : : int bGotScale;
1576 : 0 : double myScale = GDALGetRasterScale( myGdalBand, &bGotScale );
1577 : 0 : if ( bGotScale && qgsDoubleNear( myScale, 0.0 ) )
1578 : 0 : return 0.0;
1579 : :
1580 : : int bGotOffset;
1581 : 0 : double myOffset = GDALGetRasterOffset( myGdalBand, &bGotOffset );
1582 : 0 : if ( bGotOffset )
1583 : 0 : return myOffset;
1584 : : else
1585 : 0 : return 0.0;
1586 : 0 : }
1587 : :
1588 : 0 : int QgsGdalProvider::bandCount() const
1589 : : {
1590 : 0 : return mBandCount;
1591 : : }
1592 : :
1593 : 0 : int QgsGdalProvider::colorInterpretation( int bandNo ) const
1594 : : {
1595 : 0 : QMutexLocker locker( mpMutex );
1596 : 0 : if ( !const_cast<QgsGdalProvider *>( this )->initIfNeeded() )
1597 : 0 : return colorInterpretationFromGdal( GCI_Undefined );
1598 : :
1599 : 0 : if ( mMaskBandExposedAsAlpha && bandNo == GDALGetRasterCount( mGdalDataset ) + 1 )
1600 : 0 : return colorInterpretationFromGdal( GCI_AlphaBand );
1601 : 0 : GDALRasterBandH myGdalBand = GDALGetRasterBand( mGdalDataset, bandNo );
1602 : 0 : return colorInterpretationFromGdal( GDALGetRasterColorInterpretation( myGdalBand ) );
1603 : 0 : }
1604 : :
1605 : 0 : bool QgsGdalProvider::isValid() const
1606 : : {
1607 : 0 : QgsDebugMsgLevel( QStringLiteral( "valid = %1" ).arg( mValid ), 4 );
1608 : 0 : return mValid;
1609 : : }
1610 : :
1611 : 0 : QString QgsGdalProvider::lastErrorTitle()
1612 : : {
1613 : 0 : return QStringLiteral( "Not implemented" );
1614 : : }
1615 : :
1616 : 0 : QString QgsGdalProvider::lastError()
1617 : : {
1618 : 0 : return QStringLiteral( "Not implemented" );
1619 : : }
1620 : :
1621 : 0 : QString QgsGdalProvider::name() const
1622 : : {
1623 : 0 : return PROVIDER_KEY;
1624 : : }
1625 : :
1626 : 3 : QString QgsGdalProvider::providerKey()
1627 : : {
1628 : 6 : return PROVIDER_KEY;
1629 : : };
1630 : :
1631 : 0 : QString QgsGdalProvider::description() const
1632 : : {
1633 : 0 : return PROVIDER_DESCRIPTION;
1634 : : }
1635 : :
1636 : 0 : QgsRasterDataProvider::ProviderCapabilities QgsGdalProvider::providerCapabilities() const
1637 : : {
1638 : 0 : return ProviderCapability::ProviderHintBenefitsFromResampling |
1639 : 0 : ProviderCapability::ProviderHintCanPerformProviderResampling |
1640 : : ProviderCapability::ReloadData;
1641 : : }
1642 : :
1643 : : // This is used also by global isValidRasterFileName
1644 : 0 : QStringList QgsGdalProvider::subLayers( GDALDatasetH dataset )
1645 : : {
1646 : 0 : QStringList subLayers;
1647 : :
1648 : 0 : if ( !dataset )
1649 : : {
1650 : 0 : QgsDebugMsg( QStringLiteral( "dataset is nullptr" ) );
1651 : 0 : return subLayers;
1652 : : }
1653 : :
1654 : 0 : char **metadata = GDALGetMetadata( dataset, "SUBDATASETS" );
1655 : :
1656 : 0 : if ( metadata )
1657 : : {
1658 : 0 : const int subdatasetCount = CSLCount( metadata ) / 2;
1659 : 0 : for ( int i = 1; i <= subdatasetCount; i++ )
1660 : : {
1661 : 0 : const char *name = CSLFetchNameValue( metadata, CPLSPrintf( "SUBDATASET_%d_NAME", i ) );
1662 : 0 : const char *desc = CSLFetchNameValue( metadata, CPLSPrintf( "SUBDATASET_%d_DESC", i ) );
1663 : 0 : if ( name && desc )
1664 : : {
1665 : : // For GeoPackage, desc is often "table - identifier" where table=identifier
1666 : : // In that case, just keep one.
1667 : 0 : const char *sepPtr = strstr( desc, " - " );
1668 : 0 : if ( sepPtr && strncmp( desc, sepPtr + 3, strlen( sepPtr + 3 ) ) == 0 )
1669 : : {
1670 : 0 : desc = sepPtr + 3;
1671 : 0 : }
1672 : 0 : subLayers << QString::fromUtf8( name ) + QgsDataProvider::sublayerSeparator() + QString::fromUtf8( desc );
1673 : 0 : }
1674 : 0 : }
1675 : 0 : }
1676 : :
1677 : 0 : if ( !subLayers.isEmpty() )
1678 : : {
1679 : 0 : QgsDebugMsgLevel( "sublayers:\n " + subLayers.join( "\n " ), 3 );
1680 : 0 : }
1681 : :
1682 : 0 : return subLayers;
1683 : 0 : }
1684 : :
1685 : 0 : bool QgsGdalProvider::hasHistogram( int bandNo,
1686 : : int binCount,
1687 : : double minimum, double maximum,
1688 : : const QgsRectangle &boundingBox,
1689 : : int sampleSize,
1690 : : bool includeOutOfRange )
1691 : : {
1692 : 0 : QMutexLocker locker( mpMutex );
1693 : 0 : if ( !initIfNeeded() )
1694 : 0 : return false;
1695 : :
1696 : 0 : QgsDebugMsgLevel( QStringLiteral( "theBandNo = %1 binCount = %2 minimum = %3 maximum = %4 sampleSize = %5" ).arg( bandNo ).arg( binCount ).arg( minimum ).arg( maximum ).arg( sampleSize ), 3 );
1697 : :
1698 : : // First check if cached in mHistograms
1699 : 0 : if ( QgsRasterDataProvider::hasHistogram( bandNo, binCount, minimum, maximum, boundingBox, sampleSize, includeOutOfRange ) )
1700 : : {
1701 : 0 : return true;
1702 : : }
1703 : :
1704 : 0 : QgsRasterHistogram myHistogram;
1705 : 0 : initHistogram( myHistogram, bandNo, binCount, minimum, maximum, boundingBox, sampleSize, includeOutOfRange );
1706 : :
1707 : : // If not cached, check if supported by GDAL
1708 : 0 : if ( myHistogram.extent != extent() )
1709 : : {
1710 : 0 : QgsDebugMsgLevel( QStringLiteral( "Not supported by GDAL." ), 2 );
1711 : 0 : return false;
1712 : : }
1713 : :
1714 : 0 : if ( ( sourceHasNoDataValue( bandNo ) && !useSourceNoDataValue( bandNo ) ) ||
1715 : 0 : !userNoDataValues( bandNo ).isEmpty() )
1716 : : {
1717 : 0 : QgsDebugMsgLevel( QStringLiteral( "Custom no data values -> GDAL histogram not sufficient." ), 3 );
1718 : 0 : return false;
1719 : : }
1720 : :
1721 : 0 : QgsDebugMsgLevel( QStringLiteral( "Looking for GDAL histogram" ), 4 );
1722 : :
1723 : 0 : GDALRasterBandH myGdalBand = getBand( bandNo );
1724 : 0 : if ( ! myGdalBand )
1725 : : {
1726 : 0 : return false;
1727 : : }
1728 : :
1729 : : // get default histogram with force=false to see if there is a cached histo
1730 : : double myMinVal, myMaxVal;
1731 : : int myBinCount;
1732 : :
1733 : 0 : GUIntBig *myHistogramArray = nullptr;
1734 : 0 : CPLErr myError = GDALGetDefaultHistogramEx( myGdalBand, &myMinVal, &myMaxVal,
1735 : : &myBinCount, &myHistogramArray, false,
1736 : : nullptr, nullptr );
1737 : :
1738 : 0 : if ( myHistogramArray )
1739 : 0 : VSIFree( myHistogramArray ); // use VSIFree because allocated by GDAL
1740 : :
1741 : : // if there was any error/warning assume the histogram is not valid or non-existent
1742 : 0 : if ( myError != CE_None )
1743 : : {
1744 : 0 : QgsDebugMsgLevel( QStringLiteral( "Cannot get default GDAL histogram" ), 2 );
1745 : 0 : return false;
1746 : : }
1747 : :
1748 : : // This is fragile
1749 : 0 : double myExpectedMinVal = myHistogram.minimum;
1750 : 0 : double myExpectedMaxVal = myHistogram.maximum;
1751 : :
1752 : 0 : double dfHalfBucket = ( myExpectedMaxVal - myExpectedMinVal ) / ( 2 * myHistogram.binCount );
1753 : 0 : myExpectedMinVal -= dfHalfBucket;
1754 : 0 : myExpectedMaxVal += dfHalfBucket;
1755 : :
1756 : : // min/max are stored as text in aux file => use threshold
1757 : 0 : if ( myBinCount != myHistogram.binCount ||
1758 : 0 : std::fabs( myMinVal - myExpectedMinVal ) > std::fabs( myExpectedMinVal ) / 10e6 ||
1759 : 0 : std::fabs( myMaxVal - myExpectedMaxVal ) > std::fabs( myExpectedMaxVal ) / 10e6 )
1760 : : {
1761 : 0 : QgsDebugMsgLevel( QStringLiteral( "Params do not match binCount: %1 x %2, minVal: %3 x %4, maxVal: %5 x %6" ).arg( myBinCount ).arg( myHistogram.binCount ).arg( myMinVal ).arg( myExpectedMinVal ).arg( myMaxVal ).arg( myExpectedMaxVal ), 2 );
1762 : 0 : return false;
1763 : : }
1764 : :
1765 : 0 : QgsDebugMsgLevel( QStringLiteral( "GDAL has cached histogram" ), 2 );
1766 : :
1767 : : // This should be enough, possible call to histogram() should retrieve the histogram cached in GDAL
1768 : :
1769 : 0 : return true;
1770 : 0 : }
1771 : :
1772 : 0 : QgsRasterHistogram QgsGdalProvider::histogram( int bandNo,
1773 : : int binCount,
1774 : : double minimum, double maximum,
1775 : : const QgsRectangle &boundingBox,
1776 : : int sampleSize,
1777 : : bool includeOutOfRange, QgsRasterBlockFeedback *feedback )
1778 : : {
1779 : 0 : QMutexLocker locker( mpMutex );
1780 : 0 : if ( !initIfNeeded() )
1781 : 0 : return QgsRasterHistogram();
1782 : :
1783 : 0 : QgsDebugMsgLevel( QStringLiteral( "theBandNo = %1 binCount = %2 minimum = %3 maximum = %4 sampleSize = %5" ).arg( bandNo ).arg( binCount ).arg( minimum ).arg( maximum ).arg( sampleSize ), 2 );
1784 : :
1785 : 0 : QgsRasterHistogram myHistogram;
1786 : 0 : initHistogram( myHistogram, bandNo, binCount, minimum, maximum, boundingBox, sampleSize, includeOutOfRange );
1787 : :
1788 : : // Find cached
1789 : 0 : const auto constMHistograms = mHistograms;
1790 : 0 : for ( const QgsRasterHistogram &histogram : constMHistograms )
1791 : : {
1792 : 0 : if ( histogram == myHistogram )
1793 : : {
1794 : 0 : QgsDebugMsgLevel( QStringLiteral( "Using cached histogram." ), 2 );
1795 : 0 : return histogram;
1796 : : }
1797 : : }
1798 : :
1799 : 0 : if ( ( sourceHasNoDataValue( bandNo ) && !useSourceNoDataValue( bandNo ) ) ||
1800 : 0 : !userNoDataValues( bandNo ).isEmpty() )
1801 : : {
1802 : 0 : QgsDebugMsgLevel( QStringLiteral( "Custom no data values, using generic histogram." ), 2 );
1803 : 0 : return QgsRasterDataProvider::histogram( bandNo, binCount, minimum, maximum, boundingBox, sampleSize, includeOutOfRange, feedback );
1804 : : }
1805 : :
1806 : 0 : if ( myHistogram.extent != extent() )
1807 : : {
1808 : 0 : QgsDebugMsgLevel( QStringLiteral( "Not full extent, using generic histogram." ), 2 );
1809 : 0 : return QgsRasterDataProvider::histogram( bandNo, binCount, minimum, maximum, boundingBox, sampleSize, includeOutOfRange, feedback );
1810 : : }
1811 : :
1812 : 0 : QgsDebugMsgLevel( QStringLiteral( "Computing GDAL histogram" ), 2 );
1813 : :
1814 : 0 : GDALRasterBandH myGdalBand = getBand( bandNo );
1815 : :
1816 : 0 : int bApproxOK = false;
1817 : 0 : if ( sampleSize > 0 )
1818 : : {
1819 : : // cast to double, integer could overflow
1820 : 0 : if ( ( static_cast<double>( xSize() ) * static_cast<double>( ySize() ) / sampleSize ) > 2 ) // not perfect
1821 : : {
1822 : 0 : QgsDebugMsgLevel( QStringLiteral( "Approx" ), 2 );
1823 : 0 : bApproxOK = true;
1824 : 0 : }
1825 : 0 : }
1826 : :
1827 : 0 : QgsDebugMsgLevel( QStringLiteral( "xSize() = %1 ySize() = %2 sampleSize = %3 bApproxOK = %4" ).arg( xSize() ).arg( ySize() ).arg( sampleSize ).arg( bApproxOK ), 2 );
1828 : :
1829 : 0 : QgsGdalProgress myProg;
1830 : 0 : myProg.type = QgsRaster::ProgressHistogram;
1831 : 0 : myProg.provider = this;
1832 : 0 : myProg.feedback = feedback;
1833 : :
1834 : : #if 0 // this is the old method
1835 : :
1836 : : double myerval = ( bandStats.maximumValue - bandStats.minimumValue ) / binCount;
1837 : : GDALGetRasterHistogram( myGdalBand, bandStats.minimumValue - 0.1 * myerval,
1838 : : bandStats.maximumValue + 0.1 * myerval, binCount, myHistogramArray,
1839 : : ignoreOutOfRangeFlag, histogramEstimatedFlag, progressCallback,
1840 : : &myProg ); //this is the arg for our custom gdal progress callback
1841 : :
1842 : : #else // this is the new method, which gets a "Default" histogram
1843 : :
1844 : : // calculate min/max like in GDALRasterBand::GetDefaultHistogram, but don't call it directly
1845 : : // because there is no bApproxOK argument - that is lacking from the API
1846 : :
1847 : : // Min/max, if not specified, are set by histogramDefaults, it does not
1848 : : // set however min/max shifted to avoid rounding errors
1849 : :
1850 : 0 : double myMinVal = myHistogram.minimum;
1851 : 0 : double myMaxVal = myHistogram.maximum;
1852 : :
1853 : : // unapply scale anf offset for min and max
1854 : 0 : double myScale = bandScale( bandNo );
1855 : 0 : double myOffset = bandOffset( bandNo );
1856 : 0 : if ( myScale != 1.0 || myOffset != 0. )
1857 : : {
1858 : 0 : myMinVal = ( myHistogram.minimum - myOffset ) / myScale;
1859 : 0 : myMaxVal = ( myHistogram.maximum - myOffset ) / myScale;
1860 : 0 : }
1861 : :
1862 : 0 : double dfHalfBucket = ( myMaxVal - myMinVal ) / ( 2 * myHistogram.binCount );
1863 : 0 : myMinVal -= dfHalfBucket;
1864 : 0 : myMaxVal += dfHalfBucket;
1865 : :
1866 : : #if 0
1867 : : const char *pszPixelType = GDALGetMetadataItem( myGdalBand, "PIXELTYPE", "IMAGE_STRUCTURE" );
1868 : : int bSignedByte = ( pszPixelType && EQUAL( pszPixelType, "SIGNEDBYTE" ) );
1869 : :
1870 : : if ( GDALGetRasterDataType( myGdalBand ) == GDT_Byte && !bSignedByte )
1871 : : {
1872 : : myMinVal = -0.5;
1873 : : myMaxVal = 255.5;
1874 : : }
1875 : : else
1876 : : {
1877 : : CPLErr eErr = CE_Failure;
1878 : : double dfHalfBucket = 0;
1879 : : eErr = GDALGetRasterStatistics( myGdalBand, true, true, &myMinVal, &myMaxVal, nullptr, nullptr );
1880 : : if ( eErr != CE_None )
1881 : : {
1882 : : delete [] myHistogramArray;
1883 : : return;
1884 : : }
1885 : : dfHalfBucket = ( myMaxVal - myMinVal ) / ( 2 * binCount );
1886 : : myMinVal -= dfHalfBucket;
1887 : : myMaxVal += dfHalfBucket;
1888 : : }
1889 : : #endif
1890 : :
1891 : 0 : GUIntBig *myHistogramArray = new GUIntBig[myHistogram.binCount];
1892 : 0 : CPLErr myError = GDALGetRasterHistogramEx( myGdalBand, myMinVal, myMaxVal,
1893 : 0 : myHistogram.binCount, myHistogramArray,
1894 : 0 : includeOutOfRange, bApproxOK, progressCallback,
1895 : 0 : &myProg ); //this is the arg for our custom gdal progress callback
1896 : :
1897 : 0 : if ( myError != CE_None || ( feedback && feedback->isCanceled() ) )
1898 : : {
1899 : 0 : QgsDebugMsgLevel( QStringLiteral( "Cannot get histogram" ), 2 );
1900 : 0 : delete [] myHistogramArray;
1901 : 0 : return myHistogram;
1902 : : }
1903 : :
1904 : : #endif
1905 : :
1906 : 0 : for ( int myBin = 0; myBin < myHistogram.binCount; myBin++ )
1907 : : {
1908 : 0 : myHistogram.histogramVector.push_back( myHistogramArray[myBin] );
1909 : 0 : myHistogram.nonNullCount += myHistogramArray[myBin];
1910 : 0 : }
1911 : :
1912 : 0 : myHistogram.valid = true;
1913 : :
1914 : 0 : delete [] myHistogramArray;
1915 : :
1916 : 0 : QgsDebugMsgLevel( ">>>>> Histogram vector now contains " + QString::number( myHistogram.histogramVector.size() ) + " elements", 3 );
1917 : :
1918 : 0 : mHistograms.append( myHistogram );
1919 : 0 : return myHistogram;
1920 : 0 : }
1921 : :
1922 : : /*
1923 : : * This will speed up performance at the expense of hard drive space.
1924 : : * Also, write access to the file is required for creating internal pyramids,
1925 : : * and to the directory in which the files exists if external
1926 : : * pyramids (.ovr) are to be created. If no parameter is passed in
1927 : : * it will default to nearest neighbor resampling.
1928 : : *
1929 : : * \param tryInternalFlag - Try to make the pyramids internal if supported (e.g. geotiff). If not supported it will revert to creating external .ovr file anyway.
1930 : : * \return null string on success, otherwise a string specifying error
1931 : : */
1932 : 0 : QString QgsGdalProvider::buildPyramids( const QList<QgsRasterPyramid> &rasterPyramidList,
1933 : : const QString &resamplingMethod, QgsRaster::RasterPyramidsFormat format,
1934 : : const QStringList &configOptions, QgsRasterBlockFeedback *feedback )
1935 : : {
1936 : 0 : QMutexLocker locker( mpMutex );
1937 : :
1938 : : //TODO: Consider making rasterPyramidList modifiable by this method to indicate if the pyramid exists after build attempt
1939 : : //without requiring the user to rebuild the pyramid list to get the updated information
1940 : :
1941 : : //
1942 : : // Note: Make sure the raster is not opened in write mode
1943 : : // in order to force overviews to be written to a separate file.
1944 : : // Otherwise reoopen it in read/write mode to stick overviews
1945 : : // into the same file (if supported)
1946 : : //
1947 : :
1948 : 0 : if ( mGdalDataset != mGdalBaseDataset )
1949 : : {
1950 : 0 : QgsLogger::warning( QStringLiteral( "Pyramid building not currently supported for 'warped virtual dataset'." ) );
1951 : 0 : return QStringLiteral( "ERROR_VIRTUAL" );
1952 : : }
1953 : :
1954 : : // check if building internally
1955 : 0 : if ( format == QgsRaster::PyramidsInternal )
1956 : : {
1957 : :
1958 : : // test if the file is writable
1959 : : //QFileInfo myQFile( mDataSource );
1960 : 0 : QFileInfo myQFile( dataSourceUri( true ) );
1961 : :
1962 : 0 : if ( !myQFile.isWritable() )
1963 : : {
1964 : 0 : return QStringLiteral( "ERROR_WRITE_ACCESS" );
1965 : : }
1966 : :
1967 : : // if needed close the gdal dataset and reopen it in read / write mode
1968 : : // TODO this doesn't seem to work anymore... must fix it before 2.0!!!
1969 : : // no errors are reported, but pyramids are not present in file.
1970 : 0 : if ( GDALGetAccess( mGdalDataset ) == GA_ReadOnly )
1971 : : {
1972 : 0 : QgsDebugMsg( QStringLiteral( "re-opening the dataset in read/write mode" ) );
1973 : 0 : GDALClose( mGdalDataset );
1974 : : //mGdalBaseDataset = GDALOpen( QFile::encodeName( dataSourceUri() ).constData(), GA_Update );
1975 : :
1976 : 0 : mGdalBaseDataset = gdalOpen( dataSourceUri( true ), GDAL_OF_UPDATE );
1977 : :
1978 : : // if the dataset couldn't be opened in read / write mode, tell the user
1979 : 0 : if ( !mGdalBaseDataset )
1980 : : {
1981 : 0 : mGdalBaseDataset = gdalOpen( dataSourceUri( true ), GDAL_OF_READONLY );
1982 : : //Since we are not a virtual warped dataset, mGdalDataSet and mGdalBaseDataset are supposed to be the same
1983 : 0 : mGdalDataset = mGdalBaseDataset;
1984 : 0 : return QStringLiteral( "ERROR_WRITE_FORMAT" );
1985 : : }
1986 : 0 : }
1987 : 0 : }
1988 : :
1989 : : // are we using Erdas Imagine external overviews?
1990 : 0 : QgsStringMap myConfigOptionsOld;
1991 : 0 : myConfigOptionsOld[ QStringLiteral( "USE_RRD" )] = CPLGetConfigOption( "USE_RRD", "NO" );
1992 : 0 : myConfigOptionsOld[ QStringLiteral( "TIFF_USE_OVR" )] = CPLGetConfigOption( "TIFF_USE_OVR", "NO" );
1993 : 0 : if ( format == QgsRaster::PyramidsErdas )
1994 : 0 : CPLSetConfigOption( "USE_RRD", "YES" );
1995 : : else
1996 : : {
1997 : 0 : CPLSetConfigOption( "USE_RRD", "NO" );
1998 : 0 : if ( format == QgsRaster::PyramidsGTiff )
1999 : : {
2000 : 0 : CPLSetConfigOption( "TIFF_USE_OVR", "YES" );
2001 : 0 : }
2002 : : }
2003 : :
2004 : : // add any driver-specific configuration options, save values to be restored later
2005 : 0 : if ( format != QgsRaster::PyramidsErdas && ! configOptions.isEmpty() )
2006 : : {
2007 : 0 : const auto constConfigOptions = configOptions;
2008 : 0 : for ( const QString &option : constConfigOptions )
2009 : : {
2010 : 0 : QStringList opt = option.split( '=' );
2011 : 0 : if ( opt.size() == 2 )
2012 : : {
2013 : 0 : QByteArray key = opt[0].toLocal8Bit();
2014 : 0 : QByteArray value = opt[1].toLocal8Bit();
2015 : : // save previous value
2016 : 0 : myConfigOptionsOld[ opt[0] ] = QString( CPLGetConfigOption( key.data(), nullptr ) );
2017 : : // set temp. value
2018 : 0 : CPLSetConfigOption( key.data(), value.data() );
2019 : 0 : QgsDebugMsgLevel( QStringLiteral( "set option %1=%2" ).arg( key.data(), value.data() ), 2 );
2020 : 0 : }
2021 : : else
2022 : : {
2023 : 0 : QgsDebugMsg( QStringLiteral( "invalid pyramid option: %1" ).arg( option ) );
2024 : : }
2025 : 0 : }
2026 : 0 : }
2027 : :
2028 : : //
2029 : : // Iterate through the Raster Layer Pyramid Vector, building any pyramid
2030 : : // marked as exists in each RasterPyramid struct.
2031 : : //
2032 : : CPLErr myError; //in case anything fails
2033 : :
2034 : 0 : QVector<int> myOverviewLevelsVector;
2035 : 0 : QList<QgsRasterPyramid>::const_iterator myRasterPyramidIterator;
2036 : 0 : for ( myRasterPyramidIterator = rasterPyramidList.begin();
2037 : 0 : myRasterPyramidIterator != rasterPyramidList.end();
2038 : 0 : ++myRasterPyramidIterator )
2039 : : {
2040 : : #ifdef QGISDEBUG
2041 : : QgsDebugMsgLevel( QStringLiteral( "Build pyramids:: Level %1" ).arg( myRasterPyramidIterator->getLevel() ), 2 );
2042 : : QgsDebugMsgLevel( QStringLiteral( "x:%1" ).arg( myRasterPyramidIterator->getXDim() ), 2 );
2043 : : QgsDebugMsgLevel( QStringLiteral( "y:%1" ).arg( myRasterPyramidIterator->getYDim() ), 2 );
2044 : : QgsDebugMsgLevel( QStringLiteral( "exists : %1" ).arg( myRasterPyramidIterator->getExists() ), 2 );
2045 : : #endif
2046 : 0 : if ( myRasterPyramidIterator->getBuild() )
2047 : : {
2048 : 0 : QgsDebugMsgLevel( QStringLiteral( "adding overview at level %1 to list"
2049 : : ).arg( myRasterPyramidIterator->getLevel() ), 2 );
2050 : 0 : myOverviewLevelsVector.append( myRasterPyramidIterator->getLevel() );
2051 : 0 : }
2052 : 0 : }
2053 : : /* From : http://www.gdal.org/classGDALDataset.html#a2aa6f88b3bbc840a5696236af11dde15
2054 : : * pszResampling : one of "NEAREST", "GAUSS", "CUBIC", "CUBICSPLINE" (GDAL >= 2.0), "LANCZOS" ( GDAL >= 2.0),
2055 : : * "BILINEAR" ( GDAL >= 2.0), "AVERAGE", "MODE" or "NONE" controlling the downsampling method applied.
2056 : : * nOverviews : number of overviews to build.
2057 : : * panOverviewList : the list of overview decimation factors to build.
2058 : : * nListBands : number of bands to build overviews for in panBandList. Build for all bands if this is 0.
2059 : : * panBandList : list of band numbers.
2060 : : * pfnProgress : a function to call to report progress, or nullptr.
2061 : : * pProgressData : application data to pass to the progress function.
2062 : : */
2063 : :
2064 : : // resampling method is now passed directly, via QgsRasterDataProvider::pyramidResamplingArg()
2065 : : // average_mp and average_magphase have been removed from the gui
2066 : 0 : QByteArray ba = resamplingMethod.toLocal8Bit();
2067 : 0 : const char *method = ba.data();
2068 : :
2069 : : //build the pyramid and show progress to console
2070 : 0 : QgsDebugMsgLevel( QStringLiteral( "Building overviews at %1 levels using resampling method %2"
2071 : : ).arg( myOverviewLevelsVector.size() ).arg( method ), 2 );
2072 : : try
2073 : : {
2074 : : //build the pyramid and show progress to console
2075 : 0 : QgsGdalProgress myProg;
2076 : 0 : myProg.type = QgsRaster::ProgressPyramids;
2077 : 0 : myProg.provider = this;
2078 : 0 : myProg.feedback = feedback;
2079 : 0 : myError = GDALBuildOverviews( mGdalBaseDataset, method,
2080 : 0 : myOverviewLevelsVector.size(), myOverviewLevelsVector.data(),
2081 : : 0, nullptr,
2082 : 0 : progressCallback, &myProg ); //this is the arg for the gdal progress callback
2083 : :
2084 : 0 : if ( ( feedback && feedback->isCanceled() ) || myError == CE_Failure || CPLGetLastErrorNo() == CPLE_NotSupported )
2085 : : {
2086 : 0 : QgsDebugMsg( QStringLiteral( "Building pyramids failed using resampling method [%1]" ).arg( method ) );
2087 : : //something bad happenend
2088 : : //QString myString = QString (CPLGetLastError());
2089 : 0 : GDALClose( mGdalBaseDataset );
2090 : 0 : mGdalBaseDataset = gdalOpen( dataSourceUri( true ), mUpdate ? GDAL_OF_UPDATE : GDAL_OF_READONLY );
2091 : : //Since we are not a virtual warped dataset, mGdalDataSet and mGdalBaseDataset are supposed to be the same
2092 : 0 : mGdalDataset = mGdalBaseDataset;
2093 : :
2094 : : // restore former configOptions
2095 : 0 : for ( QgsStringMap::const_iterator it = myConfigOptionsOld.constBegin();
2096 : 0 : it != myConfigOptionsOld.constEnd(); ++it )
2097 : : {
2098 : 0 : QByteArray key = it.key().toLocal8Bit();
2099 : 0 : QByteArray value = it.value().toLocal8Bit();
2100 : 0 : CPLSetConfigOption( key.data(), value.data() );
2101 : 0 : }
2102 : :
2103 : : // TODO print exact error message
2104 : 0 : if ( feedback && feedback->isCanceled() )
2105 : 0 : return QStringLiteral( "CANCELED" );
2106 : :
2107 : 0 : return QStringLiteral( "FAILED_NOT_SUPPORTED" );
2108 : : }
2109 : : else
2110 : : {
2111 : 0 : QgsDebugMsgLevel( QStringLiteral( "Building pyramids finished OK" ), 2 );
2112 : : //make sure the raster knows it has pyramids
2113 : 0 : mHasPyramids = true;
2114 : : }
2115 : 0 : }
2116 : : catch ( CPLErr )
2117 : : {
2118 : 0 : QgsLogger::warning( QStringLiteral( "Pyramid overview building failed!" ) );
2119 : 0 : }
2120 : :
2121 : : // restore former configOptions
2122 : 0 : for ( QgsStringMap::const_iterator it = myConfigOptionsOld.constBegin();
2123 : 0 : it != myConfigOptionsOld.constEnd(); ++it )
2124 : : {
2125 : 0 : QByteArray key = it.key().toLocal8Bit();
2126 : 0 : QByteArray value = it.value().toLocal8Bit();
2127 : 0 : CPLSetConfigOption( key.data(), value.data() );
2128 : 0 : }
2129 : :
2130 : 0 : QgsDebugMsgLevel( QStringLiteral( "Pyramid overviews built" ), 2 );
2131 : :
2132 : : // Observed problem: if a *.rrd file exists and GDALBuildOverviews() is called,
2133 : : // the *.rrd is deleted and no overviews are created, if GDALBuildOverviews()
2134 : : // is called next time, it crashes somewhere in GDAL:
2135 : : // https://trac.osgeo.org/gdal/ticket/4831
2136 : : // Crash can be avoided if dataset is reopened, fixed in GDAL 1.9.2
2137 : 0 : if ( format == QgsRaster::PyramidsInternal )
2138 : : {
2139 : 0 : QgsDebugMsgLevel( QStringLiteral( "Reopening dataset ..." ), 2 );
2140 : : //close the gdal dataset and reopen it in read only mode
2141 : 0 : GDALClose( mGdalBaseDataset );
2142 : 0 : mGdalBaseDataset = gdalOpen( dataSourceUri( true ), mUpdate ? GDAL_OF_UPDATE : GDAL_OF_READONLY );
2143 : : //Since we are not a virtual warped dataset, mGdalDataSet and mGdalBaseDataset are supposed to be the same
2144 : 0 : mGdalDataset = mGdalBaseDataset;
2145 : 0 : }
2146 : :
2147 : : //emit drawingProgress( 0, 0 );
2148 : 0 : return QString(); // returning null on success
2149 : 0 : }
2150 : :
2151 : : #if 0
2152 : : QList<QgsRasterPyramid> QgsGdalProvider::buildPyramidList()
2153 : : {
2154 : : //
2155 : : // First we build up a list of potential pyramid layers
2156 : : //
2157 : : int myWidth = mWidth;
2158 : : int myHeight = mHeight;
2159 : : int myDivisor = 2;
2160 : :
2161 : : GDALRasterBandH myGDALBand = GDALGetRasterBand( mGdalDataset, 1 ); //just use the first band
2162 : :
2163 : : mPyramidList.clear();
2164 : : QgsDebugMsg( QStringLiteral( "Building initial pyramid list" ) );
2165 : : while ( ( myWidth / myDivisor > 32 ) && ( ( myHeight / myDivisor ) > 32 ) )
2166 : : {
2167 : :
2168 : : QgsRasterPyramid myRasterPyramid;
2169 : : myRasterPyramid.level = myDivisor;
2170 : : myRasterPyramid.xDim = ( int )( 0.5 + ( myWidth / static_cast<double>( myDivisor ) ) );
2171 : : myRasterPyramid.yDim = ( int )( 0.5 + ( myHeight / static_cast<double>( myDivisor ) ) );
2172 : : myRasterPyramid.exists = false;
2173 : :
2174 : : QgsDebugMsg( QStringLiteral( "Pyramid %1 xDim %2 yDim %3" ).arg( myRasterPyramid.level ).arg( myRasterPyramid.xDim ).arg( myRasterPyramid.yDim ) );
2175 : :
2176 : : //
2177 : : // Now we check if it actually exists in the raster layer
2178 : : // and also adjust the dimensions if the dimensions calculated
2179 : : // above are only a near match.
2180 : : //
2181 : : const int myNearMatchLimit = 5;
2182 : : if ( GDALGetOverviewCount( myGDALBand ) > 0 )
2183 : : {
2184 : : int myOverviewCount;
2185 : : for ( myOverviewCount = 0;
2186 : : myOverviewCount < GDALGetOverviewCount( myGDALBand );
2187 : : ++myOverviewCount )
2188 : : {
2189 : : GDALRasterBandH myOverview;
2190 : : myOverview = GDALGetOverview( myGDALBand, myOverviewCount );
2191 : : int myOverviewXDim = GDALGetRasterBandXSize( myOverview );
2192 : : int myOverviewYDim = GDALGetRasterBandYSize( myOverview );
2193 : : //
2194 : : // here is where we check if its a near match:
2195 : : // we will see if its within 5 cells either side of
2196 : : //
2197 : : QgsDebugMsg( "Checking whether " + QString::number( myRasterPyramid.xDim ) + " x " +
2198 : : QString::number( myRasterPyramid.yDim ) + " matches " +
2199 : : QString::number( myOverviewXDim ) + " x " + QString::number( myOverviewYDim ) );
2200 : :
2201 : :
2202 : : if ( ( myOverviewXDim <= ( myRasterPyramid.xDim + myNearMatchLimit ) ) &&
2203 : : ( myOverviewXDim >= ( myRasterPyramid.xDim - myNearMatchLimit ) ) &&
2204 : : ( myOverviewYDim <= ( myRasterPyramid.yDim + myNearMatchLimit ) ) &&
2205 : : ( myOverviewYDim >= ( myRasterPyramid.yDim - myNearMatchLimit ) ) )
2206 : : {
2207 : : //right we have a match so adjust the a / y before they get added to the list
2208 : : myRasterPyramid.xDim = myOverviewXDim;
2209 : : myRasterPyramid.yDim = myOverviewYDim;
2210 : : myRasterPyramid.exists = true;
2211 : : QgsDebugMsg( QStringLiteral( ".....YES!" ) );
2212 : : }
2213 : : else
2214 : : {
2215 : : //no match
2216 : : QgsDebugMsg( QStringLiteral( ".....no." ) );
2217 : : }
2218 : : }
2219 : : }
2220 : : mPyramidList.append( myRasterPyramid );
2221 : : //sqare the divisor each step
2222 : : myDivisor = ( myDivisor * 2 );
2223 : : }
2224 : :
2225 : : return mPyramidList;
2226 : : }
2227 : : #endif
2228 : :
2229 : 0 : QList<QgsRasterPyramid> QgsGdalProvider::buildPyramidList( const QList<int> &list )
2230 : : {
2231 : 0 : QList< int > overviewList = list;
2232 : 0 : QMutexLocker locker( mpMutex );
2233 : :
2234 : 0 : int myWidth = mWidth;
2235 : 0 : int myHeight = mHeight;
2236 : 0 : GDALRasterBandH myGDALBand = GDALGetRasterBand( mGdalDataset, 1 ); //just use the first band
2237 : :
2238 : 0 : mPyramidList.clear();
2239 : :
2240 : : // if overviewList is empty (default) build the pyramid list
2241 : 0 : if ( overviewList.isEmpty() )
2242 : : {
2243 : 0 : int myDivisor = 2;
2244 : :
2245 : 0 : QgsDebugMsgLevel( QStringLiteral( "Building initial pyramid list" ), 2 );
2246 : :
2247 : 0 : while ( ( myWidth / myDivisor > 32 ) && ( ( myHeight / myDivisor ) > 32 ) )
2248 : : {
2249 : 0 : overviewList.append( myDivisor );
2250 : : //sqare the divisor each step
2251 : 0 : myDivisor = ( myDivisor * 2 );
2252 : : }
2253 : 0 : }
2254 : :
2255 : : // loop over pyramid list
2256 : 0 : const auto constOverviewList = overviewList;
2257 : 0 : for ( int myDivisor : constOverviewList )
2258 : : {
2259 : : //
2260 : : // First we build up a list of potential pyramid layers
2261 : : //
2262 : :
2263 : 0 : QgsRasterPyramid myRasterPyramid;
2264 : 0 : myRasterPyramid.setLevel( myDivisor );
2265 : 0 : myRasterPyramid.setXDim( ( int )( 0.5 + ( myWidth / static_cast<double>( myDivisor ) ) ) ); // NOLINT
2266 : 0 : myRasterPyramid.setYDim( ( int )( 0.5 + ( myHeight / static_cast<double>( myDivisor ) ) ) ); // NOLINT
2267 : 0 : myRasterPyramid.setExists( false );
2268 : :
2269 : 0 : QgsDebugMsgLevel( QStringLiteral( "Pyramid %1 xDim %2 yDim %3" ).arg( myRasterPyramid.getLevel() ).arg( myRasterPyramid.getXDim() ).arg( myRasterPyramid.getYDim() ), 2 );
2270 : :
2271 : : //
2272 : : // Now we check if it actually exists in the raster layer
2273 : : // and also adjust the dimensions if the dimensions calculated
2274 : : // above are only a near match.
2275 : : //
2276 : 0 : const int myNearMatchLimit = 5;
2277 : 0 : if ( GDALGetOverviewCount( myGDALBand ) > 0 )
2278 : : {
2279 : : int myOverviewCount;
2280 : 0 : for ( myOverviewCount = 0;
2281 : 0 : myOverviewCount < GDALGetOverviewCount( myGDALBand );
2282 : 0 : ++myOverviewCount )
2283 : : {
2284 : : GDALRasterBandH myOverview;
2285 : 0 : myOverview = GDALGetOverview( myGDALBand, myOverviewCount );
2286 : 0 : int myOverviewXDim = GDALGetRasterBandXSize( myOverview );
2287 : 0 : int myOverviewYDim = GDALGetRasterBandYSize( myOverview );
2288 : : //
2289 : : // here is where we check if its a near match:
2290 : : // we will see if its within 5 cells either side of
2291 : : //
2292 : 0 : QgsDebugMsgLevel( "Checking whether " + QString::number( myRasterPyramid.getXDim() ) + " x " +
2293 : : QString::number( myRasterPyramid.getYDim() ) + " matches " +
2294 : : QString::number( myOverviewXDim ) + " x " + QString::number( myOverviewYDim ), 2 );
2295 : :
2296 : :
2297 : 0 : if ( ( myOverviewXDim <= ( myRasterPyramid.getXDim() + myNearMatchLimit ) ) &&
2298 : 0 : ( myOverviewXDim >= ( myRasterPyramid.getXDim() - myNearMatchLimit ) ) &&
2299 : 0 : ( myOverviewYDim <= ( myRasterPyramid.getYDim() + myNearMatchLimit ) ) &&
2300 : 0 : ( myOverviewYDim >= ( myRasterPyramid.getYDim() - myNearMatchLimit ) ) )
2301 : : {
2302 : : //right we have a match so adjust the a / y before they get added to the list
2303 : 0 : myRasterPyramid.setXDim( myOverviewXDim );
2304 : 0 : myRasterPyramid.setYDim( myOverviewYDim );
2305 : 0 : myRasterPyramid.setExists( true );
2306 : 0 : QgsDebugMsgLevel( QStringLiteral( ".....YES!" ), 2 );
2307 : 0 : }
2308 : : else
2309 : : {
2310 : : //no match
2311 : 0 : QgsDebugMsgLevel( QStringLiteral( ".....no." ), 2 );
2312 : : }
2313 : 0 : }
2314 : 0 : }
2315 : 0 : mPyramidList.append( myRasterPyramid );
2316 : : }
2317 : :
2318 : 0 : return mPyramidList;
2319 : 0 : }
2320 : :
2321 : 0 : QStringList QgsGdalProvider::subLayers() const
2322 : : {
2323 : 0 : return mSubLayers;
2324 : : }
2325 : :
2326 : 0 : QVariantMap QgsGdalProviderMetadata::decodeUri( const QString &uri ) const
2327 : : {
2328 : 0 : return QgsGdalProviderBase::decodeGdalUri( uri );
2329 : : }
2330 : :
2331 : 0 : QString QgsGdalProviderMetadata::encodeUri( const QVariantMap &parts ) const
2332 : : {
2333 : 0 : return QgsGdalProviderBase::encodeGdalUri( parts );
2334 : : }
2335 : :
2336 : 0 : bool QgsGdalProviderMetadata::uriIsBlocklisted( const QString &uri ) const
2337 : : {
2338 : 0 : const QVariantMap parts = decodeUri( uri );
2339 : 0 : if ( !parts.contains( QStringLiteral( "path" ) ) )
2340 : 0 : return false;
2341 : :
2342 : 0 : QFileInfo fi( parts.value( QStringLiteral( "path" ) ).toString() );
2343 : 0 : const QString suffix = fi.completeSuffix();
2344 : :
2345 : : // internal details only
2346 : 0 : if ( suffix.compare( QLatin1String( "aux.xml" ), Qt::CaseInsensitive ) == 0 || suffix.endsWith( QLatin1String( ".aux.xml" ), Qt::CaseInsensitive ) )
2347 : 0 : return true;
2348 : 0 : if ( suffix.compare( QLatin1String( "tif.xml" ), Qt::CaseInsensitive ) == 0 )
2349 : 0 : return true;
2350 : :
2351 : 0 : return false;
2352 : 0 : }
2353 : :
2354 : :
2355 : 0 : QgsGdalProvider *QgsGdalProviderMetadata::createProvider( const QString &uri, const QgsDataProvider::ProviderOptions &options, QgsDataProvider::ReadFlags flags )
2356 : : {
2357 : : Q_UNUSED( flags );
2358 : 0 : return new QgsGdalProvider( uri, options );
2359 : 0 : }
2360 : :
2361 : : /**
2362 : :
2363 : : Convenience function for readily creating file filters.
2364 : :
2365 : : Given a long name for a file filter and a regular expression, return
2366 : : a file filter string suitable for use in a QFileDialog::OpenFiles()
2367 : : call. The regular express, glob, will have both all lower and upper
2368 : : case versions added.
2369 : :
2370 : : \note
2371 : :
2372 : : Copied from qgisapp.cpp.
2373 : :
2374 : : \todo XXX This should probably be generalized and moved to a standard
2375 : : utility type thingy.
2376 : :
2377 : : */
2378 : 264 : static QString createFileFilter_( QString const &longName, QString const &glob )
2379 : : {
2380 : : // return longName + " [GDAL] (" + glob.toLower() + ' ' + glob.toUpper() + ");;";
2381 : 264 : return longName + " (" + glob.toLower() + ' ' + glob.toUpper() + ");;";
2382 : 0 : } // createFileFilter_
2383 : :
2384 : 3 : void buildSupportedRasterFileFilterAndExtensions( QString &fileFiltersString, QStringList &extensions, QStringList &wildcards )
2385 : : {
2386 : :
2387 : : // then iterate through all of the supported drivers, adding the
2388 : : // corresponding file filter
2389 : :
2390 : : GDALDriverH myGdalDriver; // current driver
2391 : :
2392 : 3 : QStringList catchallFilter; // for Any file(*.*), but also for those
2393 : : // drivers with no specific file filter
2394 : :
2395 : 3 : GDALDriverH jp2Driver = nullptr; // first JPEG2000 driver found
2396 : :
2397 : 3 : QgsGdalProviderBase::registerGdalDrivers();
2398 : :
2399 : : // Grind through all the drivers and their respective metadata.
2400 : : // We'll add a file filter for those drivers that have a file
2401 : : // extension defined for them; the others, well, even though
2402 : : // theoreticaly we can open those files because there exists a
2403 : : // driver for them, the user will have to use the "All Files" to
2404 : : // open datasets with no explicitly defined file name extension.
2405 : :
2406 : 3 : fileFiltersString.clear();
2407 : :
2408 : 3 : QgsDebugMsgLevel( QStringLiteral( "GDAL driver count: %1" ).arg( GDALGetDriverCount() ), 2 );
2409 : :
2410 : 672 : for ( int i = 0; i < GDALGetDriverCount(); ++i )
2411 : : {
2412 : 669 : myGdalDriver = GDALGetDriver( i );
2413 : :
2414 : 669 : Q_CHECK_PTR( myGdalDriver ); // NOLINT
2415 : :
2416 : 669 : if ( !myGdalDriver )
2417 : : {
2418 : 0 : QgsLogger::warning( "unable to get driver " + QString::number( i ) );
2419 : 0 : continue;
2420 : : }
2421 : :
2422 : : // in GDAL 2.0 vector and mixed drivers are returned by GDALGetDriver, so filter out non-raster drivers
2423 : 669 : if ( QString( GDALGetMetadataItem( myGdalDriver, GDAL_DCAP_RASTER, nullptr ) ) != QLatin1String( "YES" ) )
2424 : 219 : continue;
2425 : :
2426 : : // now we need to see if the driver is for something currently
2427 : : // supported; if not, we give it a miss for the next driver
2428 : :
2429 : 450 : QString myGdalDriverDescription = GDALGetDescription( myGdalDriver );
2430 : 450 : if ( myGdalDriverDescription == QLatin1String( "BIGGIF" ) )
2431 : : {
2432 : : // BIGGIF is a technical driver. The plain GIF driver will do
2433 : 3 : continue;
2434 : : }
2435 : :
2436 : : // QgsDebugMsg(QString("got driver string %1").arg(myGdalDriverDescription));
2437 : :
2438 : 447 : QString myGdalDriverExtensions = GDALGetMetadataItem( myGdalDriver, GDAL_DMD_EXTENSIONS, "" );
2439 : 447 : QString myGdalDriverLongName = GDALGetMetadataItem( myGdalDriver, GDAL_DMD_LONGNAME, "" );
2440 : : // remove any superfluous (.*) strings at the end as
2441 : : // they'll confuse QFileDialog::getOpenFileNames()
2442 : 447 : myGdalDriverLongName.remove( QRegExp( "\\(.*\\)$" ) );
2443 : :
2444 : : // if we have both the file name extension and the long name,
2445 : : // then we've all the information we need for the current
2446 : : // driver; therefore emit a file filter string and move to
2447 : : // the next driver
2448 : 447 : if ( !( myGdalDriverExtensions.isEmpty() || myGdalDriverLongName.isEmpty() ) )
2449 : : {
2450 : : #if QT_VERSION < QT_VERSION_CHECK(5, 14, 0)
2451 : : const QStringList splitExtensions = myGdalDriverExtensions.split( ' ', QString::SkipEmptyParts );
2452 : : #else
2453 : 264 : const QStringList splitExtensions = myGdalDriverExtensions.split( ' ', Qt::SkipEmptyParts );
2454 : : #endif
2455 : :
2456 : : // XXX add check for SDTS; in that case we want (*CATD.DDF)
2457 : 264 : QString glob;
2458 : :
2459 : 585 : for ( const QString &ext : splitExtensions )
2460 : : {
2461 : : // This hacking around that removes '/' is no longer necessary with GDAL 2.3
2462 : 321 : extensions << QString( ext ).remove( '/' ).remove( '*' ).remove( '.' );
2463 : 321 : if ( !glob.isEmpty() )
2464 : 57 : glob += QLatin1Char( ' ' );
2465 : 321 : glob += "*." + QString( ext ).replace( '/', QLatin1String( " *." ) );
2466 : : }
2467 : :
2468 : : // Add only the first JP2 driver found to the filter list (it's the one GDAL uses)
2469 : 525 : if ( myGdalDriverDescription == QLatin1String( "JPEG2000" ) ||
2470 : 261 : myGdalDriverDescription.startsWith( QLatin1String( "JP2" ) ) ) // JP2ECW, JP2KAK, JP2MrSID
2471 : : {
2472 : 9 : if ( jp2Driver )
2473 : 6 : continue; // skip if already found a JP2 driver
2474 : :
2475 : 3 : jp2Driver = myGdalDriver; // first JP2 driver found
2476 : 3 : if ( !glob.contains( "j2k" ) )
2477 : : {
2478 : 3 : glob += QLatin1String( " *.j2k" ); // add alternate extension
2479 : 6 : extensions << QStringLiteral( "j2k" );
2480 : 3 : }
2481 : 3 : }
2482 : 255 : else if ( myGdalDriverDescription == QLatin1String( "VRT" ) )
2483 : : {
2484 : 3 : glob += QLatin1String( " *.ovr" );
2485 : 6 : extensions << QStringLiteral( "ovr" );
2486 : 3 : }
2487 : :
2488 : 258 : fileFiltersString += createFileFilter_( myGdalDriverLongName, glob );
2489 : 264 : }
2490 : :
2491 : 441 : if ( myGdalDriverExtensions.isEmpty() && !myGdalDriverLongName.isEmpty() )
2492 : : {
2493 : : // Then what we have here is a driver with no corresponding
2494 : : // file extension; e.g., GRASS. In which case we append the
2495 : : // string to the "catch-all" which will match all file types.
2496 : : // (I.e., "*.*") We use the driver description instead of the
2497 : : // long time to prevent the catch-all line from getting too
2498 : : // large.
2499 : :
2500 : : // ... OTOH, there are some drivers with missing
2501 : : // DMD_EXTENSION; so let's check for them here and handle
2502 : : // them appropriately
2503 : :
2504 : 183 : if ( myGdalDriverDescription.startsWith( QLatin1String( "AIG" ) ) )
2505 : : {
2506 : 6 : fileFiltersString += createFileFilter_( myGdalDriverLongName, QStringLiteral( "hdr.adf" ) );
2507 : 6 : wildcards << QStringLiteral( "hdr.adf" );
2508 : 3 : }
2509 : : #if !(GDAL_VERSION_MAJOR > 2 || (GDAL_VERSION_MAJOR == 2 && GDAL_VERSION_MINOR >= 3))
2510 : : else if ( myGdalDriverDescription.startsWith( QLatin1String( "EHdr" ) ) )
2511 : : {
2512 : : // Fixed in GDAL 2.3
2513 : : fileFiltersString += createFileFilter_( myGdalDriverLongName, QStringLiteral( "*.bil" ) );
2514 : : extensions << QStringLiteral( "bil" );
2515 : : }
2516 : : else if ( myGdalDriverDescription == QLatin1String( "ERS" ) )
2517 : : {
2518 : : // Fixed in GDAL 2.3
2519 : : fileFiltersString += createFileFilter_( myGdalDriverLongName, QStringLiteral( "*.ers" ) );
2520 : : extensions << QStringLiteral( "ers" );
2521 : : }
2522 : : #endif
2523 : : else
2524 : : {
2525 : 180 : catchallFilter << QString( GDALGetDescription( myGdalDriver ) );
2526 : : }
2527 : 183 : } // each loaded GDAL driver
2528 : :
2529 : 450 : } // each loaded GDAL driver
2530 : :
2531 : : // sort file filters alphabetically
2532 : : #if QT_VERSION < QT_VERSION_CHECK(5, 15, 0)
2533 : : QStringList filters = fileFiltersString.split( QStringLiteral( ";;" ), QString::SkipEmptyParts );
2534 : : #else
2535 : 6 : QStringList filters = fileFiltersString.split( QStringLiteral( ";;" ), Qt::SkipEmptyParts );
2536 : : #endif
2537 : 3 : filters.sort();
2538 : 3 : fileFiltersString = filters.join( QLatin1String( ";;" ) ) + ";;";
2539 : :
2540 : : // VSIFileHandler (see qgsogrprovider.cpp) - second
2541 : 3 : QgsSettings settings;
2542 : 6 : if ( settings.value( QStringLiteral( "qgis/scanZipInBrowser2" ), "basic" ).toString() != QLatin1String( "no" ) )
2543 : : {
2544 : 6 : fileFiltersString.prepend( createFileFilter_( QObject::tr( "GDAL/OGR VSIFileHandler" ), QStringLiteral( "*.zip *.gz *.tar *.tar.gz *.tgz" ) ) );
2545 : 18 : extensions << QStringLiteral( "zip" ) << QStringLiteral( "gz" ) << QStringLiteral( "tar" ) << QStringLiteral( "tar.gz" ) << QStringLiteral( "tgz" );
2546 : 3 : }
2547 : :
2548 : : // can't forget the all supported case
2549 : 3 : QStringList exts;
2550 : 345 : for ( const QString &ext : std::as_const( extensions ) )
2551 : 684 : exts << QStringLiteral( "*.%1 *.%2" ).arg( ext, ext.toUpper() );
2552 : 6 : fileFiltersString.prepend( QObject::tr( "All supported files" ) + QStringLiteral( " (%1);;" ).arg( exts.join( QLatin1Char( ' ' ) ) ) );
2553 : :
2554 : : // can't forget the default case - first
2555 : 3 : fileFiltersString.prepend( QObject::tr( "All files" ) + " (*);;" );
2556 : :
2557 : : // cleanup
2558 : 3 : if ( fileFiltersString.endsWith( QLatin1String( ";;" ) ) ) fileFiltersString.chop( 2 );
2559 : :
2560 : 3 : QgsDebugMsgLevel( "Raster filter list built: " + fileFiltersString, 2 );
2561 : 3 : QgsDebugMsgLevel( "Raster extension list built: " + extensions.join( ' ' ), 2 );
2562 : 3 : } // buildSupportedRasterFileFilter_()
2563 : :
2564 : :
2565 : 0 : bool QgsGdalProvider::isValidRasterFileName( QString const &fileNameQString, QString &retErrMsg )
2566 : : {
2567 : 0 : gdal::dataset_unique_ptr myDataset;
2568 : :
2569 : 0 : QgsGdalProviderBase::registerGdalDrivers();
2570 : :
2571 : 0 : CPLErrorReset();
2572 : :
2573 : 0 : QString fileName = fileNameQString;
2574 : :
2575 : : // Try to open using VSIFileHandler (see qgsogrprovider.cpp)
2576 : : // TODO suppress error messages and report in debug, like in OGR provider
2577 : 0 : QString vsiPrefix = QgsZipItem::vsiPrefix( fileName );
2578 : 0 : if ( !vsiPrefix.isEmpty() )
2579 : : {
2580 : 0 : if ( !fileName.startsWith( vsiPrefix ) )
2581 : 0 : fileName = vsiPrefix + fileName;
2582 : 0 : QgsDebugMsgLevel( QStringLiteral( "Trying %1 syntax, fileName= %2" ).arg( vsiPrefix, fileName ), 2 );
2583 : 0 : }
2584 : :
2585 : : //open the file using gdal making sure we have handled locale properly
2586 : : //myDataset = GDALOpen( QFile::encodeName( fileNameQString ).constData(), GA_ReadOnly );
2587 : 0 : myDataset.reset( QgsGdalProviderBase::gdalOpen( fileName, GDAL_OF_READONLY ) );
2588 : 0 : if ( !myDataset )
2589 : : {
2590 : 0 : if ( CPLGetLastErrorNo() != CPLE_OpenFailed )
2591 : 0 : retErrMsg = QString::fromUtf8( CPLGetLastErrorMsg() );
2592 : 0 : return false;
2593 : : }
2594 : 0 : else if ( GDALGetRasterCount( myDataset.get() ) == 0 )
2595 : : {
2596 : 0 : QStringList layers = QgsGdalProvider::subLayers( myDataset.get() );
2597 : 0 : if ( layers.isEmpty() )
2598 : : {
2599 : 0 : retErrMsg = QObject::tr( "This raster file has no bands and is invalid as a raster layer." );
2600 : 0 : return false;
2601 : : }
2602 : 0 : return true;
2603 : 0 : }
2604 : : else
2605 : : {
2606 : 0 : return true;
2607 : : }
2608 : 0 : }
2609 : :
2610 : 0 : bool QgsGdalProvider::hasStatistics( int bandNo,
2611 : : int stats,
2612 : : const QgsRectangle &boundingBox,
2613 : : int sampleSize )
2614 : : {
2615 : 0 : QMutexLocker locker( mpMutex );
2616 : 0 : if ( !initIfNeeded() )
2617 : 0 : return false;
2618 : :
2619 : 0 : QgsDebugMsgLevel( QStringLiteral( "theBandNo = %1 sampleSize = %2" ).arg( bandNo ).arg( sampleSize ), 2 );
2620 : :
2621 : : // First check if cached in mStatistics
2622 : 0 : if ( QgsRasterDataProvider::hasStatistics( bandNo, stats, boundingBox, sampleSize ) )
2623 : : {
2624 : 0 : return true;
2625 : : }
2626 : :
2627 : 0 : QgsRasterBandStats myRasterBandStats;
2628 : 0 : initStatistics( myRasterBandStats, bandNo, stats, boundingBox, sampleSize );
2629 : :
2630 : 0 : if ( ( sourceHasNoDataValue( bandNo ) && !useSourceNoDataValue( bandNo ) ) ||
2631 : 0 : !userNoDataValues( bandNo ).isEmpty() )
2632 : : {
2633 : 0 : QgsDebugMsgLevel( QStringLiteral( "Custom no data values -> GDAL statistics not sufficient." ), 2 );
2634 : 0 : return false;
2635 : : }
2636 : :
2637 : : // If not cached, check if supported by GDAL
2638 : 0 : int supportedStats = QgsRasterBandStats::Min | QgsRasterBandStats::Max
2639 : : | QgsRasterBandStats::Range | QgsRasterBandStats::Mean
2640 : : | QgsRasterBandStats::StdDev;
2641 : :
2642 : 0 : if ( myRasterBandStats.extent != extent() ||
2643 : 0 : ( stats & ( ~supportedStats ) ) )
2644 : : {
2645 : 0 : QgsDebugMsg( QStringLiteral( "Not supported by GDAL." ) );
2646 : 0 : return false;
2647 : : }
2648 : :
2649 : 0 : QgsDebugMsgLevel( QStringLiteral( "Looking for GDAL statistics" ), 2 );
2650 : :
2651 : 0 : GDALRasterBandH myGdalBand = getBand( bandNo );
2652 : 0 : if ( ! myGdalBand )
2653 : : {
2654 : 0 : return false;
2655 : : }
2656 : :
2657 : 0 : int bApproxOK = false;
2658 : 0 : if ( sampleSize > 0 )
2659 : : {
2660 : 0 : if ( ( static_cast<double>( xSize() ) * static_cast<double>( ySize() ) / sampleSize ) > 2 ) // not perfect
2661 : : {
2662 : 0 : bApproxOK = true;
2663 : 0 : }
2664 : 0 : }
2665 : :
2666 : : // Params in GDALGetRasterStatistics must not be nullptr otherwise GDAL returns
2667 : : // without error even if stats are not cached
2668 : : double dfMin, dfMax, dfMean, dfStdDev;
2669 : 0 : double *pdfMin = &dfMin;
2670 : 0 : double *pdfMax = &dfMax;
2671 : 0 : double *pdfMean = &dfMean;
2672 : 0 : double *pdfStdDev = &dfStdDev;
2673 : :
2674 : 0 : if ( !( stats & QgsRasterBandStats::Min ) ) pdfMin = nullptr;
2675 : 0 : if ( !( stats & QgsRasterBandStats::Max ) ) pdfMax = nullptr;
2676 : 0 : if ( !( stats & QgsRasterBandStats::Mean ) ) pdfMean = nullptr;
2677 : 0 : if ( !( stats & QgsRasterBandStats::StdDev ) ) pdfStdDev = nullptr;
2678 : :
2679 : : // try to fetch the cached stats (bForce=FALSE)
2680 : : // Unfortunately GDALGetRasterStatistics() does not work as expected according to
2681 : : // API doc, if bApproxOK=false and bForce=false/true and exact statistics
2682 : : // (from all raster pixels) are not available/cached, it should return CE_Warning.
2683 : : // Instead, it is giving estimated (from sample) cached statistics and it returns CE_None.
2684 : : // see above and https://trac.osgeo.org/gdal/ticket/4857
2685 : : // -> Cannot used cached GDAL stats for exact
2686 : 0 : if ( !bApproxOK ) return false;
2687 : :
2688 : 0 : CPLErr myerval = GDALGetRasterStatistics( myGdalBand, bApproxOK, true, pdfMin, pdfMax, pdfMean, pdfStdDev );
2689 : :
2690 : 0 : if ( CE_None == myerval ) // CE_Warning if cached not found
2691 : : {
2692 : 0 : QgsDebugMsgLevel( QStringLiteral( "GDAL has cached statistics" ), 2 );
2693 : 0 : return true;
2694 : : }
2695 : :
2696 : 0 : return false;
2697 : 0 : }
2698 : :
2699 : 0 : QgsRasterBandStats QgsGdalProvider::bandStatistics( int bandNo, int stats, const QgsRectangle &boundingBox, int sampleSize, QgsRasterBlockFeedback *feedback )
2700 : : {
2701 : 0 : QMutexLocker locker( mpMutex );
2702 : 0 : if ( !initIfNeeded() )
2703 : 0 : return QgsRasterBandStats();
2704 : :
2705 : 0 : QgsDebugMsgLevel( QStringLiteral( "theBandNo = %1 sampleSize = %2" ).arg( bandNo ).arg( sampleSize ), 2 );
2706 : :
2707 : : // TODO: null values set on raster layer!!!
2708 : :
2709 : : // Currently there is no API in GDAL to collect statistics of specified extent
2710 : : // or with defined sample size. We check first if we have cached stats, if not,
2711 : : // and it is not possible to use GDAL we call generic provider method,
2712 : : // otherwise we use GDAL (faster, cache)
2713 : :
2714 : 0 : QgsRasterBandStats myRasterBandStats;
2715 : 0 : initStatistics( myRasterBandStats, bandNo, stats, boundingBox, sampleSize );
2716 : :
2717 : 0 : const auto constMStatistics = mStatistics;
2718 : 0 : for ( const QgsRasterBandStats &stats : constMStatistics )
2719 : : {
2720 : 0 : if ( stats.contains( myRasterBandStats ) )
2721 : : {
2722 : 0 : QgsDebugMsgLevel( QStringLiteral( "Using cached statistics." ), 2 );
2723 : 0 : return stats;
2724 : : }
2725 : : }
2726 : :
2727 : : // We cannot use GDAL stats if user disabled src no data value or set
2728 : : // custom no data values
2729 : 0 : if ( ( sourceHasNoDataValue( bandNo ) && !useSourceNoDataValue( bandNo ) ) ||
2730 : 0 : !userNoDataValues( bandNo ).isEmpty() )
2731 : : {
2732 : 0 : QgsDebugMsgLevel( QStringLiteral( "Custom no data values, using generic statistics." ), 2 );
2733 : 0 : return QgsRasterDataProvider::bandStatistics( bandNo, stats, boundingBox, sampleSize, feedback );
2734 : : }
2735 : :
2736 : 0 : int supportedStats = QgsRasterBandStats::Min | QgsRasterBandStats::Max
2737 : : | QgsRasterBandStats::Range | QgsRasterBandStats::Mean
2738 : : | QgsRasterBandStats::StdDev;
2739 : :
2740 : 0 : QgsDebugMsgLevel( QStringLiteral( "theStats = %1 supportedStats = %2" ).arg( stats, 0, 2 ).arg( supportedStats, 0, 2 ), 2 );
2741 : :
2742 : 0 : if ( myRasterBandStats.extent != extent() ||
2743 : 0 : ( stats & ( ~supportedStats ) ) )
2744 : : {
2745 : 0 : QgsDebugMsgLevel( QStringLiteral( "Statistics not supported by provider, using generic statistics." ), 2 );
2746 : 0 : return QgsRasterDataProvider::bandStatistics( bandNo, stats, boundingBox, sampleSize, feedback );
2747 : : }
2748 : :
2749 : 0 : QgsDebugMsgLevel( QStringLiteral( "Using GDAL statistics." ), 2 );
2750 : 0 : GDALRasterBandH myGdalBand = getBand( bandNo );
2751 : :
2752 : : //int bApproxOK = false; //as we asked for stats, don't get approx values
2753 : : // GDAL does not have sample size parameter in API, just bApproxOK or not,
2754 : : // we decide if approximation should be used according to
2755 : : // total size / sample size ratio
2756 : 0 : int bApproxOK = false;
2757 : 0 : if ( sampleSize > 0 )
2758 : : {
2759 : 0 : if ( ( static_cast<double>( xSize() ) * static_cast<double>( ySize() ) / sampleSize ) > 2 ) // not perfect
2760 : : {
2761 : 0 : bApproxOK = true;
2762 : 0 : }
2763 : 0 : }
2764 : :
2765 : 0 : QgsDebugMsgLevel( QStringLiteral( "bApproxOK = %1" ).arg( bApproxOK ), 2 );
2766 : :
2767 : : double pdfMin;
2768 : : double pdfMax;
2769 : : double pdfMean;
2770 : : double pdfStdDev;
2771 : 0 : QgsGdalProgress myProg;
2772 : 0 : myProg.type = QgsRaster::ProgressHistogram;
2773 : 0 : myProg.provider = this;
2774 : 0 : myProg.feedback = feedback;
2775 : :
2776 : : // try to fetch the cached stats (bForce=FALSE)
2777 : : // GDALGetRasterStatistics() do not work correctly with bApproxOK=false and bForce=false/true
2778 : : // see above and https://trac.osgeo.org/gdal/ticket/4857
2779 : : // -> Cannot used cached GDAL stats for exact
2780 : :
2781 : 0 : CPLErr myerval =
2782 : 0 : GDALGetRasterStatistics( myGdalBand, bApproxOK, true, &pdfMin, &pdfMax, &pdfMean, &pdfStdDev );
2783 : :
2784 : 0 : QgsDebugMsgLevel( QStringLiteral( "myerval = %1" ).arg( myerval ), 2 );
2785 : :
2786 : : // if cached stats are not found, compute them
2787 : 0 : if ( !bApproxOK || CE_None != myerval )
2788 : : {
2789 : 0 : QgsDebugMsgLevel( QStringLiteral( "Calculating statistics by GDAL" ), 2 );
2790 : 0 : myerval = GDALComputeRasterStatistics( myGdalBand, bApproxOK,
2791 : : &pdfMin, &pdfMax, &pdfMean, &pdfStdDev,
2792 : 0 : progressCallback, &myProg );
2793 : 0 : mStatisticsAreReliable = true;
2794 : 0 : }
2795 : : else
2796 : : {
2797 : 0 : QgsDebugMsgLevel( QStringLiteral( "Using GDAL cached statistics" ), 2 );
2798 : : }
2799 : :
2800 : 0 : if ( feedback && feedback->isCanceled() )
2801 : 0 : return myRasterBandStats;
2802 : :
2803 : : // if stats are found populate the QgsRasterBandStats object
2804 : 0 : if ( CE_None == myerval )
2805 : : {
2806 : 0 : myRasterBandStats.bandNumber = bandNo;
2807 : 0 : myRasterBandStats.range = pdfMax - pdfMin;
2808 : 0 : myRasterBandStats.minimumValue = pdfMin;
2809 : 0 : myRasterBandStats.maximumValue = pdfMax;
2810 : : //calculate the mean
2811 : 0 : myRasterBandStats.mean = pdfMean;
2812 : 0 : myRasterBandStats.sum = 0; //not available via gdal
2813 : : //myRasterBandStats.elementCount = mWidth * mHeight;
2814 : : // Sum of non NULL
2815 : 0 : myRasterBandStats.elementCount = 0; //not available via gdal
2816 : 0 : myRasterBandStats.sumOfSquares = 0; //not available via gdal
2817 : 0 : myRasterBandStats.stdDev = pdfStdDev;
2818 : 0 : myRasterBandStats.statsGathered = QgsRasterBandStats::Min | QgsRasterBandStats::Max
2819 : : | QgsRasterBandStats::Range | QgsRasterBandStats::Mean
2820 : : | QgsRasterBandStats::StdDev;
2821 : :
2822 : : // define if the band has scale and offset to apply
2823 : 0 : double myScale = bandScale( bandNo );
2824 : 0 : double myOffset = bandOffset( bandNo );
2825 : 0 : if ( myScale != 1.0 || myOffset != 0.0 )
2826 : : {
2827 : 0 : if ( myScale < 0.0 )
2828 : : {
2829 : : // update Min and Max value
2830 : 0 : myRasterBandStats.minimumValue = pdfMax * myScale + myOffset;
2831 : 0 : myRasterBandStats.maximumValue = pdfMin * myScale + myOffset;
2832 : : // update the range
2833 : 0 : myRasterBandStats.range = ( pdfMin - pdfMax ) * myScale;
2834 : : // update standard deviation
2835 : 0 : myRasterBandStats.stdDev = -1.0 * pdfStdDev * myScale;
2836 : 0 : }
2837 : : else
2838 : : {
2839 : : // update Min and Max value
2840 : 0 : myRasterBandStats.minimumValue = pdfMin * myScale + myOffset;
2841 : 0 : myRasterBandStats.maximumValue = pdfMax * myScale + myOffset;
2842 : : // update the range
2843 : 0 : myRasterBandStats.range = ( pdfMax - pdfMin ) * myScale;
2844 : : // update standard deviation
2845 : 0 : myRasterBandStats.stdDev = pdfStdDev * myScale;
2846 : : }
2847 : : // update the mean
2848 : 0 : myRasterBandStats.mean = pdfMean * myScale + myOffset;
2849 : 0 : }
2850 : :
2851 : : #ifdef QGISDEBUG
2852 : : QgsDebugMsgLevel( QStringLiteral( "************ STATS **************" ), 2 );
2853 : : QgsDebugMsgLevel( QStringLiteral( "MIN %1" ).arg( myRasterBandStats.minimumValue ), 2 );
2854 : : QgsDebugMsgLevel( QStringLiteral( "MAX %1" ).arg( myRasterBandStats.maximumValue ), 2 );
2855 : : QgsDebugMsgLevel( QStringLiteral( "RANGE %1" ).arg( myRasterBandStats.range ), 2 );
2856 : : QgsDebugMsgLevel( QStringLiteral( "MEAN %1" ).arg( myRasterBandStats.mean ), 2 );
2857 : : QgsDebugMsgLevel( QStringLiteral( "STDDEV %1" ).arg( myRasterBandStats.stdDev ), 2 );
2858 : : #endif
2859 : 0 : }
2860 : : else
2861 : : {
2862 : 0 : myRasterBandStats.statsGathered = QgsRasterBandStats::Stats::None;
2863 : : }
2864 : :
2865 : 0 : mStatistics.append( myRasterBandStats );
2866 : 0 : return myRasterBandStats;
2867 : :
2868 : 0 : } // QgsGdalProvider::bandStatistics
2869 : :
2870 : 0 : bool QgsGdalProvider::initIfNeeded()
2871 : : {
2872 : 0 : if ( mHasInit )
2873 : 0 : return mValid;
2874 : :
2875 : 0 : mHasInit = true;
2876 : :
2877 : 0 : QString gdalUri = dataSourceUri( true );
2878 : :
2879 : : // Try to open using VSIFileHandler (see qgsogrprovider.cpp)
2880 : 0 : QString vsiPrefix = QgsZipItem::vsiPrefix( gdalUri );
2881 : 0 : if ( !vsiPrefix.isEmpty() )
2882 : : {
2883 : 0 : if ( !gdalUri.startsWith( vsiPrefix ) )
2884 : 0 : setDataSourceUri( vsiPrefix + gdalUri );
2885 : 0 : QgsDebugMsgLevel( QStringLiteral( "Trying %1 syntax, uri= %2" ).arg( vsiPrefix, dataSourceUri() ), 2 );
2886 : 0 : }
2887 : :
2888 : 0 : gdalUri = dataSourceUri( true );
2889 : :
2890 : 0 : CPLErrorReset();
2891 : 0 : mGdalBaseDataset = gdalOpen( gdalUri, mUpdate ? GDAL_OF_UPDATE : GDAL_OF_READONLY );
2892 : :
2893 : 0 : if ( !mGdalBaseDataset )
2894 : : {
2895 : 0 : QString msg = QStringLiteral( "Cannot open GDAL dataset %1:\n%2" ).arg( dataSourceUri(), QString::fromUtf8( CPLGetLastErrorMsg() ) );
2896 : 0 : appendError( ERRMSG( msg ) );
2897 : 0 : return false;
2898 : 0 : }
2899 : :
2900 : 0 : QgsDebugMsgLevel( QStringLiteral( "GdalDataset opened" ), 2 );
2901 : :
2902 : 0 : initBaseDataset();
2903 : 0 : return mValid;
2904 : 0 : }
2905 : :
2906 : :
2907 : 0 : void QgsGdalProvider::initBaseDataset()
2908 : : {
2909 : 0 : mDriverName = GDALGetDriverShortName( GDALGetDatasetDriver( mGdalBaseDataset ) );
2910 : 0 : mHasInit = true;
2911 : 0 : mValid = true;
2912 : : #if 0
2913 : : for ( int i = 0; i < GDALGetRasterCount( mGdalBaseDataset ); i++ )
2914 : : {
2915 : : mMinMaxComputed.append( false );
2916 : : mMinimum.append( 0 );
2917 : : mMaximum.append( 0 );
2918 : : }
2919 : : #endif
2920 : : // Check if we need a warped VRT for this file.
2921 : 0 : bool hasGeoTransform = GDALGetGeoTransform( mGdalBaseDataset, mGeoTransform ) == CE_None;
2922 : 0 : if ( ( hasGeoTransform
2923 : 0 : && ( mGeoTransform[1] < 0.0
2924 : 0 : || mGeoTransform[2] != 0.0
2925 : 0 : || mGeoTransform[4] != 0.0
2926 : 0 : || mGeoTransform[5] > 0.0 ) )
2927 : 0 : || GDALGetGCPCount( mGdalBaseDataset ) > 0
2928 : 0 : || GDALGetMetadata( mGdalBaseDataset, "RPC" ) )
2929 : : {
2930 : 0 : QgsLogger::warning( QStringLiteral( "Creating Warped VRT." ) );
2931 : :
2932 : 0 : if ( GDALGetMetadata( mGdalBaseDataset, "RPC" ) )
2933 : : {
2934 : 0 : mGdalDataset =
2935 : 0 : QgsGdalUtils::rpcAwareAutoCreateWarpedVrt( mGdalBaseDataset, nullptr, nullptr,
2936 : : GRA_NearestNeighbour, 0.2, nullptr );
2937 : 0 : mGdalTransformerArg = QgsGdalUtils::rpcAwareCreateTransformer( mGdalBaseDataset );
2938 : 0 : }
2939 : : else
2940 : : {
2941 : 0 : mGdalDataset =
2942 : 0 : GDALAutoCreateWarpedVRT( mGdalBaseDataset, nullptr, nullptr,
2943 : : GRA_NearestNeighbour, 0.2, nullptr );
2944 : : }
2945 : :
2946 : 0 : if ( !mGdalDataset )
2947 : : {
2948 : 0 : QgsLogger::warning( QStringLiteral( "Warped VRT Creation failed." ) );
2949 : 0 : mGdalDataset = mGdalBaseDataset;
2950 : 0 : }
2951 : : else
2952 : : {
2953 : 0 : hasGeoTransform = GDALGetGeoTransform( mGdalDataset, mGeoTransform ) == CE_None;
2954 : : }
2955 : 0 : }
2956 : : else
2957 : : {
2958 : 0 : mGdalDataset = mGdalBaseDataset;
2959 : : }
2960 : :
2961 : 0 : if ( !mGdalTransformerArg )
2962 : : {
2963 : 0 : mGdalTransformerArg = GDALCreateGenImgProjTransformer( mGdalBaseDataset, nullptr, nullptr, nullptr, TRUE, 1.0, 0 );
2964 : 0 : }
2965 : :
2966 : 0 : if ( !hasGeoTransform )
2967 : : {
2968 : : // Initialize the affine transform matrix
2969 : 0 : mGeoTransform[0] = 0;
2970 : 0 : mGeoTransform[1] = 1;
2971 : 0 : mGeoTransform[2] = 0;
2972 : 0 : mGeoTransform[3] = 0;
2973 : 0 : mGeoTransform[4] = 0;
2974 : 0 : mGeoTransform[5] = -1;
2975 : 0 : }
2976 : :
2977 : : // get sublayers
2978 : 0 : mSubLayers = QgsGdalProvider::subLayers( mGdalDataset );
2979 : :
2980 : : // check if this file has bands or subdatasets
2981 : 0 : CPLErrorReset();
2982 : 0 : GDALRasterBandH myGDALBand = GDALGetRasterBand( mGdalDataset, 1 ); //just use the first band
2983 : 0 : if ( !myGDALBand )
2984 : : {
2985 : 0 : QString msg = QString::fromUtf8( CPLGetLastErrorMsg() );
2986 : :
2987 : : // if there are no subdatasets, then close the dataset
2988 : 0 : if ( mSubLayers.isEmpty() )
2989 : : {
2990 : 0 : appendError( ERRMSG( tr( "Cannot get GDAL raster band: %1" ).arg( msg ) ) );
2991 : :
2992 : 0 : closeDataset();
2993 : 0 : return;
2994 : : }
2995 : : // if there are subdatasets, leave the dataset open for subsequent queries
2996 : : else
2997 : : {
2998 : 0 : QgsDebugMsg( QObject::tr( "Cannot get GDAL raster band: %1" ).arg( msg ) +
2999 : : QString( " but dataset has %1 subdatasets" ).arg( mSubLayers.size() ) );
3000 : 0 : mValid = false;
3001 : 0 : return;
3002 : : }
3003 : 0 : }
3004 : :
3005 : : // check if this file has pyramids
3006 : 0 : mHasPyramids = gdalGetOverviewCount( myGDALBand ) > 0;
3007 : :
3008 : : // Get the layer's projection info and set up the
3009 : : // QgsCoordinateTransform for this layer
3010 : : // NOTE: we must do this before metadata is called
3011 : 0 : QString crsWkt;
3012 : 0 : if ( OGRSpatialReferenceH spatialRefSys = GDALGetSpatialRef( mGdalDataset ) )
3013 : : {
3014 : 0 : crsWkt = QgsOgrUtils::OGRSpatialReferenceToWkt( spatialRefSys );
3015 : 0 : }
3016 : 0 : if ( crsWkt.isEmpty() )
3017 : : {
3018 : 0 : if ( OGRSpatialReferenceH spatialRefSys = GDALGetGCPSpatialRef( mGdalDataset ) )
3019 : : {
3020 : 0 : crsWkt = QgsOgrUtils::OGRSpatialReferenceToWkt( spatialRefSys );
3021 : 0 : }
3022 : 0 : }
3023 : 0 : if ( !crsWkt.isEmpty() )
3024 : : {
3025 : 0 : mCrs = QgsCoordinateReferenceSystem::fromWkt( crsWkt );
3026 : 0 : }
3027 : : else
3028 : : {
3029 : 0 : if ( mGdalBaseDataset != mGdalDataset &&
3030 : 0 : GDALGetMetadata( mGdalBaseDataset, "RPC" ) )
3031 : : {
3032 : : // Warped VRT of RPC is in EPSG:4326
3033 : 0 : mCrs = QgsCoordinateReferenceSystem::fromOgcWmsCrs( QStringLiteral( "EPSG:4326" ) );
3034 : 0 : }
3035 : : else
3036 : : {
3037 : 0 : QgsDebugMsgLevel( QStringLiteral( "No valid CRS identified" ), 2 );
3038 : : }
3039 : : }
3040 : :
3041 : : //set up the coordinat transform - in the case of raster this is mainly used to convert
3042 : : //the inverese projection of the map extents of the canvas when zooming in etc. so
3043 : : //that they match the coordinate system of this layer
3044 : :
3045 : : //metadata();
3046 : :
3047 : : // Use the affine transform to get geo coordinates for
3048 : : // the corners of the raster
3049 : 0 : double myXMax = mGeoTransform[0] +
3050 : 0 : GDALGetRasterXSize( mGdalDataset ) * mGeoTransform[1] +
3051 : 0 : GDALGetRasterYSize( mGdalDataset ) * mGeoTransform[2];
3052 : 0 : double myYMin = mGeoTransform[3] +
3053 : 0 : GDALGetRasterXSize( mGdalDataset ) * mGeoTransform[4] +
3054 : 0 : GDALGetRasterYSize( mGdalDataset ) * mGeoTransform[5];
3055 : :
3056 : 0 : mExtent.setXMaximum( myXMax );
3057 : : // The affine transform reduces to these values at the
3058 : : // top-left corner of the raster
3059 : 0 : mExtent.setXMinimum( mGeoTransform[0] );
3060 : 0 : mExtent.setYMaximum( mGeoTransform[3] );
3061 : 0 : mExtent.setYMinimum( myYMin );
3062 : :
3063 : : //
3064 : : // Set up the x and y dimensions of this raster layer
3065 : : //
3066 : 0 : mWidth = GDALGetRasterXSize( mGdalDataset );
3067 : 0 : mHeight = GDALGetRasterYSize( mGdalDataset );
3068 : :
3069 : : // Check if the dataset has a mask band, that applies to the whole dataset
3070 : : // If so then expose it as an alpha band.
3071 : 0 : int nMaskFlags = GDALGetMaskFlags( myGDALBand );
3072 : 0 : const int bandCount = GDALGetRasterCount( mGdalDataset );
3073 : 0 : if ( ( nMaskFlags == 0 && bandCount == 1 ) || nMaskFlags == GMF_PER_DATASET )
3074 : : {
3075 : 0 : mMaskBandExposedAsAlpha = true;
3076 : 0 : }
3077 : :
3078 : 0 : mBandCount = bandCount + ( mMaskBandExposedAsAlpha ? 1 : 0 );
3079 : :
3080 : 0 : GDALGetBlockSize( GDALGetRasterBand( mGdalDataset, 1 ), &mXBlockSize, &mYBlockSize );
3081 : : //
3082 : : // Determine the nodata value and data type
3083 : : //
3084 : : //mValidNoDataValue = true;
3085 : 0 : for ( int i = 1; i <= bandCount; i++ )
3086 : : {
3087 : 0 : GDALRasterBandH myGdalBand = GDALGetRasterBand( mGdalDataset, i );
3088 : 0 : GDALDataType myGdalDataType = GDALGetRasterDataType( myGdalBand );
3089 : :
3090 : 0 : int isValid = false;
3091 : 0 : double myNoDataValue = GDALGetRasterNoDataValue( myGdalBand, &isValid );
3092 : : // We check that the double value we just got is representable in the
3093 : : // data type. In normal situations this should not be needed, but it happens
3094 : : // to have 8bit TIFFs with nan as the nodata value. If not checking against
3095 : : // the min/max bounds, it would be cast to 0 by representableValue().
3096 : 0 : if ( isValid && !QgsRaster::isRepresentableValue( myNoDataValue, dataTypeFromGdal( myGdalDataType ) ) )
3097 : : {
3098 : 0 : QgsDebugMsgLevel( QStringLiteral( "GDALGetRasterNoDataValue = %1 is not representable in data type, so ignoring it" ).arg( myNoDataValue ), 2 );
3099 : 0 : isValid = false;
3100 : 0 : }
3101 : 0 : if ( isValid )
3102 : : {
3103 : 0 : QgsDebugMsgLevel( QStringLiteral( "GDALGetRasterNoDataValue = %1" ).arg( myNoDataValue ), 2 );
3104 : : // The no data value double may be non representable by data type, it can result
3105 : : // in problems if that value is used to represent additional user defined no data
3106 : : // see #3840
3107 : 0 : myNoDataValue = QgsRaster::representableValue( myNoDataValue, dataTypeFromGdal( myGdalDataType ) );
3108 : 0 : mSrcNoDataValue.append( myNoDataValue );
3109 : 0 : mSrcHasNoDataValue.append( true );
3110 : 0 : mUseSrcNoDataValue.append( true );
3111 : 0 : }
3112 : : else
3113 : : {
3114 : 0 : mSrcNoDataValue.append( std::numeric_limits<double>::quiet_NaN() );
3115 : 0 : mSrcHasNoDataValue.append( false );
3116 : 0 : mUseSrcNoDataValue.append( false );
3117 : : }
3118 : : // It may happen that nodata value given by GDAL is wrong and it has to be
3119 : : // disabled by user, in that case we need another value to be used for nodata
3120 : : // (for reprojection for example) -> always internally represent as wider type
3121 : : // with mInternalNoDataValue in reserve.
3122 : : // Not used
3123 : : #if 0
3124 : : int myInternalGdalDataType = myGdalDataType;
3125 : : double myInternalNoDataValue = 123;
3126 : : switch ( srcDataType( i ) )
3127 : : {
3128 : : case Qgis::Byte:
3129 : : myInternalNoDataValue = -32768.0;
3130 : : myInternalGdalDataType = GDT_Int16;
3131 : : break;
3132 : : case Qgis::Int16:
3133 : : myInternalNoDataValue = -2147483648.0;
3134 : : myInternalGdalDataType = GDT_Int32;
3135 : : break;
3136 : : case Qgis::UInt16:
3137 : : myInternalNoDataValue = -2147483648.0;
3138 : : myInternalGdalDataType = GDT_Int32;
3139 : : break;
3140 : : case Qgis::Int32:
3141 : : // We believe that such values is no used in real data
3142 : : myInternalNoDataValue = -2147483648.0;
3143 : : break;
3144 : : case Qgis::UInt32:
3145 : : // We believe that such values is no used in real data
3146 : : myInternalNoDataValue = 4294967295.0;
3147 : : break;
3148 : : default: // Float32, Float64
3149 : : //myNoDataValue = std::numeric_limits<int>::max();
3150 : : // NaN should work well
3151 : : myInternalNoDataValue = std::numeric_limits<double>::quiet_NaN();
3152 : : }
3153 : : #endif
3154 : : //mGdalDataType.append( myInternalGdalDataType );
3155 : :
3156 : : // define if the band has scale and offset to apply
3157 : 0 : double myScale = bandScale( i );
3158 : 0 : double myOffset = bandOffset( i );
3159 : 0 : if ( !qgsDoubleNear( myScale, 1.0 ) || !qgsDoubleNear( myOffset, 0.0 ) )
3160 : : {
3161 : : // if the band has scale or offset to apply change dataType
3162 : 0 : switch ( myGdalDataType )
3163 : : {
3164 : : case GDT_Unknown:
3165 : : case GDT_TypeCount:
3166 : 0 : break;
3167 : : case GDT_Byte:
3168 : : case GDT_UInt16:
3169 : : case GDT_Int16:
3170 : : case GDT_UInt32:
3171 : : case GDT_Int32:
3172 : : case GDT_Float32:
3173 : : case GDT_CInt16:
3174 : 0 : myGdalDataType = GDT_Float32;
3175 : 0 : break;
3176 : : case GDT_Float64:
3177 : : case GDT_CInt32:
3178 : : case GDT_CFloat32:
3179 : 0 : myGdalDataType = GDT_Float64;
3180 : 0 : break;
3181 : : case GDT_CFloat64:
3182 : 0 : break;
3183 : : }
3184 : 0 : }
3185 : :
3186 : 0 : mGdalDataType.append( myGdalDataType );
3187 : : //mInternalNoDataValue.append( myInternalNoDataValue );
3188 : 0 : }
3189 : :
3190 : 0 : if ( mMaskBandExposedAsAlpha )
3191 : : {
3192 : 0 : mSrcNoDataValue.append( std::numeric_limits<double>::quiet_NaN() );
3193 : 0 : mSrcHasNoDataValue.append( false );
3194 : 0 : mUseSrcNoDataValue.append( false );
3195 : 0 : mGdalDataType.append( GDT_Byte );
3196 : 0 : }
3197 : 0 : }
3198 : :
3199 : 0 : QgsGdalProvider *QgsGdalProviderMetadata::createRasterDataProvider(
3200 : : const QString &uri,
3201 : : const QString &format,
3202 : : int nBands,
3203 : : Qgis::DataType type,
3204 : : int width, int height,
3205 : : double *geoTransform,
3206 : : const QgsCoordinateReferenceSystem &crs,
3207 : : const QStringList &createOptions )
3208 : : {
3209 : : //get driver
3210 : 0 : GDALDriverH driver = GDALGetDriverByName( format.toLocal8Bit().data() );
3211 : 0 : if ( !driver )
3212 : : {
3213 : 0 : QgsError error( "Cannot load GDAL driver " + format, QStringLiteral( "GDAL provider" ) );
3214 : 0 : return new QgsGdalProvider( uri, error );
3215 : 0 : }
3216 : :
3217 : 0 : QgsDebugMsgLevel( "create options: " + createOptions.join( " " ), 2 );
3218 : :
3219 : : //create dataset
3220 : 0 : CPLErrorReset();
3221 : 0 : char **papszOptions = QgsGdalUtils::papszFromStringList( createOptions );
3222 : 0 : gdal::dataset_unique_ptr dataset( GDALCreate( driver, uri.toUtf8().constData(), width, height, nBands, ( GDALDataType )type, papszOptions ) );
3223 : 0 : CSLDestroy( papszOptions );
3224 : 0 : if ( !dataset )
3225 : : {
3226 : 0 : QgsError error( QStringLiteral( "Cannot create new dataset %1:\n%2" ).arg( uri, QString::fromUtf8( CPLGetLastErrorMsg() ) ), QStringLiteral( "GDAL provider" ) );
3227 : 0 : QgsDebugMsg( error.summary() );
3228 : 0 : return new QgsGdalProvider( uri, error );
3229 : 0 : }
3230 : :
3231 : 0 : GDALSetGeoTransform( dataset.get(), geoTransform );
3232 : 0 : GDALSetProjection( dataset.get(), crs.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED_GDAL ).toLocal8Bit().data() );
3233 : :
3234 : 0 : QgsDataProvider::ProviderOptions providerOptions;
3235 : 0 : return new QgsGdalProvider( uri, providerOptions, true, dataset.release() );
3236 : 0 : }
3237 : :
3238 : 0 : bool QgsGdalProvider::write( void *data, int band, int width, int height, int xOffset, int yOffset )
3239 : : {
3240 : 0 : QMutexLocker locker( mpMutex );
3241 : 0 : if ( !initIfNeeded() )
3242 : 0 : return false;
3243 : :
3244 : 0 : if ( !mGdalDataset )
3245 : : {
3246 : 0 : return false;
3247 : : }
3248 : 0 : GDALRasterBandH rasterBand = getBand( band );
3249 : 0 : if ( !rasterBand )
3250 : : {
3251 : 0 : return false;
3252 : : }
3253 : 0 : return gdalRasterIO( rasterBand, GF_Write, xOffset, yOffset, width, height, data, width, height, GDALGetRasterDataType( rasterBand ), 0, 0 ) == CE_None;
3254 : 0 : }
3255 : :
3256 : 0 : bool QgsGdalProvider::setNoDataValue( int bandNo, double noDataValue )
3257 : : {
3258 : 0 : QMutexLocker locker( mpMutex );
3259 : 0 : if ( !initIfNeeded() )
3260 : 0 : return false;
3261 : :
3262 : 0 : if ( !mGdalDataset )
3263 : : {
3264 : 0 : return false;
3265 : : }
3266 : :
3267 : 0 : GDALRasterBandH rasterBand = getBand( bandNo );
3268 : 0 : CPLErrorReset();
3269 : 0 : CPLErr err = GDALSetRasterNoDataValue( rasterBand, noDataValue );
3270 : 0 : if ( err != CPLE_None )
3271 : : {
3272 : 0 : QgsDebugMsg( QStringLiteral( "Cannot set no data value" ) );
3273 : 0 : return false;
3274 : : }
3275 : 0 : mSrcNoDataValue[bandNo - 1] = noDataValue;
3276 : 0 : mSrcHasNoDataValue[bandNo - 1] = true;
3277 : 0 : mUseSrcNoDataValue[bandNo - 1] = true;
3278 : 0 : return true;
3279 : 0 : }
3280 : :
3281 : 0 : bool QgsGdalProvider::remove()
3282 : : {
3283 : 0 : QMutexLocker locker( mpMutex );
3284 : 0 : if ( !initIfNeeded() )
3285 : 0 : return false;
3286 : :
3287 : 0 : while ( *mpRefCounter != 1 )
3288 : : {
3289 : 0 : QgsDebugMsgLevel( QStringLiteral( "Waiting for ref counter for %1 to drop to 1" ).arg( dataSourceUri() ), 2 );
3290 : 0 : QThread::msleep( 100 );
3291 : : }
3292 : :
3293 : 0 : if ( mGdalDataset )
3294 : : {
3295 : 0 : GDALDriverH driver = GDALGetDriverByName( mDriverName.toLocal8Bit().constData() );
3296 : 0 : if ( !driver )
3297 : 0 : return false;
3298 : :
3299 : 0 : closeDataset();
3300 : :
3301 : 0 : CPLErrorReset();
3302 : 0 : CPLErr err = GDALDeleteDataset( driver, dataSourceUri( true ).toUtf8().constData() );
3303 : 0 : if ( err != CPLE_None )
3304 : : {
3305 : 0 : QgsLogger::warning( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
3306 : 0 : QgsDebugMsg( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
3307 : 0 : return false;
3308 : : }
3309 : 0 : QgsDebugMsgLevel( QStringLiteral( "Raster dataset dataSourceUri() successfully deleted" ), 2 );
3310 : 0 : return true;
3311 : : }
3312 : 0 : return false;
3313 : 0 : }
3314 : :
3315 : : /**
3316 : : Builds the list of file filter strings to later be used by
3317 : : QgisApp::addRasterLayer()
3318 : :
3319 : : We query GDAL for a list of supported raster formats; we then build
3320 : : a list of file filter strings from that list. We return a string
3321 : : that contains this list that is suitable for use in a
3322 : : QFileDialog::getOpenFileNames() call.
3323 : :
3324 : : */
3325 : 15 : QString QgsGdalProviderMetadata::filters( FilterType type )
3326 : : {
3327 : 15 : switch ( type )
3328 : : {
3329 : : case QgsProviderMetadata::FilterType::FilterRaster:
3330 : : {
3331 : 3 : QString fileFiltersString;
3332 : 3 : QStringList exts;
3333 : 3 : QStringList wildcards;
3334 : 3 : buildSupportedRasterFileFilterAndExtensions( fileFiltersString, exts, wildcards );
3335 : 3 : return fileFiltersString;
3336 : 3 : }
3337 : :
3338 : : case QgsProviderMetadata::FilterType::FilterVector:
3339 : : case QgsProviderMetadata::FilterType::FilterMesh:
3340 : : case QgsProviderMetadata::FilterType::FilterMeshDataset:
3341 : : case QgsProviderMetadata::FilterType::FilterPointCloud:
3342 : 12 : return QString();
3343 : : }
3344 : 0 : return QString();
3345 : 15 : }
3346 : :
3347 : 0 : QString QgsGdalProvider::validateCreationOptions( const QStringList &createOptions, const QString &format )
3348 : : {
3349 : 0 : QString message;
3350 : :
3351 : : // first validate basic syntax with GDALValidateCreationOptions
3352 : 0 : message = QgsGdalUtils::validateCreationOptionsFormat( createOptions, format );
3353 : 0 : if ( !message.isNull() )
3354 : 0 : return message;
3355 : :
3356 : : // next do specific validations, depending on format and dataset
3357 : : // only check certain destination formats
3358 : 0 : QStringList formatsCheck;
3359 : 0 : formatsCheck << QStringLiteral( "gtiff" );
3360 : 0 : if ( ! formatsCheck.contains( format.toLower() ) )
3361 : 0 : return QString();
3362 : :
3363 : : // prepare a map for easier lookup
3364 : 0 : QMap< QString, QString > optionsMap;
3365 : 0 : const auto constCreateOptions = createOptions;
3366 : 0 : for ( const QString &option : constCreateOptions )
3367 : : {
3368 : 0 : QStringList opt = option.split( '=' );
3369 : 0 : optionsMap[ opt[0].toUpper()] = opt[1];
3370 : 0 : QgsDebugMsgLevel( "option: " + option, 2 );
3371 : 0 : }
3372 : :
3373 : : // gtiff files - validate PREDICTOR option
3374 : : // see gdal: frmts/gtiff/geotiff.cpp and libtiff: tif_predict.c)
3375 : 0 : if ( format.compare( QLatin1String( "gtiff" ), Qt::CaseInsensitive ) == 0 && optionsMap.contains( QStringLiteral( "PREDICTOR" ) ) )
3376 : : {
3377 : 0 : QString value = optionsMap.value( QStringLiteral( "PREDICTOR" ) );
3378 : 0 : GDALDataType nDataType = ( !mGdalDataType.isEmpty() ) ? ( GDALDataType ) mGdalDataType.at( 0 ) : GDT_Unknown;
3379 : 0 : int nBitsPerSample = nDataType != GDT_Unknown ? GDALGetDataTypeSize( nDataType ) : 0;
3380 : 0 : QgsDebugMsgLevel( QStringLiteral( "PREDICTOR: %1 nbits: %2 type: %3" ).arg( value ).arg( nBitsPerSample ).arg( ( GDALDataType ) mGdalDataType.at( 0 ) ), 2 );
3381 : : // PREDICTOR=2 only valid for 8/16/32 bits per sample
3382 : : // TODO check for NBITS option (see geotiff.cpp)
3383 : 0 : if ( value == QLatin1String( "2" ) )
3384 : : {
3385 : 0 : if ( nBitsPerSample != 8 && nBitsPerSample != 16 &&
3386 : 0 : nBitsPerSample != 32 )
3387 : : {
3388 : 0 : message = QStringLiteral( "PREDICTOR=%1 only valid for 8/16/32 bits per sample (using %2)" ).arg( value ).arg( nBitsPerSample );
3389 : 0 : }
3390 : 0 : }
3391 : : // PREDICTOR=3 only valid for float/double precision
3392 : 0 : else if ( value == QLatin1String( "3" ) )
3393 : : {
3394 : 0 : if ( nDataType != GDT_Float32 && nDataType != GDT_Float64 )
3395 : 0 : message = QStringLiteral( "PREDICTOR=3 only valid for float/double precision" );
3396 : 0 : }
3397 : 0 : }
3398 : :
3399 : 0 : return message;
3400 : 0 : }
3401 : :
3402 : 0 : QString QgsGdalProvider::validatePyramidsConfigOptions( QgsRaster::RasterPyramidsFormat pyramidsFormat,
3403 : : const QStringList &configOptions, const QString &fileFormat )
3404 : : {
3405 : : // Erdas Imagine format does not support config options
3406 : 0 : if ( pyramidsFormat == QgsRaster::PyramidsErdas )
3407 : : {
3408 : 0 : if ( ! configOptions.isEmpty() )
3409 : 0 : return QStringLiteral( "Erdas Imagine format does not support config options" );
3410 : : else
3411 : 0 : return QString();
3412 : : }
3413 : : // Internal pyramids format only supported for gtiff/georaster/hfa/jp2kak/mrsid/nitf files
3414 : 0 : else if ( pyramidsFormat == QgsRaster::PyramidsInternal )
3415 : : {
3416 : 0 : QStringList supportedFormats;
3417 : 0 : supportedFormats << QStringLiteral( "gtiff" ) << QStringLiteral( "georaster" ) << QStringLiteral( "hfa" ) << QStringLiteral( "gpkg" ) << QStringLiteral( "rasterlite" ) << QStringLiteral( "nitf" );
3418 : 0 : if ( ! supportedFormats.contains( fileFormat.toLower() ) )
3419 : 0 : return QStringLiteral( "Internal pyramids format only supported for gtiff/georaster/gpkg/rasterlite/nitf files (using %1)" ).arg( fileFormat );
3420 : 0 : }
3421 : : else
3422 : : {
3423 : : // for gtiff external pyramids, validate gtiff-specific values
3424 : : // PHOTOMETRIC_OVERVIEW=YCBCR requires a source raster with only 3 bands (RGB)
3425 : 0 : if ( configOptions.contains( QStringLiteral( "PHOTOMETRIC_OVERVIEW=YCBCR" ) ) )
3426 : : {
3427 : 0 : if ( GDALGetRasterCount( mGdalDataset ) != 3 )
3428 : 0 : return QStringLiteral( "PHOTOMETRIC_OVERVIEW=YCBCR requires a source raster with only 3 bands (RGB)" );
3429 : 0 : }
3430 : : }
3431 : :
3432 : 0 : return QString();
3433 : 0 : }
3434 : :
3435 : 0 : QgsPoint QgsGdalProvider::transformCoordinates( const QgsPoint &point, QgsRasterDataProvider::TransformType type )
3436 : : {
3437 : 0 : if ( !mGdalTransformerArg )
3438 : 0 : return QgsPoint();
3439 : :
3440 : : int success;
3441 : 0 : double x = point.x(), y = point.y(), z = point.is3D() ? point.z() : 0;
3442 : 0 : GDALUseTransformer( mGdalTransformerArg, type == TransformLayerToImage, 1, &x, &y, &z, &success );
3443 : 0 : if ( !success )
3444 : 0 : return QgsPoint();
3445 : :
3446 : 0 : return QgsPoint( x, y, z );
3447 : 0 : }
3448 : :
3449 : 0 : bool QgsGdalProvider::isEditable() const
3450 : : {
3451 : 0 : return mUpdate;
3452 : : }
3453 : :
3454 : 0 : bool QgsGdalProvider::setEditable( bool enabled )
3455 : : {
3456 : 0 : QMutexLocker locker( mpMutex );
3457 : 0 : if ( !initIfNeeded() )
3458 : 0 : return false;
3459 : :
3460 : 0 : if ( enabled == mUpdate )
3461 : 0 : return false;
3462 : :
3463 : 0 : if ( !mValid )
3464 : 0 : return false;
3465 : :
3466 : 0 : if ( mGdalDataset != mGdalBaseDataset )
3467 : 0 : return false; // ignore the case of warped VRT for now (more complicated setup)
3468 : :
3469 : 0 : while ( *mpRefCounter != 1 )
3470 : : {
3471 : 0 : QgsDebugMsgLevel( QStringLiteral( "Waiting for ref counter for %1 to drop to 1" ).arg( dataSourceUri() ), 2 );
3472 : 0 : QThread::msleep( 100 );
3473 : : }
3474 : :
3475 : 0 : closeDataset();
3476 : :
3477 : 0 : mUpdate = enabled;
3478 : :
3479 : : // reopen the dataset
3480 : 0 : mGdalBaseDataset = gdalOpen( dataSourceUri( true ), mUpdate ? GDAL_OF_UPDATE : GDAL_OF_READONLY );
3481 : 0 : if ( !mGdalBaseDataset )
3482 : : {
3483 : 0 : QString msg = QStringLiteral( "Cannot reopen GDAL dataset %1:\n%2" ).arg( dataSourceUri(), QString::fromUtf8( CPLGetLastErrorMsg() ) );
3484 : 0 : appendError( ERRMSG( msg ) );
3485 : 0 : return false;
3486 : 0 : }
3487 : :
3488 : : //Since we are not a virtual warped dataset, mGdalDataSet and mGdalBaseDataset are supposed to be the same
3489 : 0 : mGdalDataset = mGdalBaseDataset;
3490 : 0 : mValid = true;
3491 : 0 : return true;
3492 : 0 : }
3493 : :
3494 : 0 : GDALRasterBandH QgsGdalProvider::getBand( int bandNo ) const
3495 : : {
3496 : 0 : QMutexLocker locker( mpMutex );
3497 : 0 : if ( !const_cast<QgsGdalProvider *>( this )->initIfNeeded() )
3498 : 0 : return nullptr;
3499 : :
3500 : 0 : if ( mMaskBandExposedAsAlpha && bandNo == GDALGetRasterCount( mGdalDataset ) + 1 )
3501 : 0 : return GDALGetMaskBand( GDALGetRasterBand( mGdalDataset, 1 ) );
3502 : : else
3503 : 0 : return GDALGetRasterBand( mGdalDataset, bandNo );
3504 : 0 : }
3505 : :
3506 : : // pyramids resampling
3507 : :
3508 : : // see http://www.gdal.org/gdaladdo.html
3509 : : // http://www.gdal.org/classGDALDataset.html#a2aa6f88b3bbc840a5696236af11dde15
3510 : : // http://www.gdal.org/classGDALRasterBand.html#afaea945b13ec9c86c2d783b883c68432
3511 : :
3512 : : // from http://www.gdal.org/gdaladdo.html
3513 : : // average_mp is unsuitable for use thus not included
3514 : :
3515 : : // from qgsgdalprovider.cpp (removed)
3516 : : // NOTE magphase is disabled in the gui since it tends
3517 : : // to create corrupted images. The images can be repaired
3518 : : // by running one of the other resampling strategies below.
3519 : : // see ticket #284
3520 : 0 : QList<QPair<QString, QString> > QgsGdalProviderMetadata::pyramidResamplingMethods()
3521 : : {
3522 : 0 : static QList<QPair<QString, QString> > methods;
3523 : :
3524 : 0 : if ( methods.isEmpty() )
3525 : : {
3526 : 0 : methods.append( QPair<QString, QString>( QStringLiteral( "NEAREST" ), QObject::tr( "Nearest Neighbour" ) ) );
3527 : 0 : methods.append( QPair<QString, QString>( QStringLiteral( "AVERAGE" ), QObject::tr( "Average" ) ) );
3528 : 0 : methods.append( QPair<QString, QString>( QStringLiteral( "GAUSS" ), QObject::tr( "Gauss" ) ) );
3529 : 0 : methods.append( QPair<QString, QString>( QStringLiteral( "CUBIC" ), QObject::tr( "Cubic" ) ) );
3530 : 0 : methods.append( QPair<QString, QString>( QStringLiteral( "CUBICSPLINE" ), QObject::tr( "Cubic Spline" ) ) );
3531 : 0 : methods.append( QPair<QString, QString>( QStringLiteral( "LANCZOS" ), QObject::tr( "Lanczos" ) ) );
3532 : 0 : methods.append( QPair<QString, QString>( QStringLiteral( "BILINEAR" ), QObject::tr( "Bilinear" ) ) );
3533 : 0 : methods.append( QPair<QString, QString>( QStringLiteral( "MODE" ), QObject::tr( "Mode" ) ) );
3534 : 0 : methods.append( QPair<QString, QString>( QStringLiteral( "NONE" ), QObject::tr( "None" ) ) );
3535 : 0 : }
3536 : :
3537 : 0 : return methods;
3538 : 0 : }
3539 : :
3540 : 0 : QgsProviderMetadata::ProviderCapabilities QgsGdalProviderMetadata::providerCapabilities() const
3541 : : {
3542 : 0 : return FileBasedUris;
3543 : : }
3544 : :
3545 : 3 : QList<QgsDataItemProvider *> QgsGdalProviderMetadata::dataItemProviders() const
3546 : : {
3547 : 3 : QList< QgsDataItemProvider * > providers;
3548 : 3 : providers << new QgsGdalDataItemProvider;
3549 : 3 : return providers;
3550 : 3 : }
3551 : :
3552 : 3 : QgsGdalProviderMetadata::QgsGdalProviderMetadata():
3553 : 9 : QgsProviderMetadata( PROVIDER_KEY, PROVIDER_DESCRIPTION )
3554 : 3 : {
3555 : 3 : }
3556 : :
3557 : : ///@endcond
|
