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